An Eulerian multi-material scheme for elastic–plastic impact and penetration problems involving large material deformations

Accepted Manuscript An Eulerian multi-material scheme for elastic-plastic impact and penetration problems involving large material deformations C.D. Sijoy, Shashank Chaturvedi PII: DOI: Reference:

S0997-7546(15)00052-7 http://dx.doi.org/10.1016/j.euromechflu.2015.04.004 EJMFLU 2888

To appear in:

European Journal of Mechanics B/Fluids

Received date: 27 November 2014 Revised date: 11 March 2015 Accepted date: 10 April 2015 Please cite this article as: C.D. Sijoy, S. Chaturvedi, An Eulerian multi-material scheme for elastic-plastic impact and penetration problems involving large material deformations, European Journal of Mechanics B/Fluids (2015), http://dx.doi.org/10.1016/j.euromechflu.2015.04.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Click here to view linked References

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

An Eulerian multi-material scheme for elastic-plastic impact and penetration problems involving large material deformations Sijoy C. D. a,∗ Shashank Chaturvedi a a Computational

Analysis Division, Bhabha Atomic Research Centre, Visakhapatnam, India

Abstract An Eulerian multi-material scheme is developed for studying elastic-plastic impact and penetration problems involving large material deformations and shock waves. It consists of a Lagrangian plus remap (onto the original mesh) strategy at each time-step. The stress components of each material in a mixed cell are computed independently by assuming a common strain-rate to all materials. Similarly, the energy of each material in a mixed cell is updated independently and the mean pressure is determined using volume weighted average. This simplified approach eliminates commonly used iteration based pressure relaxation procedures and equation of state (EOS) closure models to find an equilibrium pressure in a mixed cell. A similar procedure is used to find equivalent stress in a mixed cell. The volume fraction of each material in a mixed cell is updated using a volume-of-fluid (VOF) interface tracking scheme. A new predictor-corrector (PC) scheme is developed to integrate the governing equations of the proposed multi-material model. The remapping of the stress components of the individual materials is achieved by a combination of a second-order monotonic upwind scheme and the VOF method. It is demonstrated that, despite using a ‘simple’ volume weighted scheme for calculating the total mean stress in mixed cells, the proposed scheme yields reasonable agreement with known analytical and experimental results of various impact and penetration problems. Key words: VOF, multi-material, hydrodynamics, compressible flow, volume-of-fluid, material-strength, elastic-plastic flow, impact and penetration

∗ Corresponding author, Tel: 91+891-2892156, Fax: 91+891-2742460 Email address: [email protected] (Sijoy C. D.).

Preprint submitted to Elsevier

11 March 2015

1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Introduction

Numerical simulation of hydrodynamic phenomena make use of two classical descriptions: the Lagrangian and Eulerian descriptions. In Lagrangian simulations the computational cell move with fluid velocity. There is no mass flow crossing the cell boundaries, i.e. it allows an easy tracking of material interfaces. The Lagrangian methods, therefore, capture the contact discontinuity sharply in multi-material flows. However, it fails when the material deformations are large enough to give severe distortions to the computational mesh. On the other hand, the Eulerian methods are characterized by a fixed computational mesh through which the material flows. Therefore, the Eulerian methods are capable of handling large material deformations. However, they are less accurate and diffusive compared to the Lagrangian methods. Also, the multi-material simulations using pure Eulerian methods demand mixed equation of state (EOS) evaluation or multi-material EOS closure models in mixed cells (computational cells containing more than one material). The continuous rezoning of computational mesh (e.g. ALE type algorithm) can overcome the limitations of pure Lagrangian method. This method provides the accuracy of the Lagrangian and the robustness of the Eulerian approaches. The rezoning can be performed continuously or only when the Lagrangian simulation starts failing due to non-convex or negative-volume cell (large deformations). Another approach for multi-material hydrodynamic problems is to track the material interface on a fixed (Eulerian) or arbitrary mesh (ALE) using a specialized technique called interface reconstruction or interface tracking [1–8]. The terminology of material interface-tracking in multi-material hydrodynamics simulations stands for an advanced numerical scheme to explicitly track the material interfaces within a computational cell and to treat them distinctly so that the materials do not interpenetrate. The Eulerian and ALE formulations permit the material interfaces to be independent of the computational mesh. The authors are aware of the existence of other methods in the context of the simulation of fluid/elastic-plastic flows using a Lagrangian or ALE scheme, such as the classical cell-centered schemes [9–15]. The cell-centered scheme seems to be a promising alternative to the usual staggered scheme and it does not require the artificial viscosity terms since the viscosities are inherent in the scheme. Moreover, in the context of ALE strategy, it allows straightforward implementation of conservative remapping schemes due to the fact that the physical variables are located at the cell centers, see for instance [16–18]. The higher order extension of these schemes can be found in Refs. [12,19,20]. In the present work, however, we have followed an Eulerian scheme (La2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

grangian plus remap strategy at each time-step) in which the material interfaces are explicitly tracked by using a volume-of-fluid (VOF) method. The VOF scheme is also used to remap/advect various physical quantities at each time-step. In general, multi-material Eulerian or ALE schemes employ a pressure relaxation algorithm based on an iteration procedure so that the individual material pressure get relaxed to a common equilibrium pressure [21–23]. In each iteration step, the volume/mass fraction of each material is adjusted until the pressures of those materials are converged to a common equilibrium pressure. However, in reality, the pressure equilibrium in a computational cell cannot take place in a time shorter than that required for a sound wave to propagate through the cell at least once. For numerical stability, most of the hydrodynamic calculations are performed with a numerical time-step much smaller than this time. Therefore, such assumptions are not strictly valid while numerically solving most of the real-life problems. Apart from this, for materials having highly non-linear EOS, the iteration procedure may diverge. Similarly, the assumption of same temperature and/or pressure to all materials in a mixed cell may lead to unrealistic energy flow from a hot material to a cold material after the iteration procedure. For more details on various multi-material mixture theories and related issues, see Refs. [23–25]. Also, a detailed discussion on different multi-material relaxation methods for ALE hydrodynamics and their analysis can be found in Ref. [26]. Similarly, a comparative study of multi-material Lagrangian and Eulerian methods with different multi-material closure models (pressure relaxation schemes) is described in Ref. [27]. The objective of the present work, therefore, is to develop a ‘simple’ and robust multi-material Eulerian algorithm for elastic-plastic problems involving large material deformations (e.g. high velocity impact and penetration problems) by using a relatively ‘simple’ multi-material model to evaluate mean stress in mixed cells. In the present algorithm we do not attempt to relax individual material pressure to a common equilibrium pressure by adjusting the material volume fractions. We calculate the effective cell pressure by using a ‘simple’ volume weighted average. Similarly, the equivalent stress in a mixed cell is evaluated using volume weighted average, where the individual stress components for each material in a mixed cell are updated with an assumption of common strain rates to all materials present in a cell. Thus, the iteration procedure for pressure relaxation in a mixed cell can be avoided. Such an approach for fluid problems (without material strength) can be found in Ref. [28], which is based on the theoretical formulation given in Ref. [22]. In Ref. [22], the cell volume changes are partitioned among the materials in a mixed cell in proportion to their bulk modulus. This avoids unphysical material compression when materials with different compressibility are mixed in a cell. Even though the present approach assumes common strain rates to all materials present in a mixed cell, irrespective to their compressibility, it is ro3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

bust and yields satisfactory results, with the proposed PC integration scheme, for various complex impact and penetration problems presented in this work. Also, see Ref. [28] for various other numerical tests without material strength model. A brief description of the present scheme and the preliminary results of various impact and penetration problems with material strength model are given in our previous work [31]. It is worth to mention here that the present scheme, however, demands the evaluation of M different energy equations in a mixed cell (discussed in detail in later sections), where M is the total number of materials contained in a mixed cell [22,28,31]. In the remaining sections, we have described this ‘simple’ Eulerian algorithm with multi-material strength model for the simulation of impact and penetration problems involving large material deformations and its validations. As mentioned earlier, we have also proposed a new PC integration scheme to solve the governing equations. The remainder of the paper is organized as follows. The basic physical model (with multi-material strength model) and governing equations are given in Sec. 2. In Sec. 3, we have discussed the proposed predictor-corrector integration scheme. The details of the advection/remap step based on a combination of a second order accurate monotonic upwind scheme (MUSCL) [32] and VOF method are discussed in Sec. 3.2. The VOF scheme used for material interface tracking is briefly discussed in Sec. 6. The details of this scheme can be found in Refs. [28–30]. In Sec. 4, we have discussed the details of flux limited scalar artificial viscosity [33,34]. Several test problems are provided for the validation of the code in Sec. 7. Finally, the conclusions and future scope of this work are given in Sec. 8.

2

Physical model and governing equations

The conservation equations in two-dimensional cylindrical Eulerian coordinates can be written as Dρ + ρ [ǫzz + ǫrr + ǫθθ ] = 0 Dt

(1)

∂σzz ∂σrz σrz Du = + + Dt ∂z ∂r r

(2)

∂σrr ∂σrz σrr − σθθ Dv = + + Dt ∂r ∂z r

(3)

ρ

ρ

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

ρ

DE = σij ǫij Dt

(4)

where, ρ is the density, u and v are the velocity components in axial and radial directions respectively and E is the specific internal energy. Also, the following definitions are used.

ǫzz

∂u = ; ∂z

σij = Sij − pδij

(5)

D ∂ ∂ ∂ = +u +v Dt ∂t ∂z ∂r

(6)

∂v ǫrr = ; ∂r

ǫθθ

v = ; r

ǫrz

1 = 2

∂u ∂v + ∂r ∂z

!

(7)

Here, Sij and ǫij are the stress components and the strain rates respectively. The hydrodynamic pressure, p, is evaluated using EOS of the material: p = f (ρ, ǫ). The multi-material hydrodynamic formulation used in this paper can be seen in Refs. [22,28], where the strength of the material is neglected. We have included the material strength model in a similar manner. The mass conservation equation for individual materials in terms of its volume fraction in a mixed cell is given by the following expression [22,28]. ∂ (f m ρm ) ∂ (f m ρm u) 1 ∂ (rf m ρm v) + + = 0; ∂t ∂z r ∂r

