The development of modern industry has prompted the exploration of material properties under extremely complex service conditions. For materials ranging from traditional ones that have been used for centuries to recently designed ones relying on new manufacturing technologies, studies of performance under extreme service conditions provide an indispensable basis for their applications in a wide variety of fields, such as defense-related and astronautics industries,1–5 the nuclear industry,6–13 and nano-manufacturing.14,15 In this context, the mechanical responses of ductile metals under high-strain-rate loading conditions are of great interest, including their dynamic yielding behavior and dynamic fracture process. The fundamental mechanism governing the dynamic fracture of ductile metals has been widely studied, especially in recent years, since it not only provides a basis for further theoretical studies of the mechanical behavior of these materials, but is also very different from the mechanism of fracture in quasistatic conditions.

Spalling is one of the most important dynamic fracture modes in metals, and observations of spalling using relatively simple experimental methods can reveal fundamental mechanical aspects of the dynamic fracture process.1 The spalling process can be described as follows. A sufficiently energetic strike resulting from plate impact, explosive detonation, or laser shock will first produce shock compression pulses, which are then transformed into rarefaction waves after free-surface reflection. The interference between these rarefaction waves produces a state of high tensile stress in the target, causing planar separation of material. Spalling in brittle metals is characterized by cracks with sharp tips around which little plastic deformation takes place, whereas in ductile metals, it is characterized by voids that tend to be spherical up to a certain size.1 From a microscopic point of view, the spalling in ductile metals results from nucleation, growth, and coalescence of these voids.16 In the past few decades, an abundance of experimental and computational studies have helped to reveal a wealth of information about the nucleation and growth process from microscopic to macroscopic aspects.1–4 On the one hand, these studies have provided increasingly clear physical pictures, on the basis of which we have come to understand the underlying mechanism of the damage initiation process and describe this process in precise mathematical terms. On the other hand, the development of experimental and simulation techniques in this field has led to a need for a powerful and effective theory to describe current observations and to provide a reliable reference for further investigations. For instance, a theoretical model that could accurately predict the critical nucleation stress would be very helpful if the appropriate velocity of the flyer needed to be set in a plane impact experiment to cause only incipient damage to the target. To sum up, the study of theoretical models helps to reveal the mechanism of void nucleation and growth under dynamic loading and makes an important contribution to the design of damage-resistant materials. Therefore, in this review, we summarize classical work and recent research progress on theoretical models related to void nucleation and growth under dynamic loading.

The remainder of the paper is organized as follows. In Secs. II and III, theoretical works on void nucleation and growth, respectively, are introduced, with the former covering models of void nucleation from common crystallographic defects and the latter dealing with growth models based on classical and on newly developed constitutive relations. A summary and outlook are presented in Sec. IV.

Since the stages of void nucleation and growth are usually coupled with each other, the definition of nucleation remains controversial among researchers owing to different perspectives that they adopt. Specific definitions of void nucleation are basically of two types (1) size-based and (2) process-based.

From a size-based perspective, Curran et al.17 defined the nucleation of microscopic fracture as the appearance of microscopic flaws (or voids) whose sizes are comparable to the size of the continuum limit of the material. However, this definition inevitably also includes submicroscopic growth in the vast majority of cases. From a process-based perspective, nucleation is defined as the creation of a flaw (or a void) large enough to grow under an applied stress field.18,19

For nucleation occurring at the microscale, the above two definitions are often equivalent to each other. For nucleation at the nanoscale, however, a size-based definition may lead to an inaccurate division of the stages involved. Although they recognized the theoretical attraction of the process-based definition, Curran et al.17 held the view that the difficulty of applying continuum mechanics to describe events at a scale much smaller than the continuum limit hinders its adoption. This point of view was reasonable, considering the experimental techniques available at that time, but subsequent technical developments have made it somewhat outdated. Recently, extensive experiments and simulations have focused on nanoscale nucleation, and the abundant information that these have provided has now revealed the physical mechanisms of nanoscale nucleation, making it feasible to model some of these mechanisms in the framework of continuum mechanics.20,21

In this section, we overview the microscopic mechanisms of void nucleation under dynamic loading, which can be classified into two categories: heterogeneous and homogeneous nucleation.1,22 Heterogeneous nucleation usually occurs at micrometer-sized defect sites such as grain boundaries, triple junctions, and second-phase particles. Homogeneous nucleation, as suggested by previous studies,1,19,22 refers to the nucleation process originating at perfect lattice structures or submicroscopic heterogeneities within grains. Besides the nucleation mechanism, the size distribution and statistical evolution of the nucleated voids are also important issues in this field. The associated models are covered in previous comprehensive reviews3,23 and are not discussed in this paper.

Deformation around inclusions and second-phase particles is an important mechanism for void nucleation under dynamic loading. As experimental evidence for this, Pedrazas et al.25 presented the strong correlation between second-phase particles and the dimples produced in spallation at very high strain rates [∼(2–4) × 10⁶ s⁻¹], even for commercially pure aluminum (99%). The nucleation mechanism of inclusions and second-phase particles includes cracking of inclusions, debonding at interfaces, and matrix fracture near inclusions.1,17 Models of void nucleation in material containing a dispersion of inclusions or second-phase particles have been proposed in the past few decades to predict the onset of damage initiation. For void nucleation at the matrix–particle interface, there are two necessary conditions for the void nucleation process, namely, the mechanical debonding condition and the energetically favorable condition.17

