Micromechanics of Void Nucleation and Early Growth at Incoherent Precipitates: Lattice-Trapped and Dislocation-Mediated Delamination Modes

Abstract

The initial stages of debonding at hard-particle interfaces during rupture is relevant to the fracture of most structural alloys, yet details of the mechanistic process for rupture at the atomic scale are poorly understood. In this study, we employ molecular dynamics simulation of a spherical Al₂Cu θ precipitate in an aluminum matrix to examine the earliest stages of void formation and nanocrack growth at the particle-matrix interface, at temperatures ranging from 200–400 K and stresses ranging from 5.7–7.2 GPa. The simulations revealed a three-stage process involving (1) stochastic instantaneous or delayed nucleation of excess free volume at the particle-matrix interface involving only tens of atoms, followed by (2) steady time-dependent crack growth in the absence of dislocation activity, followed by (3) dramatically accelerated crack growth facilitated by crack-tip dislocation emission. While not all three stages were present for all stresses and temperatures, the second stage, termed lattice-trapped delamination, was consistently the rate-limiting process. This lattice-trapped delamination process was determined to be a thermally activated brittle fracture mode with an unambiguous Arrhenius activation energy of 1.37 eV and an activation area of 1.17 Ų. The role of lattice-trapped delamination in the early stages of particle delamination is not only relevant at the high strain-rates and stresses associated with shock spallation, but Arrhenius extrapolation suggests that the mechanism also operates during quasi-static rupture at micrometer-scale particles.

1. Introduction

Void nucleation is the first step towards fracture in many different contexts, including quasi-static tearing, dynamic spall, creep rupture, irradiation creep, and wear debris generation. Voids are predominantly thought to nucleate at second phase particles, either when the particles crack or when the interface between the particle and matrix debonds [1]. Here, we focus on the earliest stages of void nucleation via particle delamination. Subsequent to nucleation, these voids grow until they induce fracture. While many studies have employed continuum models, such as finite element modeling [2,3,4,5], to evaluate the process of void nucleation, the atomistic mechanisms governing nucleation are less well studied. And yet, since void nucleation begins at the nanoscale, it is intrinsically an atomistic process [6]. Our goal in this work is to evaluate void nucleation in a model system with an incoherent, second-phase particle (θ-particle in Al) in an effort to reveal the underlying micromechanics and kinetically limiting processes.

Given its central role in fracture, continuum damage models commonly invoke void nucleation in their underlying formalisms. Perhaps the most popular approach is the porous plasticity model of Gurson, Tvergaard, and Needleman [7] which utilizes a yield criterion that is a function of the void volume fraction, $f$. In this model, the void volume fraction is governed by a differential equation of the form:

$$\dot{f}={\dot{f}}_{\mathrm{nuc}}+{\dot{f}}_{\mathrm{growth}}$$

where ${\dot{f}}_{\mathrm{nuc}}$ and ${\dot{f}}_{\mathrm{growth}}$ account for the contributions from nucleation and growth events, respectively. Ideally, the term ${\dot{f}}_{\mathrm{nuc}}$ would be derived from a fundamental, micromechanics-based understanding of void nucleation. Given that this understanding is lacking, a more phenomenological approach is commonly utilized. For example, it is commonly assumed that void nucleation occurs at particles when a critical plastic strain, ${\epsilon }_c^P$, is reached, and that ${\epsilon }_c^P$ varies from particle to particle according to a probability density function $F\left({\epsilon }^P/{\epsilon }_c^P\right)$ where ${\epsilon }^P$ is the equivalent plastic strain [8,9]. The void volume fraction then increases in time as a result of void nucleation according to the expression:

$${\dot{f}}_{\mathrm{nuc}}=F\left({\epsilon }^P/{\epsilon }_c^P\right){\dot{\epsilon }}^P$$

where ${\dot{\epsilon }}^P$ is the equivalent plastic strain rate. Usually, $F\left({\epsilon }^P/{\epsilon }_c^P\right)$ is assumed to be a normal distribution with mean ${\overline{\epsilon }}_c^P$ and standard deviation $s_{{\epsilon }_c}$ [8], however there is no direct evidence which justifies this choice. Furthermore, ${\overline{\epsilon }}_c^P$ and $s_{{\epsilon }_c}$ are treated as empirical parameters that are fitted against experimental data (e.g., stress-strain curves). While this and other similar phenomenological approaches have been applied pervasively across the literature, the Sandia Fracture Challenges have recently shown that these models often fare poorly when making blind fracture predictions [10,11,12]. This motivates a deeper look at the micromechanics of void nucleation, so that strong assumptions about what governs fracture (a critical strain?) and how the propensity for fracture varies across the population of particles (normally distributed?) can be lifted.