m = 1, .., M

(8)

Here, f m (0 ≤ f m ≤ 1) is the volume fraction of material m and M is the total number of materials present in a mixed cell. The total density is given by P ρ = m f m ρm , where the summation is for all the materials present in a mixed cell. The assumption of locally adiabatic evolution of internal energy leads to the following equation for internal energy update for each material [22,28].

∂E m ∂E m ∂E m +u +v ∂t ∂z ∂r " # ∂ρm ∂ρm ∂Esm pm ∂ρm +u +v ; m = 1, .., M + = m 2 (ρ ) ∂t ∂z ∂r ∂t

(9)

The details of the derivation of Eq.(9) and the physical assumptions used are 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

given in [22,28]. Here, an additional term ∂Esm /∂t due to stress components is added to energy equation. ∂Esm 1 = m (Szz ǫzz + Srr ǫrr + Sθθ ǫθθ + 2Srz ǫrz ) ; ∂t ρ

m = 1, .., M

(10)

The stress components in Eq. (10) are the average values evaluated using mixture theory (discussed below in this section). Individual materials can now possess different energies and temperatures. This leads to an evaluation of M different energy equations in each mixed cell. One more equation, apart from the equation of state (EOS) of each material, is required to close the system which is given in the following. ∂f m ∂f m ∂f m +u +v = 0; ∂t ∂z ∂r

m = 1, .., M

(11)

The average pressure in a mixed cell is calculated using the following equation. p=

X

f m pm

(12)

m

Assuming common strain rates to all materials present in a cell, the material stress components for each material can be evaluated using the following equations. 







m m ∂Szz ∂Srr Sv Sv m m + δzz , + δrr , = 2µm ǫzz − = 2µm ǫrr − ∂t 3 ∂t 3   m m ∂Srz ∂Sθθ Sv m , = 2µm ǫθθ − = 2µm ǫrz + δrz ; m = 1, .., M ∂t 3 ∂t

(13)

where, µm is the shear modulus for material m and Sv = tr [ǫij ]. Also, δijm ’s are the correction for the rigid body rotation for each material as given by the following expressions [35].

m m Szz − Srr m (cos 2ω − 1) − Srz sin 2ω 2 m S m − Srr m = Srz (cos 2ω − 1) + zz sin 2ω 2 m m δrr = −δzz

m δzz = m δrz

Here, ω is the rotation angle defined as in the following. 6

(14)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

∆t sin ω = 2

∂v ∂u − ∂z ∂r

!

(15)

q √ The Von-Mises Yield condition ( J m − 23 Y m ≤ 0, where J m is the second invariant of the stress tensor and Y m is the flow stress of the material m [35]) is applied to the stress components of each material to model elastic-plastic flow. The quantity J m is determined by using Eq. (16).

J

m

=

m 2 (Szz )

+

m 2 (Srr )

+

m 2 (Sθθ )

+

m 2 2(Szr );

Sijm = Sijm ; Sijm = Sijm Jvm ;

s

2 Ym √ 3 Jm if (Jvm ≥ 1) if (Jvm < 1)

Jvm =

(16)

The mean stress is approximated by using Eq. (17), which is similar to the equation used for calculating average pressure in a mixed cell [23,36].

Sij =

X

f m Sijm

(17)

m

3

Numerical method and governing equations

A fixed Eulerian mesh is used with a virtual Lagrangian calculation. At each time-step, the effects of Lagrangian deformations are remapped onto the original mesh during the Eulerian advection and remap step. We have used a spatially staggered mesh, where velocities are defined at cell faces and density, internal energy, pressure, stress etc. are defined at the mid point of the cell, see Fig. 1. All the strain-rates are defined at the mid point of the cell except the shear strain rate which is first calculated at the cell corners and then averaged to the cell centers using its four corner values. The physical variables are updated using a ‘directional split’ method. In each time-step, first a Lagrangian calculation is performed in one direction (say radial) followed by a radial remap step. After these sets of calculations, a Lagrangian step is performed in the other direction (axial) followed by an axial remap step. To avoid the biasing effects the sequence of steps in the radial and axial directions are alternated during subsequent time-steps. In the following, we have described the details of this numerical algorithm. 7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

j+1

v

i+1/2,j+1

u

u

i+1/2, j+1/2

i,j+1/2

j

i+1,j+1/2

v

i+1/2,j

i+1

i

Fig. 1. Schematic of staggered computational grid.

3.1 Lagrangian step The Lagrangian phase calculation is divided into radial and axial directions using a directional split strategy. The radial equations can be obtained by dropping axial terms and vice versa. The stress update equations are further simplified using the continuity equation. The radial part of the equations are more complicated than the axial part. Therefore, the detailed discussions are given only for the calculations in the radial direction. 3.1.1 Radial direction update The equations (only the final equations are given) to be solved in the Lagrangian phase (radial direction) can be summarized as given below: ∂ρm ρm ∂ (rv) + = 0; ∂t r ∂r

m = 1, .., M

(18)

#

(19)

"

∂v 1 ∂σrr σrr − σθθ = + ∂t ρ ∂r r "

∂u 1 ∂σrz σrz = + ∂t ρ ∂r r "

#

#

∂u 1 ∂Esm = m σrr ǫrr + σθθ ǫθθ + σrz ; ∂t ρ ∂r 



m ∂Srr 2 1 = 2µm ǫrr − ǫθθ ; ∂t 3 3

8

(20)

m = 1, .., M

m = 1, .., M

(21)

(22)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65





m ∂Szz 1 = 2µm − (ǫrr + ǫθθ ) ; ∂t 3 m ∂u ∂Srz = µm ; ∂t ∂r

m = 1, .., M

m = 1, .., M

m Sθθ = −(Srr + Szz );

(23)

(24)

m = 1, .., M

(25)

The Eq. (18) is obtained by considering the fact that the Lagrangian step conserve the mass and volume fraction of each material in a cell. The average stress components, Sij , in mixed cells are obtained using Eq. (17). The above equations are solved using a predictor-corrector (PC) integration scheme as described in the following. The right-hand side of the equations are approximated at half time-step (predictor step). The half time-step pressure is evaluated using Eq. (26) based on the assumption of isentropic evolution of material pressure [28,37,38]. dρm dpm = (cm )2 ; dt dt

m = 1, .., M

(26)

In Eq. (26), cm is the sound velocity of the material and dρm /dt is the Lagrangian rate of change of density. The Lagrangian step conserves the mass and volume fraction each material in a cell. Therefore, using Eq. (18), the Eq. (26) can be discretized for a half-time step as below:

(pm )n+1/2





∆tρm (cm )2 rj+1 vi+1/2,j+1 − rj vi+1/2,j ; = pm − 2 rj+1/2 (rj+1 − rj )

m = 1, .., M (27)

Here, rj and rj+1/2 are the radial position coordinates at the nodes and cellcenter respectively. Various terms in the right hand side of Eq. (27) are at time-step ‘n’. Also, common subscripts (i + 21 , j + 12 ) to pm , ρm and cm are omitted in the Eq. (27). The predicted values of stress components are obtained by using an explicit time integration of Eq. (22) to Eq. (25) using a time-step of ∆t/2. The governing equations are given below.

m n+1/2 m n (Srr ) = (Srr ) +

9

∆t (µm )n n [2ǫrr − ǫnθθ ] 3

(28)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

m n+1/2 m n (Szz ) = (Szz ) −

m n+1/2 (Srz )

=

m n (Srz )

∆t (µm )n n [ǫrr + ǫnθθ ] 3

∆t (µm )n + 2

∂u ∂r

h

(29)

!

m n+1/2 m n+1/2 m n+1/2 (Sθθ ) = − (Srr ) + (Szz )

(30)

i

(31)

The variables in the above equations are defined atthecell center (i + 21 , j + 21 ). The cell-centered shear strain rate for radial step ∂u is calculated by taking ∂r the average of the shear strain rates defined at the corners of the cell, see Eq. (32). The shear strain rate at the corner of the cell is calculated by using Eq. (33). ∂u ∂r

!



1 ∂u =  4 ∂r

!

i,j

∂u + ∂r

∂u ∂r

!

i,j

!

i+1,j



∂u + ∂r

!

∂u + ∂r i+1,j+1

ui,j+1/2 − ui,j−1/2

= 

rj+1/2 − rj−1/2



!

i,j+1



 

(32)

(33)

The strain-rate components ǫrr and ǫθθ are calculated using the following equations.

ǫrr =



vi+1/2,j+1 − vi+1/2,j (rj+1 − rj )



;

ǫθθ =



vi+1/2,j+1 + vi+1/2,j 2rj+1/2



(34)

Once the calculation of predicted stress components of individual materials in a mixed cell are completed, the average predicted stress values (Sij ) are estimated with the help of mixture theory (volume weighted average). For pure cells, where only one type of material is present, the average predicted stress values (Sij ) are equal to the material stress values (Sijm ). These predicted values of pressure and stress components are then used to calculate tentative full time-step velocities (˜ u, v˜) by an explicit time integration of Eqs. (19) and (20).

n v˜i+1/2,j = vi+1/2,j +

∆t ρn+1/2

"

σrr(i+1/2,j+1/2) − σrr(i+1/2,j−1/2) (σrr,θθ )a + rj+1/2 − rj−1/2 rj

10

#n+1/2

(35)

n+1/2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Here, σij = [Sij − (p + q)δij ]n+1/2 and ρn+1/2 = Σm f m (ρm )n+1/2 (half-time step value of ρ). The half-time step values of ρm for each material in a mixed cell are calculated using Eq. (18), since the Lagrangian step conserve the mass and volume fraction of each material in a cell. The discretization of this equation is already discussed (see, Eq. (27)). Also, (σrr,θθ )a is the average of the quantity (σrr − σθθ ) from the cells (i + 12 , j + 21 ) and (i + 12 , j − 12 ). It is worth to note, the singularities at r = 0 in the Eq. (22) and Eq. (23) will not effect the discretization described above due to the staggered mesh formulation employed here. The radial velocities at r = 0, i.e. vi+1/2,j for j = 1 is calculated using boundary conditions. For axisymmetric case those velocities are equal to zero. Similarly, the axial velocities are defined at x-edges (i, j + 21 ) of the cell, see Fig. 1, where rj+1/2 6= 0. The tentative half time-step velocities, (˜ un+1/2 , v˜n+1/2 ) are defined as an average of tentative velocities at full time-step, (˜ u, v˜) and velocities at the beginning of this time-step, (u, v).

