Simulation study on cavity growth in ductile metal materials under dynamic loading

Abstract

Cavity growth in ductile metal materials under dynamic loading is investigated via the material point method. Two typical cavity effects in the region subjected to rarefaction wave are identified: (i) part of material particles flow away from the cavity in comparison to the initial loading velocity, (ii) local regions show weaker negative or even positive pressures. Neighboring cavities interact via coalescence of isobaric contours. The growth of cavity under tension shows staged behaviors. After the initial slow stage, the volume and the dimensions in both the tensile and transverse directions show linear growth rate with time until the global tensile wave arrives at the upper free surface. It is interesting that the growth rate in the transverse direction is faster than that in the tensile direction. The volume growth rate linearly increases with the initial tensile velocity. After the global tensile wave passed the cavity, both the maximum particle velocity in the tensile direction and the maximum particle velocity in the opposite direction increase logarithmically with the initial tensile speed. The shock wave reflected back from the cavity and compression wave from the free surface induce the initial behavior of interfacial instabilities such as the Richtmyer-Meshkov instability, which is mainly responsible for the irregularity in the morphology of deformed cavity. The local temperatures and distribution of hot spots are determined by the plastic work. Compared with the dynamical process, the heat conduction is much slower.

Graphical abstract

Keywords

material point method / cavity growth / dynamic loading / interfacial instability

1 Introduction

Failure of ductile metal materials under dynamic loading is an important and fundamental issue in the fields of science and technology. The failure process is complicated because it couples various physical and mechanical mechanisms in the microscopic, mesoscopic and macroscopic scales. Globally speaking, spallation or fragmentation of metal material is mainly composed of the following typical stages, nucleation, growth and coalescence of microscopic voids and/or larger scale cavities. There have been extensive studies on the quasistatic growths of voids and cavities. The dynamical growth behaviors are much more complicated and far from being well understood [1].

In 1972 Carroll and Holtz [2] studied the static and dynamic cavity-collapse relations for ductile porous materials and found that the compression effect on the cavity growth is not pronounced when the material is not sensitive to the loading rate. The research was extended to the visco-plastic materials by Johnson [3] in 1981. In 1987 Becker [4] numerically analyzed the effect of a nonuniform distribution of porosity on flow localization and failure in a porous material. The void density distribution and properties used to characterize the material behavior were obtained from measurements on partially consolidated and sintered iron powder. The calculations were carried out using an elastic viscoplastic constitutive relation for porous plastic solids. Local material failure is incorporated into the model through the dependence of the flow potential on void volume fraction. The region modeled is a small portion of a larger body under various stress conditions. Both plane strain and axisymmetric deformations are considered with imposed periodic boundary conditions. It was found that interactions between regions with higher void fractions promote plastic flow localization into a band, and that local failure occurs via void growth and coalescence within the band. The results of this study suggested a failure criterion based on a critical void volume fraction that is only weakly dependent on stress history. The critical void fraction depends on the initial void distribution and material hardening characteristics. In 1992 Ortiz and Molinari [5] studied the effect of strain hardening and rate sensitivity on the dynamic growth of a void in a plastic material and pointed out that the inertial effect, hardening effect, loading rate effect can significantly influence the void growth. The studies of Benson [6] in 1993 and of Ramesh and Wright [7] in 2003 showed that the inertia effect is responsible for stable growth of the cavity. In 1998 Pardoen et al. [8] investigated the ductile fracture of round copper bars within the scope of the local approach methodology. Two damage models, the Rice–Tracey model and the Gurson–Leblond–Perrin model, were analyzed. Four coalescence criteria, (i) a critical value of the damage parameter, (ii) the Brown and Embury criterion, (iii) the Thomason criterion and (iv) a criterion based on the reaching of the maximum von Mises equivalent stress in a Gurson type simulation, were comparatively studied. Ellipsoidal void growth and void interaction were accounted for. As far as possible, all the parameters of the models were identified from experiments and physical observations. The effect of stress triaxiality was studied using specimens presenting a wide range of notch radii. The effect of strain-hardening was analyzed by comparing the behaviour of the material in the cold drawn state and in the annealed state. In 2000 Pardoen et al. [9] proposed an extended model for void growth and coalescence. This model integrated two existing contributions, the Gologanu–Leblond–Devaux model extending the Gurson model to void shape effects and the Thomason scheme for the onset of void coalescence. Each of these was extended heuristically to account for strain hardening. In addition, a micromechanically-based simple constitutive model for the void coalescence stage was proposed to supplement the criterion for the onset of coalescence. The fully enhanced Gurson model depends on the flow properties of the material and the dimensional ratios of the void-cell representative volume element. It incorporates the effect of void shape, relative void spacing, strain hardening, and porosity. In 2001 Orsini et al. [10] developed an inelastic rate-dependent crystalline constitutive formulation and specialized computational schemes to obtain a detailed understanding of the interrelated physical mechanisms which can result in ductile material failure in rate-dependent porous crystalline materials subjected to finite inelastic deformations. Results of this study are consistent with experimental observations that ductile failure can occur either due to void growth parallel to the stress axis, which results in void coalescence normal to the stress axis, or void interaction along bands, which are characterized by intense shear-strain localization and that intersect the free surface at regions of extensive specimen necking. In 2002 Tvergaard and and Hutchinson[11] discussed two mechanisms of ductile fracture, void by void growth versus multiple void interaction, Zohdi et al. [12] discussed the plastic flow in porous material.