Unfortunately, the critical strain ${\epsilon }_c^P$ is not easily studied using micromechanical simulations because plastic strain is really a homogenized, macroscale concept; at the microscale where discrete dislocations interact with particles, plastic strain is not a very relevant concept. On the other hand, some damage mechanics models employ a critical stress ${\sigma }_c$ at which nucleation occurs [9], which is more consistent with micromechanics modeling (the stress state can be specified in molecular dynamics, for example). Here, we argue, however, that rather than focusing on a “critical stress” at which void nucleation occurs, it makes more sense to consider how the nucleation rate varies with stress state, temperature, etc. In other words, void nucleation can occur over a range of stresses, with the nucleation rate increasing as the stress is increased. This view is more consistent with other works on crack nucleation, which focus on the nucleation rate [13,14]. Within this view, the critical stress is the stress at which the nucleation rate goes to infinity, meaning that nucleation occurs instantaneously. We note that the nucleation rate for a given state is likely only well-defined in an average sense, because nucleation is a stochastic phenomenon. This means that at each state, there is a distribution of nucleation rates (which could be interpreted in probabilistic terms). We argue that the possibility of “subcritical” nucleation, i.e., with $\sigma <{\sigma }_c$, and the statistical aspects of nucleation could be important to the overall nucleation process. For these reasons, our focus here is on the stress and temperature dependence of the void nucleation rate.

An important nuance to the study of void nucleation is deciding when exactly a void is said to “nucleate.” As soon as a crack appears within the particle-matrix interface? Or after a significant fraction of the interface has delaminated? We may expect that a clear nucleation event occurs whereby a crack “suddenly” appears along the interface, allowing us to disentangle this terminological ambiguity, although the appearance of a crack is often limited by the spatial and temporal resolution of the techniques employed. In the present approach, with atomic-scale and picosecond resolution, the initial emergence of a crack is still difficult to define: we observed steady growth of an interfacial crack starting from a vacancy-sized nucleus. We were unable to determine the precise mechanism by which the vacancy-sized nucleus appeared, however. Furthermore, the appearance of the vacancy-sized nucleus did not control the kinetics of void nucleation. Instead, we found that it was the subsequent growth of the crack that governed the overall delamination (e.g., void nucleation) process. Hence, we find that it is the delamination rate, controlled by the growth of a crack, which governs the void nucleation rate. For this reason, we refer to our simulations as studying void nucleation and “early growth.”

Void nucleation has been studied in perfect crystals [15,16,17,18,19,20,21,22], at grain boundaries [23,24], ahead of crack tips [25], and at second-phase particles [26,27,28,29] using molecular dynamics. In most cases, void nucleation results from interactions between several crystallographic defects, such as grain boundaries and twins/dislocations [23,24], pairs of intersecting stacking faults [21], and particles and dislocations [26,27,30]. The previous work on particle-mediated void nucleation is most relevant here. Coffman et al. [31] studied void nucleation in Si under uniaxial tension with a cubic nanograin “particle” that delaminated from the surrounding matrix. They first performed atomistic simulations to calibrate a continuum fracture model (a cohesive zone model), and then compared atomistic and continuum predictions of void nucleation. In general, they found poor agreement between the models, motivating the need for further studies of void nucleation with atomistic resolution. Pogorelko and Mayer [26,27,28,29] and Cui and Chen [30] studied void nucleation at spherical particles in a variety of material systems, considering the influence of strain rate, temperature, simulation box size, and particle volume fraction on the delamination behavior under a fixed uniaxial strain rate. In simulations with face-centered cubic (FCC) matrices [26,27,29,30], nucleation was observed to occur in two stages: first a crack nucleated at the top and bottom poles along the loading axis (similar to the behavior predicted by continuum models [2]), and then after some subsequent growth dislocations were emitted from the crack tips. On the other hand, nucleation with body-centered cubic and hexagonal close-packed (HCP) matrices seemed to initiate from defects in the matrix rather than at the particle-matrix interface [28]. While the tensile strength of these systems has been characterized extensively using these simulation results, the nucleation rate could not be estimated because of the fixed-strain-rate boundary conditions.

Our study here had two goals: (1) to assess the stress and temperature dependence of the void nucleation and early growth rate with MD and (2) identify the micromechanical processes which govern the kinetics. In contrast to previous work [26,27,28,29,30], we perform simulations here with a fixed stress state (and temperature), so that our results can be used to estimate the stress and temperature-dependent rates. Our findings indicate that void nucleation may be rate limited by the kinetics of crack growth processes rather than the kinetics of crack nucleation (e.g., the time it takes for a crack to appear). Furthermore, we show two distinct delamination modes with drastically different growth kinetics. Finally, we conclude that a brittle delamination mode which we term lattice-trapped delamination may be an important contributor to void nucleation. While the stress range employed here (5.7 to 7.2 GPa) is high relative to quasi-static loading, it is in the range where shock spallation is observed under high loading rates [32,33]. Furthermore, through a thermal activation analysis of our data we are able to extrapolate our results down to lower stress conditions. The remainder of the work is organized as follows. In Section 2 we discuss our simulation setup and analysis methods, in Section 3 we present our results, in Section 4 we discuss the implications of our findings and compare results with existing theories, and finally conclude the manuscript in Section 5.