The mechanical debonding condition of the nucleation process was analyzed by Ashby24 based on a stress calculation, and it was found that nucleation occurs when the tensile stress along the tensile axis (which may result from the applied tensile stress along x′1 in Fig. 1 or the applied compression stress along x′2) exceeds the assumed fracture stress, which can be the normal fracture stress of the particle–matrix interface. However, this condition itself is insufficient, since the energetically favorable condition is also essential for nucleation. Although the energy expressions may be different from model to model, the core of this condition can be expressed as saying that the work done should be no less than the energy required to form the new surfaces of the void. In this review, these two conditions are referred to more concisely as the “stress condition” and the “energy condition,” as in previous studies.26,27 Two typical models are introduced subsequently, which take both of these conditions into account.

Tanaka et al.26 analyzed the critical strains that satisfy the stress condition and the energy condition, respectively. The critical strain that corresponds to the stress condition26P33I≥σs was obtained asεc≥δ,α<1,δ1α,α≥1,(1)whereδ=(7−5ν)(1+ν∗)+(1+ν)(8−10ν)α10(7−5ν).(2)Here, α = E∗/E, where E∗ and E are the Young’s moduli of the particle and the matrix, respectively. ν∗ and ν are the Poisson’s ratios of the particle and the matrix, respectively. σ_s is the interface strength and takes a value of σ_s = max{α, 1} · E/10 in their work. From a energetic perspective, Tanaka et al.26 calculated the total energy of the specimen based on calculation of the strain energy in and around a spherical particle. The energy condition is therefore expressed as G₁ ≥ G₂, where G denotes the total energy of the specimen, with the subscripts 1 and 2 referring to the cases before and after nucleation, respectively. On this basis, the critical strain is derived asεc≥β12r,α<1,β12rα,α≥1,(3)whereβ=l48[(7−5ν)(1+ν*)+(1+ν)(8−10ν)α][(7−5ν)(1−ν*)+5(1−ν2)α](7−5ν)2[2(1−2ν*)+(1+ν)α],(4)with l = 1 × 10⁻⁹ cm, so the unit of β is cm. r is the radius of the spherical particle. Based on Eq. (3), the critical strain is inversely dependent on the square root of the particle size. Combining the above two conditions, Tanaka et al.26 concluded that the satisfaction of the energetically favorable one, Eq. (3), is a sufficient and necessary condition for void nucleation at particles of radius smaller than 10–15 nm. For particles of radius larger than 10–15 nm, the void will gradually grow, with further plastic deformation and under control by the stress condition, Eq. (1). However, Eq. (1) takes no account of the influence of particle size and underestimates the critical strain required for large inclusions compared with experimental results (Fig. 2).

Based on the deformation incompatibility between the hard particles and the matrix, a modified model was proposed from a different perspective.27,32 The stress condition of the model in Goods and Brown27 accounted for the influence of particle size through the work-hardening effect of the inclusion. The local flow stress was derived asσl=κμbρl,(5)where ρ_l is the local dislocation density, μ is the shear modulus of the matrix, b is the magnitude of the Burgers vector, and κ is a material constant. The local dislocation density is related to the particle radius and the plastic strain asρl=5εp3rb.(6)Considering the plastic constraint resulting from the particle, the relation between the local stress at the interface σ_c and the local flow stress σ_l was determined by Brown and Stobbs32 as σ_c ≃ 4.2σ_l. Therefore, the critical strain can be derived asεc≥130σcκμ2rb.(7)According to Eq. (7), the critical strain increases with the particle size, which is in keeping with the experimental results (Fig. 2). In the energy condition of their model, the particle is assumed to be purely elastic, while the matrix material undergoes both elastic and plastic deformation. The deformation incompatibility between the hard particles and the matrix is defined as εp∗. For situations where stress relaxation does not occur, εp∗ is equal to the plastic strain ε_p. For situations where stress relaxation does occur, Brown and Stobbs32 adopted the approximate formεp∗=bεpr.(8)The energy condition states that void nucleation occurs only when the energy released by the separation between the particle and the matrix is no less than the energy needed to form the new surface. The energy condition in Brown and Stobbs32 is expressed asΔEelas+ΔEs<0,(9)where ΔEelas=−43πμ∗r3εp∗∗2 and ΔE_s = 4πr²γ, with μ∗ being the shear modulus of the particle and γ the surface energy. Based on Eqs. (7) and (9), Brown and Stobbs32 concluded that the energy condition is always satisfied before the stress condition, except for very small particles. This conclusion is in keeping with that obtained by Tanaka et al.26 However, Eq. (9) means that nucleation needs separation of the whole surface of the particle. This kind of setting results in an overestimate of the required energy, and consequently an overestimate of the critical strain compared with the experimental results. An approximate correction was proposed by reconsidering the proportion of the particle–matrix interface that separated during the nucleation process.27 In the improved version of the model, the separation of the particle–matrix interface was assumed to occur only at the particle poles in the direction of the maximum tensile stress, as shown in Fig. 3. This assumption had already been adopted in previous theoretical work24,26 and verified by experiments.26

In summary, there are two necessary conditions for void nucleation at inclusions, namely, the mechanical debonding condition and the energetically favorable condition. The dependence on particle size varies from r^−1/2 through no dependence to r owing to the different assumptions adopted in the models. It is worth mentioning that although the above models are based on quasistatic analyses, the feasibility of applying them to the dynamic loading condition has been proved.17

