Numerical simulations of the damage evolution for plastic-bonded explosives subjected to complex stress states

Mechanics of Materials 139 (2019) 103179

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Research paper

Numerical simulations of the damage evolution for plastic-bonded explosives subjected to complex stress states

T

Ming Liua,b, Xicheng Huangb, , Yanqing Wua, , Chengjun Chenb, Fenglei Huanga ⁎

a b

⁎

State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, PR China Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang, Sichuan 621900, PR China

ARTICLE INFO

ABSTRACT

Keywords: PBX explosives Numerical simulation Viscoplasticity Damage evolution Tension–compression asymmetry

A viscoelastic–viscoplastic damage model is developed to describe the inelastic stress–strain responses and fracture processes of plastic-bonded explosives (PBXs) under complex stress states. This model improves the viscoelastic cracking constitutive model called ViscoSCRAM in two respects. (1) The growth rate factor of pressure-dependent tensile damage is incorporated into microcrack evolution equations, and the tension–compression asymmetry feature of quasi-brittle solids can be more reasonably described than the previous one. (2) The viscoplastic property, which is essential to PBXs, especially 1,3,5-triamino-2,4,6-trinitrobenzene-based explosives, is introduced into the framework through an additive decomposition of the total strain rate based on Bodner–Partom viscoplasticity. An efficient numerical scheme in the framework of large inelastic strains is developed. Several single-element tests under various loading paths and a specimen containing a cavity under compression are used to illustrate the main features and predictive capabilities of the proposed model. The numerical simulation can accurately capture the tensile cracks and diffuse shear-dominated material damage of the PBX 9502 perforated plate in comparison with the test conducted by Los Alamos National Laboratory (Liu and Thompson, 2014).

1. Introduction Plastic-bonded explosives (PBXs), such as PBX 9502, are considered particulate composite materials containing majority of solid-phase high explosive crystals embedded into a viscoelastic polymer binder (Thompson et al., 2010; Xiao et al., 2017). When these materials are subjected to mechanical and/or thermal loads, their internal structure may be altered, thus making them essentially a different material of unknown mechanical performance and ignition sensitivity (Drodge and Williamson, 2016). Therefore, understanding the mechanical behavior and damage evolution mechanisms under complex stress states is vital for safely using PBXs. The two classes of PBXs in accordance with the difference in manufacturing procedure are cast and pressed. The latter, especially 1,3,5triamino-2,4,6-trinitrobenzene (TATB)-based explosive formulations, is focused on in this work. Numerous experiments have been conducted to address the mechanical properties of pressed PBXs. The results indicated the following features of the explosives. The first feature is the pressure-dependent mechanical strength (Wiegand et al., 2011; Wiegand and Reddingius, 2005). Slow-crack failure processes are important in a low-pressure range, whereas plastic flow dominates in a

⁎

high-pressure range. The second feature is the tension–compression asymmetry feature (Thompson et al., 2012). The maximum stress and its strain under tension are several times less than those recorded during compressive tests. The plasticity of PBX 9502 is evident, and an approximately 60° fracture line appears under compression at high strains. Brittle rupture with a fracture surface perpendicular to the loading axis is observed for tensile specimens. Thus, quasi-brittle behavior is exhibited. The third feature includes strain rates and temperature sensitivity. The principles of time–temperature superposition have been applied to the compressive failure stress of EDC37 (Williamson et al., 2008). Furthermore, Thompson et al. (2012) applied the principles to analyze the interrelationships of stress, strain, and modulus values for PBX 9501 and PBX 9502, and found that the principles are valid for the two materials for temperatures lower than approximately Tg + 30 K. Many studies have been conducted to address the constitutive model and damage/failure simulation for PBXs based on experimental knowledge. Several researchers have used mesoscale simulations to investigate the damage behavior of PBXs by explicitly considering the microstructure, constituent properties, and interfacial responses (Barua and Zhou, 2013; LaBarbera and Zikry, 2015; Wang et al., 2016;

Corresponding authors. E-mail addresses: [email protected] (X. Huang), [email protected] (Y. Wu).

https://doi.org/10.1016/j.mechmat.2019.103179 Received 17 May 2019; Received in revised form 16 September 2019; Accepted 17 September 2019 Available online 18 September 2019 0167-6636/ © 2019 Elsevier Ltd. All rights reserved.

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Wu and Huang, 2009). Several different constitutive models have been adopted for explosive crystals, binders, and interfaces. Mesoscale modeling is an excellent tool for capturing the partial failure process. However, the extension from mesoscale to full-scale tests is costly and unsuitable in practical applications. Generally, PBXs are regarded as a quasi-brittle concrete-like material in terms of pressure dependence. Several researchers have borrowed the strength models that are commonly used for rocks, concretes, and soils to describe the yielding and hardening behavior of explosives (Gratton et al., 2009; Huang et al., 2017; Reaugh and Jones, 2010). Reasonable results can be obtained using these phenomenological models, and the pressure-dependent feature of PBXs can be effectively characterized. However, the microcrack-based damage mechanism, which is important for PBXs under tensile/shear stresses, is seldom considered in these phenomenological models (Zuo et al., 2010). The viscoelastic cracking constitutive model ViscoSCRAM is a viscoelastic continuum damage model, and a micromechanism, namely, microcracking, is used to consider nonlinear behavior and softening observed in experiments. The model has been used to study the thermal and mechanical behavior and non-shock ignition of PBX 9501 (Bennett et al., 1998; Hackett and Bennett, 2000). The damage model is based on the work of Addessio and Johnson (1990) combined with the statistical crack mechanics (SCRAM) approach of Dienes (1996), and several assumptions have been made to provide a simplified formulation for the characterization of brittle materials with isotropic deformation. ViscoSCRAM model involves the following assumptions: (1) the distribution and size of microcracks are assumed to be isotropic, and the growth of an average crack size is adopted to characterize the damage evolution of PBXs; (2) the distance between cracks is assumed to be larger than the crack size until complete material failure occurs, thus, the effects of crack bifurcations and interactions are ignored (Addessio and Johnson, 1990). Considerable efforts have been exerted to validate and improve the performance of the ViscoSCRAM model for PBXs. In the aspect of mechanical characterization, Rangaswamy et al. (2010) exhibited and discussed the validity of the ViscoSCRAM model using the examples on three-point bending, split Hopkinson pressure bar, and Brazilian disk tests. The results showed that the material model can be used to capture the mechanical behavior of PBX 9501 up to or near failure. Buechler and Luscher (2014) presented a semi-implicit update scheme that is suitable for simulations using relatively large time steps; moreover, the crack evolution equations were modified to address the diverse effect of pressure for crack growth under tension and compression. Yang et al. (2018) incorporated the generalized Griffith instability criterion with the dominant crack algorithm (Zuo et al., 2006) and the rate-dependent plasticity of linear hardening into the framework of ViscoSCRAM; the mechanical response and damage evolution of mildimpacted PBXs were also quantified. The ViscoSCRAM model is used to simulate high explosives in various scenarios, such as the ignition of PBX 9501 in a projectile during penetration (Sun et al., 2015), damage responses of RDX-based PBXs under uniaxial and multi-axial impact loadings (Xiao et al., 2017), and Steven test for HMX-based PBXs (Liu and Chen, 2018). Despite these successes, several issues still require further explorations to characterize the mechanical behavior of PBXs. These issues can be summarized as follows: (1) except for the viscoelastic and damage effects of the ViscoSCRAM model, PBXs must be characterized accurately in terms of viscoplastic behavior (Gratton et al., 2009; Le et al., 2010). (2) PBXs deserve considerable attention to study the microcracking evolution equations for reflecting the tension–compression asymmetry feature, which is important for quasi-brittle solids. (3) The parameters for PBXs, such as PBX 9501, PBX 9502, and EDC37, must be recalibrated. The viscoplasticity of PBXs may originate from crystal plastic slip or irreversible reorientation of the binder under high-stress states. The Bodner–Partom model (B–P model hereinafter) is a state variable-based

viscoplastic constitutive model and originally designed to simulate the behavior of metal alloy that undergoes creep at high temperatures (Bodner and Partom, 1975). The B–P equations can model strain rates, temperatures, and loading history effects. Such equations exhibit some apparent advantages for applications in hydrocodes, especially under high strain rate loading conditions. The B–P model presents no explicit yield surface and can thus reproduce the experimentally observed smooth transition from elasticity to plasticity while the loads are increased. The benefits of the B–P model have been adopted recently to simulate the mechanical response of certain nonmetallic materials, such as polymers (Frank and Brockman, 2001; Zaïri et al., 2008) and solid propellants (Pyrz and Zalewski, 2015). In the current work, a constitutive model that combines the ViscoSCRAM model and B–P viscoplasticity is presented to characterize the viscoelasticity, viscoplasticity, and microcracking-based damage behavior of pressed PBXs, especially TATB-based explosive formulations. The microcracking evolution equations of the ViscoSCRAM model are modified to describe the tension–compression asymmetry feature of PBXs effectively. This work aims to provide an enhanced understanding of the mechanical behavior and damage evolution under complicated stress states for PBXs. The remainder of the paper is organized as follows. In Section 2, the theoretical formulation of the proposed model is introduced, and the reasons for improving damage evolution equations are discussed in detail. The 3D numerical algorithm for the model is described in Section 3. After the model parameters for PBX 9502 (95 wt. % TATB and 5 wt.% Kel-F 800 binder) and another TATB-based explosive are identified, the numerical results of several single-element tests and a structural application are presented in Section 4. The features and predictive capabilities of the model under various loading paths and complex stress states are also discussed. The conclusions drawn from this work are provided in Section 5. Standard notations are used throughout this work. All tensor components are written with respect to a fixed Cartesian coordinate system, and the summation convention is used for repeated Latin indices, unless otherwise indicated. The following definitions are used in this work: 3

3

aij bij =

aij bij , aii = i, j = 1

3

aii , ai, j = ai / xj , Aijkl bkl = i, j = 1

Aijkl bkl , k, l= 1

and ai = ai / t . aiis a vector, aijis a second-order tensor, xj is a material point coordinate, Aijklis a fourth-order tensor, and t is time. A dot over a variable indicates the derivative of the variable with respect to time. 2. Model formulation 2.1. Decomposition of stress and strain The proposed model exhibits viscoelastic, viscoplastic, and damage properties. A schematic of the model is depicted in Fig. 1. The bottom branch of the model is a generalized viscoelastic Maxwell model with N parallel elements. The middle and top branches represent the viscoplastic properties derived from the B–P model and the isotropic damage model of ViscoSCRAM, respectively. These branches are arranged in series; thus, the stress is common, and the strain rate is additive for each individual component. The Cauchy stress tensor σij is decomposed into the hydrostatic and deviatoric components as follows: ij

= Sij