2. Materials and Methods

As a model system, we consider void nucleation at θ-particles in an FCC Al matrix. θ is the thermodynamically stable intermetallic phase of the Al-Cu system and is commonly observed in Al-Cu-copper alloys (e.g., 2xxx series) in the overaged condition [34]. θ-particles have a composition of Al₂Cu and a body-centered tetragonal C16 crystal structure. They are incoherent with the matrix and typically adopt plate-like geometries [35]. However, for simplicity in this work we will use a spherical precipitate geometry. While there is evidence that voids may nucleate at θ-particles [36], we emphasize that we are using the θ-Al system as a model incoherent precipitate system with the goal of gaining general insight into the micromechanics and kinetics of void nucleation.

MD simulations were performed using LAMMPS [37] with the Al-Cu angular-dependent interatomic potential of Apostol and Mishin [38]. The angular-dependent potential framework is an enriched version of the embedded atom method that enables incorporation of angular-dependent interactions. These interactions enable the potential to capture the lattice constants, anisotropic elastic moduli, and formation energy of the θ-phase with reasonable accuracy. Our simulation cell geometry is shown in Figure 1a; we initially inserted an incoherent spherical particle with a radius of $R=50$ Å into a pure Al lattice using the zero-temperature lattice parameters predicted by the potentials. The Al lattice and θ-particle are oriented so that their unit cell axes are aligned with the simulation box. The c-axis of the θ-lattice is oriented in the z-direction of the simulation cell. Periodic boundary conditions were used in all directions with a 200 Å cubic simulation cell. The sequence of each simulation is shown in Figure 1b. During the relaxation stage, we used an NPT (constant number of atoms N, pressure P, and temperature T) ensemble and simulated 2 ps at the chosen temperature and zero hydrostatic stress. Subsequently, during the ramping stage we ramped the hydrostatic stress up to the target value ${\sigma }^H$ over a duration of 23 ps. This duration was chosen empirically; if the stress was ramped too quickly we observed “premature” void nucleation, likely because of stress spikes resulting from imperfect performance of the barostat. See Appendix A for additional information. Finally, during the holding stage, the hydrostatic stress was held constant for the duration of the simulation until the particle completely debonded from the matrix or the simulation terminated after 1-week of wall time. All simulations used a thermostat damping parameter value of 0.01 ps, barostat damping parameter value of 1 ps, and a time step size of 0.001 ps.

We note that while the MD barostat controls the average (virial) stress state in the simulation cell, the local stress state may vary. In fact, we expect there to be variation because the elastic constants between the matrix and particle differ, i.e., this is an Eshelby inhomogeneity problem [39]. Furthermore, the particle images resulting from periodic boundary conditions will interact with each other, further complicating the stress field. These effects are quite complex, especially given the anisotropic nature of the C16 θ-phase. For simplicity, in our analysis of the data we assume that the applied hydrostatic stress ${\sigma }^H$ dominates the delamination behavior; our successful thermal activation analysis below justifies this assumption. We note that Pogorelko and Mayer have analyzed the spatially varying stress field near a second-phase particle under uniaxial loading and its influence on the delamination process [26,27,29].

Results from a total of 290 MD simulations are reported in this study at temperatures ranging from 200 to 400 K and stresses in the range of 5.7 to 7.2 GPa (the precise stress range differed for each temperature). In most cases, 10 simulations with different thermalization histories (initial atomic velocities) were performed at each stress-temperature condition. In many simulations, we observed nucleation of dislocations at the particle interface. To enable efficient detection of the appearance of Shockley partial dislocations in our simulation cell, we exploited the fact that atoms situated in stacking faults (e.g., produced by a Shockley partial dislocation) appear as HCP atoms when analyzed via common neighbor analysis (CNA) [40]. Hence, by simply monitoring the number of “HCP” atoms in the simulation cell $N_{HCP}$, we could identify when a dislocation appeared. We note that there was always a small, non-zero number of HCP atoms detected, due to thermal noise in the lattice and the imperfect detection capacity of CNA. To prevent false detection of a dislocation, we established a threshold value for the appearance of a dislocation, $N_{HCP}^d$, based on empirical analysis of our data. We set this threshold at $N_{HCP}^d=40$ HCP atoms, so if $N_{HCP}>40$ we “detected” appearance of a dislocation. Furthermore, we applied a moving average to the raw $N_{HCP}$ vs. time data using a window of width 0.005 ps. This served to smooth out the data a bit and remove spurious spikes in $N_{HCP}$ which did not lead to a sustained increase in $N_{HCP}$ over time (as was expected if a dislocation had nucleated and remained in the system). This approach to dislocation detection was validated by manually analyzing several datasets in OVITO Pro [41]. All simulation snapshots were produced using OVITO Pro [41].

3. Results

We begin the Results section by discussing the overall behaviors observed in our simulations, showing that two modes of delamination were observed which we call lattice-trapped delamination and dislocation-mediated delamination.