Besides inclusions and second-phase particles, grain boundaries are also important heterogeneous nucleation sites. For void nucleation on the grain boundary, only the energetically favorable condition needs to be satisfied, i.e., the external work done on an element containing a virtually growing void must exceed the energy required to create the surface of the new void.17,18 The critical condition is expressed asσmdv>γda,(10)where σ_m is the mean tensile stress, γ is the surface energy, and dv and da are the increases in void volume and void area, respectively. For a spherical void, the critical condition can be derived according to Eq. (10) asrc=2γσm.(11)Equation (11) reveals that the size of the nucleated void and the external loading are closely related to each other. The minimum size of a void that can spontaneously grow during the current loading can be obtained according to Eq. (11), and vice versa. Equation (10) provides a simple but basic way to represent the energy relation for void nucleation on a grain boundary. The influence of other factors such as void shape and dislocation pile-up has also been considered on the basis of the original energy relation, Eq. (10).18,33,34 However, some important properties of the grain boundary void nucleation revealed by experiments and simulations under dynamic loading conditions cannot be explained by Eq. (10), such as the influence of the misorientation angle35–40 and the loading direction.41 Theoretical models that can precisely describe the dynamic response of the grain boundary are urgently required.

In contrast to alloys and high-purity polycrystalline metals, few pre-existing nucleation sites are available for high-purity monocrystalline metals under high-strain-rate loading. Dislocation networks and vacancy clusters generated during the deformation process serve as potential nucleation sites for the homogeneous nucleation process in this kind of material. In this section, models related to the vacancy-mediated nanovoid nucleation mechanism are introduced in detail.

The void nucleation stage begins after the shock wave has propagated in the interior of the sample. Vacancy clusters generated by diffusion-mediated aggregation are potential nucleation sites that may finally evolve into macroscopic voids and lead to fracture. As suggested by Reina et al.,19 the “already nucleated” state of a potential nucleation site corresponds to the onset of plastic cavitation, i.e., the formation of vacancy clusters of a size such that subsequent growth can happen by plasticity. Therefore, the triggering of plastic cavitation can be regarded as the criterion for void nucleation.19 It should be emphasized that the emission of dislocation loops is the mechanism of plastic cavitation under high-strain-rate loading, which is different from the quasistatic case and has been confirmed by experiments and simulations.42–49 In the case of quasistatic loading, plastic cavitation can occur by a vacancy-diffusion-based mechanism, i.e., vacancy condensation under an external stress field or elevated temperature.23,50 However, the time available for vacancy diffusion under shock loading is very limited and is far less than the time required for a diffusion-based cavitation mechanism.20,51 Previous studies have focused on the external stress required to trigger dislocation emission on void surfaces in the interior of ductile metals. The so-called critical emission stress has been simulated in different metals and under different loading conditions.43,46–49 Theoretical models have also been proposed to analyze the emission process and derive a stress-based emission criterion.

Two-dimensional (2D) elasticity theory is adopted to develop the emission criteria from scratch. An explicit expression for the critical emission stress was first proposed by Lubarda et al.20 The first step in establishing the emission criterion is to derive the total glide force actingon the dislocation. This is obtained asF(ξ)=2σbξξ2+122−μb2πa(1−ν)ξξ4+14ξ2+122ξ4−14,(12)where μ is the shear modulus, ν is Poisson’s ratio, a is the void radius, b is the magnitude of the Burgers vector, and ξ = x/a is the position of the dislocation. Considering the existence of an unstable equilibrium position of the dislocation when σ takes the values achieved in shock loading experiments (∼GPa), it is assumed that dislocation emission from the void surface will be triggered if the equilibrium distance of the dislocation from the void surface is equal to the dislocation core cutoff radius. Thus, the criterion is derived asσcr=μb/a2π(1−ν)1+2w/a4+11+2w/a4−1,(13)where w is the dislocation core cutoff radius, which is a measure of the width of the dislocation. This criterion is based on the 2D configuration where an edge dislocation is emitted from an infinitely extended cylindrical void in an infinite isotropic medium (Fig. 4). The influence of the void radius and the dislocation width, as well as the dislocation interaction, was also investigated in the work of Lubarda et al.20 The larger the void radius and the wider the dislocation core, the smaller is the critical nucleation stress, as shown in Fig. 5. Besides, the interaction between dislocations impedes the emission of the dislocation and result in higher value of the critical emission stress compared with the situation where the dislocation interaction is neglected.

However, there are several drawbacks in their model. First, it assumes that the dislocation is emitted from the void surface along the direction of the maximum shear stress under biaxial equal tension.20 The emission direction (i.e., the rotation angle of the slip plane relative to the radial direction) is not taken as variable. Therefore, their model cannot account for the finite number of slip systems available in real materials, which limits its range of application. Second, effects resulting from the existence of other voids cannot be investigated. The void in their model is assumed to be in an ideal elastic medium, whereas actual void nucleation and growth process take place in solids with dispersed voids.