p ij,

(1)

where Sij is the deviatoric stress tensor, p = kk /3 is the pressure, and δij is the second-order identity tensor. The rate-of-deformation tensor Dij is similarly split into dilatational (volumetric) and deviatoric components, which are defined by 1

Dij = Dijdev + 3 Dkk ij, 2

(2)

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

to define the viscoplastic deformation. The relationship between the deviatoric stress and the plastic portion of the deviatoric strain rate is

Dijp = Sij,

(8)

where λ (λ ≥ 0) is a scalar function of current state quantities. Squaring Eq. (8) yields

D2p =

2J , 2

and

=

D2p / J2

(9)

= where and J2 = are the second scalar invariants of the deviatoric plastic strain rate and deviatoric stress, respectively. The growth law of D2p is 1 p p D D 2 ij ij

D2p

1 S S 2 ij ij

Z2

D2p = D02 exp

( S 2 )h ,

(10)

eff

D02

Seff = 3J2 is where (the square of D0) is the limiting value of the effective stress, and h (h > 0) is a parameter that is mainly related to strain rate sensitivity. The loading history-dependent scalar “hardening” function Z is a measure of resistance to inelastic deformation. The growth law used in Eq. (10) assumes that plastic strain rate is very small at low stresses and is limited at high stresses (Bodner, 2002). No evidence can indicate that PBX exhibits a directional hardening effect. Thus, isotropic hardening is adopted, and a simple function for Z is used.

Fig. 1. Schematic of the proposed model. Sij is the deviatoric stress tensor. Dijve , Dijp , and Dijc correspond to the deviatoric strain rate tensors of viscoelasticity, viscoplasticity, and microcracks, respectively. G(k) is the shear modulus, and η(k) is the damping constant that corresponds to the k th (of N) Maxwell element.

Z (t ) = Z1

1

is the velocity of the material point, and Dijdev and Dkk correspond to the deviatoric and volumetric parts of Dij. Dijdev can be additively decomposed into viscoelastic, viscoplastic, and damaging parts on the basis of superimposing theory of strain rates given various physical processes used in the SCRAM model (Dienes et al., 2006):

p=

Dijp ,

Various damage evolution equations of ViscoSCRAM have been proposed (Bennett et al., 1998; Hackett and Bennett, 2000; Buechler and Luscher, 2014); these equations exhibit different expressions in the growth rate of microcracks. In this section, the damage model is further investigated in terms of the pressure dependence and tension–compression asymmetry of microcracking evolution based on the recent work of Buechler and Luscher (2014).

(3)

Dijc

2.4.1. Original damage model of ViscoSCRAM In the framework of ViscoSCRAM, only shear cracks are considered, and the effect of damage to volumetric strain is neglected. The relationship between the damage strain and the deviatoric stress is presented in Bennett et al. (1998) as

(4)

K ·Dkk ,

where K is the bulk modulus. This relationship can be conveniently replaced by a general pressure–volume–temperature equation-of-state. 2.2. Viscoelastic evolution

Sij = 2G0 Dijve (k )

Sij = 2G (k ) Dijve

Sij(k ) , (k )

N k=1

G (k ),

c 3 ij a

(12)

and its rate form is expressed as

2G0 Dijc =

(5)

( )S

c 3 ij a

+3

()

c 2c S, a a ij

(13)

where c is the average crack size that stands for the damage degree of the material in a phenomenological sense, a is a normalizing parameter of c, and eijc is the total damage strain that arises from crack evolution. The evolution equation that defines the average crack growth rate is discussed on the basis of the work of Buechler and Luscher (2014):

(6)

where G0 is the instantaneous shear modulus and is expressed as

G0 = G +

( )S,

2G0 eijc =

On the basis of a generalized Maxwell model, the deviatoric stress (k ) rate Sij and deviatoric stress rate of the kth (of N) component Sij are provided using Eqs. (5) and (6), respectively. Sij(k ) N , k = 1 (k )

(11)

2.4. Damage evolution

where and are the deviatoric strain rate tensor of viscoelasticity, viscoplasticity, and damage, respectively. Here, we use a linear bulk response for the pressure and volumetric strain rate, that is,

Dijve ,

Z0 )exp( m1 W p (t )),

(Z1

where Z0 (Z0 > 0) and Z1 (Z1 > 0) are the initial and saturation values of Z, respectively. m1 (m1 > 0) controls the hardening rate, and t W p (t ) = 0 ij (t ) Dijp (t )dt is the plastic work per unit volume.

where Dij = 2 (vi, j + vj, i ) is the symmetric part of the velocity gradient, vi

Dijdev = Dijve + Dijp + Dijc ,

D2p ,

K

(7)

c=

G∞ is the steady-state shear modulus, and G(k) and τ(k) are the shear modulus and relaxation time for the kth (of N) Maxwell element, respectively.

vres ( KI )m

for KI

1

vres 1

K 0µ ( K )2 I

K

otherwise

, (14)

where

2.3. Viscoplastic evolution

KI =

The equations developed by Bodner and Partom (1975) are adopted 3

3 c S S 2 ij ij

for p

3 c 2 ij ij

otherwise

0

, (15)

Mechanics of Materials 139 (2019) 103179

M. Liu, et al. 2 m

K = K 0µ 1 +

K1 = K 0µ (1 +

2(3

+

µc

m 1/ m ) , 2

c

K0 1 2

45 2µc2)

c¯ (K¯ ) in the following parts of this work in consideration of the monotonicity of c¯ (K¯ ) and irrelevance of vres with stress states. Combining Eqs. (15) and (16)–(18) with the definition of K¯ yields

(16)

2 1/2 ) (1 m

K 0µ = K 0 1 +

µc =

,

(1 +

µc ,

(17) µc

K0

c

),

3 c S S 2 ij ij

(18)

K¯ =

(20)

.

(21)

K¯ ( 1,

m is a model parameter, K0 is the threshold value of stress intensity, µc is the static coefficient of friction, va and vb are the crack velocity parameters, and < · > is the Macaulay bracket. The initial value of crack size c in Eq. (14) is c0, which represents the initial flaw size of the material. vres is limited by the theoretical terminal crack speed vmax (the Raleigh wave speed in principle). The frictional threshold stress intensity K0μ in Eq. (18) depends on pressure p as a way of retarding the damage evolution rate under compressive loading. Such behavior is physically consistent in terms of internal frictional resistance to crack evolution. Eq. (20) establishes the relationship between the strain and crack growth rates and can be used to reflect the effect of loading rate on microcrack propagation. However, the dimension between the left and the right of Eq. (20) is non-identical. The following expression is adopted:

log10 (vres) = va log10 (¯/¯0 ) + vb,

c¯ (K¯ ) =

1

K¯ < 1 K¯

1

,

otherwise (24)

2)

2 2 1 + 2

3 cm

=

K 0 2(m + 2)

,

(25)

¯ is that is, K¯ ( 1, 2 ) = K¯ ( 1, 2 ) . This expression indicates that K identical for biaxial tension (σ1 ≥ σ2 > 0) and tension–compression (σ1 > 0 > σ2) stress states. Such behavior is physically inconsistent with the strength asymmetry observed in experiments and can be proven indirectly by the experimental strength data of TATB-based explosive from the uniaxial tension test (σ1 > 0) with Seff = 1 = 8.56 MPa and Brazilian test ( 2 = 3 1) with Seff = 13 1 = 17.11 MPa (Tang, 2016). Rangaswamy et al. (2010) emphasized that the parameters of ViscoSCRAM are obtained from uniaxial compression measurements; the extra tensile damage growth rate factor is required to simulate the uniaxial tension response. However, the physical meaning is ambiguous when a constant coefficient is used. 2.4.2. Improvement of the damage evolution model To overcome the shortcoming discussed above, the normalized stress intensity K¯ is modified to K¯ M as follows:

(22)

where ¯0 is the reference value of the effective strain rate. Eq. (14) can be transformed into a dimensionless format as follows: 2 K¯ m m+2 m ¯ 2 K m+2

.

3 c 2 ij ij 2 K0 1 + m

Under the plane stress state, the principal stresses are σ1, σ2, and 0 (σ1 ≥ σ2). The normalized stress intensity K¯ for 2 [ 1, 1] (σ1 > 0) becomes

and 3 dev dev D Dij 2 ij

for p > 0

µ cp c µ cp c (1 + ) K0 K0

1+

(19)

log10 (vres) = va log10 (¯) + vb,

¯=

2 m

K0 1 +

(26)

K¯ M = KI Kt / K , where

Kt = 1 +

(23)

µt< m> K0

c

(1 +

µt < m> K0

c

).

(27)

Kt depends on the mean stress m = p and average crack radius c as a way of accelerating the damage evolution rate under negative pressure region, that is, tensile state ( m < 0 ). The definition of the tensile coefficient µt is similar to that in Eq. (19).

where c¯ (K¯ ) = c (K¯ )/ vres is the normalized crack velocity and K¯ = KI/K is the normalized stress intensity. c¯ (K¯ ) is a monotonically increasing function of K¯ , as stated in Fig. 2. K¯ = 1 is the transition point that corresponds to stable (K¯ < 1) and unstable (K¯ > 1) crack growth in Eq. (23). K¯ is a function of the current stress state and average crack radius. The damage evolution model is discussed using K¯ rather than

µt =

1 2

45 2(3

2µt2)

µt ,

(28)

where µt is the model parameter. K¯ M degenerates to K¯ when µt = 0 . Eq. (26) can be interpreted as an effect caused by hydrostatic pressure. The damage evolution rate is retarded under compressive load given internal frictional resistance. Such a rate is accelerated under tensile load due to the open microcracks. The normalized stress intensities K¯ and K¯ M are discussed under the Haigh–Westergaard (H–W) coordinate system (ξ, ρ, θ) to illustrate the 1 effectiveness of the modified damage model, where = 3 I1 is the volumetric effective stress,

=

1 3

cos 1 (

3 3 J3 ) 2 J 3/2 2

=

2J2 is the deviatoric effective stress, and

is the Lode angle. I1 =

ij ij

is the first invariant of

stress tensor σij; and J2 = and J3 = are the second and third invariants of the deviatoric stress tensor Sij, respectively. The principal stresses σ1, σ2, and σ3 (σ1 ≥ σ2 ≥ σ3) can be denoted using the H–W system as

Fig. 2. Relationship between the normalized crack rate c¯ (K¯ ) = c (K¯ )/ vres and the normalized stress intensity K¯ = KI/ K in Eq. (23), where parameter m is 10. Point T illustrates the transition between stable and unstable crack growth.

1(

, , )=

2 3

2(

, , )=

2 3

, , )=

2 3