3.1. Void Nucleation Process

In our simulations, we generally observed the sequence of events depicted in Figure 2. At the start of the load holding step there was no evidence of cracking or voiding at the particle interface. After some time, one or more clusters of atoms began exhibiting relatively large atomic volumes, indicating the nucleation of a crack with the size of a small vacancy cluster. For example, in Figure 2a we show atoms in blue whose atomic volume exceeds 30 ų based on Voronoi analysis. For reference, the atomic volume of aluminum is 16.7 ų. The specific mechanism by which this vacancy clustered appeared could not be determined. Over time this crack grew, leading to a larger patch of atoms with volumes exceeding 30 ų as shown in Figure 2b. Importantly, this crack growth was not accompanied by any dislocation activity. Instead, the crack grew steadily in time; in the Discussion we demonstrate that this delamination rate is governed by the lattice trapping phenomenon [42], and hence refer to this as lattice-trapped delamination. Eventually, once the crack reached a critical size, Shockley partial dislocations nucleated at the tip of the crack approximately in the plane of the crack, as shown in Figure 2c. These dislocations then rapidly glided away from the particle into the bulk and began to multiply, joined by additional dislocations nucleating from the crack tip. The crack growth rate increased rapidly upon appearance of the dislocations, leading to total delamination of the particle from the matrix. We call this process dislocation-mediated delamination.

Figure 3 shows atomic displacements associated with the two delamination modes. Figure 3a shows the atomic displacement vectors as black arrows over an 18 ps time period from the snapshot in Figure 2b to the moment just before dislocation nucleation Figure 2c. As shown, most of the crack growth is accommodated by displacements of atoms at the crack tip along the circumference of the particle. These displacements under lattice-trapped delamination are rather incoherent, in the sense that their direction and magnitude vary along the crack front in an uncoordinated manner. For example, several atoms experience large displacements while their immediate neighbors do not displace much at all (see yellow circled atoms). Displacements associated with dislocation-mediated crack growth over a 2 ps time period between the snapshots shown in Figure 2c,d are shown in Figure 3b. Once again, the displacements largely occur at the crack tip and in the circumferential direction, however these displacements are more coordinated in their direction and magnitude. These displacements appear to be due to nucleation and glide of the dislocations visible in the figure, i.e., they are generally aligned with the Burgers vectors.

To determine whether there were preferential nucleation sites on the particle’s surface, we extracted the approximate crack nucleation location from 40 simulations and plot these locations in Figure 4 as a point cloud projected onto the x-z plane. This figure shows that while there may be a slight preference for nucleation at the negative z-axis pole of the precipitate (since there is a small cluster there), nucleation was also common at other points around the surface. Similarly, there is a lack of data points at the positive z-axis pole, indicating nucleation there is unfavorable. These results imply that our boundary conditions, simulation cell size, and precipitate orientation did not significantly influence the simulation behaviors (i.e., they did not introduce strong preferential sites).

To assess the delamination behavior in the absence of thermal fluctuations, we performed molecular statics simulations of hydrostatic straining. We progressively increased the hydrostatic strain of the box by increasing the volume in 0.003% increments and minimizing the energy of the system after each strain increment. Each minimization step iterated until the change in energy during a minimization step was less than 10⁻⁶% or the norm of the global force vector was less than 10⁻⁸ eV/Å. The resulting stress-strain curve is shown in Figure 5. We observe that at a hydrostatic stress of around 10 GPa, the particle catastrophically delaminates from the matrix. Hence, 10 GPa can be regarded as the athermal critical stress for void nucleation. The first peak in Figure 5 corresponds to the nucleation of dislocations at the poles of the particle along the z-axis. The subsequent peaks correspond to nucleation of dislocations around the entire circumference of the particle. Hence, it seems that athermal void nucleation is governed by the athermal nucleation of dislocations.

To get information about the crack growth, we estimate the crack volume as $V_{crack}\left(t\right)=V\left(t\right)-V_0$, where $V\left(t\right)$ is the volume of the simulation cell at time t and $V_0$ is the volume at the start of the holding phase. In Figure 6 we show a few examples of how the crack volume evolves over time during the load holding phase. In some cases, there appears to be an “incubation period” where the volume does not increase at all, followed by a gradual increase over time indicating the nucleation and growth of a crack. This gradual increase corresponds to the lattice-trapped delamination mode. In other cases, the volume appears to increase from the start of the load holding phase, with no obvious incubation period. In most cases it was difficult to unambiguously identify a clear “nucleation” event which correlated with a local increase in atomic volume at the void’s surface. For this reason, we were unable to analyze any sort of “nucleation” rate directly from the incubation time. Regardless, lattice-trapped delamination was always observed in our simulations. We also mark in Figure 6 the time where the first dislocation nucleated. Upon nucleation of one or more dislocations, the system volume increases precipitously as the crack growth rate accelerates to complete delamination.

3.2. Kinetics of Lattice-Trapped Delamination