A revised dislocation emission model taking account of the emission direction was proposed by Lubarda.52 Two related sets of polar coordinates (ρ, θ) and (r, φ) were adopted to locate the position of the dislocation relative to a point on the void surface and the center of the void, respectively. The total glide force under biaxial equal tension was derived asF=σba2r2sin⁡2(θ−φ)+μb2⁡cos(θ−φ)2π(1−ν)a2r32r2−a2r2sin2⁡θ−r2a2r2r2−a2.(14)Using the condition for emission proposed by Lubarda et al.,20 the critical emission stress is derived asσcr=μb/a4π(1−ν)r02a21f0−2f0,(15)wherer02=a2+w2+2aw⁡cos⁡θ,f0=r02−a2r02sin⁡θ.(16)The angle θ_cr specifies the emission direction in Lubarda’s revised model, which corresponds to the minimum of the far-field applied stress,52 i.e.,∂σcr∂θθ=θcr=0.(17)Besides taking account of the emission direction, Lubarda52 also extended the model to incorporate a general biaxial loading condition and analyzed dislocation emission under several typical biaxial loading conditions. It was found that remote pure shear loading requires a lower value of the critical emission stress compared with the equal biaxial loading condition and the uniaxial loading condition. The stress concentration at the void surface is considered to be the factor underlying this phenomenon.

To consider the influence of other voids, a criterion for dislocation emission in porous solids was proposed by Wilkerson and Ramesh.21 In their work, an auxiliary problem that is simple enough for closed-form solutions to be derived was constructed from the original one. The dislocation emission from a spherical void in a porous material under a general 3D macroscopic stress state was simplified to a 2D problem in which dislocation emission from an infinitely extended cylindrical void occurs in an infinite isotropic linear-elastic medium. Dislocation emission was proved to occur more easily with increasing porosity in their model.

The 2D criteria lay the foundation for a theoretical description of dislocation emission. However, these models based on a 2D configuration introduce an excess of assumptions. These assumptions help in the derivation of closed-form solutions that are as simple as possible, but they also leave out some important information. First, the real physical process of dislocation emission, namely, the emission of 3D dislocation loops,43,54–59 is overly simplified in 2D configurations. Second, the influence of some important factors in the actual 3D process, such as the actual porosity or the actual stress triaxiality, cannot be investigated directly in 2D configurations. Owing to the dimensionality reduction, further hypotheses must be introduced in 2D models. For example, the original 3D macroscopic stress is approximated as an effective 2D stress state,21 while dispersed spherical voids are approximated as arrays of cylindrical voids to derive an effective porosity.21

Pioneering work on the mechanical description of 3D dislocation loops was done by Willis and Bullough,60 who analyzed the interaction between a spherical void and a circular interstitial dislocation loop using spherical harmonics. On this basis, Ahn et al.61,62 proposed the concept of a threshold applied stress for the emission of prismatic dislocation loop (PDL) in an ideal elastic material under hydrostatic loading. However, the accuracy of their model is limited by the hypotheses required and the use of artificially introduced parameters. Besides, some important factors were not been included in their model, such as emission position, porous softening, and stress triaxiality.

Recently, a 3D criterion in porous media has been proposed to determine the critical emission stress of a PDL, σ_cr, and the corresponding emission angle θ_cr specifying the emission position on the spherical void surface,53 as illustrated in Fig. 6. Based on this model, Sui et al.53 systematically studied the factors influencing PDL emission in a 3D configuration, such as the geometric parameters, the stress triaxiality, and the porosity. The criterion for dislocation was derived as−1+2η2cos⁡2θcr+(1−η)37−5ν(sin4θcr−3⁡sin2θcr⁡cos2θcr)−3⁡cos⁡2θcr+3γaσcrcos⁡2θcr=0,(18a)τgimage(r0,θ0)+τgexternal(r0,θ0,σcr,f)+τgsurface(r0,θ0)+τgfriction=0,(18b)where a is the void radius, γ is the surface energy, η is the loading parameter, f is the porosity, τ_g is the glide shear stress on the glide cylinder, and the superscripts “image,” “external,” “surface,” and “friction” refer to the image stress field, the remote stress field, the surface energy, and the lattice fraction, respectively. Other variables in Eq. (18) are introduced in Fig. 6. It is revealed that the critical emission stress for a PDL decreases with the dislocation width and the void size. Emission of a PDL emission from a void in a porous medium will be facilitated by the nucleation and growth of other PDLs, since they result in an increase in porosity (Fig. 7). The above conclusions are in keeping with those obtained from 2D models. In addition, the critical emission stress is notably affected by the stress triaxiality (Fig. 8). When the loading state varies from uniaxial tension (i.e., η = 0 in Fig. 8) to hydrostatic tension (i.e., η = 1) with increasing stress triaxiality, the critical emission stress decreases significantly.

Void growth is vital for the evolution of ductile damage and dramatically influences the service life and performance of materials. In this section, phenomenological and physics-based growth models are introduced. Phenomenological models adopt the traditional constitutive relationships to describe the dynamic growth process of a void, while physics-based models describe the same process by focusing on the laws governing the evolution of microstructures under high-strain-rate loading conditions.

It is worth mentioning that quasistatic void growth models provide the basis for the dynamic ones. The earliest models describing the growth of isolated voids in ductile metals under quasistatic loading conditions were first developed in the last century.63–68 Since these models cannot capture the important effect of an evolving void volume fraction in actual situations, the micromechanics-based Gurson model69 and its modifications were subsequently proposed to consider the influence of void interactions. The development of dynamic growth models is closely related to that of quasi-static models. In this section, phenomenological dynamic growth models are systematically introduced. We start from the dynamic growth of an isolated void in an infinite medium as the basis for the following models. Subsequently, statistically based models and micromechanics-based models that describe the macroscopic mechanical response of the damaged material are introduced.

