Why are some Interfaces in Materials Stronger than others?

Grain boundaries (GBs) are often the preferred sites for void nucleation in ductile metals. However, it has been observed that all boundaries do not contribute equally to this process. We present a mechanistic rationale for the role of GBs in damage nucleation in copper, along with a quantitative map for predicting preferred void nucleation at GBs based on molecular dynamics simulations in copper. Simulations show a direct correlation between the void nucleation stress and the ability of a grain boundary to plastically deform by emitting dislocations, during shock compression. Plastic response of a GB, affects the development of stress concentrations believed to be responsible for void nucleation by acting as a dissipation mechanism for the applied stress.

The microstructure of metals and alloys consists of a network of single crystals (called grains) and the boundaries between these crystals (called grain boundaries). While microstructure represents the overall uniform description of a material, it is the imperfections or defects like grain boundaries, the weak links, which dictate some material properties, especially its propensity to fail. In the case of ductile metals, like Cu, Au, or Al, it is known that all grain boundaries are not equal in their propensity for either strengthening or weakening a material. However, the reason for this selectiveness is not well understood. To develop next generation predictive capabilities for dynamic failure in engineering materials and to design future damage-tolerant materials, it is imperative to not just qualitatively understand the reasons behind why all grain boundaries are not equal but also to quantify this selectiveness in a manner that can be useful for both design and prediction. In this work, we focus on developing a qualitative and quantitative understanding of failure at grain boundaries in ductile materials.

Fracture in ductile materials is typically characterized by void nucleation, growth and coalescence. A broad-base fundamental problem in materials science is where, when, and why voids nucleate and grow in materials as an initiation of damage and failure. Previous experimental research has shown that microstructural features such as grain boundaries, inclusions, vacancies and heterogeneities can act as initial void nucleation sites123456. However, we cannot accurately predict the preferred locations and stresses for void nucleation. In fact, larger scale simulations that study void growth and failure tend to begin with material microstructures where voids have been nucleated in a stochastic manner, which is not a true representation of void nucleation in real experiments7. This random nucleation of voids is tied to a lack of systematic and detailed understanding of the critical mechanisms and stresses involved in void nucleation. Direct observations of void nucleation and early stage growth are critical to understanding these underlying mechanisms. However, these in-situ observations have long remained a grand experimental challenge, largely due to the extremely fine spatial (nm) and narrow temporal (fs-ns) scales involved. Hence in this work, we utilize molecular-dynamics simulations, which can provide in-situ, time resolved insights, to quantify for the first time the relationship between void nucleation stress and grain boundary structure.

A specific type of fracture in solids, induced by shock compression-and-release termed “spall” motivates this investigation28. During spall, compression waves are generated in both the target and flyer. These compression waves are reflected from the target and flyer free surfaces and intersect each other to form a region of tension, which can lead to failure9. Although, there are differences in stress states, rates of work hardening and operative deformation mechanisms preceding damage evolution between various loading conditions, we argue that the observations from the current work on dynamic loading provide previously missing insights into the rationale for preferred void nucleation sites.

Specifically, this MD study is motivated by a plethora of experimental work with similar observations regarding void nucleation under both dynamic and quasti-static loading10111213141516. Dynamic damage in high-purity Cu demonstrates, that voids preferentially nucleate at grain boundaries10, and specifically certain special boundary types are consistently more resistant to void nucleation11. Furthermore, it has been shown that boundary structure and not just the crystallographic orientation on either side of the boundary is important with regard to the early stages of damage17. This observation is also supported by the work of Wayne et al.12 who found that grain boundaries with certain misorientations, in polycrystalline Cu, serve as preferred locations for intergranular damage. These insights regarding damage nucleation at high strain-rates apply equally to damage at lower strain-rates. Work by Mikhailovskij13 on bcc tungsten shows that under uniaxial stress conditions, special boundaries such as Σ1, Σ3, and Σ9 possess higher resistance to failure in comparison with other types of CSL and non-CSL boundaries. Experiments by Lim14 on low-cycle fatigue of polycrystalline fcc nickel samples showed that low-order CSL boundaries such as Σ3 and Σ5 did not fail during the deformation process. Evrard15 and McMurtrey16 made similar observations through both experiment and simulations showing that special boundaries were more resistant to crack nucleation in pre-irradiated austentic stainless steels. These results collectively suggest that all grain boundaries, regardless of materials or loading condition, are not equal in terms of their propensity for serving as void nucleation sites.

Similar numerical results have been ascertained by molecular-dynamics (MD) simulations. Recent MD simulations by Luo et al.1819 on a Σ3 asymmetric tilt grain boundary support the experimental findings listed above by showing that crystal anisotropy and grain boundaries both change the response of a material under shock compression, in comparison to the response a single crystal, by reducing both its spall strength and flow stress. Work by Han et al.20 on Σ3 coherent and incoherent twin GBs in copper shows that, under dynamic loading conditions, void nucleation occurs at the coherent twin boundary and not at the incoherent boundary emphasizing yet again that not all GBs are equal in terms of void nucleation.