Currently, most of the studies on cavity/void growth are focused on their relevance on macroscopic behaviors [2–12]. The quantitative relations are obtained by fitting experimental results. Those studies do not reveal or indicate the underlying idiographic physical and mechanical mechanisms of cavity/void growth. Cavity coalescence is the final stage of spallation developed from mesoscopic scale to macroscale [13]. It is also the least-known stage [14–21]. Continuous damage mechanics adopts fluid or solid description supplemented by damage modeling. The damage is generally modeled by an internal variable. The internal variable is defined by the variation of some mechanical behavior and is not dynamically relevant to the particular structures.

The molecular dynamics simulations [14, 22–25] can help understand some mechanisms from the atomic scale, but the temporal and spatial scales it can access are too small to be comparable with experiments. The Material Point (MP) method [26–29] is a newly developed mesoscopic particle method in the field of computational solid mechanics. In this method, the continuum bodies are discretized with N material particles. Each material particle carries the information of position, velocity, temperature, mass, density, Cauchy stress, strain tensor and all other internal state variables necessary for the constitutive model. At each time step, calculations consist of two parts: a Lagrangian part and an Eulerian one. Firstly, the material particles flow with the body, and is used to determine the strain increment, and the stresses in the sequel. Then, the velocity field is mapped from the particles to the Eulerian mesh nodes. The spatial derivatives are calculated and the momentum equation is integrated. The velocity and acceleration fields are mapped back to update those of the particles [30, 31]. The MP method not only takes advantages of both the Lagrangian and Eulerian algorithms but avoids their drawbacks as well. Since using Eulerian background grid, it is more stable and has a higher computational efficiency than the meshless smooth particle hydrodynamics (SPH) [29].

From the physics side, the MP method is based on continuum medium description and designing of contact force. It has been extensively used to simulate the complex dynamical behaviors of shock wave interaction on inhomogeneous materials [32–37]. The mesoscopic MP simulation can be further used to investigate the growth and evolution of defect structures such as cavities and cracks in the scales of micron and larger. Such investigations may present indicative results for improving physical modeling of fracture in larger scales. From the simulation side, in the MP method the continuous portion and the cavities are considered separately. So, it is convenient to set the particular structures according to our need and convenient to obtain the concrete information on the shapes, sizes, connectivity of relevant structures and their influences on surrounding materials. In other words, it is convenient to recover with more fidelity the physical processes of damage and failure. Simulation results may work as theoretical bases for the physical modeling of damnification. Different from the phenomenological quasistatic analysis, the MP simulation results contain intrinsically the inertial effects.