The dynamic growth of a single void in a power-law-hardening unbounded solid was studied by Ortiz and Molinari,70 who systematically investigated the effects of inertia, strain hardening, and rate dependence under conditions of high-strain-rate loading from an energetic viewpoint. They assumed the matrix to be incompressible and the void to remain spherical during the process of growth. Based on the principle of conservation of energy, the evolution laws of the void radius for both the rate-dependent [Eq. (19a)] and rate-independent [Eq. (19b)] cases were derived as follows:ddt32ρa˙24πa33+σyε˙023εy1/n×2a˙ε˙0a(m+1)/m⁡faa0,m,n4πa33=3pa˙a4πa33,(19a)32ρa˙24πa33+nσyεyn+123εyn+1/ngaa0,n4πa33=pa3−a03a34πa33,(19b)wherefaa0,m,n=∫1∞logxx−1+a03/a31/nx−(m+1)/mdx,(20)gaa0,n=∫1∞logxx−1+a03/a3(n+1)/ndx,(21)ρ is the mass density, a is the void radius, a₀ is the initial void radius, and σ_y, ε_y, ε˙0, m, and n are material constants (for details, see Ref. 70). By considering the dimensionless forms of the evolution equations (19), short- and long-term approximations of the original form are derived by leaving out the insignificant time-related term to describe the early and late stages of void growth. It is revealed that inertia, strain hardening, and rate sensitivity all have a remarkable influence on void evolution. To measure the effective inertia of the system, the following parameter is defined:D=ρa02kε˙ref2−1/m,(22)where k=σyεy−1/nε˙0−1/m, and ε˙ref is a reference strain rate that is chosen according to the problem under consideration.70 The early stage of void growth is dominated by viscous effects, whereas the late stage is dominated by inertia. The role of inertia becomes significant only when the void reaches a critical size that depends on both the mechanical properties and the rate of expansion. The rate sensitivity of the material retards the expansion of the void for the short-term solution in the work of Ortiz and Molinari,70 which is as expected. For the long-term solution, however, the rate-dependent description of the material results in faster void growth than in the rate-independent situation, which is counterintuitive, since the one would expect the rate sensitivity to influence the void growth in a consistent way.

Taking the influence of heat conduction into consideration, Wu et al.71–73 proposed a thermal–mechanical coupling model to study the dynamic growth of a single spherical void, with special attention being paid to the comprehensive effect of inertia, thermal softening, rate dependence, and plastic strain gradient. The evolution law of the void in their work was derived asp(t)−∫02⁡ln(a/a0)σeexp(3ε/2)−1dε=ρaa⃜+3a˙22,(23)where p(t) is the applied far-field stress, ρ is the material density, a₀ is the initial size of the void, and σ_e is the von Mises equivalent stress. Similar to Eq. (22), an inertial scaling factor was defined by Wu et al.71 according to the dimensionless form of Eq. (23) to investigate the influence of inertia, i.e.,Iinertia=ρa02p˙2σyps2,(24)where p˙ is the loading rate, σ_y is the yield strength, and p_s is the steady value of the applied far-field stress. Compared with Eq. (22), Eq. (24) decouples the influence of rate sensitivity and inertia and quantifies the effect of inertia in a more direct way. It is found that inertial effects are dependent on the loading rate and the initial size of the void. The influence of inertia on void growth becomes more obvious for situations with a larger initial void and a higher loading rate. In addition, Wu et al.71 also analyzed the influence of thermal conduction. The constitutive relationship of the matrix, includes thermal softening and strain hardening, was assumed to take the formσeσy=f(ε)h(T∗),(25)where f(ε) is the power-law strain hardening function. h(T∗) is the thermal softening function, which was derived according to conservation of energy ash(T∗)=exp−cωσy(n+1)ρcp[εf(ε)−εy],(26)where c is the thermal coefficient, c_p is the specific heat at constant pressure, and ω is the fraction of plastic work transformed into heat. It is found that the critical stress for unstable growth of the void is reduced by thermal softening. By comparison, heat conduction makes a positive contribution to the stable growth of the void, which is evidently affected by the void size under dynamic loading conditions. The stabilizing influence of thermal conduction becomes remarkable for small voids because of the high heat transfer efficiency around them. The effect of rate hardening, which appears to contribute in the opposite way in Ortiz and Molinari,70 is also discussed by Wu et al.71,72 It is revealed that the rate sensitivity reduces the growth rate of the void in the early stage of dynamic growth, while it has little effect in the late stage, since the rate-hardening effect decreases with increasing void size.

