Diffuse interface model for compressible fluid – Compressible elastic–plastic solid interaction

Diffuse interface model for compressible fluid – Compressible elastic–plastic solid interaction

Journal of Computational Physics 231 (2012) 2695–2723 Contents lists available at SciVerse ScienceDirect Journal of Computational Physics journal ho...

1MB Sizes 0 Downloads 62 Views

Journal of Computational Physics 231 (2012) 2695–2723

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Diffuse interface model for compressible fluid – Compressible elastic–plastic solid interaction N. Favrie, S.L. Gavrilyuk ⇑ Aix-Marseille University, UMR CNRS 6595 IUSTI, SMASH Project, INRIA, 5 rue E. Fermi, 13453 Marseille Cedex 13, France

a r t i c l e

i n f o

Article history: Received 29 April 2011 Received in revised form 7 November 2011 Accepted 19 November 2011 Available online 8 December 2011 Keywords: Diffuse solid–fluid interfaces Shocks Elastic–plastic solids Riemann problem Godunov type methods

a b s t r a c t An Eulerian hyperbolic diffuse interface model for elastic–plastic solid–fluid interaction is constructed. The system of governing equations couples Euler equations of compressible fluids and a visco-plastic model of Maxwell type materials (the deviatoric part of the stress tensor decreases during plastic deformations) in the same manner as models of multicomponent fluids. In particular, the model is able to create interfaces which were not present initially. The model is thermodynamically compatible: it verifies the entropy inequality. However, a numerical treatment of the model is particularly challenging. Indeed, the model is nonconservative, so a special numerical splitting is proposed to overcome this difficulty. The numerical algorithm contains two relaxation procedures. One of them is physical and is related to the plastic relaxation mechanism (relaxation toward the yield surface). The second one is numerical. It consists in replacing the algebraic equation expressing a mechanical equilibrium between components by a partial differential equation with a short relaxation time. The numerical method was tested in 1D case (Wilkins’ flying plate problem), 2D plane case (impact of a projectile on a plate) and axisymmetrical case (Taylor test problem, impact with penetration effects, etc.). Numerical examples show the ability of the model to deal with real physical phenomena. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Works by Karni [17], Abgrall [1] and Saurel and Abgrall [22] showed the attractivity of the diffuse interface approach to modeling interfaces between ideal compressible fluids having different thermodynamic characteristics. The corresponding Eulerian models are reminiscent of one-velocity multiphase flow models. The method presents several advantages compared to a direct coupling of models of homogeneous fluids through a sharp interface. 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 other phase in pure phases. So, no need to use the interface tracking approach. This kind of model can describe the generation of new interfaces without having to re-mesh the domain or at least destroy and create cells. The main drawback of the Eulerian diffuse interface approach compared to a Lagrangian formulation is that the interfaces are not stiff. Indeed, the ‘‘mixture sells’’ are always present at the vicinity of moving interfaces. The thickness of ‘‘the mixture region’’ being increased in time, the method can only be used for short times. Hence, a natural application of such methods is high-velocity process (impacts, explosive phase transitions, etc.).

⇑ Corresponding author. E-mail addresses: [email protected] (N. Favrie), [email protected] (S.L. Gavrilyuk). 0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.11.027

2696

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

This approach has been extended in Petitpas et al. [20], Saurel et al. [23] to fluids exhibiting phase transitions and heat exchanges, and in Favrie et al. [8] to the description of solid–fluid interfaces. The last model was limited only by elastic solids, the plasticity effects were not taken into account. In Favrie and Gavrilyuk [9,10] an extension of the ideal plasticity model by Godunov [14] to the description of isotropic elastic–plastic solids was proposed. This model describes, in particular, the evolution of the effective elastic deformation tensor (Finger tensor) in Eulerian coordinates. The stresses are given by a hyperelastic relation. A special equation of state in a separable form was proposed: it was the sum of a hydrodynamic part depending only on the density and the entropy, and the elastic part depending on the two other invariants of the Finger tensor. The model was hyperbolic under natural convexity conditions in all domain of parameters. When plastic deformations occurred, the model was in agreement with the following requirements:    

The mass is conserved and does not evolve during the relaxation process (the plastic response is isochoric). Entropy is created (irreversibility of plastic deformations). The intensity of shear stresses decreases. The plasticity yield limit is reached at the end of the relaxation process.

The first two requirements are quite classical. The third condition comes from experimental evidence: for isotropic solids, the shear stress decreases (see Godunov and Romenskii [15] for a detail description of such continua called also ‘‘Maxwell continua’’). Finally, the last condition means that the yield surface is a global attractor for initial data which are outside the elasticity domain. Such a property (to be an attractor) is a general requirement, but the construction of the corresponding relaxation system is specific to the choice of the plasticity criteria. The difficulty to guarantee such a property is in the fact that the yield criteria is formulated in terms of stresses, while the evolution equations are written in terms of elastic deformation tensor. In our case such a construction of the relaxation system is specific to the von Mises yield criterion. However, the approach can also be applied to other yield criteria. In the case of hypoelastic models (the evolution equations are for the stress tensor which is not related, in general, to a stored energy) such a construction can also be found in Després et al. [6]. A particular case of a Maxwell continuum with the zero yield strength was developed in Godunov [14], Godunov and Romenskii [15] and Godunov and Peshkov [16]. Such a model is usually applied for very high-rate deformation of metals. The model with a non-zero yield strength we use here covers a larger field of applications. We adapt the model of visco-plastic solids and the Euler equations of compressible fluids to diffuse interface formulation. Other mathematical models for dynamical transformation in elastic–plastic solids exist (Wilkins [30], Miller and Colella [19], Kluth and Després [18] and others). However, we could not adapt them to this formulation. This paper is organized as follows. In Section 2, the elastic solid–fluid diffuse interface model is recalled. Sections 3 and 4 generalize this model to the case of elastic–plastic solids–compressible fluids diffuse interfaces. Section 5 details the numerical method to solve this model. Section 6 presents computation results for Wilkins’ flying plate problem, a 2D plane impact problem, and axisymmetric impact problems (Taylor test problem and impact of two solids with penetration effects). 2. Elastic solid–fluid diffuse interface model The diffuse interface model for solid–fluid coupling is based on the averaged formulation able to describe the dynamics of pure components (solid and fluid), as well as the interaction at the interface between components. It is necessary that in the limit of vanishing volume fractions the model gives the governing equations of pure phases. We want to derive such a model permitting to account for the sliding at the solid–fluid interfaces and plastic deformations in the solid phase. We resume first the one-velocity multiphase elastic solid–fluid model derived in Favrie et al. [8]. Let X = (Xb), b = 1, 2, 3 be the Lagrangian coordinates of the mixture particles, x = (xn), n = 1, 2, 3 be the Eulerian coordinates. The corresponding cobasis vectors Eb = rXb satisfy the equations:

@Eb þ rðEb  v Þ ¼ 0; @t rotðEb Þ ¼ 0:

ð1Þ ð2Þ

Here v(t, x) is the velocity field. In particular, the condition (2) is satisfied at t = 0, it is satisfied at all times. A non-conservative form of (1) which takes into account (2) is:

 T DEb @v Eb ¼ 0; þ Dt @x

D @ ¼ þ v  r: Dt @t

ð3Þ

Here and hereafter, the classical differential operators r, div and rot are always applied to functions defined in the Eulerian D coordinates x = (xn), n = 1, 2, 3, the material derivative is denoted by Dt . We introduce the deformation gradient F and the T 1 Finger tensor G = F F defined as:

FT ¼ ðE1 ; E2 ; E3 Þ;



3 X b¼1

Eb  Eb :

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2697

In particular, (3) implies the equation for G:

 T DG @v @v þ GþG ¼ 0: Dt @x @x The density of the solid phase qs is defined by the following algebraic relation:

as qs ¼ as0 qs0 jGj1=2 ;

ð4Þ

where index ‘‘0’’ denotes the initial values, and as > 0 is the volume fraction of the solid. This algebraic equation implies the mass conservation law for the solid phase:

@ðas qs Þ þ divðas qs v Þ ¼ 0: @t

ð5Þ

The equation of mass for the gas component is postulated in the same way:

@ðag qg Þ þ divðag qg v Þ ¼ 0; @t

ð6Þ

where ag = 1  as > 0 is the gas volume fraction, and qg is the gas density. Eqs. (5) and (6) imply the mass conservation law of the mixture

@q þ divðqv Þ ¼ 0; @t where q = asqs + agqg is the mixture density. The Hamilton principle of stationary action was used to derive the governing system (see Section 3.2 ‘‘Reversible equilibrium model’’ in Favrie et al. [8]). This principle states that the motion of a mechanical system with fixed initial and final positions is, under some constraints, a stationary point of the Hamilton action which is the time integral of the difference between the kinetic energy and the potential energy of the system. The constraints are the conservation of mass and entropy in the phases. In particular, the corresponding equilibrium elastic solid–gas diffuse interface model obtained by such an approach is:

@Eb þ rðEb  v Þ ¼ 0; @t rotEb ¼ 0; @ðag qg Þ þ divðag qg v Þ ¼ 0; @t @ðas qs Þ þ divðas qs v Þ ¼ 0; @t @ qv þ divðqv  v  ðas rs þ ag rg ÞÞ ¼ 0; @t Dgg Dgs ¼ 0; ¼ 0; Dt Dt ps ¼ pg :

ð7aÞ ð7bÞ ð7cÞ ð7dÞ ð7eÞ ð7fÞ ð7gÞ

gs and gg are the entropies, ps and pg are the pressures in the solid and in the gas phase, respectively. The last equation ps = pg coming from the variation of the Hamilton action with respect to the volume fraction describes the microstructure of the mixture region: the distribution of the volume fraction in the diffusion zone follows the pressure equilibrium. Eqs. (7) imply the energy conservation law

@ qE þ divðv qE  ðas rs þ ag rg Þv Þ ¼ 0; @t

ð8Þ

where the total energy is



qE ¼ q e þ

v  v

;

2

q ¼ as qs þ ag qg ; qe ¼ as qs es þ ag qg eg ;

and the specific energies es = es(G, gs) and eg(qg, gg) satisfy the Gibbs identities:

hs dgs ¼ des  tr



 @es dG ; @G

hg dgg ¼ deg þ pg d

1

qg

! :

ð9Þ

Here hs and hg are the corresponding temperatures. Let us recall that contrary to the pure solid case, the density qs and the deformation tensor of the mixture G are not related to each other directly: they are related through the volume fraction (see (4)). The stress tensors in the solid rs and in the gas rg are defined, respectively, as:

2698

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

rs ¼ 2qs

@es G; @G

rg ¼ pg I:

A particular form of the energy is the equation of state in separable form [12]:

g ¼ G=jGj1=3 :

es ¼ ehs ðqs ; gs Þ þ ees ðgÞ;

ð10Þ

ehs

Here is the hydrodynamic part of the energy, and allows us to write:

ees