To characterize the kinetics of lattice-trapped delamination, we extracted the crack growth rate in terms of the rate of increase of the crack volume as shown in Figure 6. Since no other defects were present in the simulation cell, it is expected that essentially all increases in volume must be due to growth of interfacial cracks. Visual inspection of simulations and identification of volume “hot spots” where atomic volumes were elevated confirmed this assumption; we never observed significant volume increases in the bulk, only at the interface near cracks. It is more conventional, of course, to characterize a crack in terms of its area or radius. However, since cracks typically adopted complex shapes and morphologies, it was difficult to directly extract their area and/or radius. For these reasons, we use changes in simulation cell volume to quantify the delamination rate.

We determined the lattice-trapped delamination rate by taking a linear regression of the volume vs. time curve from the start of the load holding phase to the time of dislocation nucleation, shown by the dashed lines in Figure 6. We find that under the same nominal conditions (stress and temperature), the delamination rate was sensitive to the initial conditions (e.g., initial random atomic velocities). For this reason, 10 replicate simulations were performed at each condition. Figure 7 shows histograms of the lattice-trapped delamination rate at a few conditions with 30 replica simulations, demonstrating the spread of the data. Interestingly, the delamination rates are systematically biased towards lower rates (rather than being symmetrical about the mean), qualitatively taking the form of exponential distributions. Figure 8 shows the mean delamination rate as a function of ${\sigma }^H$ at temperatures of 200, 250, 300, 350, and 400 K. There is clearly a strong sensitivity to both stress and temperature. We will further analyze these data in the Discussion.

3.3. Analysis of Dislocation-Mediated Delamination

While we always observed lattice-trapped delamination in our simulations, dislocation nucleation was only observed in “high stress” simulations. Furthermore, if the stress was too high then dislocations would nucleate at the beginning of the simulation, much like the behavior observed in our molecular statics simulation. One important question is: what is the interplay between lattice-trapped and dislocation-mediated delamination? What governs the transition from one mode to the other? To gain insight into this question, we extracted the crack volume at the moment of dislocation nucleation from our simulations and present the average values in Figure 9. At lower stresses, the crack volume increases to a peak value and then subsequently decays to zero, indicating immediate nucleation of dislocations. We believe that these trends result from changes in the lattice-trapped delamination and dislocation nucleation rates with stress and temperature, and changes in the driving force for dislocation nucleation as the crack grows. We defer further analysis to the Discussion.

In all simulations, as soon as the first dislocation nucleation event occurred, the delamination rate increased dramatically. The delamination rate was so large that a well-defined rate could not be determined. In some cases, dislocations nucleated from the “crack tip” near the particle interface, but not in all cases. For example, dislocation nucleation in the middle of the crack face in the Al matrix was also observed.

4. Discussion

4.1. Thermal Activation Analysis

Here we further analyze the lattice-trapped delamination rate data and demonstrate that lattice-trapped delamination is thermally activated. For a system containing a crack of length $a$ loaded by hydrostatic stress ${\sigma }^H$, the theory of thermally activated crack growth says that the growth rate is [43]

$$\dot{a}\left({\sigma }^H,a\right)={\dot{a}}_0\exp \left(\frac{-E_a+A^{*}K{\left({\sigma }^H,a\right)}^2/E}{k_{B}T}\right)$$

(1)

where $E_a$ is the activation energy for bond breaking, $K$ is the stress intensity factor, $A^{*}$ is the activation area, $E$ is the modulus of elasticity, ${\dot{a}}_0$ is the exponential prefactor related to the attempt frequency, $T$ is the absolute temperature, and $k_B$ is Boltzmann’s constant. It is important to note that Equation (1) is only valid when $K>K_c$, the critical stress intensity factor. If $K<K_c$, then the free energy of the system increases when the crack grows, violating the second law of thermodynamics [44]. The stress intensity factor can always be written in the form

$$K=Y{\sigma }^H\sqrt{\pi a}$$

(2)

where $Y$ is a geometric factor dictated by the geometry of the problem. To relate the crack volume to the crack length, we approximate the crack as an ellipsoid with two axes of radius $a$ and the other of radius $a/2$ (consistent with cracks observed in our simulations). The volumetric crack growth rate is then $\dot{V}=4\pi a^2\dot{a}$. In our simulations where lattice-trapped delamination occurs, the crack length does not increase significantly (cracks remain relatively small during lattice-trapped delamination, see Figure 2). Hence, for simplicity, we neglect changes in $a$ and assume that $K=K\left({\sigma }^H,\overline{a}\right)$ and $\dot{V}=4\pi {\overline{a}}^2\dot{a}$, where $\overline{a}$ is the average crack length during the simulation. With this assumption and using Equation (1), we obtain that the activation enthalpy (numerator in the exponential) is

$$\Delta H_a=E_a-C{k}_B{\left({\sigma }^H\right)}^2$$

(3)