n+1/2

u˜i,j+1/2 =

ui,j+1/2 + u˜i,j+1/2 ; 2

n+1/2

v˜i+1/2,j =

vi+1/2,j + v˜i+1/2,j 2

(36)

Once the tentative values of the velocity components are known, a corrector step is applied to re-calculate the pressure and stresses (corrected values of pm and Sijm ). The corrected pressure is determined using the following equation.

(pm )n+1/2





∆tρm (cm )2 rj+1 v˜i+1/2,j+1 − rj v˜i+1/2,j = pm − ; 2 rj+1/2 (rj+1 − rj )

m = 1, .., M (37)

Also, the corrected stress values for full time-step, (Sijm )n+1 are then calculated using Eq. (28) to Eq. (31) with strain-rates (ǫij ) evaluated using the velocity components (˜ un+1/2 , v˜n+1/2 ). One can now determine the half timestep values of the corrected stress components by using a simple time average, (Sijm )n+1/2 = 12 ((Sijm )n + (Sijm )n+1 ). Similar to the predictor step, the average n+1/2 stress components in a mixed cell (Sij and Sijn+1 ), are estimated with the help of mixture theory (volume weighted average). Next, the corrected Lagrangian velocities are obtained at full time-step using the following equations.

n+1 vi+1/2,j

=

n vi+1/2,j

"

∂σrr σrr − σθθ + ∂r r

∆t + ρ 11

!#n+1/2

(38)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

un+1 i,j+1/2

"

∆t = uni,j+1/2 + ρ

∂σrz σrz + ∂r r

!#n+1/2

(39)

n+1/2

Here, σij = [Sij − (p + q)δij ]n+1/2 are the corrected values at half timestep and q is the artificial viscosity. The above equations are discretized in the same manner as discussed earlier for the predictor step calculation. Finally, the strain energy and Lagrangian change in internal energy are updated using Eq. (40) and Eq. (41).

∆Esm

"

∆t ∂u = m σrr ǫrr + σθθ ǫθθ + σrz ρ ∂r

E˜ m = E m −



!#n+1/2

∗ ∗ (pm )n+1/2 rj+1 vi+1/2,j+1 − rj vi+1/2,j

(ρm )n+1/2 rj+1/2 (rj+1 − rj )



;

m = 1, .., M

(40)

i+1/2,j+1/2



∆t + ∆Esm ;

m = 1, .., M (41)



n+1 ∗ . The Eq. (41) is derived using the contiHere, vi+1/2,j = 21 vi+1/2,j + vi+1/2,j nuity and volume fraction update equations. The strain energy due to shear stress and shear strain 2σrz ǫrz = σrz ∂u in Eq. (40) is obtained by dropping ∂r axial terms from Eq. (7) for ǫrz . The other part of the shear strain energy σrz ∂v ∂z is added in the axial step of the directionally split algorithm. This completes the Lagrangian phase calculations in the radial direction.

3.1.2 Axial direction update The Lagrangian update equations in the axial direction are solved using the same PC sequence described in Sec. 3.1.1 for the radial direction. The relevant equations in the axial directions are slightly different from those described in Sec. 3.1.1. However, its derivations are straightforward. We have omitted those equations and descriptions here for brevity. 3.2 Remap or advection steps The Eulerian advection terms are treated in the remap-step of the algorithm. All the solution variables are advected in a conservative manner even though they are not necessarily governed by physical conservation laws. The advected mass, energy, stress components (for each material) and the momentum contained in the overlap region in between the deformed Lagrangian co-ordinates 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and fixed Eulerian grid positions are found with the help of VOF algorithm. The method for the evaluation of this quantity using a VOF technique is described briefly in Sec. 6 and its details can be found in Refs. [28,29].

3.2.1 Mass and energy (radial direction) The discretization of continuity equation for material m during the radial remap step can be written as

(f m ρm )n+1 = (f m ρm )n −

  ∆t m Φm − Φ j rj+1/2 (rj+1 − rj ) j+1

(42)

where the signed mass flux (Φm j ) is defined as

Φm j =

Mjm,a . ∆z∆t

(43)

The quantity Mim,a is the signed advected mass, evaluated using VOF method [28]. The same quantity is used to get the advected energy. The energy flux is obtained by multiplying the mass flux and the internal energy after Lagrangian step in the donor cell [28]. The following equation is used for updating individual energies in the remap step.

∆t rj+1/2 (rj+1 − rj ) h    i m m n+1/2 ∗ ∗ m m m ˜ × (f p ) rj+1 vj+1 − rj vj + E˜j+1 Φm − E Φ j+1 j j (f m ρm E m )n+1 = (f m ρm E m )n −

(44)

The equation includes the effect of the Lagrangian deformation and the Eulerian advection in a single expression. The upwind values appearing in the advective fluxes are given by

m E˜jm = E˜j−1/2 , m m E˜j = E˜j+1/2 ,

vj∗ > 0 vj∗ ≤ 0

(45)

The new density for each material is calculated using the volume fraction f m updated by the volume-of-fluid method. 13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

 

n+1 ρm j+1/2

em j+1/2

n+1

=

=

(f m ρm )n+1 j+1/2 m,n+1 fj+1/2

(46)

(f m ρm em )n+1 j+1/2 (f m ρm )n+1 j+1/2

The above equations for the radial part of the directional split algorithm are valid for the axial direction with few more or less trivial modifications. We do not list those equations here for brevity. 3.2.2 Stress advection (radial direction) In many hydro-codes, the solution and history variables (strain, plastic strain etc.) are transported using a second order monotonic advection scheme. If the material stresses are updated independently for each material (also the materials are assumed to be separated by an interface) in a mixed cell, a general advection procedure for average stress components is not suitable. This situation is depicted in Fig. 2, where a solid material (m1 ) and gas (m2 ) are separated by a material interface in a mixed cell. For this case, the average stress of the mixed cell is mainly contributed by material m1 , due to its high material strength and volume fraction. Therefore, if one uses the average stress for the advection, the computational cell right to this mixed cell (assuming uj is positive) which contains material with no strength (gas) will be assigned unphysical stress values after the advection step. To eliminate this problem we have advected the individual stress components for each material using their advected volume evaluated using VOF algorithm. The solution update is as given by Eq. (47). 1111111111111111111111111 0000000000000000000000000 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111 0000000000000000000000000 1111111111111111111111111

m1

Uj x dt

m2 j

j−1

Fig. 2. Advection volume for a typical mixed cell to a pure cell for a typical time-step.

n+1 (φm = j+1/2 )

1 m (Vj+1/2 )n+1

h

n m n m m ˜ ˜ (φm j+1/2 ) (Vj+1/2 ) + φj ∆Vj − φj+1 ∆Vj+1

14

i

(47)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In Eq. (47), the nodes are numbered as j, j + 1 etc. and the cell defined by nodes j and j + 1 is j + 21 . Also, φm = Sijm is the stress components for material m, ∆V m is the signed volume advected across the node j for material m and φ˜j is the upwind values of stress components at node locations. The accuracy can be improved by calculating second order accurate values of φ˜j . Similarly, monotonicity can be achieved by limiting their values. We have used a second order accurate monotonic upwind scheme (MUSCL) due van Leer [32] as it is one of the simplest and most efficient algorithms that is both second order accurate and monotonic. The details of the algorithm can found in Ref. [32]. It is worth mentioning here that there is no guarantee that the individual materials will satisfy the yield criterion after the advection step. However, this discrepancy will be removed by itself in the next time-step as all material stress values are adjusted as per Von-Mises yield criterion.

3.2.3 Momentum advection (radial direction) The numerical algorithm for momentum advection follows the method given by Benson in [39]. Niem et al. [28] have given a generalized algorithm in cylindrical co-ordinates by using the formulations given in [39]. We have used these algorithms to advect the momentum in a fixed staggered Eulerian mesh. The basic theoretical model and derivations of equations used here can be found in [28,39]. The total signed Eulerian mass flux across edge j is given by

Φj =

X

Φm j .

(48)

m

The momentum remap equation in flux difference form is written as [28] 



(Mj vj )n+1 = (Mj vj )n + Φj−1/2 vˆj−1/2 − Φj+1/2 vˆj+1/2 ∆t

(49)

where, Mj vj is the radial momentum defined at cell edge, Φj+1/2 is the geometrically weighted average of face centered Φj values and vˆj+1/2 is a monotonized second order flux limited interpolated value of the velocity at cell center. Different flux limiters used are described in Section 4. More details on the interpolation function can be found in [28]. The axial momentum update during r-sweep needs more attention. Here the face centered mass fluxes are first geometrically averaged to node locations. Finally, the axial momentum is updated using these node located mass fluxes. 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

 1 Φi+1/2,j + Φi−1/2,j 2 (50) = (Mi,j+1/2 ui,j+1/2 )n + (Φi,j uˆi,j − Φi,j+1 uˆi,j+1 ) ∆t

Φi,j =

(Mi,j+1/2 ui,j+1/2 )n+1

Here, uˆi,j is a monotonized second order flux limited interpolated value of the velocity at node locations. A geometrical weight factor consistent to cylindrical geometry given in [28] is used for averaging and interpolation. Also, it is claimed in [28] that such a choice will conserve the total energy and momentum in cylindrical co-ordinate system. The equations for the momentum advection during z-sweep follow the same procedures described above, which we are omitted here for brevity.

3.2.4 Stress correction for rigid body rotation The new stresses, after material rotation or deformation, have to be expressed in terms of the original co-ordinate system. One possibility is to include this correction into the governing differential equations used for the calculations in the radial and axial directions. However, we have followed the method given in Ref. [35]. The rotational correction to the stress components is applied after completing the Lagrangian and remap step in both directions (radial and axial) before they are used for next time-step calculations. These corrections are performed by using the equations given in Sec. 2. The equations are discretized using a formulation given in Ref. [35].

4

Flux limited artificial viscosity