2 Theoretical model of the material

We assume that the material follows an associative von Mises plasticity model with linear kinematic and isotropic hardening [38]. The stress and strain tensors, \sigmaand \epsilon, reads \sigma =s-PI, P=-\mathrm{Tr}(\sigma )/3, \epsilon =e+\theta I/3,\theta =\mathrm{Tr}(\epsilon )/3, where P is the pressure scalar, s the deviatoric stress tensor, and e the deviatoric strain. The strain e is generally decomposed as e=e^e+e^p, where e^e and e^p are the traceless elastic and plastic components, respectively. The material shows a linear elastic response until the von Mises yield criterion is reached. The yield {\sigma }_Y increases linearly with the second invariant of the plastic strain tensor e^p, i.e., {\sigma }_Y={\sigma }_{Y0}+E_{\mathrm{Harden}}\Vert e^p\Vert, where {\sigma }_{Y0} is the initial yield stress and E_{\mathrm{Harden}} the hardening coefficient. The deviatoric stress s is related to the Young’s module Y, the Poisson’s ratio \nu and the traceless elastic stress tensor by s=e^{e}Y/(1+\nu ). The shock speed U_s and the particle speed U_p after the shock follows a linear relation, U_s=c_0+\lambda U_p, where c_0 is the sound speed and \lambda a characteristic coefficient of material. The pressure P is calculated by the following Mie–Grüneisen state of equation [39]:

\begin{array}{}\text{(1)}& P-P_{\mathrm{H}}={\displaystyle \frac{\gamma (V)}{V}}[E-E_{\mathrm{H}}(V_{\mathrm{H}})]\end{array}

where P_{\mathrm{H}}, V_{\mathrm{H}} and E_{\mathrm{H}} are pressure, specific volume and energy on the Rankine–Hugoniot curve, respectively. The relation between P_{\mathrm{H}} and V_{\mathrm{H}} can be written as

\begin{array}{}\text{(2)}& P_{\text{H}}=\{\begin{array}{l}\frac{{\rho }_0{c}_0^2\left(1-\frac{V_{\text{H}}}{V_0}\right)}{{\left(\text{λ}-1\right)}^2{\left(\frac{\text{λ}}{\text{λ}-1}\times \frac{V_{\text{H}}}{V_0}-1\right)}^2},\begin{array}{cc}& V_{\text{H}}\le V_0\end{array}\\ {\rho }_0{c}_0^2\left(\frac{V_{\text{H}}}{V_0}-1\right),\begin{array}{cccc}& & & V_{\text{H}}>V_0\end{array}\end{array}\end{array}

The increment of specific internal energy E-E_{\mathrm{H}}(V_{\mathrm{H}}) is taken as the plastic energy. Both the shock compression and the plastic work contribute to the increasing of temperature. The temperature increment from shock compression is calculated by

\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{T}_{\mathrm{H}}}{\mathrm{d}{V}_{\mathrm{H}}}}={\displaystyle \frac{c_0^2\cdot \lambda {(V_0-V_{\mathrm{H}})}^2}{c_v{\left[(\lambda -1)V_0-\lambda V_{\mathrm{H}}\right]}^3}}-{\displaystyle \frac{\gamma (V)}{V_{\mathrm{H}}}}{T}_{\mathrm{H}}\end{array}

where c_v is the specific heat. The increasing of temperature from plastic work is \mathrm{d}{T}_p=\mathrm{d}{W}_p/c_v.

In this work we choose aluminum as the sample material. The corresponding parameters are as below: initial material density in the solid portion {\rho }_0=2700 kg/m³, Y=69 GPa, \nu =0.33, {\sigma }_{Y0}=120 MPa, E_{\mathrm{Harden}}=384 MPa, c_0=5.35 km/s, \lambda =1.34, c_v=880 J/(kg\cdotK), k=237 W/(m\cdotK) and {\gamma }_0=1.96 when the pressure is below 270 GPa. The initial temperature of the material is 300 K.