where $C=A^{*}{Y}^2\pi \overline{a}/\left(E{k}_B\right)$ and simple algebraic manipulation further shows that

$$\ln \dot{V}=\ln {\dot{V}}_0-\frac{E_a}{k_{B}T}+C\frac{{\left({\sigma }^H\right)}^2}{T}$$

(4)

where ${\dot{V}}_0=4\pi {\overline{a}}^2{\dot{a}}_0.$ Hence, if we plot ln$\langle \dot{V}\rangle$ from our simulations as a function of ${\left({\sigma }^H\right)}^2/T$, the dataset for each temperature should form a straight line with slope $C$ if growth is thermally activated. In Figure 10 we plot the data in this way and see a consistent linear trend across all datasets. Specifically, we find that the same slope fits all datasets, indicating that $C=258$ K/GPa². Next, we extract the y-intercept from each of these linear fits, and Equation (4) indicates that these intercepts should scale with $1/T$. Figure 11 plots the y-intercepts from Figure 10 as a function of $1/T$ and once again a linear behavior is recovered as expected. The slope of Figure 10 is $-E_a/k_B$ and the y-intercept is $\ln {\dot{V}}_0$; we obtain values of $E_a=1.37$ eV and ${\dot{V}}_0=1.23\times {10}^9$ ų/ps. The fact that our data so strongly reproduces the behaviors predicted by Equation (4) indicates that the lattice-trapped delamination observed in our MD simulations is indeed thermally activated, and that our neglect of the crack length dependence of $K$ does not introduce any significant errors into our analysis.

The thermally activated delamination mode that we observe is likely governed by the so-called lattice trapping phenomenon [42]. Under lattice trapping, a crack which is loaded supercritically (i.e., with $K>K_c$) grows in a step-wise manner as atomic bonds at the crack tip are sequentially broken, often by a kink-pair mechanism [43,44,45,46]. As the crack grows, it experiences an oscillating potential energy landscape due to the periodicity of the lattice, and the height of these oscillations dictates the activation energy for growth. To our knowledge this crack growth mode has only been observed in brittle materials like Si [45] and glass [47], but not in ductile materials like Al considered here. This is nonetheless reasonable, since in effect the Al system is acting in a brittle manner in the absence of dislocation activity.

To further analyze our extracted parameters, we need to estimate the stress intensity factor. Tan and Gao [48] numerically determined the stress intensity factor for an axisymmetric interfacial crack on a sphere which forms an angle $\varphi$ from the radial axis, as shown in the inset in Figure 12a, with various modulus ratios $E_p/E_m$ where m stands for matrix and p for precipitate. We can express their results in the form

$$K_0=Y\left(\varphi ,E_p/E_m\right){\sigma }^H\sqrt{\pi R}$$

(5)

where $K_0=\sqrt{K_I^2+K_{II}^2}$ is the “overall” stress intensity factor for the mixed-mode loading (the crack is generally mixed mode) and $Y\left(\varphi ,E_p/E_m\right)$ was determined by Tan and Gao via numerical boundary integral methods. Using Equation (5) enables us to obtain the activation area as

$$A^{*}=C \frac{k_B\overline{E}}{\pi RY{\left(\varphi ,E_p/E_m\right)}^2}$$

(6)

where $\overline{E}=2E_m{E}_p/\left(E_m+E_p\right)$ is the bimaterial modulus for interfacial fracture [49]. For approximation purposes, we employ the typical Young’s modulus of untextured polycrystals instead of the anisotropic single crystal elastic constants. Estimating $E_m=70$ GPa (experimental value for pure Al) and $E_p\approx 120$ GPa [50] gives $\overline{E}\approx 88$ GPa and $E_p/E_m\approx 1.7$. Tan and Gao found that for cracks varying from $\varphi =22.5$ to $67.5$ degrees with a modulus ratio of $E_p/E_m=2$, $Y\left(\varphi ,E_p/E_m\right)$ varied from 1.31 to 1.79. Unfortunately, cracks observed in our simulations typically had angles $\varphi <22.5°$, so it is difficult to apply Tan and Gao’s solution. Note that according to Equation (5), $Y\to 0$ as $\varphi \to 0$ since $K_0$ must go to zero when the crack length is zero (i.e., $\varphi =0$). Hence, we expect the $Y$ values in our simulations to be less than 1.31. Nonetheless, to gain insight into orders of magnitude for the thermal activation parameters we assume $Y\approx 1.31$ in the analysis below. Using these parameter values with $R=50$ Å in Equation (6) leads to $A^{*}\approx 1.17 {\AA }^2$. According to Schoeck [43]

$$A^{*}=\left(1+\beta \right)\Delta A$$

(7)

where $\Delta A$ is the atomistic area of crack advance between the equilibrium and saddle position of the crack front and $\beta$ is a factor in the range between 0 and 2. Schoeck argued that for crack advance by breaking of individual bonds, $\Delta A\approx 1 {\AA }^2$; hence, our $A^{*}$ value gives the correct order of magnitude, further bolstering our conclusion that lattice-trapped delamination is thermally activated.