In previous MD simulation work, we have attempted to understand and predict void nucleation at grain boundaries by using average GB properties212223. To help frame the relationship between these properties and void nucleation21, we considered a standard model for the fracture toughness of a material2425: where γ and γ are the fracture energy, the surface energy, the grain boundary energy, and the energy associated with plastic work, respectively. We have examined the importance of γ, γ, and related average properties, such as excess volume to predict the stress required to void nucleation at a grain boundary in a ductile metal. Those results suggest that there is no direct correlation between grain boundary energy21, excess volume and the void nucleation stress for a ductile material. In a brittle material, where γ = 0, the susceptibility of a boundary to fracture can be predicted using the grain boundary and the free surface energies. However, the plastic energy term, γ can be dominant in ductile materials101112. In this work, we determine if the plastic work term, γ is a better determinant of the void nucleation stress for an interface in ductile bi-crystals.

This paper presents a correlation of the susceptibility of a boundary to nucleate a void under tensile unloading to its ability to undergo plastic deformation under shock compression, using MD simulations of copper bi-crystals. Specifically, we simulate the impact of a dissipative process such as plastic deformation by dislocation emission on the development of the stress concentrations believed to drive void nucleation. The effort outlined above is the first step towards examining and predicting void nucleation at grain boundaries under high strain-rate, uniaxial strain loading. This paper presents the simulation methodologies, the substance of the simulation results, and finally conclusions concerning the effect of dissipative plastic deformation on the stress to initiate void nucleation at a boundary.

Results

GB structures are 0K

Four boundaries in Cu were chosen for this study - representing a range of GBs encountered in “real” polycrystalline materials:

Σ11 <101> {545}-181} 50.48° asymmetric tilt

Σ11 <101> {545}-{181} 50.48° asymmetric tilt boundary whose structure was disordered by annealing and subsequently quenching the Σ11 asymmetric tilt boundary

Σ5 <100> {310} 36.9° symmetric tilt

Σ3 <110> {112} 36.9° symmetric tilt

The relaxed structures for these boundaries at 0 K are shown in Fig. 1 and compare well with other GB structures found in other studies262728. The atoms in Fig. 1 are colored by the centrosymmetry parameter29 where blue sphere denotes an atom in an environment corresponding to the fcc structure and the red sphere indicates an environment that is furthest away the from fcc structure. The Σ3 and both the Σ11 boundaries have pre-existing stacking faults protruding from the grain boundary into one of the grains. These stacking faults are connected with Shockley partial dislocations that are generated due to the dissociation of grain boundary dislocations that reduce the stress concentrations at the boundary grain boundary. The average calculated properties for these boundaries are available in Ref. 30.

Figure 1

GB structures at 0 K calculated using an EAM potential for Cu33 for the boundaries utilized in this work.

These structures are viewed along tilt axis. The atoms are colored by the centrosymmetry29 parameter where blue represents atoms close to an FCC environment and green represents atoms no longer in FCC configurations. The range of the centrosymmetry parameter in the Σ3 and Σ5 boundaries is 0–3 and 0–13 in the Σ11 boundaries.

GB Plasticity under Mechanical Loading

Under shock compression, dislocation emission is observed from grain boundaries in the perfect bi-crystals. The simulations cells are designed such that other complexities like the presence of pre-existing dislocations that can lead to pile-ups, triple points and point defects are avoided. Hence, in this work, as a starting point, we investigate one particular dissipative process, plastic deformation by dislocation emission at a grain boundary.

The extent of plastic deformation by the emission of Shockley partials varies in each boundary. The Σ3 boundary emits regularly spaced Shockley partials in both the boundaries as shown in Fig. 2. This difference in plastic response is attributed to grain boundary structures. This is best highlighted by the Σ11 asymmetric tilt boundaries, in Fig. 2, that have different local boundary structures, but the same average grain orientations across the boundary, and very different emission responses.

Figure 2

Snapshots showing emission of Shockley partials from the boundaries studied in this work.

The atoms are coloured by the centrosymmetry29 parameter where blue represents atoms close to an FCC environment and green represents atoms no longer in FCC configurations.

We quantify the difference in the plastic response of these four boundaries by using the total excess energy (as defined in the Methods Section). Figure 3 shows the calculated total excess energy (Eq. 2) as a function of time. Before the shock compression wave reaches the GB, this excess energy is equal to the grain boundary free energy. The arrival times of the shock compression wave at the GB vary by 1–2 picoseconds at each boundary because of slight differences in the sound speed, in each bi-crystal, due to the varying elastic constants (as shown in the inset of Fig. 3).