is the energy related to shear stresses. One can prove that such a form

rs ¼ ps I þ Ss ; trðSs Þ ¼ 0 with

@ehs ðqs ; gs Þ ; @ qs e @e Ss ¼ 2qs s G: @G ps ¼ q2s

Let us remark that the pressure ps is determined here only by the hydrodynamic part of the energy ehs . For isotropic solids, the energy responsible for shearing effects can be written as a function of the invariants of g.

ees ðgÞ ¼ ees ðn1 ; n2 Þ; where

n1 ¼ trðgÞ ¼

J1

n2 ¼ trðg2 Þ ¼

;

jGj1=3

J2 jGj2=3

;

J i ¼ trðGi Þ;

i ¼ 1; 2:

In terms of eigenvalues j1, j2 and j3 of G, the variables ni, i = 1, 2 can be expressed as

n1 ¼

j1 þ j2 þ j3 j2 þ j22 þ j23 ; n2 ¼ 1 : 1=3 ðj1 j2 j3 Þ ðj1 j2 j3 Þ2=3

Obviously, the tensor

Ss ¼ 2qs

   !  e    @ees 1 @ees J1 2 @ees J2 @es n1 @ees n2 2 2 G  G  q g  g  G ¼ 2qs I þ I ¼ 2 I þ 2 I s @G 3 3 @n1 3 @n2 3 jGj1=3 @n1 jGj2=3 @n2

has zero trace. Indeed, this general fact is a consequence of the fact that ees is a homogeneous function of degree zero of G. Particular forms of the equations of state used in applications are:

ehs ð

qs ; g s Þ ¼

ees ðgÞ ¼

As exp

ls 4qs0

0 tr@





qcs þ ðcs  1Þp1s ; ðcs  1Þqs

gs gs0 cv s

G jGj1=3

!2 1

I



ls 4qs0

trððg  IÞ2 Þ:

ð11Þ ð12Þ

Here As, gs0, p1s, cs > 1, cvs and ls are constants. In the limit of small deformations, the Hooke law is recovered. In the limit of large deformations, the solid behaves as a fluid. In particular, the corresponding stress tensor in the solid is:

rs ¼ ps I þ Ss ;    ! qs 1 J2 1 J1 2 Ss ¼ ls G  I  1=3 G  I ; trðSs Þ ¼ 0: qs0 jGj2=3 3 3 jGj

ð13Þ

Model (7) with equations of state (11) and (12) is hyperbolic at least in 1D case [8]. 3. Elastic–plastic solid–fluid diffuse interface model 3.1. Introduction of a plastic dissipation The introduction of dissipation is always a phenomenological procedure. However, the dissipative terms should be compatible with the mass conservation law and the entropy inequality. Also, they should be in agreement with additional physical characteristics of the dissipation process. For example, it is experimentally observed in isotropic solids that the shear stresses decrease during relaxation. The model should take into account this phenomena (Maxwell type model).

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2699

The dissipation due to plastic deformations in solids will be introduced in the same manner as in pure solids [14,15] in the case of zero yield strength, and Favrie and Gavrilyuk [9,10] in the case of finite yield strength). Let us take the governing equations in the form:

 T DG @v @v þG G ¼ U; þ Dt @x @x @ðas qs Þ þ divðas qs v Þ ¼ 0; @t @ðag qg Þ þ divðag qg v Þ ¼ 0; @t @ qv þ divðqv  v  rÞ ¼ 0; @t @ðqEÞ þ divðqv E  rv Þ ¼ 0: @t ps ¼ pg :

with U ¼

GR þ RG

srel

;

ð14Þ

The tensor G plays the role of the effective elastic deformation tensor (effective elastic Finger tensor) [14,15]). We abuse the notation by using the same notation G as in the case of pure elastic deformations. The relaxation term U is a symmetric second order tensor depending on an unknown second order tensor R, and srel is a relaxation time. They will be chosen later in accordance with the principles stated in the introduction. The formulation of governing equations in terms of the effective Finger tensor G is useful only for computation of smooth solutions, but it is useless for computation of discontinuous solutions (this equation can not be put in conservative form). So, a conservative formulation is needed. Let us try to find an equivalent form of the equation for G in terms of the effective cobasis vectors Eb. As in the elastic case, we suppose that



3 X

Eb  Eb :

b¼1

The equation on Eb can be taken in the form:

 T 3 X D b @v 1 E þ Eb þ Aba Ea ¼  REb Dt @x srel a¼1

ð15Þ

with Aab ¼ Aba being an antisymmetric matrix. Let us show that (15) implies the evolution equation for G. Indeed,

! 3 3 X DG X D b DEb DEb ¼ ðE  Eb Þ ¼  Eb þ Eb  Dt Dt Dt Dt b¼1 b¼1 ! !    T T 3 3 3 3 X b a X X X @v 1 @v 1 b a b b b b b b E þ Aa E þ RE  E  E  E þ Aa E þ RE ¼ 2srel 2srel @x @x b¼1 b¼1 a¼1 a¼1  T  T 3 X 3   X @v @v 1 @v @v 1 GG ðRG þ GRÞ þ Aba þ Aab Eb  Eb ¼  GG ðRG þ GRÞ: ¼   @x @x srel @x @x srel b¼1 a¼1

ð16Þ

The symmetric effective elastic Finger tensor has six independent components while the effective basis Eb, b = 1, 2, 3 has nine independent components. So, it is natural to have three additional degrees of freedom defined by the matrix Aba having three independent components. To determine the stress tensor, we only need to know the Finger tensor, so the choice of the matrix A will have no influence. However, a specific procedure is needed to model the plastic relaxation in terms of Eb. This procedure, based on the singular value decomposition, will be explained below. One can also calculate observable Finger tenb which is not related to stresses. It is defined in a classical way. We calculate first the particle trajectories: sor G

dx ¼ v; dt

xjt¼0 ¼ X:

Then we calculate the observable deformation gradient:

@x b F¼ @X and the observable Finger tensor:

b¼b F 1 : G F T b Both of them, effective elastic Finger tensor and observable Finger tensor are important for the description of plastic transformations. As noted in Godunov and Peshkov [16], the effective elastic Finger tensor contains information about stresses, but

2700

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

not about the observable deformations. The observable Finger tensor describes the observable deformations but does not contain any information about stresses. The hyperbolicity of this diffuse interface model in all domain of parameters was already proved (in the 1D case) in Favrie et al. [8] for the energy given in separable form (10). The source term U added here does not change the hyperbolicity property, however its form is important. Indeed, the symmetric tensor R = RT determining U and the relaxation time srel > 0 should be determined in such a way that:    

They should be compatible with the mass conservation law. They should be compatible with the entropy inequality. The intensity of shear stresses decreases during the relaxation process (Maxwell type model). They should be compatible with a yield criterion (in our case, it will be the von Mises yield criterion).

The proof of the hypotheses formulated here is reminiscent of that made in the case of pure elastic–plastic solids [9]. The difference is in the proof of the entropy inequality which will depend on the volume fraction evolution. 3.2. Compatibility with the mass conservation law Let us calculate the material derivative of the average density:

!!    T Dq Dðq0 jGj1=2 Þ DjGj1=2 q0 1=2 @jGj DG q0 1=2 @v @v 1 1 ¼ tr tr jGjG G G ðGR þ RGÞ ¼ q0 ¼ ¼ jGj jGj  Dt @G Dt Dt Dt 2 2 @x @x srel !!  T q @v @v 1 q ¼ tr G ðR þ G1 RGÞ ¼ qdivv  trðRÞ:   G1 2 @x @x srel srel Hence, for the mass conservation it is necessary that

trðRÞ ¼ 0: 3.3. Compatibility with the entropy inequality We transform the energy Eq. (8) to the entropy equation:

   @   v  v v  v q eþ þ div qv e þ  ðas rs þ ag rg Þv @t 2 2  e     De @v p p Das @es DG @es Dgs ¼q  tr ðas rs þ ag rg Þ þ Ys ¼ q Y s s divv  Y s s þ Y s tr Dt @x qs as qs Dt @G Dt @ gs Dt !   pg pg Dag @eg Dgg @v þ Yg  tr ðas rs  ag pg IÞ Y g divv  Y g qg ag qg Dt @ gg Dt @x



Das Dag ¼ ðas ps þ ag pg Þdivv  ps  pg Dt Dt  e   ! @es DG @es Dgs @eg Dgg @v þ q Y s tr þ Ys þ Yg  tr ðas rs  ag pg IÞ @G Dt @ gs Dt @ gg Dt @x ¼ ðas ps þ ag pg Þdivv     ! @ees @v R @es Dgs @eg Dgg @v þ Ys þ as tr 2qs G þ þ Yg  tr ðas rs  ag pg IÞ @G @x srel @ gs Dt @ gg Dt @x ¼ ðas ps þ ag pg Þdivv  !    @v R @es Dgs @eg Dgg @v þ Ys þ tr as Ss þ þ Yg  tr ðas rs  ag pg IÞ @x srel @ gs Dt @ gg Dt @x ! 1 @es Dgs @eg Dgg trðas Ss RÞ þ Y s ¼ þ Yg srel @ gs Dt @ gg Dt ¼ a s qs hs

Dgg as Dgs trðSs RÞ: þ ag qg hg þ Dt Dt srel

Finally, we obtain the entropy balance law:

as qs hs

Dgg as Dgs trðSs RÞ ¼ 0: þ ag qg hg þ srel Dt Dt

ð17Þ

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2701

Eq. (17) should be separated into two equations for the entropy evolution of each component. Such a separation is also phenomenological. However, it should respect the entropy inequality and basic physical principles. Since there is no reason that plasticity deformations have influence on the entropy of the gas, we can use the following separation:

as qs hs

Dgs as trðSs RÞ; ¼ Dt srel

ag qg hg

Dgg ¼ 0: Dt

Obviously, they are compatible with (17): the sum of these equations gives us Eq. (17). We consider a simplest possibility to satisfy the entropy inequality by taking

R ¼ aSs ; where a is a positive constant. It implies the positivity of the term



as 2srel

trðSs RÞ:

Finally, the entropy inequality will be satisfied:

@ðqgÞ þ divðqgv Þ P 0; @t where g is the mixture entropy:

g ¼ Y s gs þ Y g gg : 3.4. Compatibility with the von Mises yield criterion The von Mises yield criterion says that the material starts yielding when the corresponding yield function

2 f ðSs Þ ¼ Ss : Ss  Y 2 3 becomes positive. The surface f(Ss) = 0 is called yield surface. Here Y is the yield strength, and Ss is the deviatoric part of the stress tensor:

rs ¼ ps I þ Ss : If the yield function is negative, the material is elastic, no plastic deformations will occur. In the elasticity domain we use srel = 1. Other plasticity criteria could be used but we will only deal here with the von Mises yield criterion. The relaxation procedure we describe here is in somewhat a numerical one. As we already told before, the resolution of Eq. (15) is preferable than the resolution of Eq. (16) for G. Indeed, one can use the splitting procedure for solving these equations: the relaxation step will follow the hyperbolic step. If during the hyperbolic step the material state is in the plasticity domain, after the relaxation step it should be at the yield surface. At the hyperbolic step, we solve the equation (the procedure will be described later):