Finally, we note that an important aspect of thermally activated crack growth is the athermal stress intensity factor, $K_{\mathrm{ath}}$, at which the activation enthalpy goes to zero. At and above this load, the crack growth rate is no longer governed by thermally activated bond breaking. According to Equation (1) and accounting for the modulus mismatch, the athermal stress intensity factor is

$$K_{\mathrm{ath}}=\sqrt{\frac{E_a\overline{E}}{A^{*}}}.$$

(8)

Using our estimates for parameters above, we obtain that $K_{\mathrm{ath}}\approx 1.29$ MPa·m^1/2. According to Equation (5), $K_{\mathrm{ath}}$ is reached when

$${\sigma }_{\mathrm{ath}}^H=\frac{K_{\mathrm{ath}}}{Y\sqrt{\pi R}}$$

(9)

and using the estimates above we obtain ${\sigma }_{\mathrm{ath}}^H\approx 7.8$ GPa. This value is consistent with our data since we did not observe lattice-trapped delamination above 7.2 GPa (although we did not attempt to obtain a maximum stress where lattice-trapped delamination rates could be obtained).

4.2. Dislocation Nucleation

Delamination and crack growth via dislocation nucleation have been observed by other researchers in the past [30,51]. Similar to the thermally activated lattice-trapped delamination mode discussed above, dislocation nucleation is thermally activated. Above some critical load, the activation barrier for dislocation nucleation goes to zero and spontaneous nucleation occurs [51]. Hence, we have a race between two thermally activated processes—lattice-trapped delamination and dislocation nucleation—each characterized by their own activation enthalpy which varies with the stress intensity factor. This race gives rise to the crack volume trends exhibited in Figure 9. At lower stress, lattice-trapped delamination is slow so that the crack does not grow much before a dislocation is nucleated. As the stress is increased, the lattice-trapped delamination rate increases and there is a time lag before the dislocation nucleates, causing the crack volume at dislocation nucleation to increase. Finally, at very high stresses the athermal load is reached for dislocation nucleation so that a dislocation nucleates immediately during the simulations, causing the crack volume at dislocation nucleation to go to zero. This qualitatively explains the trends observed in Figure 9. A quantitative explanation would require a detailed understanding of the nucleation parameters for a dislocation at the crack tip, which are non-trivial to compute and beyond the scope of this work [51].

Unfortunately, the conditions of the simulations here were such that as soon as a dislocation nucleated, the particle delaminated almost instantaneously. This indicates that our systems were overdriven with the respect to dislocation-mediated delamination. It is possible that quantitatively useful delamination rates could be obtained under dislocation-mediated delamination from simulations at lower stresses if a large enough pre-crack is manually inserted into the system. Here, if our applied stress was too low, the lattice-trapped delamination rate was so small that a crack never became large enough to enable dislocation nucleation. Future work should focus on quantifying the dislocation-mediated delamination rate.

4.3. Implications for Damage and Rupture of Materials

The present simulations bear direct relevance to spall formation under shock loading conditions; for example, spall by microvoid coalescence has been observed in 1100-O aluminum at shock stresses as high as 6 GPa [33]. However, it is reasonable to ask if the present observations also relate to quasi-static ductile rupture of aluminum alloys. The main gaps between our simulation conditions and those of typical experiments are that our applied stresses are much higher (for example, the yield strength of 2xxx series aluminum alloys is in the range of 150–350 MPa depending on alloy and heat treatment) and our particles sizes are smaller.

First, we address dislocation-mediated growth. The reality is that aside from high rate loading scenarios where high stresses are attained [52] and highly localized loading (e.g., under a nano-indenter), dislocation nucleation is typically far too slow to meaningfully impact material behaviors [51,53]. Hence, we do not expect that dislocation-mediated growth, as observed here through nucleation of dislocations, will be relevant under most loading scenarios. Indeed, even within the present spectrum of simulations, the dislocation-mediated growth only occurred at the highest stresses approaching the athermal critical stress. On the other hand, Sills and Boyce recently showed that dislocation-mediated growth processes are greatly accelerated when dislocations are already present in the material, and thus do not require nucleation [54]. Kinematically speaking, it is equivalent to nucleate a dislocation at a crack tip or to adsorb an existing dislocation of opposite sign. We believe that dislocation-mediated growth via adsorption of existing dislocations is likely to play an important role in particle delamination, just as Sills and Boyce argue it does for void growth.