3 Simulation results and physical interpretation

3.1 Global scenario

In our simulations the body of aluminum material with cavity is connected with a rigid wall fixed at the bottom with the coordinate z=0. The simulated body is located within the volume, [-20,20]\times [-20,20]\times [0,50] with the length unit μm. Initially, a spherical cavity with radius r = 5 μm is located at the position (0,0,z) within the material body. At the time t=0 the material body with cavity starts to move upward at the velocity v_{z0}. Thus, the rarefaction wave or tensile wave occurs at the plane with z=0. The rarefaction wave propagates upwards within the material body. In our MP simulations, the mesh size is 1 μm and the diameter of the material particle is 0.5 μm. Periodic boundary conditions are used in the horizontal directions and free boundary condition is used in the upper surface of the material body. The rigid wall is assumed to be the same kind of material with the material body. Therefore, no special treatment of contact surfaces is needed in the MP simulations.

Figure 1 shows the snapshots of configurations with v_z field at four different times for the case with z=10 μm and initial v_{z0}=100 m/s. Figures 1(a)–(d) correspond to the times t=0.8 ns, 1.2 ns, 2.0 ns and 3.0 ns, respectively. The contours for v_z=0 are shown in the plots. Since no material particles are located within the cavity, the velocities at the nodes within the cavity are equal zero. Before the arrival of the global rarefaction wave, the upper contour with v_z=0 presents the initial morphology of the cavity. In Figs. 1(a)–(d) the moving upwards of the lower v_z=0 contour shows the propagation of rarefaction wave. It is clear from Fig. 1(a) that the lower v_z=0 contour is approaching the lower boundary of the cavity at the time t=0.8 ns. The velocities of particles below the cavity had begun to decrease before t=0.8 ns. Below the lower v_z=0 contour some material particles show negative velocities. With propagating of the rarefaction wave, the lower v_z=0 contour begins to get connection with that corresponding to the cavity. When the rarefaction wave arrives at the cavity, compression wave is reflected back. Under the action of the reflected compression wave, more material particles show negative velocities and their amplitudes continue to increase.[See Figs. 1(b) and (c).] The deformation rate of the cavity is slower compared with the propagation speed of rarefaction wave surrounding the cavity. After passing the cavity the surrounding rarefaction waves begin to converge. Therefore, stronger negative pressure appears on the top of the cavity. Material particles on the top of the cavity are accelerated by the upward stresses. At the time t=3 ns some material particles show velocities larger than 100 m/s. See Fig. 1(d).

Figures 2(a) and (c) show the configurations with pressure and velocity fields within the plane x=0 at times t=2 ns and t=3 ns, respectively. To investigate the amplitudes of particle velocities we show the distribution of v_z along the vertical direction in the other plots of Fig. 2. Figures 2(b) and (d) correspond to Figs. 2(a) and (c), respectively. At the times t=2 ns and t=3 ns, the maximum downward or minimum particle velocities are -230 m/s and -300 m/s, respectively. From Fig. 2(d) we can observe the vertical distribution of material particles with velocities larger than 100 m/s. In Figs. 2(a) and (c) the color from blue to red corresponds to the increase of pressure. From Figs. 2(a) and (c), one can observe the deformation of the cavity under tensile loading. The irregularities in the morphology of the cavity result from the following three aspects. (i) The initial cavity represented by the placed particles is not strictly spherical. (ii) The shock waves reflected back from the cavity induce the well-known Richtmyer–Meshkov (RM) hydrodynamic interfacial instability. The RM instability is the main mechanism for the initial irregularities of the deformed cavity. (iii) Compared with the dimension of the cavity, the mesh size is not small. For point (iii), it should be commented that if decrease the mesh unit, the body size can be simulated becomes smaller. We have to make compromise between the simulated body size and the mesh unit. It should also be pointed out that the practical cavities in materials are generally not strict spherical, which is qualitatively accordance with the simulated one.