D b @ vT E þ ¼ 0; Dt @x while the relaxation step corresponds to the solution of the system 3 @ b X a E þ Aba Ea ¼ S Eb : @t srel s a¼1

ð18Þ

These equations need a closure because the coefficients Aba are unknown. We will now detail some hypotheses to reduce the system of nine differential Eq. (18) to a simpler system of only three equations (Section 3.4.1). Also, we will propose a simplification reducing the resolution of (18) to the resolution of an algebraic system of equations (Section 3.4.4). 3.4.1. Singular value decomposition Eq. (18) can be rewritten as an equation for the effective deformation gradient FT = (E1, E2, E3) (the vectors Ea are the columns of the matrix FT):

@ T a S FT F þB¼ @t srel s

ð19Þ

  with matrix B defined as B ¼ A12 E2 þ A13 E3 ; A12 E1 þ A23 E3 ; A13 E1  A23 E2 . Consider the singular value decomposition of FT: FT = VKUT, where U and V are orthogonal matrices, and K is the diagonal matrix (see, for example, [15]). The singular values ja of the matrix G = FT F1 (which are also the eigenvalues of G) are related to the singular values ka of the matrix FT by the following obvious formula:

2702

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

ja ¼ k2a : Indeed, it follows from the identity

G ¼ VK2 VT ¼ VK2 V1 : The deviatoric part of the stress tensor depends only on G, so it can also be rewritten under the form Ss = VDVT, where D is a diagonal matrix with the eigenvalues Sa, a = 1, 2, 3. Eq. (19) can be rewritten under the following form:

! @K @UT aDK T @V T þ V KþK : U þ V BU ¼ @t @t @t srel

The following hypotheses will determine the ‘‘relaxation path’’ to simplify the computation of this equation:  The matrices V and U do not change during the relaxation process.  The matrix B vanishes. The first hypothes is was proposed by Godunov [14], Godunov and Romenskii [15] and Godunov and Peshkov [16]. These hypotheses are equivalent to the choice of the ‘‘plastic relaxation path’’ needed to avoid the indeterminacy of the relaxation Eq. (15) for Eb. In general, it does not correspond to the orthogonal projection onto the yield surface. With these hypothesis Eq. (19) becomes a coupled system of three equations:

dka a k S ; ¼ dt srel a a

ð20Þ

where ka are the singular values of FT. Here we replaced the partial derivative with respect to time by dtd because during the relaxation all variables depend only on time. 3.4.2. Lyapunov function System (20) coming from the singular value decomposition of FT has nice mathematical properties. First, it admits the mass conservation law:

@k1 k2 k3 ¼0 @t because S1 + S2 + S3 = 0. Second, it admits the Lyapunov function S:S for our specific choice of the energy (12) and (13). Indeed, in terms of the eigenvalues ja of G, the eigenvalues Sa of the deviatoric part of the stress tensor Ss are:





!

 q 1 j2 þ j22 þ j23 1 j þ j2 þ j3  Sa ¼ ls s j2a  1 ja  1  : 2=3 1=3 qs0 ðj1 j2 j3 Þ 3 3 ðj1 j2 j3 Þ In the following, we will denote

x ¼ j1 j2 j3 ¼ ðk1 k2 k3 Þ2 ¼



q q0

2

¼ jGj ¼ const:

We introduce new variables

ba ¼

ja ðj1 j2 j3 Þ1=3

:

In these variables the relaxation Eq. (20) are

dba 2a b S ; ¼ dt srel a a and the deviatoric components of the stress tensor are:

Sa ¼ ls x

! b21  b1 þ b22  b2 þ b23  b3 ba  ba  3    2  2 !  2 2 b  1 þ b2  12 þ b3  12 1 : ba   1 2 2 3 2

1=2

¼ ls x1=2

One can be proved that the function

L ¼ Ss : Ss ¼

X a

2

2 s

Sa ¼ l

0 1 4  2 !2 X 1 1 X 1 @ A x ba   ba  2 3 a 2 a

ð21Þ

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2703

is a Lyapunov function [9]:

dL 6 0: dt The existence of the Lyapunov allows us to guarantee the decreasing the tangential stresses during the plastic deformations. However, it does not allow us to guarantee that the material state will be at the yield surface at the end of the plastic relaxation process. 3.4.3. Relaxation of the stress tothe yield limit The invariant Ss : Ss ¼ tr S2s is decreasing during the relaxation step, so we have to chose the relaxation time function in a proper way to ‘‘stop’’ this decreasing. The choice we will propose is specific to the von Mises criterion. A modification is needed, if other plasticity criteria are used. However, the general idea is very simple: the relaxation time should vanish at the yield surface. For the von Mises yield criterion it can be done, for example, by using the following relaxation time function:

1

srel

¼

8 P n P 2 2 2 S 2  2Y 2 > > a a 3 < s1 ; if Sa  3 Y > 0; l2 0 a

> > : 0;

if

P a

ð22Þ

S2a  23 Y 2 6 0:

Here s0 is a characteristic constant time, and n > 0. With such a choice of the relaxation time the trajectories are attracted to the yield surface in finite time if 0 < n < 1, and in infinite time if n P 1. 3.4.4. Simplified relaxation equation for small deformations Even if the ordinary differential equations for ba can easily be solved numerically, it would be useful to have simple algebraic relations (exact solutions) of the time dependence of ba, to make numerical codes more efficient. Here we present such a simplification in the case of small deformations. Let ba = 1 + ca where ca is small. We can estimate now the principal eigenvalues of the deviatoric part of the stress tensor:

Sa ¼ ls x

b2  b1 þ b22  b2 þ b23  b3 ba  ba  1 3 2

1=2

!

¼ ls

  q c þ O c2a q0 a

Or, equivalently, q

2als q Sa dSa 0 ¼ : dt srel

ð23Þ

We choose n = 1 in (22) to have:

1

srel

¼

P

1

2

a Sa

s0

 23 Y 2

!

l2

¼

1

s0

X S2 a a

n 2

l

! 2

;

s0 ¼ const:; n2 ¼

2 Y2 : 3 l2

ð24Þ

Also, one can choose a such that

2als x1=2 ¼ 1: Then (23) and (24) admit the following explicit solution [9]:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 t=s 0 C 0  1e ffi; Sa ¼ Sa jt¼0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 C 0 e2n t=s0  1

ð25Þ

where

P C0 ¼ P

S2a j

t¼0

a l2

S2a j

t¼0

a l2

> 1:

 n2

This explicit formula reflects quite well the qualitative behavior of the deviatoric part of the stress tensor. In particular, if the characteristic time scale Dt is large with respect to the relaxation scale s0, the final values of Sa can be estimated as:

Saf

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi C0  1  pffiffiffiffiffiffi Sa jt¼0 : C0

ð26Þ

This asymptotic result (26) can be used for construction of efficient relaxation solvers in the case where the temporal time step is much larger than characteristic relaxation time.

2704

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

4. Non-equilibrium model 4.1. Modified governing equations In Section 3, we have explained how to solve the equation for the effective elastic tensor G for ‘‘master’’ system (14). However, system (14) contains another equation which poses serious computational challenges, those are even present in the context of interfaces between two fluids. This is the last algebraic equation (‘‘microstructure’’ equation) expressing the equality of pressures: pg = ps. This equation determines the volume fraction in the mixture zone. Instead of solving the non-linear algebraic equation, we replace it by an ordinary differential equation as in Baer and Nunziato [3]:

Dag ¼ l0 ðpg  ps Þ Dt with a large parameter l0. When l0 goes to infinity, we recover the equilibrium condition pg = ps. We will follow the same strategy presented in Saurel et al. [23] in the context of fluid–fluid interfaces and Favrie et al. [8] for elastic solid– compressible fluid interfaces. A general non-equilibrium model is:

 T DG @v @v 2a þG G¼ GS ; þ Dt @x @x srel s @ðas qs Þ þ divðas qs v Þ ¼ 0; @t @ðag qg Þ þ divðag qg v Þ ¼ 0; @t @ qv þ divðqv  v  rÞ ¼ 0; @t @ðqEÞ þ divðqv E  rv Þ ¼ 0: @t Dag ¼ l0 ðpg  ps Þ: Dt

ð27Þ

Compared to ‘‘master’’ system (14), the only modification in (27) is the last equation for the volume fraction. Let us check that the model is thermodynamically compatible, i.e. the entropy inequality is satisfied. 4.2. Compatibility with the entropy inequality We transform the energy Eq. (8) to the entropy equation:

     @   v  v  v  v De @v  tr ðas rs þ ag rg Þ q eþ þ div qv e þ  ðas rs þ ag rg Þv ¼ q @t Dt 2 2 @x !  e   pg pg Dag p p Das @es DG @es Dgs @eg Dgg þ Ys Y g divv  Y g ¼ q Y s s divv  Y s s þ Y s tr þ Yg qs as qs Dt @G Dt @ gs Dt qg ag qg Dt @ gg Dt   @v  tr ðas rs  ag pg IÞ @x  e   ! Das Dag @es DG @es Dgs @eg Dgg @v ¼ ðas ps þ ag pg Þdivv  ps þ Ys  pg þ q Y s tr þ Yg  tr ðas rs  ag pg IÞ Dt Dt @G Dt @ gs Dt @ gg Dt @x



Das Dag  pg Dt Dt  !    @ees @ v aSs @es Dgs @eg Dgg @v þ Ys tr 2qs G  þ Yg  tr ðas rs  ag pg IÞ @G @x srel @ gs Dt @ gg Dt @x

¼ ðas ps þ ag pg Þdivv  ps þ as

Das Dag  pg Dt Dt  !    @ v aSs @es Dgs @eg Dgg @v þ Ys þ tr as Ss  þ Yg  tr ðas rs  ag pg IÞ @x srel @ gs Dt @ gg Dt @x !   Das Dag aas @es Dgs @eg Dgg ¼ ps tr S2s þ Y s  pg þ  þ Yg Dt Dt srel @ gs Dt @ gg Dt Dgg Dgs Das aas  2  ¼ as qs hs tr Ss : þ ag qg hg  ðps  pg Þ  Dt Dt Dt srel ¼ ðas ps þ ag pg Þdivv  ps

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2705

Finally, we get

as qs hs

Dgg Dgs Das aas  2  tr Ss ¼ 0: þ ag qg hg  ðps  pg Þ  Dt Dt Dt srel

ð28Þ

Eq. (28) should be separated into two equations for the entropy evolution of each component. Such a separation is also phenomenological but it should respect the entropy inequality. We use here the following separation:

as qs hs

Dgs Das aas trðS2s Þ; ¼ ðps  pI Þ þ Dt Dt srel