We have used a cell centered flux limited scalar artificial viscosity [33,34]. A detailed discussion on flux limited artificial viscosity and its advantages are given in [33,34]. Due to the directional split algorithm a one-dimensional formulation is sufficient to evaluate the artificial viscosity. The general form of flux limited artificial viscosity is given by Eq. (51) 



ql,i+1/2 = qg,i+1/2 1 − Ψi+1/2 ;



if ∆ui+1/2 < 0



(51)

where, Ψi+1/2 is a higher order limiter function, ui+1/2 is the velocity difference (ui+1 − ui ) and qg,i+1/2 is the normal scalar artificial viscosity. Currently, three different artificial viscosity qg,i+1/2 options are implemented: VonNeumann quadratic [40], VonNeumann quadratic plus linear and Wilkins strain 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

rate (along the direction of shock wave) dependent artificial viscosity [35]. A typical flux limiter given in [34] is given by Eq. (52).

Ψi+1/2 = max [0, min (0.5 (rl + rr ) , 2rl , 2rr , 1)] (ui+2 − ui+1 )/(zi+2 − zi+1 ) rl = (ui+1 − ui )/(zi+1 − zi ) (ui − ui−1 )/(zi − zi−1 ) rr = (ui+1 − ui )/(zi+1 − zi )

(52)

Various other flux limiter functions are implemented for flux limited artificial viscosity and momentum advection calculations. There are currently five possible choices for the flux limiter function: Minmod, Superbee, monotonized centered, van Leer and Osher [41].

5

Time-step estimate

The time-step calculation is based on the ideas given in [42]. Presently the following form given in Ref. [42] is used. cCF L ∆L

∆t = q+

rh

(C + u)2 + q 2

i

(53)

where, cCF L is the CFL constant, q ∆L grid length, q is the artificial viscosity, C is the acoustic wave speed ( (B + 4/3G) /ρ), u is the particle velocity, B is the Bulk modulus and G is the Shear modulus.

6

Volume-of-fluid (VOF) method

Detailed description of VOF method can be found in Refs. [2,3,43–45] and a review on VOF methods can be found in [46]. Present formulation follows our earlier work [29]. The algorithm consists of three parts: Interface reconstruction, Lagrangian deformation of material interfaces and finally the Eulerian transport (advection). A ‘directional-split’ advection procedure is used to extend this method into two-dimensional case. More details can be found in Ref. [29]. However, here, we have listed few important aspects of the algorithm. 17

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