With increasing of upward stress on the top of the cavity, the accelerations and velocities of particles within this region become larger. At the time t=7.2 ns, the global rarefactive wave arrives at the upper free surface. The maximum velocity of particles on the top of the cavity is about 430 m/s. At this moment, there exists a region where the particles have large downward velocities below the cavity. The largest downward velocity is about -325 m/s. In the plot of v_z versus z, there is a valley between the peak and the rarefaction wave front. The smallest particle velocity is about 6 m/s. When the rarefaction wave arrives at the upper free surface and compression wave is reflected back. Within the region scanned by the reflected compression wave, material particles obtain downward accelerations. Several characteristics are typical for the unloading of rarefaction wave and reflecting back of compression wave. The first is the decreasing of velocities of material particles representing the upper free surface. The second is that the valley continues to move toward the upper free surface. The third is that the maximum velocity between the valley and the cavity continues to increase. At the same time, the region with maximum downward particle velocity moves toward the bottom. Since we use periodic boundary conditions in the horizontal directions, the simulation results for the case with single cavity are also indicative for interaction of neighboring cavities. From the pressure field, it is clear that the negative pressures within regions among the neighboring cavities are weaker. The contours of negative pressure with small amplitudes get connection. The strength of compression wave reflected back from the cavity increases with increasing the strength of rarefaction wave. Consequently, local positive pressures occur among the neighboring cavities. The occurrence of positive pressures within the region scanned by the global rarefaction wave is a typical cavity effect.

Before the reflected compression wave arrive at the cavity, the deformation of the cavity is still controlled by the tensile loading. Below, we discuss the pressure distributions within the material at two times, 9 ns and 11 ns. Figure 3 shows the configurations with pressure field at the time 9 ns. The pressure contours in Figs. 3(a)–(f)are for -300 MPa, -350 MPa, -400 MPa, -450 MPa, -500 MPa and -550 MPa, respectively. Figure 3 shows that the contours around the cavities for pressure lower than -300 MPa are connected. The neighboring cavities get interaction via the connection of pressure contours. At the time 9 ns, there is still no positive pressure occur among the neighboring cavities. Figure 4 shows various pressure contours at the time 11 ns. The contours in Fig. 3(a)–(f) are for 0 MPa, -50 MPa, -100 MPa, -150 MPa, -200 MPa and -250 MPa, respectively. Pressure distribution around the cavity is as below. (i) The pressure surrounding the cavity is zero. (ii) With increasing the pressure, the corresponding contour moves away from the cavity and its surface area becomes larger. (iii) Among cases shown in the figure the contour for -150 MPa has the maximum area. If further increase the pressure, the contour area becomes smaller. Pressure distribution between the cavity and rigid wall is as below. There are four regions around the cavity show positive pressure. The pressure contours for -100 MPa between the nearest cavities are connected. The contours for -150 MPa have a higher connectivity. All contours for -200 MPa, -250 MPa, etc. are connected. The pressure distribution on the top of the cavity is as below. The region with the highest pressure does not locate above the cavity but above the middle of neighboring cavities. Since the rarefaction wave propagates more quickly within the solid region, the wave firstly arrives at the upper free surface and get reflected. The weaker the negative pressure, the closer to the upper free surface the corresponding contour, and the more planar the corresponding contour.

3.2 Morphology versus tensile strength

Figure 5 shows the configurations with temperature field at the time 6 ns. The contours in (a)–(d) are for 310 K, 320 K, 330 K and 340 K, respectively. Compared with the dynamical process, the thermal process is much slower. The temperature and distribution of hot-spots are mainly determined by the corresponding plastic work.