ag qg hg

Dgg Dag ¼ ðpg  pI Þ : Dt Dt

Obviously, they are compatible with (28): the sum of these equations gives us (28). Here pI is the interfacial pressure. In general, the interfacial pressure is supposed to be a linear combination of average pressures in each phase [13]:

pI ¼ bg ps þ bs pg ;

bs > 0;

bg > 0;

bs þ bg ¼ 1:

Hence, it is always between ps and pg. The coefficients bg and bs should also be determined somehow. For example, in Saurel et al. [24] the pressure pI was estimated as the pressure obtained from the solution of the linearized Riemann problem for which the left and right states have pressures pg and ps and equal velocities (see Eq. (34) in that paper for the entropies in the limit of equal velocities). This limit corresponds to the choice:

bg ¼

Zg ; Zg þ Zs

bs ¼

Zs ; Zg þ Zs

where Zk = qkck, k = s, g are acoustical impedances. The choice of pI does not really matter. Indeed, as shown in Saurel et al. [23] in the case of pure fluids, the choice of this variable has no influence on the numerical results. Finally, the entropy inequality will be satisfied:

@ðqgÞ þ divðqgv Þ P 0; @t where g is the mixture entropy:

g ¼ Y s gs þ Y g gg : For numerical purposes, the equations for the entropies are usually replaced by the evolution equations for the specific energies of each phase:

  @ @v Das ¼ pI ðas qs es Þ þ divðas qs es v Þ  as tr rs ; @t @x Dt   @ @v Das ¼ pI ðag qg eg Þ þ divðag qg eg v Þ  ag tr rg : @t @x Dt

ð29aÞ ð29bÞ

For the equation of state in separate form (10), they can be written as:

   @  Da a as qs ehs þ div as qs ehs v þ as ps divv ¼ pI s  s trðSs RÞ; @t Dt 2srel @ Das ðag qg eg Þ þ divðag qg eg v Þ þ ag pg divv ¼ pI : @t Dt

ð30aÞ ð30bÞ

5. Numerical approximation of the elastic–plastic diffuse interface model In this Section we propose a Godunov type scheme for the resolution of the diffuse interface model for elastic–plastic solids. We will consider the one-dimensional case where all the variables depend only on (x, t). Moreover, the velocity field has only two non trivial components: v = (u, v, 0)T. We suppose that A(3) = B(3) = C(1) = C(2) = 0. Indeed, one can prove that if these variables vanish initially, they will vanish at any time. This model can be rewritten under the following form:

@ as @ as þu ¼ l0 ðps  pg Þ; @t @x @ ag qg eg @ ag qg eg u @u  ag r11g ¼ pI l0 ðpg  ps Þ; þ @x @x @t @ as qs es @ as qs es u @u @v  as r11s  as r12s þ ¼ pI l0 ðps  pg Þ; @x @x @t @x ð1Þ ð1Þ ð1Þ ð1Þ @A @A u @v S11s A þ S12s B þ Bð1Þ ; þ ¼a @x @t @x srel

ð31aÞ ð31bÞ ð31cÞ ð31dÞ

@Að2Þ @Að2Þ u @v S11s Að2Þ þ S12s Bð2Þ þ Bð2Þ ; þ ¼a @x @t @x srel

ð31eÞ

@Bð1Þ @Bð1Þ S12s Að1Þ þ S22s Bð1Þ ; þu ¼a @t @x srel

ð31fÞ

2706

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

@Bð2Þ @Bð2Þ S12s Að2Þ þ S22s Bð2Þ ; þu ¼a @t @x srel

ð31gÞ

@C ð3Þ @C ð3Þ S33s C ð3Þ ; þu ¼a @t @x srel @ as qs @ as qs u ¼ 0; þ @x @t @ ag qg @ ag qg u ¼ 0; þ @x @t @ qu @ qu2  ðag r11g þ as r11s Þ þ ¼ 0; @t @x @ qv @ quv  ðag r12g þ as r12s Þ ¼ 0; þ @x @t @ðqEÞ @ðuðqE  ðag r11g þ as r11s ÞÞ  ðag r12g þ as r12s Þv Þ þ ¼ 0: @t @x

ð31hÞ ð31iÞ ð31jÞ ð31kÞ ð31lÞ ð31mÞ

Here

1 1 E ¼ Y s es þ Y g eg þ u2 þ v 2 ; 2 2

eg ¼ ehg ;

es ¼ ehs þ ees :

The thermodynamic closure is achieved by taking the following equations of state for the gas phase:

qg ehg ¼

pg þ cg p1g ; cg  1

ð32Þ

and for the solid phase:

qs ehs ¼

ps þ cs p1s ; cs  1

qs ees ¼

ls q s G trððg  IÞ2 Þ; g ¼ 1=3 4qs0 jGj

ð33Þ

with G given by:

0

ðAð1Þ Þ2 þ ðAð2Þ Þ2 B ð1Þ ð1Þ G ¼ @ A B þ Að2Þ Bð2Þ

Að1Þ Bð1Þ þ Að2Þ Bð2Þ ðBð1Þ Þ2 þ ðBð2Þ Þ2

0

0

0

1

C 0 A C ð3Þ

and

det G ¼ C ð3Þ ðAð1Þ Bð2Þ  Að2Þ Bð1Þ Þ2 : The solid stress tensor is:

rs ¼ ps I þ Ss ; Ss ¼ 





ls qs 2 trðg2 Þ trðgÞ I g I g  3 3 qs0

 :

The stress tensor for the gas phase is:

rg ¼ pg I: This resolution is done in three main steps:  Resolution of the hyperbolic part.  Relaxation of the stresses (or strains) from the plasticity region to the yield surface.  Relaxation of the pressures to the equilibrium pressure. The first and the third step are similar to the method presented in Favrie et al. [8]. 5.1. Hyperbolic step We present here only the first order Godunov method. A second order extension could be formulated in the same manner as in Favrie et al. [8] using a MUSCL procedure. We want to solve the following part of the model:

@ as @ as þu ¼ 0; @t @x @ ag qg eg @ ag qg eg u @u þ ag pg ¼ 0; þ @x @x @t

ð34aÞ ð34bÞ

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2707

@ as qs es @ as qs es u @u @v  as r11s  as r12s þ ¼ 0; @x @x @t @x

ð34cÞ

@Að1Þ @Að1Þ u @v þ Bð1Þ þ ¼ 0; @x @t @x

ð34dÞ

@Að2Þ @Að2Þ u @v þ Bð2Þ þ ¼ 0; @x @t @x

ð34eÞ

@Bð1Þ @Bð1Þ þu ¼ 0; @t @x

ð34fÞ

@Bð2Þ @Bð2Þ þu ¼ 0; @t @x

ð34gÞ

@C ð3Þ @C ð3Þ þu ¼ 0; @t @x

ð34hÞ

@ as qs @ as qs u ¼ 0; þ @x @t

ð34iÞ

@ ag qg @ ag qg u ¼ 0; þ @x @t

ð34jÞ

@ qu @ qu2  ðag r11g þ as r11s Þ þ ¼ 0; @t @x

ð34kÞ

@ qv @ quv  ðag r12g þ as r12s Þ ¼ 0; þ @x @t

ð34lÞ

@ðqEÞ @ðuðqE  ðag r11g þ as r11s ÞÞ  ðag r12g þ as r12s Þv Þ þ ¼ 0: @t @x

ð34mÞ

5.1.1. Characteristic speeds This system can alternatively be rewritten in primitive variables

qs ; qg ; u; v ; ps ; pg ; Að1Þ ; Að2Þ ; Bð1Þ ; Bð2Þ ; C ð3Þ ÞT

W ¼ ðas ; under the form:

@W @W þ AðWÞ ¼ 0; @t @x

ð35Þ

where the matrix A is given by:

0

u

B B 0 B B B 0 B B B ps pg B q B B B 0 B B B B 0 B B B 0 B B B B 0 B B B 0 B B B 0 B B B 0 @ 0

0 0

0

0

0

0

0

0

0

0

u 0

qs

0

0

0

0

0

0

0

0 u

qg

0

0

0

0

0

0

0

 aqs @@Ar11s ð1Þ

 aqs @@Ar11s ð2Þ

 aqs @@Br11s ð1Þ

 aqs @@Br11s ð2Þ

0 0

u

0

as q

ag q

0 0

0

u

0

0

 aqs @@Ar12s ð1Þ

 aqs @@Ar12s ð2Þ

 aqs @@Br12s ð1Þ

 aqs @@Br12s ð2Þ

0 0

qs c2s

0

u

0

0

0

0

0

0 0 q

2 g cg

0

0

u

0

0

0

0

0 0

A

ð1Þ

ð1Þ

0

0

u

0

0

0

A

ð2Þ

ð2Þ

0

0

0

u

0

0

0 0

B B

0 0

0

0

0

0

0

0

u

0

0 0

0

0

0

0

0

0

0

u

0 0

0

0

0

0

0

0

0

0

0

1

C C C C C 0 C C as @ r11s C  q @Cð3Þ C C C C  aqs @@Cr12s ð3Þ C C C C 0 C: C C 0 C C C C 0 C C C 0 C C C 0 C C C 0 A 0

u

Here the sound velocities cs and cg are determined only by the hydrodynamic parts of energy:

c2s ¼ cs

ps þ p1s

qs

;

c2g ¼ cg

pg þ p1g

qg

:

2708

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

The characteristic wave speeds are solutions of the following polynomial:

  1 0  ð1Þ @r ð2Þ @ r11s ð2Þ @ r12s as A @A11s þ Bð1Þ @@Ar12s  as qs c2s þ ag qg c2g ð1Þ þ A ð1Þ þ B @Að2Þ @Að2Þ A ðu  mÞ þ ðu  mÞ 4

2@

0 þ as @

q

  @ r12s @ r12s ðAð1Þ Bð2Þ  Bð1Þ Að2Þ Þ @@Ar11s  @@Ar11s as ð1Þ ð2Þ @Að2Þ @Að1Þ

q2

  1 ð2Þ @ r12s Bð1Þ @@Ar12s ag qg c2g þ as qs c2s ð1Þ þ B @Að2Þ A ¼ 0: 

q2

In the limit of small deformations (A(1)  1, A(2)  0, B(1)  0, B(2)  1, C(3)  1) we obtain approximate eigenvalues

u;

rffiffiffiffiffiffiffiffiffiffiffi l u  Y s s;

qs

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   4 ls þ Y g c2g : u  Y s c2s þ 3 qs