3(

4

( sin ( sin ( sin

1 S S S 3 ij jk ki

1 S S 2 ij ij

+

2 6 5 6

)+ )+ )+

1 3 1 3 1 3

. (29)

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Then, K¯ M is 3

K¯ M (

, , c) =

2m + 4 3K 02

mc 3 µ c K0

2 c + µ t 2c

3K 02 + 3 µ t K 0

2 c + µ c 2c

mc

for < 0

3. Numerical implementation

2+ 2

otherwise

m + 2 K 03

6

Fig. 3. The subsequent sections will further demonstrate the efficiency of the modified damage model using several detailed examples.

(30)

In this section, an efficient numerical scheme that couples explicit and implicit time integration for the proposed model is developed and implemented in the parallel 3D explicit finite element method (FEM) software called PANDA-Impact. This software was developed by the Institute of Systems Engineering of CAEP. The numerical algorithms are based on an explicit viscoelastic damage predictor step followed by an implicit plastic corrector step. The resulting nonlinear systems of equations are solved through a Newton–Raphson iterative method. [tn, tn + 1], t = tn + 1 tn , and Generic time intervals ( ) = ( )(tn + 1) ( )(tn) are considered. Values at tn and tn + 1 are designated with superscripts n and n + 1, respectively, to simplify the algorithmic expressions. The algorithm is strain driven. That is, the solutions of Cauchy stress n+1 and internal variables y n + 1 are required based on a known state of ij n Cauchy stress ijn , internal variables y n = {W p, n, c n, c n, pn , Sij , Sij(k ), n} , and a new incremental strain Δɛij. A suitable objective stress rate must be selected for hypoelasticbased models to evaluate the large-deformation, materially nonlinear response of solids. Several internal variables, such as the deviatoric stress tensors of the Maxwell elements, are tensors in the proposed model. The Jaumann rate commonly used in the transient dynamics codes is difficult to utilize as the objective stress rate for tensors. An appropriate choice of the frame of reference is the unrotated configuration, and Green–Naghdi rate is adopted (Flanagan and Taylor, 1987).

when µt = 0 , andK¯ is obtained using Eq. (30). 3

K¯ ( , , c ) =

2m + 4 3K 02

mc 3 µ c K0

3 mc

2+ 2

2m + 4 K 0

2 c + µ c 2c

for < 0 , otherwise

(31)

where is a function of ξ, ρ, and c, except for θ. The concept of the yield surface used in plastic mechanics is adopted for the following explanation. The surface of K¯ M = A , where A is constant and A > 0 is assumed, denotes an isosurface in which the crack growth rate is identical. If c is constant, then the shape of the isopleth for K¯ M = A is a circle in the deviatoric plane when ξ is fixed. Thus, the evolution of the damage model in 3D stress space can be described by the meridian, which is obtained by the intersection of the surface of K¯ M = A with a plane containing the hydrostatic axis. The shape of the meridian is determined by the average crack radius c and parameters K0, µc , µt , and m. Fig. 3 demonstrates that (1) the crack growth rate is pressure-dependent before and after improving the damage model. The rate is low with a large pressure under the same deviatoric effective stress. (2) The crack growth rate is dependent on the average crack size c. The stress level is low with a large crack size to achieve the same stress intensity (for instance, K¯ = K¯ M = 1.0 in Fig. 3). (3) After the damage model is improved, the crack is easy to expand under tension with the same crack size (for instance, c = 0.03 mm in Fig. 3), and the difference in the crack growth rate between tension and compression is increasingly evident. The damage evolution model exhibits no concept of the damage surface, where damage starts to accumulate or the material is completely failed; thus, the damage will accumulate when stress is present. The average crack size c is adopted to indicate the damage degree of the material. A rapid crack growth rate in each loading path indicates a low material strength. It can be used to illustrate the trend of material strength (peak stress) over the entire stress-space through the relationships of the meridians with various loading paths exhibited in

K¯ M

3.1. Viscoelastic damage predictor The crack growth kinetics is difficult to solve in a reliable and robust manner within a fully implicit scheme (Buechler and Luscher, 2014); thus, the equations of viscoelastic stress in conjunction with the damage evolution are solved through a fourth-order explicit Runge–Kutta scheme. Combining Eqs. (3), (5), (6), and (13) yields the rate form of the total deviatoric stress

Sij =

(k )

Dijp )