Since the rarefaction wave propagates in the sound speed, all rarefaction waves reach the upper free surface at the same time. With increasing the tensile strength, the growth rate of the cavity increases. Figure 6 shows the evolution of the cavity morphology. Figure 6(a) shows the cavity volume versus time. The points are simulation results and the lines are plotted to guide the eyes. The sizes of the initial tensile velocity v_{z0}, 100, 200, 400 and 1000, are shown in the legend. The unit is m/s. An enlarged portion of the curve for v_{z0}=100 m/s is shown in the inset. The growth of cavity can be described by the following stages: (i) initial slow growth stage, (ii) linear growth stage which ends when the global rarefaction wave arrives at the upper free surface, (iii) slower growth stage which ends when the reflected compression wave arrives at the cavity, (iv) quicker growth stage and (v) linear growth stage. Figure 6(b) shows the evolutions of the cavity dimensions in Horizontal (H) and Vertical (V) directions. The points are simulation results and the lines are plotted to guide the eyes. It is interesting to observe that the growth in horizontal direction is quicker than that in vertical direction. Such a mechanism is equivalent to the “necking effect” in macroscale. There exists also a linear stage in the growths of cavity dimensions. The growth rates increase with increasing the strength of tensile loading. Figure 6(c) shows the initial linear growth rate of cavity volume versus initial strength of tensile loading v_{z0}. The points are for the slopes of fitting lines in (a) for the first linear growth stage, and the line are linear fitting result for the points. It is clear that within the checked range the volume growth rate linearly increases with the initial tensile velocity v_{z0}.

Figures 7(a) and (b) show the density fields of the material at two times, 7.2 ns and 12 ns. Figures 7(c) and (d) show the corresponding pressure fields. The lower boundary of the cavity gradually becomes planar and parallel to the rigid wall.

3.3 Energy transformation versus tensile strength

For the case with uniform material, during the tensile loading, kinetic energy of the material transforms gradually to elastic potential energy and plastic work. Those energies distribute uniformly in planes parallel to the rigid wall. Although the material is three-dimensional, the dynamical and thermodynamical process is in fact one-dimensional. But for the case with cavities, the situation becomes much more complex. Figures 8(a) and (b) show the configurations with temperature field at the same two times as in Fig. 7, from which, besides the cavity morphology, we can understand better the energy transformation from kinetic to thermal. There is a high temperature layer surrounding the deformed cavity. That is because the plastic work by the stresses is pronounced in that region. Figures 8(c) and (d) show the configurations with v_z field at the same two times. With the reflecting back of compression wave from the upper free surface, the distribution range of high particle velocity becomes narrower.

Figure 9 shows the maximum upward particle velocity, v_{z\max }, above the cavity and maximum downward particle velocity, v_{z\min }, below the cavity versus the initial tensile velocity v_{z0}. The points are simulation results and the lines are fitting results. It is interesting to observe that both v_{z\max } and |v_{z\min }| logarithmically increase with the initial tensile velocity v_{z0}.

4 Conclusions

A three-dimensional material point simulation study on cavity growth in metal materials subjected to dynamic loading is conducted. Interactions of rarefaction wave with an existing cavity and the ultimate interactions of the cavity with its periodic images are carefully investigated. During the tensile loading procedure, some material particles below the cavity show high speeds in the opposite direction. Within the region subjected to the global rarefaction wave some local regions may show positive pressures. Neighboring cavities get interaction via the coalescence of isobaric contours. The deformation of cavity under tensile loading shows staged behaviors. After the initial slow growth stage, the volume and the dimensions in both the tensile and transverse directions show linear growth rate with time until the global tensile wave reaches the upper free surface. The growth rate in the tensile direction is slower than that in the transverse direction. The volume growth rate linearly increases with the initial tensile velocity. After the global tensile wave passed the cavity, the maximum particle velocity in the tensile direction and the maximum particle velocity in the opposite direction increase logarithmically with the initial tensile speed. The shock wave reflected back from the cavity and compression wave from the free surface induce the initial behavior of interfacial instabilities such as the Richtmyer–Meshkov instability, which is mainly responsible for the irregularity of the cavity morphology. Temperature and distribution of hot spots are determined by the plastic work. Compared with the dynamical process, the heat conduction is much slower.

References