We use the HLLC solver [26,27]. It preserves positivity of the volume fraction and the density, and is able to deal with strong shocks. It was shown in Favrie et al. [8] that, even if the HLLC Riemann solver has only 3 waves instead of 5, it is able to capture all five waves. With this solver, each wave is considered as a discontinuity and, consequently, jump relations are needed. The system being non-conservative, shock relations are non-conventional. 5.1.2. Shock relations In Favrie et al. [8] approximate shock relations have been derived for a non-equilibrium diffuse elastic solid–fluid interface model. These Rankine–Hugoniot relations verify the following properties:  Single phase limit is recovered for which jump relations are unambiguously known.  Conservation of the mass, momentum and energy of the mixture are satisfied.  Second law of thermodynamics for the mixture is satisfied. Below, we remind these Rankine–Hugoniot relations which fulfill these requirements. The relations are accurate enough to reach convergence for solid–fluid interface problems. 5.1.2.1. Jump relations for conservative variables. The jump relations for the mass equations are:

as qs ðu  DÞ ¼ a0s q0s ðu0  DÞ ¼ ms ; ag qg ðu  DÞ ¼ a0g q0g ðu0  DÞ ¼ mg :

ð36Þ ð37Þ

In particular, the mass fraction is continuous through the jump:

Y k ¼ Y 0k ;

k ¼ s; g:

The variables with superscript ‘‘0’’ denote the unshocked state. The speed of the discontinuity is denoted by D. Let us denote by

r11 ¼ as r11s þ ag r11g ; r12 ¼ as r12s þ ag r12g the mixture normal and tangential stresses, and

m ¼ ms þ mg the mixture mass flux. With these notations, the momentum jump relations read:

1

r11  r011  m2 ðs  s0 Þ ¼ 0; s ¼ ; q r12  r012  mðv  v 0 Þ ¼ 0:

ð38Þ ð39Þ

The mixture energy relation is:

e  e0 ¼

  1   1 ðs  s0 Þ r11 þ r011 þ ðs0 ðAð2Þ Þ0  sAð2Þ Þ r12 þ r012 2 2

with e = Yses + Ygeg and s = Ygsg + Ysss, sk = 1/qk.

ð40Þ

2709

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

5.1.2.2. Jump relations for geometrical variables. In absence of relaxation effects, the geometric equations read:

@Að1Þ @Að1Þ u @v þ Bð1Þ þ ¼ 0; @x @t @x

ð41Þ

@Að2Þ @Að2Þ u @v þ Bð2Þ þ ¼ 0; @x @t @x

ð42Þ

@Bð1Þ @Bð1Þ þu ¼ 0; @t @x

ð43Þ

@Bð2Þ @Bð2Þ þu ¼ 0; @t @x

ð44Þ

@C ð3Þ @C ð3Þ þu ¼ 0: @t @x

ð45Þ

It follows from (43)–(45) that we have the following invariants across left- and right-facing waves:

Bð1Þ ¼ ðBð1Þ Þ0 ; Bð2Þ ¼ ðBð2Þ Þ0 ; C ð3Þ ¼ ðC ð3Þ Þ0 : Hence the jump relations for (41) and (42) are:

Að1Þ ðu  DÞ þ Bð1Þ v ¼ ðAð1Þ ðu  DÞ þ Bð1Þ v Þ0 ; Að2Þ ðu  DÞ þ Bð2Þ v ¼ ðAð2Þ ðu  DÞ þ Bð2Þ v Þ0 : 5.1.2.3. Jump relations for non-conservative variables. The volume fraction is continuous through the jump:

ak ¼ a0k ; k ¼ g; s: The internal energy equations being in non-conservative form, they are not adapted to the determination of jump relations. A path has to be defined to circumvent this difficulty, since the results will depend on the path these shock relations will be used only to predict the state then a correction will be used to correct the energy of each part. The convergence of this method has been shown in Saurel et al. [23]. The following jump relations can be proposed:

  1   1 ðs  s0 Þ r11 þ r011 þ ðs0 ðAð2Þ Þ0  sAð2Þ Þ r12 þ r012 ; 2 2    1  0 ð2Þ 0  1 ek  e0k ¼ sk  s0k r11k þ r011k þ s ðA Þ  sk ðAð2Þ Þ r12k þ r012k ; 2 2 k

e  e0 ¼

ð46Þ k ¼ g; s:

ð47Þ

This choice generalizes classical additivity principle for heterogeneous mixtures (see, for example, [21]) and is compatible with the following requirements  Single phase limit When the solid volume fraction vanishes, the preceding jump relations correspond to the jump relations for the Euler equations. This condition is fulfilled because r12g is zero is the gas and r11g = pg. When the gas volume fraction vanishes, the mixture jump condition (46) reduces to the conditions in solids. Let us show that (47) implies the energy conservation law (46).  Conservation of the total energy Multiplying each Eq. (47) by Yk = akqk/q and summing them we obtain:

  r11g þ r0     r11s þ r011s   11g Y s es  e0s  Y s ss  Y s s0s þ Y g eg  e0g  Y g sg  Y g s0g 2 2   0   0 r þ r 12g 12g ðr12s þ r12s Þ  ððAð2Þ Þ0 Y s s0s  Að2Þ Y s ss Þ  ðAð2Þ Þ0 Y g s0g  Að2Þ Y g sg ¼ 0: 2 2 Or 0

ee 

as r11s þ ag r11g þ as r011s þ ag r011g 2

0

ð2Þ 0 0

ðs  s Þ  ððA Þ

ð2Þ

s  A sÞ

as r12s þ ag r12g þ as r012s þ ag r012g 2

! ¼ 0:

This is equivalent to (46). With the help of these relations, approximate Riemann solvers can be built as developed hereafter.

2710

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

5.1.3. Interface relations Consider system (31) in absence of relaxation terms. Eqs. (41) and (42) show that if B(1) and B(2) are discontinuous across contact discontinuities, then necessarily the tangential velocity should be continuous, [v] = 0. If not, the non-conservative product is not well defined. Thus, the resulting set of interface relations is:

½u ¼ 0;

ð48aÞ

½v  ¼ 0 if

ð48bÞ

as – 0; ½r11  ¼ ½ag r11g þ as r11s  ¼ 0; ½r12  ¼ ½as r12s þ ag r12g  ¼ 0:

ð48cÞ ð48dÞ

A difficulty appears with system (48). When an interface appears between two solid states, the jump condition [v] = 0 holds. However, if an interface separates two fluid states or solid and fluid, the sliding is possible: [v] – 0. Since the Lagrangian coordinates are discontinuous when the sliding is present, the cobasis equations have no sense. Hence, it will necessitate a specific numerical treatment. This issue will be addressed later. 5.1.4. HLLC type Riemann solver Consider a cell boundary separating a left state (L) and a right state (R). The left - and right facing wave speeds are obtained following Davis [5] estimates:

SR ¼ maxðuL þ cL ; uR þ cR Þ; SL ¼ minðuL  cL ; uR  cR Þ; where cL,R are the longitudinal sound speed of the mixture: c2 ¼ Y s c2ls þ Y g c2g with cls the longitudinal wave speed in the pure solid. As it was shown in Gavrilyuk et al. [12] and Favrie et al. [8], the Godunov method with HLLC solver (formulated in terms only three waves!) was able to capture all the fives waves. The speed of the intermediate wave (or contact discontinuity) is estimated under HLL approximation:

SM ¼

ðqu2  r11 ÞL  ðqu2  r11 ÞR  SL ðquÞL þ SR ðquÞR ¼ u ; ðquÞL  ðquÞR  SL qL þ SR qR

with the mixture normal stress r11 = asr11(s) + agr11(g) and mixture density q = asqs + agqg defined previously. From these wave speeds conservative state variables in the star region (see Fig. 1) are determined:

ðak qk Þ L;R ¼ ðak qk ÞL;R

SL;R  uL;R ; SL;R  u

r 11 ¼

ðuR  SR ÞqR r11L  ðuL  SL ÞqL r11R þ ðuL  SL ÞqL ðuR  SR ÞqR ðuR  uL Þ ; ðuR  SR ÞqR  ðuL  SL ÞqL

r 12 ¼

ðuR  SR ÞqR r12L  ðuL  SL ÞqL r12R þ ðuL  SL ÞqL ðuR  SR ÞqR ðv R  v L Þ ; ðuR  SR ÞqR  ðuL  SL ÞqL

v L;R ¼ v L;R þ

ðr12 Þ  ðr12 ÞL;R ðas r12s Þ  ðas r12s ÞL;R ¼ v L;R þ : ðuL;R  SL;R ÞqL;R ðuL;R  SL;R ÞqL;R

The components v L and v R are equal for the case of a solid–solid contact discontinuity. But we prefer to make a distinction between them, as they will be different in the case of a solid–fluid contact. The right and the left mixture total energies are obtained as:

E L;R ¼

qL;R EL;R ðuL;R  SL;R Þ  r11L;R uL;R  r12L;R v L;R þ r 11 u þ r 12 v L;R q L;R ðu  SL;R Þ

with

1 1 E ¼ Y s es þ Y g eg þ u2 þ v 2 : 2 2 Then, the geometrical variables are determined as: ð1Þ

ðAð1Þ Þ L;R ¼

ð1Þ

AL;R ðuL;R  SL;R Þ þ BL;R

u  SL;R



v L;R  v L;R

 ;

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2711

Fig. 1. HLLC approximate solver. Solution in the ‘‘star’’ region consists of two constant states separated from each other by a middle wave of speed SM.

ð2Þ

ðAð2Þ Þ L;R

¼

ð2Þ

AL;R ðuL;R  SL;R Þ þ BL;R



v L;R  v L;R

u  SL;R ð1Þ

ðBð1Þ Þ L;R ¼ BL;R ;

ð2Þ

ðBð2Þ Þ L;R ¼ BL;R ;

 ; ð3Þ

ðC ð3Þ Þ L;R ¼ C L;R

The non-conservative variables are finally determined:

a kL;R ¼ akL;R ; a kL;R ¼ akL;R ; q kL;R ¼ q0kL;R

uL;R  SL;R ; SM  SL;R

k ¼ s; g:

The jumps of the internal energies are determined with the help of the Hugoniots (47). Indeed, for example, with equations of state (32) and (33), the phase pressures p kL;R are calculated as functions of the phases densities q kL;R and (for k = s) geometrical variables ðAðbÞ Þ L;R , ðBðbÞ Þ L;R along their respective Hugoniot curves defined by (47). Then the internal energies are determined as:

  e gL;R ¼ eg p gL;R ; q gL;R ;   e sL;R ¼ es p sL;R ; q sL;R ; ðAð1Þ Þ L;R ; ðAð2Þ Þ L;R ; ðBð1Þ Þ L;R ; ðBð2Þ Þ L;R ; ðC ð3Þ Þ L;R : With the help of this approximate solver, it is possible to derive a Godunov type scheme. 5.1.5. Godunov type scheme The method used to solve interface problems with the diffuse interface formulation (31) consists in a sequence of steps. The system contains conservative and non-conservative equations, as well as relaxation terms. The algorithm we use follows the one presented in Saurel et al. [23] in a simplified situation involving two fluids. The conservative part of system (31), in absence of relaxation terms is updated with the conventional Godunov scheme:

Unþ1 ¼ Uni  i

  Dt   n n  F1 Ui ; Uiþ1  F 1 Uni1 ; Uni ; Dx