A piece-wise linear interface construction (PLIC) method is used. The interface line is represented by ~n · ~x − α = 0; where, ~n = ˆinx + ˆjny is the exterior normal to the line and α is the perpendicular distance from a local origin in a cell. The local origin in a cell for each material is decided by checking the slope of the interface line for that material [29]. The slope of a material interface line ~ m /|▽f ~ m |) are determined using LVIRA and hence the normal vector (~n = ▽f method [47]. Having obtained the interface normals, the interface parameter α has to be determined so that volume behind the interface line is equal to the material volume. The volume behind an interface line is calculated using an algorithm in Ref. [29]. This method explicitly constructs the points of polygon enclosing the fluid in a cell by checking intersection between material interface line and cell boundaries. The interface parameter α is found when the function f (α) = V (α) − Vactual , becomes zero (root-finding). Here, V (α) is the material volume in the cell bounded by the interface line with line parameter α. We have used Brents method [48] to find α with αmin = 0 and αmax = d as initial bracket; where d is the cell diagonal distance. The material interfaces in a mixed cell are constructed using ‘onion-skin’ or ‘layered’ model with the help of a ‘material order list’ determined dynamically. The ‘material order list’ is also used for specifying the order of transport for each material in a mixed cell to its neighboring cells. We have used a simple combination of centroid check algorithm [49] and Benson’s least squares fit to centroid algorithm [50] for dynamic material ordering (see Ref. [29] for details). The authors are aware of more advanced interface construction and tracking schemes [5–8] in which an automatic material ordering in a mixed cell has been achieved. They are more accurate than the present scheme and are capable of handling the ‘T’ and ‘Y’ junctions in a mixed cell. This possible improvement in our algorithm will be explored in future. The material volume advected (Vim ) is evaluated using a generalized algorithm given in [29]. The material volume fraction for each material is updated using the formulations given in [28] for cylindrical co-ordinates. The advected mass (Mim ) is calculated using a characteristic trace back method [28] by assuming a linear interpolation of velocity along the streamline in a computational cell. Authors in Ref. [28] have demonstrated that this method conserve total mass in cylindrical co-ordinate system.

7

Results and discussions

The numerical experiments described in following sections are performed with Courant limited time-step and a CFL coefficient ∼0.6 is used. 18

7.1 Sedov problem in cylindrical geometry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

First, since the present algorithm employs a directional splitting strategy, we assess the ability of the scheme to preserve the radial symmetry of the solutions for spherically symmetric flows in 2D. For this analysis, we have chosen the classical Sedov point-blast problem in 2D axisymmetric geometry. The Sedov point-blast problem in cylindrical geometry is often used to deal with strong shock-waves in curvilinear co-ordinates. In this case, the resulting shock wave evolves in a self-similar fashion in which the shock radius, shock velocity, peak density, etc. can be found analytically.

The material is treated as an ideal gas with γ = 1.4. The initial density is equal to unity and the initial velocity is zero. At t = 0, the total blast energy E0 is purely internal and deposited at the origin. The analytical solution gives an expanding shock of radius, rs = (E0 /(αd ρ0 ))1/(2+d) t2/(2+d) , where d = 3 for spherical geometry and αd = 0.850937 with a peak shock density, ρs = 6. The total energy E0 is defined such that rd = 1 at time t = 1. Due to the symmetry of the problem only one-quarter of the system in 2D cylindrical coordinates is simulated. We have used a square domain of size 1.2 divided into 200×200 square zones. The initial energy E0 is deposited in a computational cell (i,j)=(1,1) near the origin.

The results obtained at a time t=1, is shown in Fig. 3. The ρs ∼ 5.0 obtained in the simulation (with 200×200 mesh) is close to the analytical value 6. Similarly, inspecting the 2D density profile shown in the left side panel of Fig. 3, it is clear that the spherical symmetry of the solution is maintained reasonably well. In order to assess this property in detail, the radial (R = √ x2 + y 2 ) density profile along the 45 degree line w.r.t the z-axis (the red line in the left side panel of Fig. 3) is compared with the density profile along the z-axis (along the cells with j = 1, where ‘j’ is the radial cell indices). This is plotted in the right side panel of Fig. 3. It is clear that these two density profiles agree well (within ∼3.5% error in ρs ), i.e., the expected symmetry is preserved reasonably well. Note that no significant difference in the density profile along the z-axis (along j = 1 cells) and r-axis (along i = 1 cells) is observed. The symmetry error, defined as the ratio of the difference in ρs obtained in these directions w.r.t. the average ρs value, with 100×100 and 400×400 grids are ∼8.2 and 2.2%, respectively. The symmetry error in the shock radius rs is found to be negligible when the number of cells in each direction is greater than typically 80, below which the spherical symmetry of solution is not well maintained. For example with a 40×40 mesh the symmetry error in rs is ∼14%. 19

6

Analytical o along 45 line (red) axial variation (blue)

Mesh = 200x200

5

4 Density

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3

2

1

0 0

0.2

0.4

0.6

0.8 1 Radial distance

1.2

1.4

1.6

1.8

p Fig. 3. The left and right side plots show the 2D RZ & 1D radial R = x2 + y 2 density profiles, respectively at a time t = 1 for the Sedov problem. The red line in the left side plot shows the 45 degree line w.r.t the z-axis and the corresponding 1D radial density profile is shown in the right side panel in red color.

7.2 Two material Sod problem In this section, we have considered a 1D Sod shock tube problem consists of two different materials [24]. The computational domain is [0,1] with a material interface at x = 0.5. The initial conditions for the left and right materials are (ρL , PL , γL ) = (1, 2, 2) and (ρR , PR , γR ) = (0.125, 0.1, 1.4), respectively. The values mentioned above are in M.K.S units. The initial velocity is assumed to be zero in the entire computational domain. An ideal gas EOS is assumed. The initial grid is made of 100 (case 1) and 200 (case 2) cells. In each case, we have considered two different limiters (upwind and Vanleer) in the advection procedures. The numerical results against the analytical solutions are plotted in Fig. 4. The results are in good agreement with the analytical solutions. Moreover, no spurious pressure oscillations are created at the contact discontinuity. We note that the pressure rarefaction fan and the shock discontinuities are best reproduced by the use of Vanleer limiter in the advection calculations. However, an overshoot in velocity and a corresponding density dip near the tail of the rarefaction fan are observed with the use of Vanleer limiter. The use of upwind limiter in the advection calculations produces slightly diffusive results. Note that the overshoot in velocity is absent with the upwind limiter. Finally, it is worth mentioning here that our simulation results are consistent with the results presented in Refs. [16,24].

7.3 1-D impact problem (Wilkins’s test problem) Next, a 1-D impact test problem due to Wilkins [51] is used to study the performance of the proposed predictor-corrector algorithm. A piece of aluminum flyer with a thickness of 4 mm and an impact velocity of 2 km/s, impacts an 20

1

200 cells

0.8 Density (kg/m3)

0.7 0.6 0.5 0.4

0.7 0.6 0.5 0.4

0.3

0.3

0.2

0.2

0.1

0.1 0

0.1

0.2

0.3

0.4 0.5 0.6 Distance (m)

2

0.7

0.8

0.9

1

0

1.6

0.1

0.2

0.3

0.4 0.5 0.6 Distance (m)

2

Analytical upwind limiter VanLeer limiter

1.8

100 cells

0.7

0.8

0.9

1

Analytical upwind limiter VanLeer limiter

1.8 1.6

200 cells

1.4 Pressure (Pa)

1.4 Pressure (Pa)

Analytical upwind limiter VanLeer limiter

0.9

100 cells

0.8 Density (kg/m3)

1

Analytical upwind limiter VanLeer limiter

0.9

1.2 1 0.8

1.2 1 0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0 0

1.4

0.2

0.4 0.6 Distance (m)

0.8

1

0

1.4

Analytical upwind limiter VanLeer limiter

1.2

0.2

0.4 0.6 Distance (m)

0.8

1

0.4 0.6 Distance (m)

0.8

1

Analytical upwind limiter VanLeer limiter

1.2

100 cells

200 cells 1 Velocity (m/s)

1 Velocity (m/s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.8 0.6

0.8 0.6

0.4

0.4

0.2

0.2

0

0 0

0.2

0.4 0.6 Distance (m)

0.8

1

0

0.2

Fig. 4. The two material Sod problem density, pressure and velocity with two different limiter options (upwind and Vanleer) versus the corresponding analytical solutions (solid line). The left and right side plots show the results with 100 & 200 cells, respectively.

initially stationary aluminum half space (slab), see Fig. 5. The rear (left) surface of the flyer is a free space (vacuum). In our computation, we have used a 1-D problem domain of 50 mm wide (along the dashed line shown in Fig. 5), in which the aluminum half space initially occupies from 12.5 mm to 45 mm and the aluminum flyer occupies from 8.5 mm to 12.5 mm, i.e. at t = 0, the flyer and the target are in contact. The rest of the domain is filled with vacuum. The vacuum material is modeled with an extremely low density ideal gas with ρ0 = 10−12 kg/m3 and γ = 1.4. The purpose of this material is to fill-up the initial and generated vacuum space without effecting the system dynamics. Mie-Gruneisen equation of state is used with parameters provided for alu21

minum in [22,35,51]. A flux limited HEMP artificial viscosity [35] is used with a typical limiter function given in [34]. The numerical results are compared with the solutions provided in [22], where a high order Godunov method is used. In order to be consistent with the results given in [22], the strength of the material is set equal to zero. Therefore, this test case is basically to study the performance of the proposed predictor-corrector algorithm. Flyer Target v = 0 v>0 1−D simulation line

Vacuum

Fig. 5. Impact problem: a flyer with velocity > 0, impacts a target half-space (initially at rest). Rest of the space is a vacuum. The simulation is carried out in one-dimension along the dashed line shown in the figure. We start our simulation when the flyer just hit the target, i.e. at t = 0, the flyer and the target are in contact.

The propagation of shock and rarefaction waves at different times are shown in Fig. 6. Initially, left and right traveling shocks propagate outward from the point of impact of the flyer with the half space. When the left traveling shock wave reaches the free surface of the flyer, a right traveling rarefaction is created. This rarefaction wave has attenuated the right traveling shock and reaches the shock front at ∼ 5.0 µs after the impact, see Fig. 6. In the following, we examine these results quantitatively.

Pressure (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

20 18 16 14 12 10 8 6 4 2 0 -2

t=1.5 µs t=2.5 µs t=4.0 µs t=5.0 µs

0

5

10

15

20 25 30 35 Distance (mm)

40

45

50

Fig. 6. The propagation of shock and rarefaction waves at different times for Wilkins’s impact problem described in Sec. 7.3. A mesh-size of 0.25 mm is used.

The pressure and the normalized density (normalized to ρ0 = 2780 kg/m3 ) profiles at 2.5 µs are shown in Fig. 7. Our results are in good agreement with the corresponding results given in Ref. [22], especially when the mesh22

20 18 16 14 12 10 8 6 4 2 0 -2

1.2

Godunov ∆x = 1 mm ∆x = 0.5 mm ∆x = 0.25 mm ∆x = 0.125 mm ∆x = 0.0625 mm

Normalized density

Pressure (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

1.15 1.1 1.05 1

10

15

20 25 Distance (mm)

18

30

20

22

24

26

28

30

32

Distance (mm)

Fig. 7. A zoomed-in plot of pressure (left) and density (right) profiles at 2.5 µs after impact for the Wilkins’s problem described in Sec. 7.3. The results of Godunov method [22] is obtained with a mesh-size of 0.1 mm. Table 1 The L1 norm error in the density and convergence rate for Wilkins’s problem described in Sec. 7.3. No. of cells ∆x (mm) L1 norm Rate 50

1.0

7.288e-2

-

100

0.5

1.952e-2

1.9

200

0.25

6.161e-3

1.66

400

0.125

1.811e-3

1.77

800

0.0625

5.201e-4

1.8

resolution approaches to 0.1 mm (a mesh-size used in [22]). However, with a coarse-mesh simulation, small amount of density fluctuations behind the rarefaction fan are observed in our simulation. Consequently, small amplitude oscillations at similar points are present in the pressure profile. Note that, with increased spatial resolution, the numerical solutions are free from such oscillations, see Fig. 7. Next, we have calculated theL1 norm of the  error in the density. The L1 norm c e error is calculated as L1 = ΣN |f − f | /N ; where fic,e are the computed i=1 i i and exact solutions respectively, and N is the total number of cells. The L1 norm error is calculated by comparing the computed solution with the solution obtained with a mesh-resolution of 1500 cells for the same problem. Also, the rate of convergence pξ = log(L12h /L1h )/log(2) is derived from the L1 error values obtained for different cases, where L1h and L12h are the L1 norm errors with grid-size ‘h’ and ‘2h’ respectively. The results are summarized in Table 1. It is clear from the table, for this 1D test problem, the proposed predictorcorrector algorithm yielded an average order convergence rate ∼ 1.78. It is worth to see the spatial variation of numerical error ξi = |fic − fie |. The 23

variation of ξi at 2.5 µs for a typical case with ∆x = 0.25 mm is shown in fig. 8. The numerical error is comparatively higher at the rear end of the flyer. Also, the ξi calculated at the impact plane where the multi-material treatment occurs is less than the error estimated at the shock front located inside the target. A similar trend for numerical error is seen in Ref. [22] and also for the cases with different ∆x in our simulation. 0.25

∆x = 0.25 mm

0.2 Numerical error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.15 0.1 0.05 0 0

5

10 15 20 25 30 35 40 45 50 Distance (mm)

Fig. 8. Spatial variation of numerical error ξi in density at 2.5 µs after the impact with ∆x = 0.25 mm for Wilkins’s problem described in Sec. 7.3.

Next, we have studied the effect of artificial viscosity on the accuracy of numerical solution. The pressure profiles obtained by using various flux limited artificial viscosities at 2.5 µs (a mesh-size of 0.25 mm is used) are shown in Fig. 9. Note that the multi-material theoretical formulation [28,22] followed here (without material strength model) use an artificial viscosity term to resolve the shock front. In Ref. [28], only a quadratic viscosity term is used. However, we have observed that, for impact driven shock problems, both the linear and quadratic viscosity terms are necessary in our algorithm to yield an oscillation (post shock oscillation) free solution especially when the material strength model is enabled. Inspecting Fig. 9, it is clear that the Wilkins’s strain rate (along the shock direction) dependent viscosity yields slightly better results for this problem. 7.4 1-D Wilkins’s impact problem with elastic precursor Next, we have examined Wilkins’s flying plate problem [51] with material strength model included. This involves the propagation of an elastic precursor. The problem involves a 5 mm thick aluminum flyer plate impacting an initially stationary half-space, see Fig. 5. A one-dimensional calculation is used. We have used Wilkins’s rate model [35] to model EOS of the aluminum. The shear modulus (µ0 ) and yield strength (σY ) used were µ0 = 27.6 GPa and σY = 0.3 GPa, respectively. The work hardening is neglected assuming the flow is perfectly plastic (a rate dependent material strength model is not included in the present algorithm). The Von Mises yield criterion is used to model 24

30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

quadratic viscosity alone linear and quadratic HEMP viscosity

25 20 15 10 5 0 -5 10

15

20

25

30

35

Fig. 9. Spatial variation of pressure at 2.5 µs after the impact with ∆x = 0.25 mm for Wilkins’s problem described in Sec. 7.3 obtained with various flux limited viscosity options.

plasticity. A Cartesian grid is used with 500 cells along 50 mm computational domain in which the aluminum half space initially occupies from 12.5 mm to 45 mm and the aluminum flyer occupies from 7.5 mm to 12.5 mm. The rest of the domain is filled with vacuum. The vacuum material is modeled with an extremely low density ideal gas with ρ0 = 10−12 kg/m3 and γ = 1.4. Two cases are considered: one with an impact velocity of 800 m/s and the second one with an impact velocity of 2.5 km/s. Fig. 10 shows the propagation of longitudinal stress σzz for the case with an impact velocity 800 m/s. For this case the plastic shock is behind the leading elastic precursor. When the shock wave reaches the left free surface of the flyer right traveling elastic and trailing plastic rarefaction waves begin to overtake the initial right traveling shocks, see Fig. 10. For the case with impact velocity 2.5 km/s, the shock stress is above the elastic limit. Therefore, only plastic shocks are formed. Once the shock wave reaches the left free surface, a right traveling elastic rarefaction is formed followed by the plastic wave. Fig. 11 shows the propagation of longitudinal stress σzz for this case. The results discussed above are in good agreement with those described in Ref. [51] and therefore validate the correctness of the implementation of the PC integration scheme with material strength model.

7.5 1-D impact problem with elastic precursor - comparison with analytical solution Next, we present the results of impact on elasto-plastic material with material strength model to verify the ability of the algorithm to correctly model the discontinuities in an elastic-plastic medium. For materials described with MieGruneisen EOS, it is possible to determine the exact solution to the uniaxial strain equations under shock loading by a velocity discontinuity at impact 25

0

t=0.5 µs

-1

-1

-2

-2

Stress σzz (GPa)

Stress σzz (GPa)

0

-3 -4 -5 -6 -7 -10

0

-3 -4 -5

0

10 20 30 Distance (mm)

40

-7 -10

50

0

t=3.0 µs

-1

-1

-2

-2

-3 -4 -5 -6 -7 -10

t=1.0 µs

-6

Stress σzz (GPa)

Stress σzz (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0

10 20 30 Distance (mm)

40

50

10 20 30 Distance (mm)

40

50

t=5.0 µs

-3 -4 -5 -6

0

10 20 30 Distance (mm)

40

-7 -10

50

0

Fig. 10. Stress (σzz ) profile at different times for impact problem with impact velocity 800 m/s.

face, see Fig. 12 for the description. The analytical solution for this problem is described in Ref. [52]. The computational domain extends from x = 0 to x = 2 m. The impact at x = 1 is modeled by a velocity change jump in initial velocity profile at x = 1 m (40 m/s). The impact at a velocity of 40 m/s at a stationary material results in a 20 m/s final particle velocity. In order to be consistent with the analytical solutions provided in Ref. [52], the Gruneisen EOS parameters for copper is taken from Ref. [52]. The Gruneisen equation of state parameters are Γ = 2, c0 = 3.94 km/s, s = 1.49, with initial density ρ0 = 8.93 g/cc, shear modulus G = 45 GPa, and yield stress Y = 90 MPa. The mesh-size used is varied from 1.25 mm to 5 mm. A flux limited artificial viscosity with both quadratic and linear term is used. The plots in Fig. 13 show the density, pressure, stress (Sxx ) and velocity profiles obtained at 170 µs after the impact. Only the zoomed-in results in the target domain (1 to 2 m) are given in the figures, since the solutions are symmetric in the other side (0 to 1 m) of the impact plane located at 1 m. The analytical solutions given in Ref. [52] are also shown in those figures. It is clear from those figures that the numerical results obtained using our algorithm show reasonable agreement with the analytical solutions of various physical quantities. The numerical diffusion experienced near the elastic precursor can be minimized with a high resolution simulation (mesh-size of the order of 0.1 26

0

0

t=0.9 µs

-10 -15 -20 -25 -10

0

-10 -15 -20

0

10 20 30 Distance (mm)

40

-25 -10

50

0

t=2.7 µs

0

10 20 30 Distance (mm)

40

50

10 20 30 Distance (mm)

40

50

t=4.0 µs

-5 Stress σzz (GPa)

-5 -10 -15 -20 -25 -10

t=1.8 µs

-5 Stress σzz (GPa)

Stress σzz (GPa)

-5

Stress σzz (GPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-10 -15 -20

0

10 20 30 Distance (mm)

40

-25 -10

50

0

Fig. 11. Stress (σzz ) profile at different times for impact problem with impact velocity 2.5 km/s.

v>0

v=0

1D simulation line

0

x

2x

Fig. 12. Geometry setup for impact problem generated by a velocity discontinuity. One dimensional calculations are performed along the solid line parallel to x-axis.

mm). Our simulation results, however, are in reasonable agreement with the Eulerian results given in Ref. [52]. We note that, a higher-order cell-centered Lagrangian simulation of the same problem, described in Ref. [14], yields lesser diffusion near the elastic precursor. The final particle velocity obtained is 20 m/s which is in good agreement with analytical value (half the impact velocity). A slight overheating near impactorslab interface is experienced, which was also reported in [52]. The increase in internal energy is compensated by a decrease in density leading to a correct estimate for pressure. This overheating errors in energy and density can be 27

0.7

8970

3

0.5

Analytical ∆x = 5 mm ∆x = 2.5 mm ∆x = 1.25 mm

8980 Density (kg/m )

Pressure (GPa)

8990

Analytical ∆x = 5 mm ∆x = 2.5 mm ∆x = 1.25 mm

0.6

0.4 0.3 0.2 0.1

8960 8950 8940 8930 8920

0

8910 1.6

0

1.65

1.7 1.75 1.8 Distance (m)

1.85

1.9

1.6

1.65

20

Analytical ∆x = 5 mm ∆x = 2.5 mm ∆x = 1.25 mm

-20 -30 -40

1.7 1.75 1.8 Distance (m)

1.85

1.9

Analytical ∆x = 5 mm ∆x = 2.5 mm ∆x = 1.25 mm

15 Velocity (m/s)

-10 Sxx (MPa)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

10

5 -50 -60

0 1.4

1.5

1.6 1.7 1.8 Distance (m)

1.9

2

1.6

1.65

1.7 1.75 1.8 Distance (m)

1.85

1.9

Fig. 13. A zoomed-in pressure, density, stress (Sxx ) and velocity profiles for the impact problem [52] at 170 µs after the impact.

fixed by using an isobaric fix as proposed by Fedkiw et al. [53] or by adding a heat diffusion term at these points as proposed in [54]. Present formulation lack these models. However, the pressure, velocity, density and stress match the analytical results, despite a small error in internal energy near impact plane due to overheating. Similarly, the axial stress solution shows some very small oscillations near the impact plane. These oscillations are found to be very small and does not grow in time. Also, we have noted that these oscillations were able to suppress by increasing the value of the coefficient of linear artificial viscosity. 7.6 2-D axisymmetric Taylor impact problem with material strength (Cylindrical geometry) The problem of an axisymmetric impact of a long cylindrical rod on a rigid surface (Taylor impact problem [35,52,55]) is often used as a benchmark problem to test code capability to deal two-dimensional impact and deformation in cylindrical geometry. The problem setup is shown in Fig. 14. The cylindrical rod is assumed to be made of copper and initially has a radius of 3.2 mm and a length of 32.4 mm. The impact velocity is 227 m/s. The 28

Computational domain

111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111

Rigid surface V −−> Cylindrical rod axis

10 9 8 7

v=0

Vacuum

R (mm)

r

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

6 5

Rigid surface

4

z

3 2 Cylindrical rod 1 0 0

5

10

15

20

25

30

35

Z (mm)

Fig. 14. Geometry setup for the Taylor impact problem. The right side plot shows the initial computational mesh. The grids shown in green color indicate the cylindrical rod and the background grids (red) are filled with a low density gas. A fixed boundary condition is used at the right side of the mesh.

deformation is assumed to be symmetric about the z-axis. Therefore, the computations are performed only on the symmetric half of the problem domain. A uniform mesh covering a 2D domain of 35.2 mm along the z-direction and 10 mm along the radial direction is used. A mesh-size of 0.4 mm is used in each direction. The cylindrical rod initially fills up a region of 2.8 mm to 35.2 mm along the axial direction with a radius of 3.2 mm. The rest of the domain is filled with vacuum. The vacuum material is modeled with an extremely low density ideal gas with ρ0 = 10−12 kg/m3 and γ = 1.4. To model the rigid surface (right axial end of the domain) a fixed boundary condition is used. The material properties chosen were: Youngs modulus = 117 GPa, Poissons ratio = 0.35, yield stress = 0.4 GPa and density of copper = 8930 kg/m3 . Flux limited HEMP artificial viscosity is used. The calculations were carried out up to a time of 80 µs and the results are compared with those given in Refs. [52,55]. Fig. 15, shows the material deformation at different times (8, 24, 56 and 80µs) after the impact. The copper rod come to rest at a time ∼ 80.5 µs. This value for time of rest agrees well with the time of rest predicted (80 µs) in Ref. [55]. Thus our algorithm predicts fairly accurately the rate of conversion of initial kinetic energy to plastic work. We have compared the final rod length to that given in Refs. [52,55]. The final rod length obtained from our simulation is 21.6 mm, which is in reasonable agreement with the value 21.4 mm given in [52,55]. Similarly, we have compared the maximum width of the mushroom head formed after impact. The values obtained by the various authors fall in the range of 6.8 to 7.24 mm. The width of the mushroom head obtained from our simulation was 6.82 mm. Also, the shape of the rod deformation is in good qualitative agreement with the figures shown in Refs. [52,55]. 29

t = 8.0 µs

t = 24.0 µs

5 radial length (mm)

radial length (mm)

4

3

2

1

4 3 2 1

0

0 0

5

10

15

20

25

30

35

5

10

15

axial length (mm) 6

20

25

30

35

axial length (mm)

t = 56.0 µs

t = 80.0 µs 6 radial length (mm)

5 radial length (mm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4 3 2 1

5 4 3 2 1

0

0 10

15

20

25

30

35

10

axial length (mm)

15

20

25

30

35

axial length (mm)

Fig. 15. The material deformation at different times obtained for axisymmetric Taylor impact problem. Only one half of the problem domain is simulated because of symmetry about the z-axis.

7.7 2-D impact and penetration problem with material strength (Cylindrical geometry) Next, we examine the case of 2D-axisymmetric high velocity impact of two deformable materials. The initial conditions are depicted in Fig. 16, where a sphere of radius, Rc is traveling with velocity vc impacts on a thin slab. A detailed numerical analysis of the same problem have been reported in Ref. [56] by Vishal et. al using their 2-D and 3-D SPH codes. The materials and their dimensions are taken from Ref. [57]. Two different cases from [57] are considered here. The initial parameters for both cases are summarized in Table 2. The material strength is modeled using parameters given in [35] and a Mie-Gruneisen EOS [35] is used. The computational domain consists of 600 × 240 cells in z and r directions with a mesh-size of 0.1667 mm in each direction. The background material (vacuum) is modeled with an extremely low density ideal gas with ρ0 = 10−12 kg/m3 and γ = 1.4. Flux limited artificial viscosity (quadratic plus linear) is used. The sequence of material deformation (for case 1 in Table 2) at different times are depicted in Fig. 17. The Al sphere completely penetrates through the thin Al slab. Both the impactor and the target undergo severe deformation. The beginning of material breakups are clearly visible in those figures. The nature of debris cloud expansion and the deformation of both the target and the 30

50

Extremely low density gas Sphere (projectile) Thin slab (Target)

40

Background medium

V

r

R

R (mm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Axis

30

20

z Computational domain

Impactor

10

Target

0 0

t

10

20

30

40

50 60 Z (mm)

70

80

90

100

Fig. 16. Geometry for sphere-slab impact problem at t = 0. The right side plot shows the initial computational mesh. Due to the symmetry of the problem, only one half of the system in cylindrical geometry is simulated. The cells representing the projectile (sphere) and its interfaces in mixed cells are plotted in red color. Similarly, the target (thin slab) is plotted in blue color. The empty region (filled with extremely low density gas) and its interfaces in mixed cells are plotted in green color. Table 2 Parameters used for impact problem. In the table, Y and µ are the yield and shear modulus respectively (in GPa). The distances Rc and tslab are given in cm. The velocity vc is given in km/s. Case

Rc

vc

tslab

Impactor

Slab

YAl , YCu

µAl , µCu

1

0.5

6.18

0.4

Al

Al

0.29, 0.12

27.6, 47.7

2

0.5

5.75

0.15

Al

Cu

0.29, 0.12

27.6, 47.7

impactor are in good qualitative agreement with the simulation results given in Refs. [56,57]. The material interfaces calculated with and without material strength at 10 µs after the impact are shown in Fig. 18. As expected the simulation overpredicts the crater diameter by ∼22 % when the material strength model is disabled. We have compared few experimental data mentioned in [57] with our corresponding simulation results for the validation of the code. Also, we have included the simulation results obtained in Ref. [56] for a comparison with 3-D SPH simulations. The Comparison of simulated crater diameter (dcr ) and length to width (l/w) ratio of the debris cloud against experimental results are summarized in Table 3 and Table 4. In order to have a consistency with the experimental values given in Ref. [57] the comparisons are made at 20 µs after the impact. The simulation results show reasonable agreement with experimentally observed values. However, the l/w ratio obtained in our simulation is higher than the simulation results reported in [57]. The reason for this difference (also mentioned in [57]) may be due to the planar symmetry used in their simulations. Inspecting the values given in Table 4, it is clear 31

0.06 (a) t=0

0.05

Radial distance (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(b) t = 6 µs

(c) t = 10 µs

0.04 0.03 0.02 0.01 0 0

0.02

0.04 0.06 0.08 Axial length (m)

0.1

Fig. 17. Material interfaces at different times after the impact for sphere impact problem. Different sub-plots Fig. (b) and Fig. (c) in the figure are conveniently shifted along the axis, 2 cm and 5 cm respectively. Only one half of the actual system is simulated.

Table 3 Comparison of crater diameter, dcr , with experimental [57] and simulation values given in [57,56] with our results for impact problem at 20 µs. The values are given in cm. Set Expt. [57] Sim. [57] Our Sim. 2D-Sim. [56] 3D-Sim. [56] 1

2.75-3.45

3.5

2.7

3.2−3.3

2.4−3.2

2

2.12

2.3

2.0

2.1−2.5

2.0−2.3

that our simulation results are close to the 3-D SPH simulation results. Similarly, the reason for the difference in our simulation results with experimental values may be due to the assumption of a constant yield value for materials in our simulations. Perhaps, a rate dependent material strength model may give better match with experimental results. However, this test demonstrates that the present algorithm is capable of handling large material deformations with a fairly accurate prediction of overall dynamics. The maximum energy conservation error for this simulation was ∼2.2 % and the linear momentum conservation error was negligible. This increase in the energy conservation error (compared to other problems) is due to the material breakups occurred in the final stage of the simulation (during debris cloud expansion) and the associated loss of small amount of mass. This can be easily minimized by using a finer mesh. 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.05 with strength

without strength

0.03

0.01 0 0

0.04

0.08

0.12

Fig. 18. Material interfaces at 10µs after impact for sphere impact problem. The plot obtained without using material strength model is shifted ∼ 7 cm along the axial direction. Table 4 Comparison of debris length to width ratio, l/w, with experimental [57] and simulation values given in [57,56] with our results for impact problem. Set

Expt. [57]

Sim. [57]

Our Sim.

2D-Sim. [56]

3D-Sim. [56]

1

1.39

1.33

1.58

1.31−1.44

1.53−1.9

2

1.39

1.11

1.51

1.41−1.51

1.7−2.3

7.8 Penetration of metal targets by ogive nose metal projectile

Finally, we have studied the penetration of metal targets by ogive nose metal projectile. The initial configuration of the problem is depicted in Fig. 19. The authors in Ref. [58] have performed a series depth-of-penetration experiments using an ogive nose steel projectile and aluminum targets. The experiments were conducted for a range of velocities varies form 0.5 to 3 km/s. We have compared the depth-of-penetration for a particular case reported in Ref. [58] with our simulation result. The dimensions of the projectile and target are taken from Ref. [58]. The projectile diameter and length (including ogive nose) are 7.11 mm and 71.1 mm respectively. A steel projectile is used. The projectile nose is ogive shaped with an axial length of 11.8 mm. The aluminum target is ∼ 250 mm in diameter. The calculations were performed with a zoning of 550 × 75 cells. A variable meshing is used. The r-direction meshing is started at r = 0 and extended to a radius of 8.88 mm using a uniform mesh size of 0.355 mm. This was followed by a non-uniform meshing consists of 50 cells and extending to a radius of 15.0 cm. Similarly, in the z-direction, a uniform meshing consists of 500 cells (starting at z = 0) with a mesh-size of 0.355 mm 33

0.1

Extremely low density gas Projectile -- ogive Target

target

r

0.08

0.06 R (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.04

0.02

impactor 0 0

axis

0.02

z

0.04

0.06

0.08

0.1

Z (m)

Fig. 19. Penetration of metal targets by an ogive nose metal projectile: The left and right side plots show the initial geometry and the corresponding computational mesh (zoomed-in). Due to the symmetry of the problem, only one half of the system in r-z geometry is simulated. The cells representing the projectile (ogive) and its interfaces in mixed cells are plotted in red color. Similarly, the target is plotted in blue color. The empty region (filled with extremely low density gas) and its interfaces in mixed cells are plotted in green color.

followed by a non-uniform meshing consists of 50 cells extending to an axial length of 30.0 cm is used. The material parameters are taken from Ref. [58]. Also, both the projectile and target were modeled with Mie-Gruneisen EOS. The background material (vacuum) is modeled with an extremely low density ideal gas with ρ0 = 10−12 kg/m3 and γ = 1.4. The initial projectile velocity was 570 m/s. In our simulation the ogive shape is approximated with an outer surface defined by Eq. (54).

f (x, y) =

q





[d (C 2 + 0.25)]2 − x2 − d C 2 − 0.25 − y

(54)

Here, C = L/d; where L and d are the length and diameter of the ogive nose. Fig. 20 shows the shape of the projectile and target at different times. The penetration has stopped at about 390 µs. Note that the projectile shows no permanent deformation (globally), which is consistent with the experiment [58]. However, a small amount of local material erosion is observed (towards the end of penetration, see Fig. 20) for both the projectile and target, which was not reported in [58]. A contact algorithm such as boundary-layer or free-slip algorithm [59] may mitigate this problem. We have compared the depth-ofpenetration with the experimental value (55 mm for 570 m/s projectile velocity) given in Ref. [58]. The penetration depth obtained from the simulation was 52 mm, shows reasonable agreement despite small erosion happened to both the projectile and target in our simulation. The energy conservation error was less than 1.6%. 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.04

0.04

0.03

0.03

0.02

0.02

0.01

0.01

0

0

-0.01

-0.01

-0.02

-0.02

-0.03

-0.03

-0.04 0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.04

0.1

0.06

0.08

0.1

0.12

Fig. 20. A zoomed-in plot shows the shape of the projectile (red), target (blue) and the depth-of-penetration at different times for the penetration problem described in Sec. 7.8. The second figure is at 390 µs. The lower symmetric part (the negative y-axis) is duplicated from the upper symmetric part.

8

Summary and conclusions

We have presented a ’simple’ multi-material Eulerian algorithm with material strength model. The material interfaces are tracked by using a classical volume-of-fluid (VOF) method. The algorithm consists of an Eulerian (Lagrangian plus remap onto the initial mesh at each time-step) formulation with a new predictor-corrector (PC) integration scheme. The multi-material formulation used in this work avoids commonly used iteration based mixed EOS evaluation or pressure equilibration schemes in a mixed cell. The mean stress in a mixed cell is determined by taking a ‘simple’ volume weighted average of individual material stress components. These individual stress components are updated with an assumption of common strain rates to all materials present in a mixed cell. The advection of stress components for individual materials and the advection of many other solution and history variables are performed by using a combination of VOF method and a second order accurate monotonic upwind scheme (MUSCL). One and two dimensional examples show that the algorithm reasonably approximates (both qualitatively and quantitatively) the physical solutions for the shock and rarefaction waves arises from material impacts. That is the approximation of equivalent stress and pressure in a multi-material mixed cell using a ‘simple’ volume weighted average works well to reasonably predict the physical solutions of various impact and penetration problems. Also, it is demonstrated that the algorithm is capable of handling multi-material impact and penetration problems involving large material deformations. Further, for these problems, the algorithm yields results that are in good agreement with the reported experimental values. The algorithm exhibits an overall convergence rate of ∼1.8 for one dimensional problems. Moreover, for spherically symmetric problems in 2D, we have demonstrated that the algorithm preserves 35

the radial symmetry of the solution reasonably well. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In ongoing work several extensions of the formulations presented in this work are in progress. They are: rate-dependent material strength model, void and crack growth model and material breakup model. Similarly, the implementation of a contact algorithm such as boundary-layer and free-slip algorithm [59] to predict rigid-body response of the projectile is in progress.

9

Acknowledgments

One of the author (CDS) is grateful to Vishal Mehra (CAD, Bhabha Atomic Research Centre, Visakhapatnam, India) for many fruitful and stimulating discussions. We thank our Referees and Editor for their valuable comments and helpful suggestions.

References [1] S. E. Hieber and P. Koumoutsakos, A Lagrangian particle level set method, J. Comp. Phys., vol 210, pp 342-367, (2005) [2] Hirt C. W and Nicholas B. D., Volume of Fluid (VOF) method for the dynamics of free boundaries, J. Comp. Phys., vol 39, pp 201-225, (1981) [3] W. J. Rider and D. B. Kothe, Reconstructing volume tracking, J. Comp. Phys., vol 141, pp 112-152, (1998) [4] D. L. Youngs, Time dependent multi-material flow with large fluid distortions, Numerical methods for Fluid dynamics, K. W. Morton and M. J. Baines (eds), Academic Press, pp 273-285, (1982) [5] S. P. Schofield, R. V. Garimella, M. M. Francois and R. Loubere, Material order independent interface reconstruction using power diagrams, Int. J. for Numer. Meth. in Fluids, vol 56(6), pp. 643-659, (2008) [6] S. P. Schofield, R. V. Garimella, M. M. Francois and R. Loubere, A second order accurate material-order-independent interface reconstruction technique for multi-material flow simulations J. Comp. Phys., vol 228, pp 731-745, (2009) [7] V. Dyadechko and M. Shashkov, Reconstruction of multi-material interfaces from moment data, J. Comp. Phys., vol 227, pp 5361-5384, (2008) [8] H. T. Ahn and M. Shashkov, Multi-material interface reconstruction on generalized polyhedral meshes, J. Comp. Phys., vol 226, pp 2096-2132, (2007) [9] P. H. Maire and B. Nkonga, Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics, J. Comp. Phys., vol 228, pp 799-821, (2009)

36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[10] G. Carre, S. Del Pino, B. Despres and E. Labourasse, A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension, J. Comp. Phys., vol 228, pp 5160-5183, (2009) [11] P. H. Maire, A higher order cell-centered Lagrangian scheme for twodimensional compressible fluid flows on unstructured meshes, J. Comp. Phys., vol 228, pp 2391-2425, (2009) [12] P. H. Maire, A higher order cell-centered Lagrangian scheme for compressible fluid flows in two-dimensional cylindrical geometry, J. Comp. Phys., vol 228, pp 6882-6915, (2009) [13] S. K. Sambasivan, M. J. Shashkov and D. Burton, A cell-centered Lagrangian finite volume approach for computing elastic-plastic response of solids in cylindrical axisymmetric geometries, J. Comp. Phys., vol 237, pp 251-288, (2013) [14] P. H. Maire, R. Abgrall, J. Breil, R. Loubere and B. Rebourcet, A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two-dimensional unstructured grids, J. Comp. Phys., vol 235, pp 626-665, (2013) [15] D. E. Burton, T. C. Carney, N. R. Morgan, S. K. Sambasivan and M. J. Shashkov, A cell-centered Lagrangian Godunov-like method for solid dynamics, Computers & Fluids, vol 83, pp 33-47, (2013) [16] S. Galera, P. H. Maire and J. Breil, A two-dimensional unstructured cellcentered multi-material ALE scheme using VOF interface reconstruction, J. Comp. Phys., vol 229, pp 5755-5787, (2010) [17] S. Galera, J. Breil and P. H. Maire, A 2D unstructured multi-material cellcentered arbitrary Lagrangian-Eulerian (CCALE) scheme using MOF interface reconstruction, Computers & Fluids, vol 46, pp 237-244, (2011) [18] J. Breil, S. Galera and P. H. Maire, Multi-material ALE computation in inertial confinement fusion code CHIC, Computers & Fluids, vol 46, pp 161-167, (2011) [19] P. H. Maire, A higher order one-step sub-cell force-based discretization for cellcentered Lagrangian hydrodynamics on polygonal grids, Computers & Fluids, vol 46, pp 341-347, (2011) [20] P. H. Maire and J. Breil, A second-order cell-centered Lagrangian scheme for two-dimensional compressible flow problems, Int. J. for Num. Meth. in Fluids, vol 56, pp 1417-1423, (2008) [21] D. J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Comput. Meth. Appl. Mech. Eng., vol 99, pp 235-394, (1992) [22] G. H. Miller and E. G. Puckett, High-order Godunov method for multiple condensed phases, J. Comp. Phys., vol 128, pp 134-164, (1996) [23] D. J. Benson, A mixture theory for contact in multi-material Eulerian formulations, Comput. Meth. Appl. Mech. Eng., vol 140, pp 59-86, (1997)

37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[24] M. Shashkov, Closure models for multimaterial cells in arbitrary LagrangianEulerian hydrocodes, Int. J. for Num. Meth. in Fluids, vol 56, pp 1497-1504, (2008) [25] A. Barlow, A new Lagrangian scheme for multimaterial cells, Proceedings of European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS Computational Fluid Dynamics Conference, Swansea, Wales, U.K., (2001) [26] Y. V. Yanilkin, E. A. Goncharov, V. Yu. Kolobyanin, V. V. Sadchikov, J. R. Kamm, M. J. Shashkov and W. J. Rider, Multi-material pressure relaxation methods for Lagrangian hydrodynamics, Computers & Fluids, vol 83, pp 137143, (2013) [27] M. M. Francois, M. J. Shashkov, T. O. Masser and E. D. Dendy, A comparative study of multimaterial Lagrangian and Eulerian methods with pressure relaxation, Computers & Fluids, vol 83, pp 126-136, (2013) [28] D. de Niem, E. Kuhrt and U. Motschmann, A volume-of-fluid method for simulation of compressible axisymmetric multi material flow, Comp. Phys. Comm., vol 176, pp 170-190, (2007) [29] C. D. Sijoy and Shashank Chaturvedi, Volume of fluid algorithm with different modified dynamic material-ordering methods and their comparisons, J. Comp. Phys.,vol 229, pp 3848-3863, (2010) [30] C. D. Sijoy and Shashank Chaturvedi, An Eulerian MHD model for the analysis of magnetic flux compression by expanding diamagnetic fusion plasma sphere, Fusion Engineering and Design, vol 87, pp 104-117, (2012) [31] C. D. Sijoy, V. Mehra, V. Mishra and S. Chaturvedi, A VOF based multimaterial method to study impact and penetration problems, Journal of Physics: Conference Series, vol 377, pp 012106, (2012) [32] B. Van Leer, Toward the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme, J. Comput. Phys.,14, pp 361-370, (1974) [33] D. J. Benson, A new two-dimensional flux-limited shock viscosity for impact calculations, Comput. Meth. Appl. Mech. Eng., vol 93, pp 39-95, (1991) [34] E. J. Caramana, M. J. Shashkov and P. P. Whalen, Formulations of artificial viscosity for multi-dimensional shock wave computations, J. Comp. Phys., vol 144, pp 70-97, (1998) [35] M. L. Wilkins, Computer simulation of dynamic phenomena, Scientific Computation, Springer-Verlag, Berlin, New York, (1998) [36] P. M. Suquet, Elements of homogenization for inelastic solid mechanics, Homogenization Techniques for Composite Media, Lecture Notes in Physics, Springer-Verlag, vol 272, pp 194-280, (1985)

38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[37] V. V. Shuvalov, Multimaterial hydrodynamic code SOVA for interfacial flows, application to thermal layer effect, Shockwaves, vol 9, pp 381, (1999) [38] V. V. Shuvalov, N. A. Artem’eva and I. B. Kosarev, 3D hydrodynamic code SOVA for multimaterial flows, application to Shoemaker-Levy 9 comet impact problems, Int. J. Impact Eng., vol 23, pp 847, (1999) [39] D. J. Benson, Momentum advection on staggered mesh, J. Comp. Phys., vol 100, pp 143-162, (1992) [40] J. VonNeumann and R. D. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys., vol 21, pp 232, (1950) [41] R. J. Leveque, Finite volume methods for hyperbolic problems, Cambridge University Press, 2006 [42] D. J. Benson, Stable times step estimation for multi-material Eulerian hydrocodes, Comput. Meth. Appl. Mech. Eng., vol 167, pp 191-205, (1998) [43] D. Gueyffier, J. Li, R. Scardovelli and S. Zaleski, Volume-of-Fluid interface tracking with smoothed surface stress methods for three-dimensional flow, J. Comp. Phys., vol 152, pp 423, (1999) [44] J. E. Pilliod and E. G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, J. Comp. Phys., vol 199, pp 465, (2004) [45] E. Aulisa, S. Manservisi, R. Scardovelli and S. Zaleski, A geometrical areapreserving volume-of-fluid advection method, J. Comp. Phys., vol 192, pp 355, (2003) [46] D. J. Benson, Volume of fluid interface reconstruction methods for multimaterial problems, Appl. Mech. Rev., vol 55 (2), pp 151, (2002) [47] E. G. Puckett, A Volume-of-Fluid interface tracking algorithm with applications to computing shock wave refraction, Proceedings of the Fourth International Symposium on Computational Fluid dynamics, H. Dwyer (Ed.), Davis, CA, pp 933-938, (1991) [48] R. P. Brent, An Algorithm with Guaranteed Convergence for Finding a Zero of a Function, Algorithms for Minimization without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, ISBN 0-13-022335-2, (1973) [49] S. Mosso and S. Clancy, A geometrically derived priority system for Youngs interface reconstruction, Technical Report LA-CP-95-81, LANL, (1994) [50] D. J. Benson, Eulerian finite element methods for the micromechanics of the hetrogeneous materials: Dynamic prioritization of material interfaces, Comput. Meth. Appl. Mech. Eng., vol 151, pp 343-360, (1998) [51] M. L. Wilkins, in Methods in computational physics, Edited by B. Alder, S. Fernbach and M. Rotenberg, vol 3, pp 211, Academic press, New York, (1964)

39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[52] H. S. Udaykumar, L. Tran, D. M. Belk and K. J. Vanden, An Eulerian method for computation of multimaterial impact with ENO shock-capturing and sharp interfaces, J. Comp. Phys., vol 186, pp 136-177, (2003) [53] R. P. Fedkiw, A. Marquina, B. Merriman, An isobaric fix for the overheating problem in multi-material compressible flows, J. Comp. Phys., vol 148, pp 545578, (1999) [54] R. Donat, A. Marquina, Capturing shock reflections: an improved flux formula, J. Comp. Phys., vol 125, pp 42-58, (1996) [55] G. T. Camacho, Adaptive Lagrangian modelling of ballistic penetration of metallic targets, Comput. Methods. Appl. Mech. Engrg., vol 142, pp 269-301, (1997) [56] Vishal Mehra and Shashank Chaturvedi, High velocity impact of metal sphere on thin metallic plates: A comparative smooth particle hydrodynamics study, J. Comp. Phys., vol 212, pp 318-337, (2006) [57] S. Hiermaier, D. Konke, A. J. Stilp and K. Thoma, Computational simulation of the hypervelocity impact of Al-spheres on thin plates of different materials, Int. J. Impact Engng., vol 20, pp 363-374, (1997) [58] A.J. Piekutowski, M.J. Forrestal, K.L.Poormon and T.L. Warren, Penetration of 6061-T6511 Aluminum targets by ogive-nose steel projectiles with striking velocities between 0.5 and 3.0 km/s, Int. J. Impact Engng., vol 23, pp 723-724, (1999) [59] J. M. McGlaun, S. L. Thompson and M. G. Elrick, CTH: A three dimensional shock wave physics code, Int. J. Impact Eng., vol 10, pp 351, (1990)

40