The evolution law for an isolated void in an infinite plastic matrix71 was adopted by Molinari and Wright74 to propose a statistical description of damage evolution. Equation (23) is rewritten asp(t)−pc=ρaa⃜+3a˙22,(27)for the case of p(t) ≥ p_c. The local critical nucleation stress for a potential nucleation site is defined aspc=lima/a0→∞∫02⁡ln(a/a0)σeexp(3ε/2)−1dε=23σy+∫εy∞σeexp(3ε/2)−1dε.(28)The value of p_c is well determined by the tensile yield stress σ_y, the tensile yield strain ε_y, and the von Mises equivalent stress σ_e, since they contain all the constitutive information of the local material around the potential nucleation site. The fact that p_c has the same value when either a → ∞ or a₀ → 0 indicates that void nucleation (i.e., the initial size of the void tending to zero) can be regarded as a bifurcation (i.e., the void radius growing to infinity) in the homogeneous solution.74 Besides, Molinari and Wright74 proved that voids of all different sizes grow at the same equilibrium rate for a given applied stress above the critical stress. Therefore, they ignored the initial sizes of the voids (i.e., the initial porosity) in their work and assumed that all voids in a local area nucleate at the local critical stress. Considering the material heterogeneity due to microstructures, the local critical stress varies form site to site. Molinari and Wright74 introduced a statistical theory to describe these variations. The probability density function of p_c, denoted by g(p_c), is given according to the Weibull law or the Gaussian law in their work. For instance, the probability density function of p_c based on the Weibull law can be expressed asg(pc)=β2β1〈pc−p0c〉β1β2−1⁡exp−〈pc−p0c〉β1β2,(29)with the conventionx=12x+x.β₁ and β₂ are statistical parameters, and p_0c is the lower limit for having a positive probability of void nucleation. When the loading rate is a positive constant, i.e., p(t)=p˙t, a simple similarity solution of Eq. (27) can be derived asa=833p˙ρt−pcp˙3/2.(30)It should be emphasized that potential nucleation sites are equivalent if they have the same value of the critical nucleation stress. As a result, the relationship between the void volume generated per unit initial volume V_void and the critical nucleation stress p_c can be derived asdVvoiddpc=43πa3Ng(pc),(31)where N is the number of potential nucleation sites per unit volume. Therefore, the porosity can be derived asf=Vvoid1+Vvoid,(32)withVvoid=∫pc43πa3Ng(pc)dpc.(33)On combining Eqs. (29) and (30) with Eqs. (32) and (33), the time evolution of the porosity can be obtained. It is worth noting that the probability density function g(p_c) is the bridge between the microscopic information and the macroscopic damage description. A physical-based modeling of g(p_c) will help to describe the damage process as precisely as possible. However, the governing equation (27) is for a configuration of an infinite body containing a single spherical void and subjected to external hydrostatic tensile stress. Since the influence of void interaction is not taken into consideration, the porosity evolution function proposed by Molinari and Wright74 is only applicable to the early stage of void growth. Even so, the pioneering idea of using a probability density function to represent the variation of critical nucleation stress from site to site has been widely adopted in subsequent work.75–81 In Sec. III A 3, we will see that the combination of this method and the dynamic Gurson-type model is able to describe the damage evolution process in a rather precise way.

The studies described in Sec. III A 1 focused on the dynamic response of a single void in an elastoplastic matrix. These are highly conducive to investigations of the influence of different factors during the process of void growth, but they cannot be used as a basis for studying the macroscopic mechanical response of a material under dynamic loading, which requires a different perspective.

Gurson69 was the first to propose approximate yield criteria and flow rules for ductile porous materials in a micromechanics-based framework, with the matrix material being idealized as rigid–perfectly plastic and obeying the von Mises yield criterion. The yield function was derived asΦ=Σeσ02+2f⁡cosh3Σm2σ0−(1+f2)=0,(34)where σ₀ is the microscopic equivalent tensile yield stress, and Σ_e and Σ_m are the macroscopic equivalent stress and the macroscopic mean stress, respectively. Besides its ability to describe the progressive failure of porous materials, the original form of the Gurson model possesses nice properties for a variety of special conditions, such as the zero-porosity condition, the purely deviatoric loading condition, and the purely hydrostatic loading condition (see also Ref. 3).

Because the rigid–perfectly plastic assumption69 is too strict to be satisfied by the majority of materials, extensions have been proposed to study more complex constitutive behaviors and geometric configurations, such as viscoplastic effects,82,83 plastic anisotropy,84–86 crystal plasticity,87,88 and void shape effects.89 The most important modification of the Gurson model is the Gurson–Tvergaard–Needleman (GTN) model,90–92 which provides a mixed phenomenological and micromechanics-based framework to precisely describe the whole fracture process of a porous medium. The yield condition is expressed asΦ=Σeσ¯e2+2q1f∗⁡cosh3q2Σm2σ¯e−(1+q3f∗2)=0,(35)where σ¯e is the microscopic equivalent tensile flow stress, which accounts for the strain hardening effect, and f^∗ is the equivalent porosity, which accounts for the interaction between the microscopic shear bands and the voids.

The Gurson model has also been used to deal with the dynamic failure of ductile materials. A Gurson-type model has been directly applied to dynamic ductile fracture93–95 and dynamic crack growth.96–98 Since the inertial effects are not considered in the original form of the Gurson model, further modification is necessary for the situation of dynamic loading.99 Similar to the original form in Eq. (34), a dynamic correction of the Gurson model was derived by Wang and Jiang99 based on the principle of virtual work:Φ=Σe−Σedσ02+2f⁡cosh32Σm−Σmdσ0−(1+f2)=0,(36)where Σ_e and Σ_m are the macroscopic effective and mean stresses. Both of these can be separated into two parts, i.e., Σe=Σes+Σed and Σm=Σms+Σed, where the superscripts s and d indicate the quasistatic and dynamic parts, respectively. Detailed expressions for these variables can be found in Ref. 99. An extension of this work was proposed by Wang100 to consider rate sensitivity and thermal effects in addition to the influence of inertia.