[1] G. T. Gray, P. J. Maudlin, L. M. Hull, Q. K. Zuo, and S. R. Chen, Predicting material strength, damage, and fracture the synergy between experiment and modeling, J. Fail. Anal. Prev., 2005, 5(3): 7

[2] M. M. Carroll, and A. C. Holt, Static and dynamic Pore-Collapse relations for ductile porous materials, J. Appl. Phys., 1972, 27(3): 1626

[3] J. N. Johnson, Dynamic fracture and spallation in ductile solids, J. Appl. Phys., 1981, 52(4): 2812

[4] R. Becker, The effect of porosity distribution on ductile failure, J. Mech. Phys. Solids, 1987, 35(5): 577

[5] M. Ortiz and A.  Molinari, Effect of strain hardening and rate sensitivity on the dynamic growth of a void in a plastic material, J. Appl. Mech., 1992, 59(1): 48

[6] D. J. Benson, An analysis of void distribution effects on the dynamic growth and coalescence of voids in ductile metals, J. Mech. Phys. Solids, 1993, 41(8): 1285

[7] X. Y. Wu, K. T. Ramesh, and T. W. Wright, The dynamic growth of a single void in a viscoplastic material under transient hydrostatic loading, J. Mech. Phys. Solids, 2003, 51(1): 1

[8] T. Pardoen, I. Doghri, and F. Delannay, Experimental and numerical comparison of void growth models and void coalescence criteria for the prediction of ductile fracture in copper bars, Acta Mater., 1998, 46(2): 541

[9] T. Pardoen and J. W.  Hutchinson, An extended model for void growth and coalescence, J. Mech. Phys. Solids, 2000, 48(12): 2467

[10] V. C. Orsini and M. A.  Zikry, Void growth and interaction in crystalline materials, Int. J. Plast., 2001, 17(10): 1393

[11] V. Tvergaard and J. W.  Hutchinson, Two mechanisms of ductile fracture: Void by void growth versus multiple void interaction, Int. J. Solids Struct., 2002, 39(13–14): 3581

[12] T. I. Zohdi, M. Kachanov, and I. Sevostianov, On perfectly plastic flow in porous material, Int. J. Plast., 2002, 18(12): 1649

[13] D. R. Curran, L. Seaman, and D. A. Shockey, Dynamic failure of solids, Phys. Rep., 1987, 147(5–6): 253

[14] E. T. Seppala, J. Belak, and R. Rudd, Onset of void coalescence during dynamic fracture of ductile metals, Phys. Rev. Lett., 2004, 93(24): 245503

[15] A. K. Zurek, W. R. Thissell, J. N. Johnson, D. L. Tonks, and R. Hixson, Micromechanics of spall and damage in tantalum, J. Mater. Process. Technol., 1996, 60(1–4): 261

[16] A. K. Zurek, J. D. Embury, A. Kelly, W. R. Thissell, R. L. Gustavsen, J. E. Vorthman, and R. S. Hixson, Microstructure of depleted uranium under uniaxial strain conditions, AIP Conf. Proc., 1998, 429: 423

[17] D. L. Tonks, A. K. Zurek, and W. R. Thissell, Void coalescence model for ductile damage, AIP Conf. Proc., 2002, 620: 611

[18] J. P. Bandstra, D. M. Goto, and D. A. Koss, Ductile failure as a result of a void-sheet instability: experiment and computational modeling, Mater. Sci. Eng. A, 1998, 249(1–2): 46

[19] J. P. Bandstra, and D. A. Koss, Modeling the ductile fracture process of void coalescence by void-sheet formation, Mater. Sci. Eng. A, 2001, 319–321: 490

[20] J. P. Bandstra, D. A. Koss, A. Geltmacher, P. Matic, and R. K. Everett, Modeling void coalescence during ductile fracture of a steel, Mater. Sci. Eng. A, 2004, 366(2): 269

[21] M. F. Horstemeyer, M. M. Matalanis, A. M. Sieber, and M. L. Botos, Micromechanical finite element calculations of temperature and void configuration effects on void growth and coalescence, Int. J. Plasticity, 2000, 16(7): 979