2G0 (Dijdev

1+

Sij N k = 1 (k ) c 3

3

()

()

c 2c S a a ij

, (32)

a

and deviatoric stress of the kth Maxwell element (k )

Sij = 2G (k ) (Dijdev

Sij(k )

Dijp)

(k )

G (k ) G0

( )S

c 3 ij a

+3

()

c 2c S a a ij

.

(33)

With the assumption that the plastic strain rate Dijp is constant (a known variable), Eqs. (4), (14), (32), and (33) are solved simultaneously. With (k ) Y {Sij(k), Sij, p , c } and Y {Sij , Sij, p , c} , the fourth-order explicit Runge–Kutta scheme can be described as

Y n +1 = Yn +

1 6

t (K1 + 2K2 + 2K3 + K 4),

(34)

where

K1 = Y (tn, Y n)

Fig. 3. Relationships of the meridians with various loading paths. Parameters K 0 = 15.81 MPa mm , m = 10 , and µc = 0.5 are adopted in accordance with Hackett and Bennett (2000). µt = 0.5 is presumed. Loading path: biaxial compression (OA, 1 = 0 > 2 = 3 ), uniaxial compression (OB, 1 = 2 = 0 > 3 ), biaxial shear (OC, 1 = 3 > 0 = 2 ), uniaxial tension (OD, 1 > 0 = 2 = 3 ), and biaxial tension (OE, 1 = 2 > 0 = 3 ).

( = Y (t

K2 = Y tn +

1 2

t, Y n +

1 2

K3

1 2

Yn

1 2

n

+

t,

+

). tK ) t K1 2

K 4 = Y (tn + t , Y n + t K3) 5

(35)

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

3.2. Plastic corrector and Newton–Raphson iteration

Original model, and that with the improved damage model discussed in Section 2.4.2 is denoted as VEVPD-Modified model. The unique distinction between the two models is the obtained normalized stress intensity in the negative-pressure region. The VEVPD-Original model without plasticity is called ViscoSCRAM model.

The adopted iterative radial corrector method is mainly based on Cook et al. (1992) with necessary modifications. The plastic corrector step is deduced in an implicit backward Euler method (Belytschko et al., 2014). As derived in Appendix A, the increments in plastic strain and other internal variables are calculated at the end of each step, and the yield condition is also enforced at the end of the step. The implicit backward Euler scheme for plastic corrector step avoids nonphysical effects, such as spurious unloading driven by nonconverged values of plastic strain. Eqs. (8) and (10) are rearranged as follows: 2 3

p Deff =

D0 exp

(

1 Z 2h ( ) 2 Seff

),

4.1. Identification of model parameters The parameters of the ViscoSCRAM model provided in Hackett and Bennett (2000) are extensively used by researchers to describe the mechanical behavior of PBX 9501. However, these parameters cannot be used to characterize TATB-based explosives because the modulus and strength are different from those of PBX 9501 (Thompson et al., 2012). In addition, the original ViscoSCRAM model contains no viscoplasticity property. The proposed model exhibits numerous parameters, including the Maxwell model of viscoelasticity G∞, G(k), and τ(k); the B–P model of viscoplasticity Z0, Z1, n, and m1; and the microcracking evolution of damage a, c0, vmax, K0, m, μt, µc , va, and vb. These parameters are difficult to identify through an analytical method given the mathematical complexity of the formulation. Thus, a numerical method is adopted to determine the parameters for PBX 9502 and PBX1. The calibration procedure is conducted as an optimization problem following Pyrz and Zairi (2007). The objective function for the optimization problem is formulated as the sum of the errors between the experimental and the numerical stress–strain curves at the selected points of the investigated domain, and a minimization problem is represented for the curve-fitting problem. The corresponding objective function is

(36)

2

p = 3 Dijp Dijp is the effective plastic strain rate. An expression where Deff for Seff obtained from the inverted form of Eq. (36) is given by

Seff = Z

2 ln

1 2h

p 3 Deff

2D0

.

(37)

With the expression derived in Appendix A, an estimation of the new R is computed by effective stress Seff R Seff =

tr S¯eff

p 3G 0 eeff . 1 + c, n + 1

(38)

The values of and are obtained by the fourth-order explicit p p = t ·Deff Runge–Kutta scheme as shown in Section 3.1, and eeff . Seff in p R Eq. (37) and Seff in Eq. (38) are functions of Deff . Thus, the nonlinear equations can be written in a residual form as follows: tr S¯eff

p p R (Deff ) = Seff (Deff )

c, n + 1

p R Seff (Deff ).

J (c1,

(39)

R

1

p,(k) Deff

p,(k) R (Deff ).

1 M

M i=1

(

exp i

num (c1, i

, cN )) 2

min,

(41)

where c1, ⋅⋅⋅, cN represent the N material constants that must be determined, inum is the numerical stress that corresponds to strain inum at each time step i (i = 1, ⋅⋅⋅, M), and iexp is the experimental stress by piecewise linear interpolation at each strain inum . A stochastic search method based on evolutionary algorithms is adopted to solve the unconstrained optimization problems considering its efficiency in difficult, large size, high cardinality, and multimodal problems (Pyrz and Zairi, 2007). The method is implemented in PANDA-Impact. Before the calibration procedure, several parameters are fixed in accordance with references, such as density = 1.89 × 10 3 g/mm3 and = 0.4 cited in Huang et al. (2017), parameter Poisson ratio D0 = 105 ms 1 in Bodner and Partom (1975), and initial crack size c0 = 0.003 mm and theoretical terminal crack speed vmax = 300 mm/ms in Hackett and Bennett (2000). The experimental curves selected to calibrate the model parameters for PBX 9502 are the stress–strain curves labeled as “Charge 01″ for 50 °C tension and compression tests under quasi-static loading in Thompson et al. (2010). During the calibration procedure, the correspondence between the experimental and the simulated curves under uniaxial compression is first guaranteed. The peak stress of uniaxial tension is then calibrated by parameter μt for the VEVPD-Modified model. The fitting curves are plotted in Fig. 4, and the corresponding parameters are listed in Tables 1–3. The sensitivity analysis of these parameters is given in Appendix C. If the proposed model exhibits no plastic property, namely, only the ViscoSCRAM model is adopted, then the fitting result of the compressive curve is poor (pink curve in Fig. 4). The softening stage drops steeply with the constraint of the same peak compressive stress in experiment. The result of the VEVPD-Original model under uniaxial tension is also included in Fig. 4 for comparison. The quasi-static longitudinal stress–strain curve (Ambos et al., 2015) is chosen to calibrate the model parameters of PBX1, and the detailed loading and unloading processes are described in Section 4.2.3. The tensile curve of PBX1 is not given in the literature; thus, parameter μt cannot be fixed. The fitting result of the longitudinal stress–strain

p ) through the The Newton–Raphson method is applied to solve R (Deff linearized formulation

p,(k + 1) p,(k) Deff = Deff

, cN ) =

(40)

The iteration procedure is repeated until convergence, that is, p,(k) p,(k + 1) p,(k) )| < tol with a given tolerance |Deff Deff | < tol and |R (Deff tol = 1 × 10 10 . After convergence is achieved, the internal variables are updated, and true stresses at the end of the step are calculated. The entire algorithmic procedure is summarized in Appendix B. 4. Results and discussion A series of numerical examples is provided in this section to demonstrate the effectiveness of the proposed model for pressed TATBbased high explosive materials. The model parameters for PBX 9502 (Thompson et al., 2010) and another TATB-based explosive (Ambos et al., 2015) are initially identified through a numerical method. In the following sections, the material studied in Ambos et al. (2015) is called PBX1 for convenience. Several singleelement tests, which are intended as an elementary verification of the basic modeling capabilities, are considered. A structural application, that is, the perforated plate under compression, is used to validate the proposed model. The predictive ability of the proposed model under various loading paths and complex stress states is discussed in this section. All examples are performed using eight-node solid elements with one-point Gauss integration scheme, except the case for mesh dependence in Section 4.3.3. In addition, the loading strain rates are set to 1 s−1 for PBX 9502 and 5 s−1 for PBX1, and the same loading conditions are adopted for the parameter identification as well. The use of high loading rate calculation, that is, explicit algorithm, for quasistatic problems is discussed in Section 4.3.1. For ease of reference, the viscoelastic–viscoplastic model with the original damage model presented in Section 2.4.1 is called VEVPD6

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 4. Fitting results of uniaxial tension and compression with tests (Thompson et al., 2010) for the VEVPD-Modified model. The results of the VEVPD-Original model under uniaxial tension and the ViscoSCRAM model under uniaxial compression are included for comparison.

Fig. 5. Experiment results (Ambos et al., 2015) versus curves by VEVPDModified model for material response to cyclic uniaxial compression. Table 4 Model parameters of viscoelasticity for PBX1.

Table 1 Model parameters of viscoelasticity for PBX 9502 at 50 ∘C. ρ[g/mm3] 0.00189

0.4

G∞[MPa]

G(k)[MPa]

τ(k)[ms]

296.065

237.751, 1200, 1063.24, 30.013

10, 1, 0.1, 0.01

ρ[g/mm3] 0.00189

Table 2 Model parameters of viscoplasticity for PBX 9502 at 50 ∘C. Z0 [MPa] 257.02

Z1 [MPa] 335.497

D0 [ms−1] 10

5

0.4

G∞[MPa]

G(k)[MPa]

τ(k)[ms]

319.558

750.149, 1486.52, 150.095, 1948.72

100, 1, 0.01, 0.0001

Table 5 Model parameters of viscoplasticity for PBX1.

h

m1

0.578442

16.4986

Z0 [MPa] 38.2188

Z1 [MPa] 62.2017

D0 [ms−1] 10

5

h

m1

2.5

16.3802

curve is plotted in Fig. 5 (the red solid curve with strain < 0), and the corresponding parameters are listed in Tables 4–6. The transverse stress–strain curve obtained by calculation (the red solid curve with strain > 0) and experiment is also depicted in Fig. 5 for comparison.

Sections 4.2.1 and 4.2.2, and the parameters of PBX1 listed in Tables 4–6 are adopted in Section 4.2.3. µt = 0 is used for the VEVPD-Original model, and G = 2400 MPa and G (k) = 0 MPa are used for the EVP and ED models.

4.2. Model features under various loading paths

4.2.1. Uniaxial loading The simulated results under uniaxial tension and compression are depicted in Figs. 6–8. The results of the VEVPD-Original model under compression are the same as those of the VEVPD-Modified model, whereas they are not illustrated in Fig. 7. The tension–compression asymmetry features, such as the distinctness of peak stress between compression and tension, are evident for PBX 9502. In Figs. 6a and 7a, the ratio of compressive to tensile peak stress is 2.72 using the VEVPD-Modified model. This ratio is close to the experimental result, as stated in Section 4.1. The ratio calculated using the VEVPD-Original model is 1.22, which is far from the experimental result. The damage parameter d = c 3/(a3 + c 3) is adopted in accordance with Rangaswamy et al. (2010). As described in Section 2.4, c is the average crack size with an initial value of c0, and a is a normalizing parameter. The nonzero value c0 stands for the initial flaw size inside or between crystals during manufacturing (Le et al., 2010). Thus, d is a nonzero value at the beginning of the simulation, as depicted in Figs. 6b

A series of single-element tests, including uniaxial tension ( 1 > 0 = 2 = 3 ), uniaxial compression ( 1 = 2 = 0 > 3), biaxial tension ( 1 = 2 > 0 = 3 ), biaxial compression ( 1 = 0 > 2 = 3 ), biaxial shear ( 1 = 3 > 0 = 2 ), and cyclic uniaxial compression, is provided to illustrate the main features of the proposed model. The stress response, damage evolution, and plastic accumulation under various loading paths are discussed in detail. The performance of the improved damage evolution equations is demonstrated by comparing the numerical results of the VEVPD-Original and VEVPD-Modified models under uniaxial tension. The role of plasticity is demonstrated by comparing the numerical results of the VEVPD-Modified and ViscoSCRAM models under uniaxial compression. Furthermore, the performance of viscoelasticity (VE), elastic viscoplasticity (EVP), elastic damage (ED), and viscoelastic damage (ViscoSCRAM) under loading/ unloading is demonstrated by the cyclic uniaxial compression test. The parameters of PBX 9502 listed in Tables 1–3 are adopted in Table 3 Model parameters of damage for PBX 9502 at 50 ∘C. a[mm]

c0[mm]

vmax[mm/ms]

K0[MPa mm1/2]

m

µc

μt

va

vb

¯0 [ms−1]

0.009

0.003

300

2.0079

6

0.339447

1.16

0.155984

−0.5

1

7

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Table 6 Model parameters of damage for PBX1. a[mm]

c0[mm]

vmax[mm/ms]

K0[MPa mm1/2]

m

µc

μt

va

vb

¯0 [ms−1]

0.009

0.003

300

4.51838

6

0.257213

–

0.123789

−0.5

1

and 7b. The average crack size c denotes the degree of damage accumulation for materials and is a monotonically increasing variable. The normalized stress intensity K¯ M represents the growth rate of microcracks. In Figs. 6c and 7c, the value of K¯ M simulated by the VEVPD-Modified model first increases with the development of the stress level (O–B), continues to grow until Point C after the peak stress (B), and then slowly decreases with the stress level (C–D). According to the theory of fracture mechanics, the crack speed is correlated with the critical stress intensity factor of crack tips. The critical stress intensity factor is a function of ¯ c , where ¯ is the effective stress at crack tips, and c is the effective length of cracks. Thus, K¯ continues to grow until Point C given the increment in the value of ¯ c . The pressure dependence of the damage evolution rate is illustrated

by the relationship between the uniaxial loading path and the meridians of Points A, B, C, and D in Figs. 6d and 7d. The deviatoric effective stress ρ is a monotonic function of the volumetric effective stress ξ in 1 each meridian, and = 3 I1 is proportional to pressure. For conciseness, the meridian of Point B is not presented in Fig. 6d. Furthermore, the meridian is calculated using Eqs. (30) and (31) with the given values of stresses and crack radius; the comparison result of Figs. 6c and 7c exhibits that the value of K¯ is identical for simulation and theory; thus the validity of code implementation can be partially demonstrated. The proposed model is based on the additive decomposition of total strain rate, and the magnitude of peak stress depends on the growth rate of microcracks. In Fig. 6, the tendency of the calculated results using the VEVPD-Original model under tension is similar to those of the

Fig. 6. Comparison of uniaxial tension responses using the VEVPD-Original and VEVPD-Modified models: (a) stress–strain responses, (b) evolution of the damage with strain, (c) evolution of the normalized stress intensities K¯ and K¯ M with strain, and (d) relationship between the uniaxial loading path and the meridians of K¯ and K¯ M . The selected moments labeled from A to D are along the VEVPD-Modified curve, and those from E to H are along the VEVPD-Original curve. Point O is the origin of simulation, A and E are the starting points of rapidly developing damage, B and F are the peak points of stress. In addition, C and G are the peak points of normalized stress intensity, and D and H are the end points of simulation with d = 0.9 . 8

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 7. Comparison between the uniaxial compression responses using the VEVPD-Modified and ViscoSCRAM models: (a) stress–strain responses, (b) evolution of the damage with strain, (c) evolution of the normalized stress intensity with strain, and (d) relationship between the uniaxial loading path and the meridians of K¯ M . The selected moments labeled from A to D are along the VEVPD-Modified curve, and those from F to H are along the ViscoSCRAM curve. Point O is the origin of simulation, A is the starting point of rapidly developing damage, and B and F are the peak points of stress. C and G are the peak points of normalized stress intensity, and D and H are the end points of simulation with d = 0.97 .

VEVPD-Modified model, except that the peak stress is large, and the damage evolution rate is slow. The different results are attributed to the damage evolution equations that are distinct, as demonstrated Fig. 6d and Section 2.4.2. In Dienes et al. (2006) and Zuo et al. (2010), PBXs are known to undergo plastic deformation even under normal application conditions (of less confinement and low-pressure environment). Neglecting plastic deformation can lead to an over-prediction of the actual stress in the materials. If the model exhibits no plastic property, as the results simulated using the ViscoSCRAM model, then the peak stress is high, and the growth rate of microcracks is rapid (Fig. 7). The evolution of damage and plasticity interacts with one another through the link of the stress level. The B–P model shows no dependence on an explicit yield criterion or loading and unloading conditions. If a nonzero stress is present, then a plastic accumulation will exist in theory. However, the plastic work slightly increases when the stress level is low. A cutoff value of 10 10 , as depicted in Section 3.2, is used for the effective plastic strain rate to ensure the robustness and efficiency of the program. That is, no plastic work will accumulate if the effective plastic strain rate at each time step

Fig. 8. Responses of stress (left) and plastic work per unit volume (right) versus strain under compression. Point A is the starting point of plastic accumulation, B is the peak point of stress, C is the end of plastic accumulation, and D is the end of simulation. 9

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 9. Comparison of the responses predicted using the VEVPD-Modified model under uniaxial tension, uniaxial compression, biaxial tension, biaxial compression, and biaxial shear: (a) effective stress–strain response, (b) evolution of the damage with strain, (c) evolution of the normalized stress intensity with strain, and (d) evolution of the plastic work per unit volume with strain.

is less than the cutoff value. This implementation gives the B–P model an initial yielding value. Fig. 8 plots the stress and plastic work per unit volume versus strain curves under compression using double y-coordinates. The material starts to deform plastically at Point A, which is the initial yield point. The plastic work gradually increases with the development of the stress level (A–B), continues to grow after the peak stress (B), and then ends at Point C. Given the plastic-hardening effect, the absolute value of compressive stress is larger at Point C than at Point A (8.34 MPa vs. 7.19 MPa). The strength of 7.19 MPa at Point A (the effective stress or strength value is the same as the uniaxial stress for uniaxial loading) is the initial yielding value. Thus, no plastic deformation exists under tension for the VEVPD-Modified model because the peak stress is only 4.59 MPa (Fig. 5a), which is consistent with the quasi-brittle feature of the experimental observation for PBXs.

arrangement sequences for peak stresses from low to high are listed as follows: biaxial tension, uniaxial tension, biaxial shear, uniaxial compression, and biaxial compression. The damage evolution rate and normalized stress intensity K¯ M are opposite, as displayed in Figs. 9b and c. As previously illustrated, the magnitude of peak stress depends on the growth rate of microcracks, and the damage evolution equation is pressure dependent. A proper interpretation for the sequences of peak stresses under various loading paths can be obtained by combining the meridians of K¯ M illustrated in Figs. 6d and 7d. The relationships between meridians and various loading paths are plotted in Fig. 3; thus, the meridians of K¯ M can be used to predict the features of the proposed model subject to other loading cases. The results demonstrated in Fig. 9d show that the plastic work per unit volume under biaxial compression exhibits a large value with 0.24 MPa. The cases of uniaxial compression and biaxial shear are small, whereas those of biaxial and uniaxial tensions are zero. The plastic work is an effective stress-induced variable, as depicted in Section 2.3. If the effective stress is high and sustains above the yield value, then the plastic work accumulates considerably.

4.2.2. Biaxial loading The simulated results using the VEVPD-Modified model under biaxial tension, biaxial compression, and biaxial shear are depicted in Fig. 9. The results of uniaxial loading are also presented for comparison. The results exhibited in Fig. 9a show that the peak effective stress and the corresponding effective strain are the least under biaxial tension and largest under biaxial compression (3.11 MPa vs. 14.82 MPa and 0.06% vs. 2.52% for stresses and strains, respectively). The

4.2.3. Cyclic uniaxial compression In the cyclic uniaxial compression, four cycles are applied in one experiment (Ambos et al., 2015). In each cycle, several different phases 10

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 10. Simulation of the cyclic uniaxial compressive test: (a) prescribed displacement–time curve for the platen, and longitudinal stress–strain curves predicted by (b) VE, (c) EVP, (d) ED, (e) ViscoSCRAM, and (f) VEVPD-Modified models. The red line in each figure stands for the second cycle, and the selected moments labeled from A to E are along the red line. Points A and E are the starting points of the loading phase, B and C are the starting and end points of the holding phase, and D is the starting point of the recovery phase.

11

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

are observed in the stress–strain responses (Fig. 5). First, stress and strain increase monotonically during the loading phase. Second, stress drops while strain holds constant (holding phase). Third, stress and strain decrease monotonically during the unloading phase. Fourth, strain recovers while stress remains zero (recovery phase). The drop in stress in the holding phase is consistent with viscoelasticity or viscoplasticity theory. Strain recovery at zero stress in the recovery phase indicates viscoelastic property of PBXs, and the significant residual strain at the end of each cycle refers to plastic strain. Furthermore, strain softening in the fourth cycle indicates damage behavior. Thus, the viscoelastic–viscoplastic damage model is suggested to describe the material response to cyclic uniaxial compression. The results by the model are shown in Fig. 5, and they illustrate the validity of the proposed model. The behavior of the viscoelasticity, viscoplasticity, and damage in the VEVPD-Modified model under cyclic uniaxial compression is given in detail in this section. The following analyses mainly focus on the inelastic deformation in the holding and recovery phases. In numerical simulation, a single element (sample) is loaded by a rigid platen, of which the displacement is prescribed, as shown in Fig. 10a, and the platen holds for 1 ms before returning to zero displacement. Therefore, the strain rate of the sample during the loading (unloading) phase is 5 s−1 (−5 s−1). The interaction between the

Fig. 11. Responses of plastic work per unit volume (left) and damage (right) versus longitudinal strain under cyclic uniaxial compression. The red solid line and red dash line stand for the second cycle, and the selected moments labeled from A to E (and A` to E`) correspond to the points in Fig. 10f.

Fig. 12. Perforated plate specimen of PBX 9502 for damage failure simulation subjected to compression: (a) Physical dimensions, where characteristic points A1, B1, C1, and D1 are located on the surface of the plate, C1 is the midpoint of the right side of the plate, and D1 is the intersection point of the circular portion and straight section. (b) Partial enlarged view of the meshes with an element size of 0.6 mm.

Fig. 13. Overall compressive load versus displacement curves: (a) comparison of the results predicted using the VEVPD-Modified model with mesh sizes of 0.3, 0.6, and 0.9 mm; and (b) comparison of the results predicted using the VEVPD-Modified model and measured by Liu and Thompson (2014). 12

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

viscoplasticity, and damage properties contribute to stress drops during the holding phase. During the recovery phase, the magnitudes of strain decrease at zero stress are 0.17% for VE, 0 for EVP and ED, and 0.17% for VEVPDModified models. Thus, the contributor to strain recovery is only viscoelastic property. However, the magnitudes of residual strain (including unrecovered viscoelastic deformation) at the end of the recovery phase are 0.26% for VE, 0.61% for EVP, 0 for ED, and 0.53% for VEVPD-Modified models. The residual strain of VE is considerably lower than that of VEVPD-Modified model. Thus, plastic deformation is evident during the entire cyclic test. According to Eq. (12), the strain of ED model changes to zero when stress is zero, as shown in Fig. 10d, and is consistent with damage mechanics. The loading/unloading modulus of ED model reduces with the increment in cycling number, whereas the loading/unloading modulus of EVP model under elastic range is constant. Compared with the stress–strain curves of VEVPD-Modified model, the curves of ViscoSCRAM model (Fig. 10e) suffer higher stress level and lower residual strain. The magnitudes of the irreversible longitudinal strain Lir and transverse strain Tir of experiment (Fig. 5) at the end of the recovery phase of the third cycle are 0.92% and 0.51%, respectively. Thus, the related irreversible volumetric strain Vir is approximately 0.1% calculated by Vir = Lir + 2 Tir . Compared with Lir , Vir is not the major factor for irreversible deformation. A linear relationship for pressure and volumetric strain and a J2-based plasticity of B–P model are adopted in this study; thus, the irreversible volumetric deformations cannot be reflected. This can be solved by the following methods: incorporating viscoelasticity for bulk response and/or incorporating pressure effect for the B–P model. The dependence on pressure of the latter method can be included parametrically through h (Eq. (10)) and/or Z (Eq. (11)) as

Fig. 14. Overall kinetic energy and external work versus displacement curves. The selected moment of point F is the same point as in Fig. 13b.

sample and the platen is modeled by contact algorithms. The longitudinal stress–strain curves predicted by VE, EVP, ED, ViscoSCRAM, and VEVPD-Modified models are depicted in Figs. 10b–f. The longitudinal stress–strain curves of the second cycle are used as the characteristic data for analysis. During the holding phase, the magnitudes of stress drops are 10.35 MPa for VE, 0.83 MPa for EVP, 5.52 MPa for ED, and 5.54 MPa for the VEVPD-Modified models. Compared with the first cycle with considerably low stress level, the major contributor to stress drops is viscoelastic property. As shown in Fig. 11, the plastic work and damage keep increasing during B-C and B`-C`. Thus, all the viscoelasticity,

Fig. 15. Comparison of the contour plots of displacement field: (a) predicted using the VEVPD-Modified model and (b) measured through the experiment (Liu and Thompson, 2014) in the horizontal direction at selected moments presented in Fig. 13b. 13

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 16. Comparison of the contour plots of displacement field: (a) predicted using the VEVPD-Modified model and (b) measured through the experiment (Liu and Thompson, 2014) in the vertical direction at selected moments exhibited in Fig. 13b.

depicted in Bodner (2002).

over the simulation. The plate is partitioned by hexahedron elements with three mesh sizes of approximately 0.3, 0.6, and 0.9 mm (element numbers: 1352,736, 178,552, and 49,336) to check the mesh convergence of FEM. Fig. 13a presents the overall compressive load histories of the specimen with the three different mesh sizes. The overall load P is normalized by the initial cross-sectional area of the specimen WT, and the displacement δ is normalized by the initial specimen height L. The peak load with a mesh size of 0.6 mm is only 0.48% larger than the peak value with a mesh size of 0.3 mm. This finding indicates that a mesh size of 0.6 mm is sufficiently small to obtain reasonable results. In the following calculations, the model partitioned with a mesh size of 0.6 mm is used, as depicted in Fig. 12b. The VEVPD-Modified and VEVPD-Original models are separately used to study the damage evolution of the perforated plate with model parameters displayed in Tables 1, 2, and 3. Notably, the strain rate of the plate is 1 s−1 in modeling against 2.8 × 10 4 s−1 in experiment. Thus, a high loading rate is adopted to model a quasi-static problem. The suitability of the methodology is as follows. First, the model parameters are identified at the same strain rate as given in Section 4.1. Thus, the parameters are appropriate. Second, the overall kinetic energy is very small over the entire simulation, and the ratio of the external work and the total kinetic energy is above 107 at the earlier stage (Fig. 14). Therefore, the inertial effect can be neglected. Third, an implicit backward Euler method is adopted in the numerical implementation as given in Section 3.2. The accuracy of the solution is guaranteed. The smoothness of the overall compressive load curve, as shown in Fig. 13b, also demonstrates the equivalence between modeling under high strain rate and experiment under quasistatic conditions.

4.3. Model features under complex stress states A compressed plate containing a cavity is simulated using the proposed model and compared with experiments conducted by Los Alamos National Laboratory (Liu and Thompson, 2014) to illustrate the efficiency of the proposed model further. The test aimed to study the process of damage initiation, accumulation, and cracking in a perforated plate specimen of PBX 9502 under overall compression by an optical digital image correlation (DIC) technique. The concrete damage plasticity model combined with the extended FEM and cohesive model was adopted to analyze the LANL experiment in Huang et al. (2017); reasonable results for the tensile cracking process are obtained, but the shear-dominated damage is not further discussed. The damage/failure processes of the specimen are studied using the proposed model, and the predictive ability of the proposed model under complex stress states is the primary focus in this section. 4.3.1. Computation model Fig. 12a displays the schematic of the computation model. The overall dimensions of the plate are L = 76.2 mm in height, W = 38.1 mm in width, and T = 12.7 mm in thickness. The cavity at the center of the plate consists of four circular arcs connected by four straight lines. The material of the specimen is PBX 9502 high explosive, and the experimental temperature is 50 °C. The bottom surface of the specimen is constrained in the vertical direction, and the top surface of the specimen is prescribed by fixed velocity in the vertical direction and toward the bottom surface to make an overall compressive load. The strain rate is approximately 1 s−1. The frictional influence is neglected 14

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 17. Comparison of the contour plots of strain field, ɛx, ɛy, and ɛxy: (a) predicted by the VEVPD-Modified model and (b) calculated by DIC results of experiment (Liu and Thompson, 2014) at moment C indicated in Fig. 13b.

4.3.2. Predictive capabilities of the VEVPD-Modified model Various experimental data, including the overall compressive load–displacement curve and contour plots of displacement field measured through DIC, as provided in Liu and Thompson (2014), are considered as validation data to the development of material constitutive models. Furthermore, the derived displacement field quantities, such as strain field, correlation coefficient field, and tensile crack growing history, are offered. The results simulated using the VEVPD-Modified model and compared with the tests are exhibited in Figs. 13 and 15–19. Fig. 13b presents the overall responses predicted using the VEVPDModified model and measured through the experiment. The selected moments labeled from A to E are provided in Liu and Thompson (2014), where D is the peak load point and B is the starting point of tensile crack initiation. The moments labeled from F to I are selected in accordance with the simulation. H is the peak load point, and F is the starting point of tensile crack initiation when the damage threshold of d = 0.9 is specified. The moments of Points C and G are identical; and moments of Points E and I are the same. In addition, the moments of A to E are used for the contour plots of the simulation demonstrated in Figs. 15–18 to compare with the experiments conveniently. In Fig. 13b, the result of the simulation significantly coincides with the experiment before the peak load. The curve drops slightly faster in modeling than in experiment at the later stage. The comparison results of the contour plots of displacement and strain fields at selected moments between the tests and modeling results, as displayed in Figs. 15–17, demonstrate that the results predicted using the VEVPD-Modified model are consistent with the experimental measurements. The contour lines at the top and bottom of the cavity tend to gather at Moment B and become dense along the center axis at Moments C, D, and E. The contour lines at the left and right of the cavity are slightly twisted along the oblique direction with an angle of approximately 45° at Moment C. This phenomenon is evident at

Moments D and E. These manifestations elucidate that the concentrated deformation region exhibits two types. These types are further indicated by the contour plots of the strain field at Moment C displayed in Fig. 17. The tensile zones of highly localized deformation along the center axis in the vertical direction and shear zones of extensive deformation on the side of the cavity arise on the surface of the specimen. The tensile zones appear earlier than the shear zones. Notably, the magnitudes of the horizontal displacement on the two sides of the plate near the cavity are different between the numerical and experimental results. At moment E, the maximum values are 0.59 and 0.8 mm in mm in Fig. 15a and b. One of the potential reasons for the distinction is represented as follows: brittle fracture will arise soon after the peak load under tension in experiment as shown in Fig. 4, and then small interaction force exists between the fracture surfaces. However, fracture criterion is not introduced over the simulation, and the VEVPD-Modified model exhibits progressive softening behavior. Therefore, the inner constraint stress is present to slow down the material fracturing process, even though the damage value d is sufficiently high; as shown in Fig. 6a, the tensile stress is 1.91 MPa when d is 0.9. Thus, the crack opening displacement along the center axis of the vertical direction is restricted, and the horizontal displacement of the plate by numerical simulation is smaller than that by experiment. Furthermore, the magnitude of the shear strain ɛxy by modeling in Fig. 17a is smaller than experimental results as well. Fig. 18a illustrates the contour plots of the damage field using the VEVPD-Modified model at selected moments from A to E. A tensile crack is already apparent at Moment A. The comparison of the curves of the tensile crack length versus displacement presented in Fig. 19 indicates that the tensile cracks rapidly expand at the early stage (F–G) and then gradually slow down in the next stage (G–I). The sheardominated damage zones start to form at Moment C (G), and the extension speed gradually accelerates. After the peak load Point H (near 15

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 18. Comparison of the contour plots of the (a) damage field predicted using the VEVPD-Modified model (d = 0.5 is shown as solid black lines), (b) correlation coefficient field calculated by DIC results of experiment (Liu and Thompson, 2014), and (c) damage field predicted using the VEVPD-Original model at selected moments depicted in Fig. 13b.

Point D), the shear-dominated damage zones rapidly spread. The overall compressive load versus displacement curve plotted in Fig. 13b exhibits that the rising speed of the overall load value slows down after Moment C (G), and the load drops after the peak load Point H. The load carrying capacity of the specimen is mainly influenced by the enlargement of shear-dominated damage regions. The damage evolution process can be further divided into three stages: tensile-dominated cracking stage (F–G), shear-dominated damage stage (H–I), and competition stage of tension and shear (G–H). Figs. 18a, 18c, and 19 demonstrate that the tensile cracks simulated using the VEVPD-Original model slowly develop. The final crack length is only 4.45 mm using the VEVPD-Original model, whereas the length reaches 20.80 mm using the VEVPD-Modified model. The same tendency is exhibited by the experimental results

displayed in Figs. 18b and 19. That is, two types of damage zones, namely, tensile cracks in a narrow and elongated shape prior to the maximum load (D) and shear-dominated extensive material damage starting to form at Moment C, are generated. The contour plots of the correlation coefficient field illustrated in Fig. 18b are derived through mathematical manipulation based on the DIC results. The correlation coefficient field demonstrates the degree of deformation concentration. However, Fig. 18a displays the damage degree of the material, that is, the magnitude of d = c 3/(a3 + c 3) . No one-to-one correspondence exists between the correlation coefficient and the damage fields given the nonlinear constitutive model. Tensile cracks evidently appear at Moment A in Fig. 18a, but Figs. 15a and 16a cannot reflect the phenomenon of deformation concentration. The growth rates of the top and bottom cracks of the cavity are 16

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

portion of the specimen is in a stationary state. Fig. 19 displays the crack length of the bottom portion only, which is accurate for comparison. Huang et al. (2017) investigated the stress distribution along the inner wall of the cavity before the onset of tensile cracks. The author indicated that the two endpoints above and below the cavity are subjected to local tensile stresses, whereas the left and right sides of the cavity are subjected to compressive and shear stresses. The stress distribution and damage evolution are further illustrated in Fig. 20 through Characteristic points A1, B1, and C1. Point A1, which initially bears compressive load (p > 0), is subjected to tensile load with the cracking procedure and finally bears an approximately stress-free state. Points B1 and C1 are subjected to compressive load during the entire loading procedure. The effective stress is larger at Point B1 than at Point C1 before softening while the pressure is low. Thus, the damage evolution is considerably faster at Point B1 than at Point C1, as depicted in Fig. 20b. As illustrated in Section 4.2.2, the growth rate of damage evolution is considerably faster under tension than under compression, and no plastic deformation is detected along the central vertical line of the specimen given the low stress level caused by cracking (Fig. 21). Considering the effect of stress concentration caused by the cavity, the pressure and effective stress are high in regions on the side of the cavity (Fig. 22). The evolution of plastic work is mainly governed by the effective stress. Thus, the distribution of the plastic work is along the shear direction of the two sides of the specimen, and the plastic work accumulates following the increment in the overall load.

Fig. 19. Comparison of the tensile crack growing history measured by the damage field of the VEVPD-Modified model, damage field of the VEVPD-Original model, and correlation coefficient field of the experiment (Liu and Thompson, 2014). The crack length is measured when d ≥ 0.9 for the simulation results.

different in the experiment. In Liu and Thompson (2014), the DIC results above the cavity were translated to justify the top loading surface, whereas those below the cavity were not translated because the bottom

Fig. 20. Comparison of curves predicted using the VEVPD-Modified model at Points A1, B1, and C1, as demonstrated in Fig. 12a: (a) effective stress versus pressure and (b) damage accumulation versus normalized displacement.

Fig. 21. Contour plots of the field of plastic work per unit volume predicted by the VEVPD-Modified model at selected moments as indicated in Fig. 13b. 17

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 22. Contour plots of pressure field p and effective stress field ¯ predicted by the VEVPD-Modified model at moment C as indicated in Fig. 13b.

the damage is progressively changed along the bar length. However, pathological mesh dependence appears after stress drops, and mesh convergence cannot be obtained by mesh refinement (Fig. 24e). The basic viscoplastic (rate-dependence) regularization mechanism is explained as follows (Niazi et al., 2013); when the deformation rate starts to increase in the softer element, the increase in strain rate stiffens the element again. Thus, the competitive result between the rate of damage evolution and viscoelasticity and viscoplasticity determines mesh dependence or not. The damage evolution rate of the proposed model over different loading paths given in Fig. 3 and Section 4.2 indicates that pathological mesh dependence will become noticeable, especially the negative pressure region, with the increment in damage growth rate from compression to tension. The effectiveness of the simulated results of the perforated plate can be illustrated as follows. On the one hand, compressive stress distributes over the entire specimen, except the vertical central line. From the discussion of mesh convergence (Section 4.3.1 and Fig. 13a) and the diffused shear-based damage with an equal width for different meshes (Fig. 25), we can deduce that the evolution of the shear-based damage is not (or seldom) influenced by pathological mesh dependence. On the other hand, the distribution of tensile cracks is parallel to the compressive load direction. Thus, tensile cracks do not contribute to the overall compressive load directly. In addition, the stress level under tension is very low after stress drops (Fig. 24d). Therefore, the influence caused by pathological mesh dependence under tension is limited.

Fig. 23. 1D bar subjected to tension and compression.

4.3.3. Discussion on mesh dependence In FEM simulations, local damage models are known to produce pathological mesh dependence. Thus, localization limiters are needed to regularize the solution. Introducing rate-dependence in the constitutive model is a technique that can regularize the mesh dependence of local damage models (Belytschko et al., 2014). Many studies, such as that of Niazi et al. (2013), have been conducted to address the effectiveness of rate-dependent (viscoplastic) regularization. The aforementioned previous study pointed out that the primary viscoplastic length scale is a function of hardening and softening parameters but does not depend upon the deformation rate. Given that the VEVPDModified model is a rate-dependent local damage model, the rate-dependent regularization of the proposed model is further checked by a test of 1D bar subjected to tensile and compressive loads in this section. Fig. 23 shows the dimensions of the 1D bar used in this study; the thickness of the bar is 1 mm. The bar is fixed at one end, and the horizontal displacement is prescribed at the other end. The width of the bar varies linearly from the fixed end to free end. The non-uniform width along the 1D bar will act as a geometric imperfection, and deformation will localize at the fixed end. Owing to the symmetrical property, one quarter of the bar is adopted for modeling to avoid buckling under compression. Coarse, medium, and fine meshes (with 5, 10, and 20 elements over the length, respectively) are used. The selective-reduced quadrature hexahedral element is adopted to eliminate volumetric locking, and the VEVPD-Modified model with parameters given in Tables 1–3 is used. The magnitude of strain rate is approximately 1 s−1 for tension and compression. The modeling results under compression are given in Figs. 24a–c. No pathological mesh dependence is observed even at moment E of the softening stage, and mesh convergence can be obtained by mesh refinement (Fig. 24b). Progressive change of damage along the bar length is exhibited, as shown in Fig. 24c. This phenomenon coincides with that under the non-uniform width of the bar. The modeling results under tension are given in Figs. 24d–f. Before stress drops, that is, moment C,

5. Conclusions In this work, an efficient viscoelastic–viscoplastic damage model is developed, and its applications for TATB-based high explosives under complex stress states are investigated. The main contributions of this work are summarized as follows: (1) The growth rate factor of pressure-dependent tensile damage is incorporated into the microcrack evolution equations of the ViscoSCRAM model, which is reasonable for characterizing the tension–compression asymmetry feature for PBXs. (2) The viscoplasticity of the B–P model is introduced into the framework of the ViscoSCRAM model to characterize the plastic behavior of PBXs. (3) An efficient numerical scheme is presented on the basis of an explicit Runge–Kutta algorithm and an implicit iterative radial corrector method to solve the equations of the proposed model. 18

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 24. Simulation of the 1D bar test: (a) Overall compressive load versus displacement curves, (compressive) damage distribution along the bar length of (b) moments C and E, and (c) different moments of the medium mesh as indicated in Fig. 24a; (d) overall tensile load versus displacement curves, (tensile) damage distribution along the bar length of (e) moments C and E, and (f) different moments of the medium mesh as indicated in Fig. 24d.

(4) The damage evolution procedure of the tensile cracks and diffused shear-based damage raised in sequence of the compressed perforated plate is accurately captured in comparison with the experiments that use the proposed model.

bulk response is assumed for simplification. Several issues still need to be addressed in the following investigation: (1) the irreversible volumetric response caused by viscoelasticity and/or viscoplasticity, which is common for granular materials; (2) the inelastic deformation induced by open cracks; and (3) the temperature effect originated from high loading rate. These issues are currently being investigated by the authors, and the related work will be reported in separate papers.

In the present model, the deviatoric deformation of viscoelastic, viscoplastic, and damage properties are major interests, and a linear 19

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Fig. 25. Contour plots of damage field of mesh sizes 0.9, 0.6, and 0.3 mm predicted by the VEVPD-Modified model at moment E indicated in Fig. 13b.

Acknowledgment

Chinese National Natural Science Foundation (Grant No. 11872119) for supporting this project. A part of this work has been supported by NSAF (Grant No. U183010034). The authors would like to thank the reviewers for their fruitful discussions and precious advices during the writing process.

The authors would like to thank the National Key R&D Program of China (2016YFB0201004), Science Challenging Program (TZ2016001), and Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmat.2019.103179. Appendix A. Derivatives for the plastic multiplier The derivatives required to set up the plastic multiplier λ as defined in Eq. (8) are shown. Based on implicit backward Euler method, the incremental integration form of the total deviatoric strain, viscoelastic strain, plastic strain, damage strain, and deviatoric stress are written as

eij = eijve + eijp + eijc,

2G0 eijve = Sij + eijp = t

=

Sijn

t

k=1

(k )

Sij(k ), n + 1,

(A.2)

n + 1S n + 1, ij

2G0 eijc =

Sijn+ 1

N

(A.1)

(A.3)

c, n + 1S n + 1 ij

c, nS n, ij

(A.4) (A.5)

+ Sij,

where eij = t· Dijdev , eijve = t ·Dijve , eijp = t ·Dijp , eijc = t· Dijc , 2G0 n + 1 tSijn + 1 + c, n+ 1Sijn+ 1 +

Sijn + 1 = Sijtr

N k= 1

c, n + 1

(c n + 1/ a)3 . Combining Eqs. (A.1)–(A.4) with (A.5) yields

t (k ), n + 1 S , (k ) ij

(A.6)

where the trial deviatoric stress Sijtr is defined as

Sijtr

Sijn + 2G0 eij +

c, nS n. ij

(A.7)

Eq. (A.6) is rearranged to give

Sijn+ 1 n+1 Seff

=

=

where tr S¯ij

1 tr S¯ , 1 + 2G 0 n + 1 t + c, n + 1 ij

(A.8)

1

tr S¯ , 1 + 2G 0 n + 1 t + c, n + 1 eff n+1 = Seff N k=1

Sijtr

3 n + 1 n + 1 ¯ tr S Sij , Seff 2 ij

(A.9)

=

3 ¯ tr ¯ tr S S , 2 ij ij

and

tr S¯ij

is defined as

t (k ), n + 1 S . (k ) ij

p = From Eq. (A.3), the incremental effective plastic strain eeff

p eeff =

2 3

t

(A.10) 2 3

eijp eijp can be written as

n + 1S n + 1. eff

(A.11)

Combining Eq. (A.9) with (A.11) gives 20

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

n+1 Seff =

tr S¯eff

p 3G0 eeff 1 + c, n + 1

.

(A.12)

With Eqs. (A.11) and (A.12), the plastic multiplier λ at time step tn + 1 can be obtained as n+ 1 =

p 3 e (1 + c, n + 1) eff p . 2 t (S¯ tr 3G 0 e ) eff eff

Inserting Eq. (A.13) into Eqs. (A.8) and (A.3), the deviatoric stress and the incremental plastic strain are given by Sijn + 1 =

p (S¯ tr 3G0 e ) tr eff eff ¯ Sij , (1 + c, n + 1) S¯ tr eff

(A.14)

and eijp =

(A.13)

p p 3 e 3 e eff ¯ tr eff n + 1 Sij = S , 2Seff ij 2S¯ tr eff

(A.15)

respectively. Appendix B. Reduced predictor-corrector solution scheme

The integration algorithms presented in Section 3 for the viscoelastic-viscoplastic damage model are summarized. Box. 1 shows the solution scheme for the viscoelastic damage predictor step. Box. 2 presents the residual computing for the plastic corrector step. Box. 3 gives the total solution procedure for the constitutive model of each time step. Box. 1 Solution scheme for the viscoelastic damage predictor step. 1. Compute K1 = Y (tn, Y n) : 1.1.

(k )

Sij

1.2.

Sij =

1.3.

p=

1.4.

c={

= 2G (k ) (Dijdev 2G0 (Dijdev

p Dij )

K · Dkk ; K v res ( I )m

v res [1

K1 K 0µ 2 ( ) ] KI

Sij(k ) (k )

Dijp)

Sij(k ) N k = 1 (k ) c 1 + ( )3 a

for KI < K otherwise

2. Compute K2 = Y (tn + 3. Compute K3 = Y (tn +

1 2 1 2

G (k ) c 3 c c [( ) Sij + 3( )2 Sij ]; G0 a a a

c c 3( )2 Sij a a

;

.

t, Yn +

t,

Yn

+

t K1) .

1 2 1 2

t K2) .

4. Compute K 4 = Y (tn + t , Y n + t K3) . 5. Compute Y n + 1 = Y n +

1 6

t (K1 + 2K2 + 2K3 + K 4) .

Box. 2 Residual computing for the plastic corrector step. 1. Deviatoric plastic strain rate Dijp = 2. Viscoelastic damage

p 3Deff

Sij .

2Seff

predictorY n + 1

= Yn +

1 6

t (K1 + 2K2 + 2K3 + K 4) , see Box. 1.

R = 3. Effective stress and its derivative by radial corrector Seff

tr S¯eff

p 3G 0 Deff t , 1 + c, n + 1

4. Effective stress and its derivative by BP model Seff = Z [ 2 ln( p,(k + 1) p,(k + 1) ) = Seff (Deff ) 5. Residual R (Deff n+1

Return: c n + 1, c n + 1, pn+ 1, pn + 1, Sijn+ 1, Sij

p 3 Deff

2D0

p,(k + 1) R Seff (Deff ).

, Sij(k ), n + 1, W p, n + 1.

21

)]

1 2h

R Seff

p Deff

,

=

Seff p . Deff

p 3G 0 Deff t

. 1 + c, n + 1

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

Box. 3 Solution scheme for the constitutive model. n

Given: t , Dijn + 1/2 , ijn , c n , c n , pn , Sij , Sij(k ), n , and W p, n . 1. Predictor and initial operations: n + 1/2 n =Rkin+ 1/2 Dkln+ 1/2 Rljn + 1/2 ; 1.1. Objective rate ^ij = Rkin kln Rljn , D^ij n n 1 1.2. Pressure and deviatoric stress pn = ^kk , Sijn = ^ij + pn ij ; 3

1.3. Viscoelastic damage predictor Y n + 1 =Y n + 6 t (K1 + 2K2 + 2K3 + K 4) , see Box. 1. p, n + 1 2. Newton-Raphson iteration for Deff : p,(0) dev = 0.5Deff 2.1. Initial guess Deff ; p,(0) p,(0) p,(0) R ) = Seff (Deff ) Seff (Deff ) , see Box. 2; 2.2. Residual R (Deff 2.3. Set k = 0 ; R p,(k + 1) p,(k) p,(k) = Deff [ p,(k) ] 1R (Deff ); 2.4. Iteration update Deff p,(k + 1) p,(k + 1) ) = Seff (Deff ) 2.5. Residual R (Deff

1

Deff

p,(k + 1) R Seff (Deff ) , see Box. 2;

p,(k) p,(k + 1) p,(k) )| < tol , then stop; Deff | < tol and |R (Deff 2.6. Convergence check if |Deff 2.7. Setk = k + 1 and go to 2.3. 3. Updates: p n+ 1 eeff 3.1. Plastic work W p, n + 1 = Seff ; n+1 n+1 n + 1 ^ p ; 3.2. ij = Sij n+ 1 3.3. Cauchy stress ijn + 1 = Rikn+ 1 ^kl Rjln+ 1 .

Return:

n+1 , ij

n+1

c n + 1, c n + 1, pn + 1, Sij

, Sij(k ), n + 1, and W p, n + 1.

Appendix C. Sensitivity analysis To get a better understanding of the effects of the parameters on the model response, a sensitivity analysis is performed by means of a perturbation method. The responses derived from uniaxial tension (UT), uniaxial compression (UC), biaxial tension (BT), biaxial compression (BC), and biaxial shear (BS) simulations of Sections 4.2.1 and 4.2.2, are analyzed, and the same loading states and parameters are adopted. Some necessary concepts are defined firstly. The parameter set p ≡ {c1, ⋅⋅⋅, ci, ⋅⋅⋅, cN} stands for the total parameters used, where N is 23. The 1%

Table C.1 Mechanism global sensitivities. Mechanisms

Viscoelasticity

Viscoplasticity

Damage

SIGlobal [%]

13.08

47.47

39.45

Table C.2 Parameter relative sensitivities for viscoelasticity. Parameters

G∞

G(1)

G(2)

G(3)

G(4)

τ(1)

τ(2)

τ(3)

τ(4)

SiVE, rel [%]

21.20

13.98

33.56

8.26

0.11

2.41

15.49

4.97

0.02

Table C.3 Parameter relative sensitivities for viscoplasticity. Parameters

Z0

Z1

D0

h

m1

SiVP, rel [%]

13.67

10.19

1.07

73.84

1.23

Table C.4 Parameter relative sensitivities for damage. Parameters

a

c0

vmax

K0

m

µc

µt

va

vb

SiD, rel [%]

5.45

9.96

0

38.78

18.85

5.25

7.65

6.63

7.43

22

Mechanics of Materials 139 (2019) 103179

M. Liu, et al.

{c1, , 1.01ci, , cN } and pi {c1, , 0.99ci, , cN } , separately. The difpositive and negative perturbations to the ith parameter are defined as p+i pert pert ref ref (p), eff (p)) and a perturbed solution (Seff ference si (p , pi) between the reference solution (Seff (pi), eff (pi)) is quantified as follows si (p , pi) =

pert ref (p, t )dt (pi, t ) eff Seff ref (p, t ) ref (p, t )dt Seff eff

ref (p, t ) Seff UT,UC,BT, BC, and BS

× 100%.

Then the differences for positive and negative perturbations are defined as si+

(C.1)

si+ (p , p+i ) and si si (p, pi ) , respectively. The deviatoric deformation mechanisms of viscoelasticity (VE), viscoplasticity (VP), and damage (D) are major interests of this study, thus the effective stress Seff (in Eq. (c.1)) is adopted for sensitivity analysis. The relative sensitivity of each parameter for each deformation mechanism is defined as SiI , rel = (si+ + si )/

si+ + I

si

× 100%,

I

where I is VE, VP, or D. Similarly, the global sensitivity of each deformation mechanism is quantified as

SIGlobal =

si+ + I

si+ +

si / I

VE,VP, and D

si

× 100%.

(C.2)

(C.3)

VE,VP, and D

The results of global and relative sensitivities are listed in Tables C.1–C.4. It can be seen that the global sensitivity of viscoplastic parameters have a major role in the model, and the damage parameters represent an important part too. The relative sensitivity of parameter h of viscoplasticity is very high, and it should be carefully calibrated. The relative sensitivity of parameter vmax is 0 for that vres (in Eq. (22)) is less than vmax over the simulation. References Addessio, F.L., Johnson, J.N., 1990. A constitutive model for the dynamic response of brittle materials. J. Appl. Phys. 67, 3275–3286. Ambos, A., Willot, F., Jeulin, D., Trumel, H., 2015. Numerical modeling of the thermal expansion of an energetic material. Int. J. Solids Struct. 60-61, 125–139. Barua, A., Zhou, M., 2013. Computational analysis of temperature rises in microstructures of HMX-Estane PBXs. Comput. Mech. 52, 151–159. Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K.I., 2014. Nonlinear Finite Elements For Continua and Structures. Bennett, J.G., Haberman, K.S., Johnson, J.N., Asay, B.W., Henson, B.F., 1998. A constitutive model for the non-shock ignition and mechanical response of high explosives. J. Mech. Phys. Solids 46, 2303–2322. Bodner, S., 2002. Unified Plasticity for Engineering Applications. Kluwer Academic/Plenum Publishers, New York. Bodner, S.R., Partom, Y., 1975. Constitutive equations for elastic-viscoplastic strain-hardening material. J. Appl. Mech. 42. Buechler, M.A., Luscher, D.J., 2014. A semi-implicit integration scheme for a combined viscoelastic-damage model of plastic bonded explosives. Int. J. Numer. Methods Eng. 99, 54–78. Cook, W.H., Rajendran, A.M., Grove, D.J., 1992. An efficient numerical implementation of the Bodner–Partom model in the EPIC-2 code. Eng. Fract. Mech. 41, 607–623. Dienes, J.K., Zuo, Q.H., Kershner, J.D., 2006. Impact initiation of explosives and propellants via statistical crack mechanics. J. Mech. Phys. Solids 54, 1237–1275. Dienes, J.K., 1996. A unified theory of flow, hot spots, and fragmentation with an application to explosive sensitivity. In: Davison, L., Grady, D.E., Shahinpoor, M. (Eds.), High-Pressure Shock Compression of Solids II. Springer Verlag, New York, pp. 366–398. Drodge, D.R., Williamson, D.M., 2016. Understanding damage in polymer-bonded explosive composites. J. Mater. Sci. 51, 668–679. Flanagan, D.P., Taylor, L.M., 1987. An accurate numerical algorithm for stress integration with finite rotations. Comput. Methods Appl. Mech. Eng. 62, 305–320. Frank, G.J., Brockman, R.A., 2001. A viscoelastic-viscoplastic constitutive model for glassy polymers. Int. J. Solids Struct. 38, 5149–5164. Gratton, M., Gontier, C., Rja Fi Allah, S., Bouchou, A., Picart, D., 2009. Mechanical characterisation of a viscoplastic material sensitive to hydrostatic pressure. Eur. J. Mech. A/Solids 28, 935–947. Hackett, R.M., Bennett, J.G., 2000. Implicit finite element material model for energetic particulate composite materials. Int. J. Numer. Methods Eng. 49, 1191–1209. Huang, X.C., Li, S.K., Wei, Q., Tian, R., Chen, C.J., Wang, L.X., Liu, M., 2017. Analysis of crack initiation and growth in PBX energetic material using XFEM-based cohesive method. Chin. J. Energ. Mater. 25, 694–700. LaBarbera, D.A., Zikry, M.A., 2015. Dynamic fracture and local failure mechanisms in heterogeneous RDX-Estane energetic aggregates. J. Mater. Sci. 50, 5549–5561. Le, V.D., Gratton, M., Caliez, M., Frachon, A., Picart, D., 2010. Experimental mechanical characterization of plastic-bonded explosives. J. Mater. Sci. 45, 5802–5813. Liu, C., Thompson, D.G., 2014. Crack initiation and growth in PBX 9502 high explosive subject to compression. J. Appl. Mech. 81, 101004. Liu, R., Chen, P.W., 2018. Modeling ignition prediction of HMX-based polymer bonded explosives under low velocity impact. Mech. Mater. 124, 106–117. Niazi, M.S., Wisselink, H.H., Meinders, T., 2013. Viscoplastic regularization of local damage models: revisited. Comput. Mech. 51, 203–216. Pyrz, M., Zairi, F., 2007. Identification of viscoplastic parameters of phenomenological constitutive equations for polymers by deterministic and evolutionary approach. Model. Simul. Mater. Sci. Eng. 15, 85–103. Pyrz, M., Zalewski, R., 2015. Modelling of the viscoplastic behaviour of homogeneous solid propellants. Cent. Eur. J. Energ. Mater. 12, 159–174. Rangaswamy, P., Thompson, D.G., Liu, C., Lewis, M.W., 2010. Modeling the mechanical response of PBX-9501. In: Proc. 14th Int. Detonation Symp. Coeur d'Alene. Idaho, USA. 11-16 April. Reaugh, J.E., Jones, A.G., 2010. Mechanical damage, ignition, and burn: experiment, model development, and computer simulations to study high-explosive violent response (HEVR). In: Proc. 14th Int. Detonation Symp. Coeur d'Alene. Idaho, USA. 11-16 April. Sun, B.P., Duan, Z.P., Wan, J.L., Liu, Y., Ou, Z.C., Huang, F.L., 2015. Investigation on ignition of an explosive charge in a projectile during penetration based on Visco-Scram model. Explosion Shock Waves 35, 689–695. Tang, W., 2016. Constitutive models and strength criterions of a TATB-based PBX under quasi-static loading. PhD thesis. Nanjing University of Science & Technology. Thompson, D.G., Brown, G.W., Olinger, B., Mang, J.T., Patterson, B., Deluca, R., Hagelberg, S., 2010. The effects of TATB ratchet growth on PBX 9502. Propellants, Explos. Pyrotech. 35, 507–513. Thompson, D.G., Deluca, R., Brown, G.W., 2012. Time-temperature analysis, tension and compression in PBXs. J. Energ. Mater. 30, 299–323. Wang, X., Wu, Y., Huang, F., Jiao, T., Clifton, R.J., 2016. Mesoscale thermal-mechanical analysis of impacted granular and polymer-bonded explosives. Mech. Mater. 99, 68–78. Wiegand, D.A., Reddingius, B., 2005. Mechanical properties of confined explosives. J. Energ. Mater. 23, 75–98. Wiegand, D.A., Redingius, B., Ellis, K., Leppard, C., 2011. Pressure and friction dependent mechanical strength - Cracks and plastic flow. Int. J. Solids Struct. 48, 1617–1629. Xiao, Y., Sun, Y., Zhen, Y., Guo, L., Yao, L., 2017. Characterization, modeling and simulation of the impact damage for polymer bonded explosives. Int. J. Impact Eng. 103, 149–158. Yang, K., Wu, Y., Huang, F., 2018. Numerical simulations of microcrack-related damage and ignition behavior of mild-impacted polymer bonded explosives. J. Hazard. Mater. 356, 34–52. Williamson, D.M., Siviour, C.R., Proud, W.G., Palmer, S.J.P., Govier, R., Ellis, K., Blackwell, P., Leppard, C., 2008. Temperature-time response of a polymer bonded explosive in compression (EDC37). J. Phys. D 41, 085404. Wu, Y.Q., Huang, F.L., 2009. A micromechanical model for predicting combined damage of particles and interface debonding in PBX explosives. Mech. Mater. 41, 27–47. Zaïri, F., Naït-Abdelaziz, M., Gloaguen, J.M., Lefebvre, J.M., 2008. Modelling of the elasto-viscoplastic damage behaviour of glassy polymers. Int. J. Plast. 24, 945–965. Zuo, Q.H., Addessio, F.L., Dienes, J.K., Lewis, M.W., 2006. A rate-dependent damage model for brittle materials based on the dominant crack. Int. J. Solids Struct. 43, 3350–3380. Zuo, Q.H., Disilvestro, D., Richter, J.D., 2010. A crack-mechanics based model for damage and plasticity of brittle materials under dynamic loading. Int. J. Solids Struct. 47, 2790–2798.

23