In terms of the relevance of lattice-trapped delamination, we are fortunate to have a fully characterized thermal activation law for the lattice-trapped delamination rate. This means that we can extrapolate the model to quasi-static conditions since the physics of thermal activation is equally valid at lower stresses and larger particle sizes. Importantly, the process of lattice-trapped delamination is bookended by two conditions on the stress intensity factor; namely, $K_c<K<K_{\mathrm{at}}$. Below the critical stress intensity factor, $K_c$, delamination is not energetically favorable. And above the athermal stress intensity factor, $K_{\mathrm{at}}$, thermal activation is not operable. Equation (9) provides an estimate for the athermal hydrostatic stress ${\sigma }_{\mathrm{ath}}^H$ at which $K_{\mathrm{at}}$ is reached. For brittle fracture, $K_c=\sqrt{\left({\gamma }_{sm}+{\gamma }_{sp}\right)\overline{E}}$ where ${\gamma }_{sm}$ and ${\gamma }_{sp}$ are the surface energies for the matrix and precipitate, respectively. The surface energy for Al is known experimentally to be ${\gamma }_{sm}=0.98$ J/m² [55]. To estimate the surface energy for θ-phase, we use the average value for the interatomic potential used here (computed for (100) and (110) surfaces), giving ${\gamma }_{sp}=1.38$ J/m² [38]. These values give a critical stress intensity factor of $K_c=0.457$ MPa·m^1/2. Using these results, we estimate the stress and particle size range where lattice-trapped delamination may operate (i.e., where $K_c<K<K_{\mathrm{at}}$) when $\varphi =22.5°$ in Figure 12a. The plot shows that in the particle size range considered here, stresses must exceed 2.8 GPa for lattice-trapped delamination to operate. On the other hand, for particles with a radius of 1 μm, lattice-trapped delamination can occur in the range of 200 to 550 MPa, which falls within the stress range relevant to quasi-static loading conditions, especially considering that inhomogeneous microstructural stresses can far exceed the homogeneous far-field applied stresses. Using this stress range, we plot the predicted lattice-trapped delamination rate $\dot{a}$ as a function of stress for particle radii of 1 and 10 μm at temperatures of 300, 500, and 700 K. Interestingly, the resulting delamination rates are large enough that they may be relevant to applications. For example, full delamination of a 1 μm particle requires ~3 μm of crack growth (half of the circumference), and for this to occur in 1 years’ time would require a delamination rate of just $\dot{a}\approx {10}^{-13}$ m/s. We emphasize that because we have made a number of approximations in our analysis, these delamination rates should be interpreted only in terms of their rough order of magnitude. Nonetheless, these results indicate that lattice-trapped delamination via thermally activated brittle crack growth may be broadly relevant to void nucleation.

While we have only considered the influence of hydrostatic loading on the delamination behavior here, shear stresses are also expected to affect delamination. For example, experimental work by Croom et al. [56] and Achouri et al. [57] has demonstrated void nucleation under shear-dominated conditions in pure Cu and particle-containing Al, respectively. In the case of particle-mediated nucleation, in addition to the applied shear stress the non-uniform plastic strain accumulation around the particle provides a driving force for delamination. In the context of lattice-trapped delamination, this driving force is difficult to quantify since it requires an elastic-plastic analysis of deformation with a significant accumulation of plastic strain. None-the-less, lattice-trapped delamination could also be operative under shear-dominated loading, but additional research is necessary to quantify the influence of shear stresses and shear deformation.

Perhaps the most confounding aspect of our simulations is that we set out to study crack nucleation and instead wound-up studying crack growth. One might have expected that the rate controlling step of void nucleation would be the appearance of an interfacial crack, in the sense that a crack spontaneously appeared along the interface. However, this was not the behavior we observed. Instead, we found that excess free volume in the form of vacancy-type defect clusters at the interface appeared quickly and then grew by lattice-trapped delamination. The majority of the simulation time was then spent growing the crack by lattice-trapped delamination until dislocation nucleation occurred. Hence, for our simulations, the rate controlling step for void nucleation was actually lattice-trapped delamination. Once the first dislocations appeared, it is true that the void fully nucleated in a catastrophic manner more consistent with a true nucleation event, but we believe that this may be an artifact of our high stresses; in our lower stress simulations, the dislocation-mediated delamination stage never occurred. We note that additional research is necessary to determine the mechanism for vacancy-type cluster nucleation so that a comprehensive picture for void nucleation can be assembled. While the present observations are directly applicable to the strain-rates associated with shock-induced spallation, extrapolation of the thermally activated process of lattice-trapped delamination suggests that the process is also relevant at quasi-static timescales for micrometer-sized particles.

5. Conclusions

We observed that early stage nanoscale delamination at a particle-matrix interface was a three-stage process: (1) a nuclei of excess free volume is formed at the particle-matrix interface, either immediately upon loading or after a brief delay; (2) subsequently, the nuclei grows in the absence of dislocation activity in a process we term lattice-trapped delamination; (3) when the crack reaches sufficient size and the crack tip stress intensity is sufficient, Shockley partial dislocations are emitted causing a rapid acceleration of crack growth and ultimately complete debonding of the particle-matrix interface. The second stage, lattice-trapped delamination, as described for the first time herein, limits the rate of early-stage growth of the nanocrack, and is not only relevant to shock spallation, but extrapolation of the kinetics of this process suggest that it is also relevant at quasi-static strain rates for micrometer-scale particles.