[22] E. T. Sepplälä, J. Belak, and R. E. Rudd, Three-dimensional molecular dynamics simulations of void coalescence during dynamic fracture of ductile metals, Phys. Rev. B, 2005, 71(6): 064112

[23] L. M. Dupuy and R. E.  Rudd, Surface identification, meshing and analysis during large molecular dynamics simulations, Model. Simul. Mater. Sci. Eng., 2006, 14(2): 229

[24] W. Pang, G. Zhang, A. G. Xu, and G. Lu, Size effect in void growth and coalescence of face-centered cubic copper crystals, Chin. J. Comp. Phys., 2011, 28: 540 (in Chinese)

[25] W. Pang, P. Zhang, G. Zhang, A. G. Xu, and X. Zhao, The nucleation and growth of nanovoids under high tensile strain rate, Sci. China – Phys. Mech. Astron., 2012, 42: 464

[26] D. Burgess, D. Sulsky, and J. U. Brackbill, Mass matrix formulation of the FLIP particle-in-cell method, J. Comput. Phys., 1992, 103(1): 1

[27] S. Bardenhagen, J. Brackbill, and D. Sulsky, The material-point method for granular materials, Comput. Methods Appl. Mech. Eng., 2000, 187(3–4): 529

[28] N. P. Daphalapurkar, H. Lu, D. Coker, and R. Komanduri, Simulation of dynamic crack growth using the generalized interpolation material point (GIMP) method, Int. J. Fract., 2007, 143(1): 79

[29] S. Ma, X. Zhang, and X. M. Qiu, Comparison study of MPM and SPH in modeling hypervelocity impact problems, Int. J. Impact Eng., 2009, 36(2): 272

[30] X. F. Pan, A. G. Xu, G. C. Zhang, P. Zhang, J. S. Zhu, S. Ma, and X. Zhang, Three-dimensional multi-mesh material point method for solving collision problems, Commun. Theor. Phys., 2008, 49(5): 1129

[31] X. F. Pan, A. G. Xu, G. C. Zhang, and J. Zhu, Generalized interpolation material point approach to high melting explosive with cavities under shock, J. Phys. D, 2008, 41(1): 015401

[32] A. G. Xu, X. F. Pan, G. C. Zhang, and J. Zhu, Material-point simulation of cavity collapse under shock, J. Phys.: Condens. Matter, 2007, 19(32): 326212

[33] A. G. Xu, G. Zhang, X. F. Pan, and J. Zhu, Simulation Study of Shock Reaction on Porous Material, Commun. Theor. Phys., 2009, 51(4): 691

[34] A. G. Xu, G. Zhang, P. Zhang, X. F. Pan, and J. Zhu, Dynamics and thermodynamics of porous HMX-like material under shock, Commun. Theor. Phys., 2009, 52(5): 901

[35] A. G. Xu, G. C. Zhang, H. Li, Y. Ying, X. Yu, and J. Zhu, Temperature pattern dynamics in shocked porous materials, Sci. China – Phys. Mech. Astron., 2010, 53(8): 1466

[36] A. G. Xu, G. C. Zhang, Y. Ying, P. Zhang, and J. Zhu, Shock wave response of porous materials: from plasticity to elasticity, Phys. Scr., 2010, 81(5): 055805

[37] A. G. Xu, G. C. Zhang, H. Li, Y. Ying, and J. Zhu, Dynamical similarity in shock wave response of porous material: From the view of pressure, Comput. Math. Appl., 2011, 61(12): 3618

[38] F. Auricchio and L. B.  da Veiga, On a new integration scheme for von-Mises plasticity with linear hardening, Int. J. Numer. Meth. Eng., 2003, 56(10): 1375

[39] B. Zhang, et al., Explosion Physics, Beijing: Ordance Industry Press of China, 1997

Rights & permissions

© Higher Education Press and Springer-Verlag Berlin Heidelberg 2013