Although the expression for the macroscopic stress includes both static and dynamic contributions, Molinari and Mercier101 held the view that the definition of the macroscopic stress in previous works99,100 was still a “static” one, i.e., the average value of the microscopic stress was given byΣ=1Ω∫ΩσdV.(37)Molinari and Mercier101 redefined the same variable asΣ=1Ω∫ΩσdV+1Ω∫Ωρx⃜⊗xdV,(38)where x⃜ is the acceleration of a particle, and ⊗ denotes the tensor product. Based on the newly defined macroscopic stress, an explicit macroscopic stress–strain-rate relationship was obtained to describe the dynamic behavior of the porous material with a rigid viscoplastic matrix:Σ−Σstatic=ρa215(f−2/3−f)D′˙+D′⋅D′−13tr(D′⋅D′)I+(f−2/3−1)DmD′+16(D′:D′)I−(f−2/3−f−1)D˙mI+3f−1−52f−2/3−12f−2Dm2I,(39)where D is the macroscopic plastic strain rate tensor, D′ is the deviatoric part of D, D_m is the trace of D, I is the second-order identity tensor, and a dot above a symbol indicates is the material time derivative. The macroscopic stress in Eq. (39) is also built from two parts: Σ^static is the viscoplastic stress tensor calculated from the quasistatic potential (for details, see Ref. 101), and the term scaled by ρa² represents the dynamic contribution to the macrostress. It is worth noting that Eq. (39) and its inverse form can be used to describe the dynamic behavior of porous material under both kinematic and stress control. However, there are still some limitations to this model. First, it is restricted to cases with low porosity, since the interactions between voids are not taken into account. Second, the spherical shape of the void in the representative volume element (RVE) requires the model to be used only in situations with large stress triaxiality.

Under the assumption that the deviatoric stress components are negligible when a large hydrostatic pressure is generated in the metallic plates during plane impact, Eq. (39) was rewritten by Czarnota et al.75 asp(t)−pstatic=ρϕ1aa⃜+ϕ23a˙22,(40)where ϕ₁ = 1 − f^1/3 and ϕ2=1−43f1/3+13f4/3. In the case of an unbounded matrix (i.e., f → 0), ϕ₁ = ϕ₂ = 1, and Eq. (40) reduces to Eq. (27) on replacing p^static with p_c, which means that the dynamic contribution of the governing equation is also suitable for the early stage of void growth, where the porosity is negligible. On this basis, Czarnota et al.75 proposed a model that can be used to describe the whole evolution process from nucleation to large size. The material resistance p^static in their model is given bypstatic=inf(pc,pvisco),(41)where p_c is the cavitation pressure above which void nucleation occurs, p^visco is the viscoplastic stress calculated from the plastic potential, Eq. (35), which accounts for the loss of stress-carrying capacity. They also derived an approximate solution of Eq. (40) for the case of linearly increasing loading p(t)=p˙t.

It should be emphasized that Eq. (40) is the governing equation of the hollow sphere configuration. Taking the hollow sphere model as the microscopic unit cell, different homogenization methods are used to link the macroscopic quantities defined at the remote boundary of the material domain to the microscopic ones defined at the unit cell level.76,78 Czarnota et al.76 extended Eq. (40) to account for an elastic–viscoplastic material response and implemented this extension in ABAQUS/Explicit software to simulate plate impact experiments, and the results reproduced experimentally measured free-surface velocity profiles with excellent accuracy (Fig. 9). Jacques et al.78 used a different homogenization scheme to simulate crack growth in a notched bar and in an edge-cracked specimen. It was revealed that the effects of micro-inertia lead to lower crack speed and higher dynamic fracture toughness (Fig. 10).

The growth models described in Sec. III A are based on traditional plastic constitutive relationships and vary from simple models (e.g., the rigid–perfectly plastic model) to complex ones (e.g., the thermal–mechanical coupled power-law-type rate-sensitivity model). However, these models may become invalid at extreme strain rates, i.e., ≥10⁸/s, since there is a dramatic transition in the strain rate dependence of the spall according to the summary presented by Reina et al.19 (Fig. 11). Therefore, phenomenological growth models are perhaps only applicable in the thermally activated glide regime, and a dislocation-based viscoplasticity constitutive relation that accounts for drag and relativistic effects is needed in the case of extreme loading conditions.83

A dislocation-based viscoplasticity model applicable at very high strain rates was first proposed by Austin and McDowell.82 Based on J₂ flow theory with isotropic hardening, the plastic deformation rate is written asDp=Nmbv¯σ¯s,(42)where N_m is the mobile dislocation density, b is the magnitude of the Burgers vector, v¯ is the mean dislocation velocity, s is the deviatoric part of the Cauchy stress tensor, and σ¯ is the von Mises equivalent stress. The mean dislocation velocity v¯ is a function of the shear stress, the glide resistance of the material, and the temperature. The initial form of v¯ is derived asv¯=L¯tw+tr,(43)where t_w is the time a dislocation spends waiting for thermal assistance to overcome an obstacle, and tr=L¯/v¯r is the time for a dislocation to move between obstacles, with v¯r being the mean velocity of dislocation gliding. Taking the influence of drag and relativistic effects into consideration, v¯r can be obtained asBv¯r=τeffb,(44)where the dislocation drag coefficient B accounts for relativistic effects and is given byB=B01−(v¯r/cs)2,(45)with c_s being the shear wave speed. τ_eff = τ_eff(τ, τ_μ) is the effective stress, which is derived from a certain form as a function of τ and τ_μ, such as τeff=τ2−τμ2 or τ_eff = τ − τ_μ, where τ is the applied shear stress, τ_μ is the athermal threshold shear stress required to overcome long-range resistance, which is dependent on the microstructure of the material. For instance, for a pure fcc metal, τ_μ can be obtained according to the Taylor hardening relation, i.e.,τμ=κμbNim+Nm,(46)where κ is the dislocation-interaction parameter, μ is the temperature-dependent shear modulus, and N_im is the immobile dislocation density. The evolution laws of the mobile and immobile dislocation densities are derived asN˙im=N˙trap−N˙rec,(47a)N˙m=N˙nuc+N˙mult−N˙ann−N˙trap,(47b)where the subscripts “trap,” “rec,” “nuc,” “mult,” and “ann” refer to dislocation trapping, recovery, nucleation, multiplication, and annihilation, respectively. The physics-based evolution equations of the terms in Eq. (47) are discussed in detail by Austin and McDowell.82 Thus, a complete constitutive relation is developed to describe the dislocation-based viscoplasticity deformation at extreme strain rates. Austin and McDowell’s model has been adopted to calculate the material velocity profiles for cases when the shock amplitude is less than 10 GPa, and the results are in good agreement with experimental measurements (Fig. 12), which indicates the validity of their model in the high-strain-rate regime.