ð49Þ

with

U ¼ ððaqÞs ;

qu; qv ; qEÞT

ðaqÞg ;

and

F1 ¼ ððaqÞs u;

ðaqÞg u;

qu2  r11 ; quv  r12 ; ðqE  r11 Þu  r12 v ÞT

The non-conservative part of system (31), in absence of relaxation terms reads:

@G @F2 @u @v þ þ K2 þ H2 ¼ 0; @t @x @x @x

ð50Þ

with

G ¼ ðas ; F2 ¼ ðas u;

Að1Þ ;

Bð1Þ ;

Að1Þ u;

Að2Þ ;

Bð1Þ u;

Bð2Þ ; Að2Þ u;

C ð3Þ ; Bð2Þ u;

as qs es ; ag qg eg ÞT ; C ð3Þ u; as qs es u;

ag qg eg uÞT ;

2712

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

H2 ¼ ðas ;

K2 ¼ ð0;

0;

Bð1Þ ;

Bð1Þ ;

0;

0;

Bð2Þ ;

0;

Bð2Þ ;

C ð3Þ ;

0;

as r12s ;

as r11s ;

ag r11g ÞT ;

ag r12g ÞT :

The non-conservative equations are solved by the following scheme, already used in Saurel et al. [23] to estimate volume fractions and internal energy equations:

Gnþ1 ¼ Gni þ i

     Dt  F2i1=2  F 2iþ1=2 þ Hn2i u i1=2  u iþ1=2 þ Kn2i v i1=2R  v iþ1=2L ; Dx

where F 2i1=2 and F 2iþ1=2 are the fluxes calculated on each cell boundary with the help of the HLLC solver. The normal velocity components u i1=2 and u iþ1=2 are also solution of the Riemann problem at the same cell boundaries, v i1=2R and v iþ1=2L are the corresponding tangential velocity components at the right and left sides of the contact discontinuity. Eqs. (49) and (50) can be rewritten under the following form:

@V @F @u @v þ þH þK ¼0 @t @x @x @x

    with VT = (UT, GT) and FT ¼ FT1 ; FT2 ; HT ¼ 0; 0; 0; 0; 0; HT2 ; KT ¼ ð0; 0; 0; 0; 0; KT2 Þ. The following scheme is used for time evolution for all the equations:

Vnþ1 ¼ Vni þ i

     Dt  Fi1=2  F iþ1=2 þ Hni u i1=2  u iþ1=2 þ Kni v i1=2R  v iþ1=2L : Dx

Regarding the non-conservative internal energy equations it is worth mentioning that the adopted jump relations used in the Riemann solver will induce entropy jump in each phase in the presence of shocks. However, there is no hope that this jump be exact as the two-phase shock is smeared over several computational cells. Moreover, the evolution scheme (50) uses cell averages of non-conservative variables and simplified approximations for non-conservative products with Hni and Kni as local constants. In order to correct this predicted solution, two extra steps are used. The first one corresponds to a pressure relaxation step, in order that mechanical equilibrium be reached everywhere and in particular at the interface. The second correction step corresponds to a reset of the internal energies with the help of the total energy conservation equation. These two correction steps are detailed in Section 5.4. 5.2. Sliding treatment A small amount of solid is always present in the air when the diffuse interface model is used. This small amount of solid will prevent sliding at the interfaces. Consequently, the tangential stress r12 is non-zero, and this creates transverse waves in the gas phase, resulting in wrong wave structure in the vicinity of the interface. This difficulty has been already mentioned when the interface solid–gas relations (48) were discussed. To circumvent this problem, a special treatment of the sliding effect has to be done. First, we have to know what type of contact (solid–solid, fluid–fluid or solid–fluid) we are considering at a given interface. The procedure is similar to the one presented in Favrie et al. [8] but the level set function is no longer useful to detect the interfaces. Indeed, the volume fraction or the mass fraction are already the ‘‘order’’ parameters to be used as a level set function. Both of them can be considered. We will use the volume fraction variable because this variable is capable to present spallation effects in solids. Indeed, the volume fraction of solid will decrease in tensile stress waves while the mass fraction will not change in such waves. Nevertheless, we will compare both choices in a specific test case. Using the volume fraction of the solid as as a ‘‘level set variable’’, and denoting alim the value corresponding to the interface, the Riemann problem solution is ‘‘modified’’ with the following correcting procedure:

r 12

8 ðuR SR Þq ðas r12s Þ ðuL SL Þq ðas r12s Þ þðuL SL Þq ðuR SR Þq ðv R v L Þ R L L R L R ; > ðuR SR ÞqR ðuSL ÞqL < ¼ if as;L > alim and as;R > alim ; > : 0; if as;L < alim or as;R < alim :

The other variables in the Riemann problem solution are unchanged. With this simple correction, the solid–solid and fluid–fluid single phase limits are preserved. However, when dealing with solid–fluid interfaces, this correction is not enough as there is a non-linear coupling with the numerical diffusion of the volume fraction. In particular, the oscillations of the tangential stress are related to the numerical diffusion of the transverse velocity. They create non-physical stress waves in the solid. To circumvent this supplementary difficulty and enforce solid–fluid interface conditions, some choice has to be made regarding the correct set of jump relations for the transverse velocity given by (48). This is done with an adaptation of the method due to Abgrall and Karni [2]. This method is only used to compute the transverse speed and not the interface flow field as in their original paper. The ‘‘Ghost fluid method’’ used here can be summarized as follows:

2713

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

1. Capture the material interface by checking the product of the values of the variable U = as  alim in the ith and i + 1th cells. A positive value indicates that the edge i + 1/2 does not contain the material interface. In this case the numerical flux and the primitive variables vector are computed without using the ‘Ghost fluid method’’. A negative value indicates that the material interface is present between the two cell centers i and i + 1. In this case two different fluxes and primitive variable vectors are computed after velocity extrapolation across the interface. This procedure can be written as follows:     If Uni Uniþ1 > 0 then W iþ1=2 ¼ WHLLC Vni ; Vniþ1 ; F iþ1=2 ¼ F W iþ1=2 . If Uni Uniþ1 < 0 then   ⁄ Extrapolate the right tangential velocity to the left and solve the Riemann problem: W iþ1=2;L ¼ WHLLC VnL;i ; Vniþ1 ,   F iþ1=2;L ¼ F W iþ1=2;L . Variables VnL;i has been changed with the tangential velocity component of the right state Vniþ1 . ⁄ Extrapolate the left tangential velocity to the right and solve the Riemann problem:     W iþ1=2;R ¼ WHLLC Vni ; VnR;iþ1 ; F iþ1=2;R ¼ F W iþ1=2;R . Variables VnR;iþ1 has been changed with the tangential velocity component of the left state Vni . 2. Following Farhat et al. [7], the prediction of the solution is necessary to detect interfaces that may change cells during two e nþ1 of the solution is determined with: consecutive time steps. A temporary value V i

~ nþ1 ¼ Vn þ Dt V i i Dx

      F i1=2;R  F iþ1=2;L þ Hni u i1=2;R  u iþ1=2;L þ Kni v i1=2;R  v iþ1=2;L :

3. For each node i, e nþ1 . If Uni Unþ1 > 0 then Vnþ1 ¼V i i i nþ1 n nþ1 e nþ1 except the transverse speed v ~ nþ1 If Ui Ui < 0 then Vi ¼ V that is taken to i i direction. 5.3. Plastic relaxation

v Ri

or

v Li .

depending on the flow

For each cell, we solve the following system:

@ as ¼ 0; @t @ ag qg eg ¼ 0; @t @ a s qs ¼ 0; @t @ qu ¼ 0; @t

ð51aÞ @ as qs es ¼ 0; @t @ ag qg ¼ 0; @t

@ qv ¼ 0; @t

@ðqEÞ ¼ 0: @t

ð51bÞ ð51cÞ ð51dÞ

The only variables that are changing are the cobasis vectors:

D b a E ¼ SEb : Dt srel This equation is a system of ordinary differential equation (ODE) and could be solved using any ODE integrator. Since this is a system of nine equations, such an approach will be ‘‘expensive’’. A possibility to reduce the ‘‘cost’’ of the integration, is detailed hereafter. 5.3.1. Singular value decomposition At the hyperbolic step we have calculated the cobasis Ea, a = 1, 2, 3. Consider the deformation gradient (FT)1 = (E1, E2, E3). Using a classical ‘‘singular value decomposition’’, we can rewrite the matrix (FT)1 under the following form: (FT)1 = VKUT where U and V are two orthogonal matrices and K is a diagonal matrix. The eigenvalues of K are singular values of (FT)1. This procedure is contained in the LAPACK libraries. The eigenvalues of the deformation tensor calculated during the elastic step are initial conditions for the relaxation Eq. (20). We also suppose that U and V are time independent during the relaxation procedure. This reduces the system of nine equations to only three equations. 5.3.2. Calculating the relaxed stress If the effective deformations are small (what is the case for most metals), we relax the cobasis equations using approximate relation (25). If the relaxation time function is more complicated, the corresponding relaxation equation should be solved directly using any ODE resolution method. If the characteristic time s0is very small compared to the time step, a good approximation could be to consider the limit solution of the ODE as t ? 1. 5.3.3. Determining new cobasis values For given values of the final stress Saf, we determine baf, a = 1, 2, 3 from the following system of algebraic equations:

2714

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

Saf ¼ ls

b21f  b1f þ b22f  b2f þ b23f  b3f qs 2 baf  baf  qs0 3

! ð52Þ

and

b1f b2f b3f ¼ 1:

ð53Þ

This system can be efficiently solved by a classical Newton–Raphson procedure. The three Eqs. (52) are not independent, the solution will be given using the first two equations of (52) and Eq. (53). Then we calculate final singular values

kaf ¼ signðkai Þ

qs qffiffiffiffiffiffiffi b : qs0 af

When kaf are computed, the final curvilinear cobasis (Ea)f can be updated using the singular value decomposition ðFT Þ1 ¼ VKf UT , where f

0

k1f B Kf ¼ @ 0 0

0 k2f 0

1 0 C 0 A: k3f

We underline that the matrices U and V do not change during the relaxation step, so they are known. 5.3.4. Summary of the plastic relaxation step We summarize here the resolution of the plastic relaxation step.  Apply the singular value decomposition to the matrix FT = VKUT. The matrices V and U will stay constant during the plastic relaxation step.  Determine the final solution after a time step of the equation @k@ta ¼ Ssa ka . Two ways are possible: rel – Use any classical numerical scheme to solve the ODE equation. – Use a simplified procedure restricted to small effective deformations and the von Mises criterion.  Compute a new cobasis matrix FT using the orthogonal matrices computed initially: FT ¼ VKf UT . f 5.4. Pressure relaxation step At the end of the hyperbolic – plastic step the two materials have different pressures. In order that mechanical equilibrium be restored in mixture cells, a new relaxation step has to be done for each cell:

@ as ¼ l0 ðps  pg Þ; @t @ ag qg eg ¼ pI l0 ðpg  ps Þ; @t @ as qs es ¼ pI l0 ðpg  ps Þ; @t

ð54aÞ ð54bÞ ð54cÞ

@AðbÞ @BðbÞ @C ð3Þ ¼ 0; ¼ 0; ¼ 0; b ¼ 1; 2; @t @t @t @ ak qk ¼ 0; k ¼ s; g; @t @ qu @ qv @ qE ¼ 0; ¼ 0: ¼ 0; @t @t @t

ð54dÞ ð54eÞ ð54fÞ

Following Saurel et al. [23] we will propose a simplified approach in the limit of stiff relaxation (l0 ? +1). The phase energy equation can be rewritten in the following form:

@ek @ sk þ pI ¼ 0; @t @t

sk ¼

1

qk

:

Integrating this equation over the time interval, we can rewrite it in the algebraic form:

^Ik efk  eik  p





sfk  sik ¼ 0;

ð55Þ

where indices f and i mean ‘‘final’’ (equilibrium) and ‘‘initial’’ (nonequilibrium) state, obtained after the hyperbolic-plastic ^Ik is a new average interface pressure. Eq. (55) should imply the conservation of the mixture energy. As a constep. Here p sequence of that, we have:

^Ig ¼ p ^I : ^Is ¼ p p

2715

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

^I is different from pI. However, if the relaxation parameter l0 is large, numerical results are insensitive to a parIn general, p ^I . It is convenient to use p ^I ¼ pf . Thus, using the algebraic energy Eq. (55) for the stiffened gas EOS, the ticular choice of p phase specific volumes become functions only of pressure:

sfk ¼ sik

pik þ ck p1;k þ pf ðck  1Þ pi þ ck p1;k þ pf ðck  1Þ i k s ¼ : k ck ðpf þ p1;k Þ pf þ ck p1;k þ pf ðck  1Þ

The closure of the model is achieved using the saturation constraint:

X

ak ¼ 1

k

Fig. 2. The initial configuration for Wilkins’s flying plate problem.

Pressure (MPa)

density (kg/m3)

7000

2950

6000

2900

5000

2850

4000

2800

3000

2750

2000

2700

1000

2650

0 0.006

0.007

0.008

0.009

0.01 0.011 abscissa (m)

0.012

0.013

0.014

0.015

2600 0.006

0.007

0.008

0.009

velocity (m/s) 1000

700

0

600

-1000

500

-2000

400

-3000

300

-4000

200

-5000

100

-6000

0.007

0.008

0.009

0.01 0.011 abscissa (m)

0.012

0.013

0.014

0.015

0.012

0.013

0.014

0.015

Stress (MPa)

800

0 0.006

0.01 0.011 abscissa (m)

0.012

0.013

0.014

0.015

-7000 0.006

0.007

0.008

0.009

0.01 0.011 abscissa (m)

Fig. 3. The results for Wilkins’s problem with the impact velocity 800 m/s are shown at time t = 0.5 ls. A double structure of the right-and-left facing waves is clearly visible: an elastic precursor of small amplitude is followed by a plastic shock wave of a large amplitude.

2716

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

which can also be rewritten as

X ai k k

ck

pik  pf pf þ p1;k

! ¼ 0:

ð56Þ

There is a unique solution of this equation which is always between pis and pig . Thus, in the limit l0 ? 1 one can use the algebraic Eq. (56) with pf as a single unknown. Once the relaxed pressure is determined, the phase volume fractions are deduced. 5.4.1. Reinitialization step As the volume fractions have been estimated previously by the relaxation method, the mixture pressure can be determined from the mixture EOS based on the mixture energy. Because the mixture total energy obeys a conservation law, its evolution is accurate in the entire flow and, in particular, at shocks. From the mixture total energy, we can extract the hydrodynamic part of the mixture internal energy:

1 1 eh ¼ Y s ehs ðps ; qs Þ þ Y g ehg ðpg ; qg Þ ¼ E  u2  v 2  Y s ees ðAð1Þ ; Að2Þ ; Bð1Þ ; Bð2Þ ; C ð3Þ Þ: 2 2 Again, considering fluids governed by the stiffened gas EOS and using the pressure equilibrium condition ps = pg = p, the mixture EOS now relates the mixture hydrodynamic internal energy, the mixture density and the volume fractions:

pðq; eh ; as ; ag Þ ¼

qeh 



ag cg p1g as cs p1s cs 1 þ cg 1 ag as cs 1 þ cg 1

 ð57Þ

:

Pressure (MPa)

density (kg/m3)

7000

2950

6000

2900

5000

2850

4000

2800

3000

2750

2000

2700

1000

2650

0

2600 0

0.01

0.02

0.03 abscissa (m)

0.04

0.05

0

0.01

0.02

velocity (m/s)

0.03 abscissa (m)

0.04

0.05

0.04

0.05

Stress (MPa)

800

1000

700

0

600

-1000

500 -2000 400 -3000 300 -4000 200 -5000

100

-6000

0 -100

0

0.01

0.02

0.03 abscissa (m)

0.04

0.05

-7000

0

0.01

0.02

0.03 abscissa (m)

Fig. 4. The evolution of flow variables for Wilkins’s flying plate problem with the impact velocity 800 m/s is shown at time instants t = 1 ls, t = 3 ls and t = 5 ls. A right facing elastic shock is followed by a right facing plastic shock. These shocks are followed by a right facing elastic rarefaction wave which starts to overtake the plastic shock. This elastic rarefaction wave is followed by a plastic rarefaction wave. Strong peaks of velocity correspond to the velocity in the air for the three instants. The wave speed in air is smaller compared to the wave speed in the solid, the wave in the air almost does not propagate during the computation time. This peak was naturally absent in Miller and Colella [19], where the same test problem was considered, since air was replaced by vacuum.

2717

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

This EOS is valid in pure phases and in the diffuse interface zone. This guarantees a correct wave dynamics on both sides of the interface. Inside the numerical diffusion zone, numerical experiments show that the method is accurate too, as the volume fractions used in the mixture EOS already came from the relaxation method. Once the mixture pressure is determined the internal energies of the phases are reinitialized with the help of their respective EOS:

ehk ¼ ehk ðp;

ðaqÞk

ak

Þ;

k ¼ s; g

ees ¼ ees ðAð1Þ ; Að2Þ ; Bð1Þ ; Bð2Þ ; C ð3Þ Þ:

t

Solid−air interface

Plastic rarefaction wave Plastic shock and elastic rarefaction waves interaction

Elastic rarefaction wave

Elastic shock x Contact between plates Fig. 5. Schematic representation of the waves present during the Wilkins’s test problem.

Fig. 6. Initial configuration for the copper impact test problem. The initial velocity of the projectile is 800 m/s, the numerical computation involves 700 700 cells.

2718

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

Fig. 7. Impact at 800 m/s of a copper projectile on a copper plate. The average density evolution is shown at time t = 83 ls, t = 280 ls, t = 495 ls and t = 1.1 ms. On the left, the results are obtained by using the volume fraction as the level set function, and on the right the classical level set function has been used.

Fig. 8. Schematic illustration of the initial and final states of Taylor impact specimen [25]. The initial state is shown on the left the final state on the right.

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2719

5.5. Summary of the numerical method The numerical method can be summarized as follows:  At each cell boundary solve the Riemann problem for system (31) without relaxation terms. The HLLC solver is recommended due to its simplicity, robustness and accuracy.  Evolve all flow variables with a higher order Godunov type method.  Use the ghost fluid method to allow sliding in the model.  Relax the deviatoric part of the stress tensor to the yield surface.  Determine the relaxed pressure and the volume fractions by solving Eq. (56). The Newton method is appropriate for this task.  Compute the mixture pressure using the Eq. (57).  Reset the internal energies with the computed pressure with the help of their respective EOS.  Start again for the next time step. This method works very well as shown in the next section. 6. Numerical results 6.1. 1D case: Wilkins’s problem We consider a 5 mm thick aluminum plate impacting a half space of the same material (see Fig. 2). The initial velocity of the plate is 800 m/s. The interface of impact is initially located at xc = 0.01 m. The calculation domain is 0.05 m long (2000 cells, CFL = 0.8). At the left of the plate the air is at rest. All the materials are at the atmospheric pressure. The air is considered as a perfect gaz with cg = 1.4 and qg = 1.2 kg/m3. The parameters of the aluminum are: cs = 3.4, P1s = 21.5 GPa, ls = 248 GPa and Ys = 297.6 MPa. In Fig. 3 we see left-and right facing shock waves propagating outward from the contact of the plate with the half space. An elastic precursor and a plastic shock are clearly present at the left and at the right of the contact discontinuity. When the left facing waves reach the free surface ‘‘air–solid’’, the right facing rarefaction waves are created. The first rarefaction wave is elastic. The second is a plastic rarefaction wave. These rarefaction waves begin to overtake the initial right facing plastic shock. The results are shown in Fig. 4 at time t = 1 ls, t = 3 ls and t = 5 ls. The drop of the density corresponds to the interface with air. All these waves are shown in the (x, t)-diagram (see Fig. 5). Other one-dimensional test cases can be found in Favrie and Gavrilyuk [9,10].

Fig. 9. Initial condition for the Taylor impact test problem. The initial velocity of the impactor is U = 227 m/s, the numerical computation involves 300 700 cells.

2720

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

Fig. 10. Axisymmetrical deformation of a copper rod impacting a rigid surface with a velocity of 227 m/s is shown. The pressure contours and the interface shapes are shown for the following time instants after impact: t = 2.5 ls, t = 20 ls, t = 40 ls, t = 80 ls and t = 100 ls. At time t = 100 ls the final state is reached (the velocity of a solid is less then 5 m/s). The interface between solid and fluid corresponds to the average density level q = 2000 kg/m3. The x-axes is scaled in mm while the y-axis is in m. At the vicinity of the triple line (rigid surface–solid–air) the interface is too diffused to present correctly the rod shape.

6.2. Impact of a copper projectile on a copper plate We consider here an impact of a solid copper projectile on a solid copper plate. The copper projectile has an initial velocity u = 800 m/s while the plate and the surrounding air are at rest. This projectile is a square of 0.1 m length. The plate is 0.5 m

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2721