Figure 3

Total excess energy increases as a function of time after interacting with the shock wave due to plastic deformation.

The red line with squares corresponds to the Σ3 symmetric tilt, green with triangles to the Σ11 asymmetric ordered tilt, blue with circles to Σ5 symmetric tilt and black with diamonds to Σ11 asymmetric disordered tilt boundaries, respectively. The inset shows the slightly different arrival times for the shock compression wave at each boundary.

As the shock wave reaches the boundary (around 7 ps), the excess energy increases sharply. The magnitude of this increase is dependent on the amount of plastic deformation at a grain boundary under shock compression. The plastic work continues to increase steadily after the compression wave is reflected from the free surfaces. Even though the reflected wave annihilates some of the dislocations emitted under the compression wave, it also generates more plasticity by emitting a different set of dislocations from the boundary.

Hence, the peak in the excess energy does not necessarily correspond to the reflection of the compression wave from the free surfaces. Eventually, the annihilation of plasticity is higher than the generation of plasticity and results in a drop in the excess energy. The increase in the excess energy at later times (after 25 ps) is associated with the nucleation of voids at the grain boundary. However, the excess energy at the boundary, immediately prior to void nucleation is more relevant to understanding and predicting its susceptibility to void nucleation.

We extract the times associated with void nucleation to be approximately 19.6, 20.2, 21.1 and 19.9 ps for the Σ11 ordered, Σ11 disordered, Σ3 and Σ5 boundaries, respectively, from our MD simulations in which the cells shown in Fig. 1 were loaded to failure. The excess energy is then calculated 2–3 ps before voids nucleate at the boundaries and is plotted along with the void nucleation stress in Fig. 4. This figure shows that there is a direct correlation between the total excess energy and the void nucleation stress. It is important to note that this is the first time a quantitative relationship has been developed and demonstrated between plasticity and the void nucleation stress.

Figure 4

The stress required for void nucleation increases linearly with the calculated excess energy.

The points in this plot correspond to the four boundaries examined in this study.

Figure 4 shows that grain boundaries that undergo higher plastic deformation under shock compression require higher stresses for void nucleation. For example, there is a 1.34 GPa difference in the void nucleation stress for the two-Σ11 asymmetric tilt boundaries. This difference can be attributed to the fact that the disordered Σ11 boundary can dissipate higher amounts of energy associated with the applied stress, by emitting more dislocations as compared to the ordered Σ11 boundary (refer Fig. 4) even though the only difference between the two is boundary structure; orientation relationship across these two is the same. These results suggest that plastic deformation by emission of dislocations from a grain boundary, controlled, at least in part, by boundary structure, acts as a principal dissipation mechanism for the applied stress and can dictate if the stress concentrations required for void nucleation at a boundary are achieved.

Discussion

In this work, we demonstrate by atomistic simulations that there is a direct correlation between the stresses required for void nucleation and the extent of plasticity at GBs. For the particular boundaries studied here, we found that the Σ3 symmetric tilt boundary was the most resistant to void nucleation due to its ability to efficiently dissipate stresses associated with mechanical loading via maximum dislocation emission.

These atomistic results are consistent with experimental findings on dynamically loaded Cu that show that Σ3 twin boundaries are more resistant to void nucleation31. Indirect experimental findings show that grain boundaries that are surrounded by “softer” grains, as defined by the Taylor factor, tend to not nucleate voids as easily as boundaries that are surrounded by “harder” grains32. Softer orientation grains promote plasticity whereas harder grains retard it. Since Taylor factor is a measure of the ease with which plastic deformation can occur, these experimental results amongst many others, indirectly verify the findings from the molecular dynamics simulations.

While these results suggest plastic deformation as a principal dissipation mechanism for the four boundaries studied here, this correlation has not been rigorously tested for other boundaries and materials. To fully quantify the operative deformation and eventually damage evolution and failure mechanisms, it will be necessary to account for potential other competing processes such as existing dislocations, in conjunction with the plastic dissipation phenomena explored here. Nonetheless the results of this study for the first time quantitatively and qualitatively highlight the importance of plasticity via dislocation emission as a dissipation mechanism during void nucleation. These findings can be utilized to design a set of criteria that can be used to predict void nucleation in a deterministic manner and improve our damage prediction capabilities. Finally, this is an important step towards modelling the complex damage processes and the competition between diverse stress accumulation and dissipation processes in polycrystalline metals.

Methods

All of the MD simulations are based on an embedded-atom method (EAM) interatomic potential model for Cu developed by Mishin et al33. This model for Cu has proven capable of properly representing a wide spectrum of material properties under these moderate conditions, specifically in reproducing the experimental shock Hugoniot34. MD simulations applied a combination of SOLVER35 for initially relaxing the grain boundaries at 0 K, and LAMMPS36 for shock-loading simulations.