Starting from the dislocation-based J₂ theory of Austin and McDowell,82 Wilkerson and Ramesh111 developed a micromechanics-based model to describe the void growth rate, taking account of the effects of micro-inertia, dislocation kinetics, and substructure evolution. The RVE in their work was still taken as the classical spherical shell model, in which a spherical void of radius a is embedded in a sphere of radius r₀ (with the inner and external radii of the spherical shell in the reference configuration being denoted by A and R, respectively). Based on the assumptions that the plastic deformation is incompressible and that the elastic deformation rate is negligible, the strain rate of a material element located at x = re_r in the current configuration of the RVE was derived asε˙=a3r33a˙a.(48)Therefore, v¯r is related to the evolution of dislocation structures by Orowan’s relation ε˙=bNmv¯r, i.e.,v¯r=a3r33a˙abNm.(49)The governing equation of void growth was derived according to the balance of radial momentum in the current configuration asp(t)−∫ar04rτdr=ρϕ1aa⃜+ϕ23a˙22,(50)where ϕ₁ = 1 − f^1/3 and ϕ2=1−43f1/3+13f4/3, as defined in Sec. III A 3, and f=(a/r0)3 is the porosity. The form of Eq. (50) is consistent with that of Eq. (40). The spatial integral on the left-hand side of Eq. (50) is decomposed into two parts, i.e.,∫ar04rτdr=Rcr+Rdd,(51)where R_cr corresponds to the quasistatic flow strength, and R_dd is associated with the effective stress τ_eff = τ − τ_μ, the rate dependence of which is obtained according to Eq. (44) asRdd=∫ar04r(τ−τμ)dr=∫ar04rB0v¯rb11−(v¯r/cs)2dr.(52)Thus, the void evolution is related to the dislocation dynamics and the structure evolution based on Eqs. (49)–(52).

It should emphasized that physics-based growth models have undergone rapid development in the past decade.83,112–115 The use of dislocation-based J₂ theory to describe the dynamic response of materials under extreme strain rates has been further extended to dislocation-based crystal plasticity by Lloyd et al.112 and Luscher et al.113 On this basis, the theory of dynamic void growth in single crystals under extreme loading has been developed by Nguyen et al.114,115 to describe the growth process under general loading states, since their previous physics-based model is only suitable for a pure hydrostatic loading state.111 Much effort has been also been spent on modifying the description of the evolution of dislocation substructures to obtain a closed-form approximation of the governing differential equations, and on implementing the constitutive relation in finite element methods to simulate experimental results.83,113–115

This review has provided a summary of theoretical research on void nucleation and growth in ductile metals under dynamic loading. Both the macroscopic mechanical response and the microscopic physical mechanism, as the two most important aspects, have been investigated over the past few decades. Much theoretical works has been done to clarify the underlying mechanisms of damage initiation and evolution and to provide increasingly precise descriptions of these mechanisms. However, there is still a long way to go to before a complete framework describing ductile damage initiation under extreme dynamic loading is established.

First, with regard to dynamic void nucleation at complex microstructures, such as grain boundaries, grain boundary triple junctions, and the dislocation cells, the currently available physics-based theoretical models suffer from deficiencies arising from the combined effects of the complex nature of the stress state and the irregular geometric configuration. With future experiments likely to provide more detailed information, related theoretical models are expected to be developed and to help in investigating the laws of evolution and the key factors involved in these process.

Second, it is important but challenging to extend the study of traditional ductile metals to newly developed metallic materials that have good performance and wide prospects for application, such as nanocrystalline materials and high-entropy alloys. The processes involved in dynamic damage of such materials under high-strain-rate loading conditions cannot straightforwardly be predicted on the basis of existing knowledge about traditional metals, since the fundamental damage mechanisms can be quite different. Therefore, combined experimental, computational, and theoretical efforts are needed to investigate the mechanism of damage initiation and evolution of these new materials.

Finally, the study of ductile damage initiation and evolution covers spatial scales from nanometers to millimeters and temporal scales from picoseconds to microseconds. Therefore, it is a great challenge to establish a multiscale theoretical framework to bridge the microscopic physical mechanism and the macroscopic mechanical response. In addition, given the spatial and temporal resolutions required to monitor the variation of microstructures, more advanced experimental and simulation techniques are needed to accurately determine the related information. Such techniques should reduce the number of phenomenological parameters that need to be used in the multiscale theoretical framework and thereby allow a more precise description of the fracture process.

Acknowledgment. Financial support for this work was provided by the Science Challenge Project (Grant No. TZ2018001) and the National Natural Science Foundation of China (Grant Nos. 11988102, 11632001, 11521202, and 12002005).