long and 0.1 m wide. Initially, a narrow space between the copper projectile and the plate is present. The other part of the computational domain contains air at atmospheric conditions. The physical domain is 0.7 m 0.7 m. The mesh contains 700 cells both in the horizontal and the vertical direction. Fig. 6 shows the initial configuration. The parameters of the copper are as follows: cs = 4.22, p1s = 34.2 GPa, the density qs = 8900 kg/m3. The shear modulus is taken as ls = 43.3 GPa, and the yield strength as Ys = 600 MPa. We took here a different value of Ys compared with the previous example. Indeed, in the present paper ‘‘copper’ is just a generic word used also for different alloys containing the majority of the copper, but having different values of yield strength. The air is considered as ideal gas with parameters cs = 1.4 and the density qg = 1 kg/m3. Initially, a small amount of air (as = 104) is present in the copper and a small amount of copper (as = 104) is present in the air. All the results presented hereafter are obtained with a high-order Godunov method using Minmod limiter. In Fig. 7, the average density evolution is shown at time t = 83 ls, t = 280 ls, t = 495 ls and t = 1.1 ms. On the left, the results are presented using the volume fraction as the level set function in the Riemann solver. On the right, the classical level set function is used. At the first time instant the results are analogous. A difference appears when cracks start to appear. However, topologically, the scenario is the same: the initial configuration (copper projectile + copper plate) is disintegrated into five pieces. The two principal pieces (higher and lower part) have exactly the same shape. It is important to note that the crack formation is dynamic: cracks were not present initially. They appear where high rarefaction occurs. 6.3. Taylor impact test problems Taylor [25] first proposed to use the normal impact of a rod against a rigid surface as a method of estimating the yield strength of materials from measurements of the initial and final lengths of the rod along with the length of the deformed section. This is a classical test for numerical modeling elastic–plastic solids [31,28,4]. The initial and final configurations are shown in Fig. 8. The Taylor theory give a simple relation between the final shape of the material (L1 is the total final length, and X is a non-deformable part), the velocity of the impact U, the yield limit Y and the density of the material q:

Y

qU 2

¼

LX 1 : 2ðL  L1 Þ lnðL=XÞ

This formula is considered to be accurate to approximately 10% [11]. We consider a copper cylinder of 32.4 mm high and 6.4 mm diameter impacting a rigid wall. The cylinder has an initial velocity of 227 m/s while the surrounding air is at rest. The contact between the cylinder and the rigid wall is already done at initial time. The other part of the domain is filled with air at atmospheric conditions. The computational domain is 35 mm height and 15 mm large. The mesh contains 700 cells in the vertical direction and 300 cells in the horizontal direction. For such an axisymmetric configuration, the equations are rewritten in cylindrical coordinates (see Miller and Colella [19] as an example of such a writing) and solved by a standard way. Fig. 9 shows the initial configuration. The parameters of the copper are as follow cs = 4.22, p1s = 34.2 GPa, the density qs = 8900 kg/m3. The shear modulus is taken as ls = 43.3 GPa and the yield strength: Ys = 490 MPa. The air is considered as ideal gas with parameters cs = 1.4 and the density qg = 1 kg/m3. To avoid the degeneracy of the model, initially a small amount of air (as = 104) is present in the copper, and a small amount of copper (as = 104) is present in the air. This test has been used in Udaykumar et al. [28], and Camacho and Ortiz [4] but they consider also hardening effects which are absent in our model. Results are presented in Fig. 10 at time instants t = 2.5 ls, t = 20 ls, t = 40 ls, t = 80 ls and t = 100 ls. At time t = 100 ls the final state is reached (the velocity of a solid is less then 5 m/s). At the final state, the total length of cylinder L1

Fig. 11. Initial configuration for the perforation of a disc.

2722

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

is 22 mm and the unplastified part X is 13 mm. From these results and using the Taylor relation we get a yield limit Y = 468 MPa. This value is 5% far from the input value (490 MPa). Results are in good agreement with the Taylor theory. 6.4. Perforation of a disc We consider here an impact of a solid copper cylinder on a solid copper disk. The projectile has an initial velocity u = 800 m/s while the disc and the surrounding air are at rest. This projectile is a cylinder of 0.1 m length and 0.1 m diameter. The disc is 1.2 m diameter and its thickness is 0.1 m. Initially, a narrow space (1 mm) is present between the projectile and the plate. The other part of the computational domain contained air at atmospheric conditions. The physical domain is 0.7 m long and 1.4 m diameter. The mesh contains 700 cells both on the longitudinal and the radial direction. The initial configuration is shown in Fig. 11. The computation is also done for the axisymmetric configuration. The parameters of the copper are as follows: cs = 4.22, p1s = 34.2 GPa, the density qs = 8900 kg/m3. The shear modulus is taken as ls = 43.3 GPa and the yield limit: Ys = 600 MPa. The air is considered as ideal gas with parameters cs = 1.4 and the density qg = 1 kg/m3. A small amount of air (as = 104) is present in the copper, and a small amount of copper (as = 104) is present in the air. All the results presented hereafter

Fig. 12. Perforation of a copper disc by a projectile of the same material. Results are presented at time instants t = 0 ms, t = 1.5 s and t = 3 ms.

N. Favrie, S.L. Gavrilyuk / Journal of Computational Physics 231 (2012) 2695–2723

2723

are obtained with a high-order Godunov method using Minmod limiter [29]. The results are shown in Fig. 12 at time instants t = 0 ms, t = 1.5 ms and t = 3 ms. The disc is perforated at the final time instant. For a better presentation of penetration phenomena, only a part of the disc has been shown. One can also see the Rayleigh waves propagating along the disc surface. 7. Conclusion An Eulerian hyperbolic diffuse interface model coupling Maxwell type compressible materials and compressible fluids has been constructed and numerically solved for different 2D cases: a disintegration of an impacted plate, the Taylor test problem and a perforation of a disc. Numerical examples show the ability of this model to deal with real physical phenomena. Numerical diffusion of solid–fluid interfaces is negligible in the case of high-speed phenomena but becomes important when the interface velocities are low (see the final stage of the Taylor test problem). The agreement between the Taylor estimation and the numerical result is quite good. Still, some work has to be done to use this model for industrial applications. In particular, the ability of the model to deal with dynamic appearance of interface (cracks) has to be studied in detail. The numerical diffusion of the interface should be reduced to be able to compute large time scale. The hardening behavior compatible with the entropy inequality should also be included in the model. Acknowledgments The authors thank Richard Saurel for fruitful discussions, and anonymous referees for important remarks and suggestions. The authors have been partially supported by grant N 07.34.048 from the DGA. References [1] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: a quasi-conservative approach, J. Comput. Phys. 125 (1996) 150– 160. [2] R. Abgrall, S. Karni, Computations of compressible multifluids, J. Comput. Phys. 169 (2001) 594–623. [3] M.R. Baer, J.W. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, Int. J. Multiphase Flow 12 (1986) 861–889. [4] G.T. Camacho, M. Ortiz, Adaptative Lagrangian modeling of ballistic penetration of metallic targets, Comput. Methods Appl. Mech. Eng. 142 (1997) 269–301. [5] S.F. Davis, Simplified second-order Godunov-type methods, SIAM J. Sci. Stat. Comput. 9 (1988) 445–473. [6] B. Després, F. Lagoutière, N. Seguin, Weak solutions to Friedrichs systems with convex constraints, Nonlinearity 4 (2011) 3055–3081. [7] C. Farhat, A. Rallu, S. Shankaran, A higher-order generalized ghost fluid method for the poor for the three-dimensional two-phase flow, computation of underwater implosions, J. Comput. Phys. 227 (2008) 7674–7700. [8] N. Favrie, S. Gavrilyuk, R. Saurel, Diffuse solid–fluid interface model in cases of extreme deformations, J. Comput. Phys. 227 (2009) 2941–2969. [9] N. Favrie, S. Gavrilyuk, Mathematical and numerical model for nonlinear viscoplasticity, Philos. Trans. R. Soc. A 369 (2011) 2864–2880. [10] N. Favrie, S. Gavrilyuk, Dynamics of shock waves in elastic–plastic solids, ESAIM Proc. 33 (2011) 50–67. [11] L.C. Forde, W.G. Proud, S.M. Walley, Symmetrical Taylor impact studies of copper, Proc. R. Soc. A 465 (2009) 769–790. [12] S. Gavrilyuk, N. Favrie, R. Saurel, Modelling wave dynamics of compressible elastic materials, J. Comput. Phys. 227 (2008) 2941–2969. [13] J. Glimm, D. Saltz, D. Sharp, Two-phase modelling of a fluid mixing layer, J. Fluid Mech. 378 (1999) 119–143. [14] S.K. Godunov, Elements of Continuum Mechanics, Nauka, Moscow, 1978 (in Russian). [15] S.K. Godunov, E.I. Romenskii, Elements of Continuum Mechanics and Conservation Laws, Kluwer Academic Plenum Publishers, NY, 2003. [16] S.K. Godunov, I.M. Peshkov, Thermodynamically consistent nonlinear model of elastoplastic maxwell medium, Comput. Math. Math. Phys. 50 (8) (2010) 1409–1426. [17] S. Karni, Multi-component flow calculations by a consistent primitive algorithm, J. Comput. Phys. 112 (1994) 31–43. [18] G. Kluth, B. Després, Perfect plasticity and hyperelastic models for isotropic materials, Continuum Mech. Thermodyn. 20 (3) (2008) 173–192. [19] G.H. Miller, P. Colella, A high-order Eulerian Godunov method for elastic–plastic flow in solids, J. Comput. Phys. 167 (2001) 131–176. [20] F. Petitpas, E. Franquet, R. Saurel, O. LeMetayer, A relaxation-projection method for compressible flows. Part II: The artificial heat exchange for multiphase shocks, J. Comput. Phys. 225 (2007) 2214–2248. [21] R. Saurel, O. Le Metayer, J. Massoni, S. Gavrilyuk, Shock jump relations for multiphase mixtures with stiff mechanical relaxation, Shock Waves 16 (2007) 209–232. [22] R. Saurel, R. Abgrall, A multiphase Godunov method for compressible multifluid and multiphase flows, J. Comput. Phys. 150 (1999) 425–467. [23] R. Saurel, F. Petitpas, R.A. Berry, Simple and efficient relaxation method for interfaces separating compressible fluids, cavitating flows and shock in multiphase mixtures, J. Comput. Phys. 228 (2009) 1678–1712. [24] R. Saurel, S. Gavrilyuk, F. Renaud, A multiphase model with internal degrees of freedom: application to shock–bubble interaction, J. Fluid Mech. 495 (2003) 283–321. [25] G.I. Taylor, The use of flat ended projectiles for determining yield stress. I. Theoretical considerations, Proc. R. Soc. A 194 (1948) 289–299. [26] E.F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves 4 (1994) 25–34. [27] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, second ed., Springer-Verlag, 1999. [28] H.S. Udaykumar, L. Tran, D.M. Belk, K.J. Vanden, An Eulerian method for computation of multimaterial impact with ENO shock-capturing and sharp interfaces, J. Comput. Phys. 186 (2003) 136–177. [29] B. Van Leer, Towards the ultimate conservative difference scheme IV. A new approach to numerical convection, J. Comput. Phys. 23 (1979) 276–299. [30] M.L. Wilkins, Calculation of elastic–plastic flow, Methods Comput. Phys. 3 (1964) 211–263. [31] M.L. Wilkins, M.W. Guinan, Impact of cylinders on a rigid boundary, J. Appl. Phys. 44 (3) (1973) 1200–1206.