GB structures at 0K

In order to obtain a relaxed GB structure at 0 K, we create an initial simulation block with two grains separated by a grain boundary. The block is periodic in the plane of the boundary (x,z directions) and non-periodic perpendicular (y-direction) to it. The grains are sandwiched between two slabs parallel to the grain-boundary plane. The atoms in each of these slabs are fixed relative to each other and can only undergo rigid body motion. The two grains are free to move and undergo displacement in all the three directions during relaxation using a viscous drag algorithm. This relaxation technique is equivalent to minimization with respect to the specific volume. A standard method of molecular dynamics that generated atomic trajectories as a function of time by integration of Newtonian equations of motion is used to relax its structures353738.

Mechanical Loading of GBs

For the shock simulations, the bi-crystal cell is divided into two regions normal to the boundary plane: flyer and target as shown in Fig. 5. The target region is twice as long as the flyer plate such that the spall plane is located at the grain boundary. The shock particle velocity is denoted by up and in this paper it was set to 0.5 km s−1. Each atom in the flyer plate is initialized with a velocity of 2up, in addition to its thermal velocity, in the velocity component along the shock direction, causing it to impact the target. The shock direction is perpendicular to the boundary plane. Shock simulations use a NVE ensemble and are performed using LAMMPS. Periodic boundary conditions are applied along the x and z-axis and the nonimpact sides of the flyer plate and target are free surfaces.

Figure 5

The schematic shows the simulation cell used during dynamic loading.

The red lines represent the grain boundary, the blue lines show examples of slip planes along which dislocations can be emitted from the boundary. γp and γs are plastic dissipation and surface free energy. The free surfaces are perpendicular to the GB (not marked).

For the bi-crystals the GB is placed approximately in the middle of the target region and the shock waves travel from grain I to grain II, respectively, before being reflected from the free surfaces. In our shock loading MD simulations, we use a time step of 1 fs and total simulation times of 100 ps. Simulations are initialized at 300 K.

Stress Required for Void Nucleation

To track the shock front and other physical properties during the simulations, the simulation cells were divided into bins of length 0.55 nm - equal to the cut-off distance for the EAM potential - along the shock direction. Properties such as particle velocity and stresses were then averaged within each bin to give the longitudinal shock profile. The atoms are coloured by the centrosymmetry parameter to characterize local deformation and structure in the simulation cell since this parameter can identify dislocation cores, stacking faults and Shockley partials29. In this work, void nucleation is defined as formation of a void that is coupled (no size or location constraints were applied) with the stress relaxation. Nucleation is determined by examining both the quantitative stress data (Fig. 6) and qualitative images (Fig. 7).

Figure 6

Averaged stress as a function of bin position at time, t.

The red rectangle represents the areas of tension.

Figure 7

Snapshots showing the Σ3 GB before and after void nucleation.

The atoms are coloured by the centrosymmetry29 parameter where blue represents atoms close to an FCC environment and green represents atoms no longer in FCC configurations.

The void nucleation stress for the boundaries is calculated as the peak stress, from the averages stress profiles, immediately before void nucleation, in the shock direction. This is a standard method for calculating stresses associated with void nucleation in shock loading studies3940. Examples of an average stress profile and qualitative images showing the GB before and after void nucleation are shown in Figs. 6 and 7.

Excess Energy

In order to quantify the amount of plastic work, excess energy was used as a first order approximation. Excess energy was calculated in each bin using the following equation: where n is the number of bins, is the total energy of atoms in bin i, is the energy of atoms in bin i if they were in a perfect FCC configuration, and AGB is area of the GB plane. Excess energy would be approximately zero in bins far away from the boundary and have a finite value in the bins that include the GB. The total excess energy would be equal to GB energy, prior to deformation. To exclude the affects from the thermal and mechanical noise, we only use a limited area around the grain boundary to calculate excess energy. This limited region is selected based on the centrosymmetry parameter and represents the region with the largest increase in the centrosymmetry parameter as compared to a reference. The reference centrosymmetry parameter is chosen as the average parameter value in bins far away from the boundary at time 0 ps.

It is important to note that the excess energy contains both the elastic and plastic energies. However, a simplifying approximation assuming the elastic component of the excess energy to be similar for all the boundaries is made in this study. This is a valid approximation since the energy per atom in the compressed region of the simulation cell, where the material is only undergoing elastic deformation, varies by ~0.01 eV amongst the boundaries.

Author Contributions

E.K.C. conceived the project and provided guidance. S.J.F. performed the theoretical computations, analysis and wrote the manuscript. S.M.V. and G.T.G. reviewed the manuscript and offered helpful suggestions during the work.