Rate-limited plastic deformation in nanocrystalline Ni

Experiments and numerical simulations explain some of the deformation mechanisms that lead to the superior mechanical properties of nanocrystalline (nc) metals,^1,2 for example, the interplay between grain boundary (GB) deformation and dislocation activity,^3–6 and twinning deformation in nc and nanolayered materials.^7–10 These studies focus primarily on length scale effects, in particular, grain size and grain size distribution. However, nanometer scale simulations are limited to high strain rate regimes and it remains an open question whether the high strain rate results from numerical simulations can be extrapolated to experimental conditions.

The advance in supercomputing has enabled recent molecular dynamics (MD) simulations to investigate the strain rate sensitivity at lower strain rates^5,11,12 than before.^3,13 Brandl et al.¹¹ investigated the influence of the strain rate in the range of 10⁷ s⁻¹–10⁹ s⁻¹ in nc Al. A 14% increase in the yield stress was observed when the strain rate was raised by one order of magnitude. This was explained by the authors by identifying two different regimes during dislocation gliding: (a) athermal and purely stress driven regime without residence time and (b) thermal regime characterized by a residence time. By reducing the strain rate, more dislocations enter the thermally activated regime. In other words, at low strain rates dislocations have enough time to overcome the energy barrier by thermal activation. Recently, Kim and Strachan¹² conducted MD simulations on Pt nc bulk and thin film samples subjected to tensile loading. They also observed that when the strain rate was elevated from 10⁸ s⁻¹ to 10⁹ s⁻¹, the yield stress increased by 31% and the flow stress (at 10% strain) increased by 18%. Fan et al.¹⁴ proposed a model to study the effect of strain rate on the flow stress in single crystal Fe based on the transition state theory.¹⁵ They calculated the escape rate of a screw dislocation using dislocation glide barriers obtained from atomistic simulations under various applied stresses. This escape rate was then used to establish a relation among the flow stress, temperature, and strain rate, which explained the sharp stress upturn¹⁶ observed in experiments when the strain rate is beyond 10³ s⁻¹.

Dislocation dynamics (DD) can probe smaller strain rates than MD. Using discrete DD simulations, Tang and Marian¹⁷ observed an overshoot in the stress-strain curve of single crystal Fe deformed at strain rates ranging from 10⁴ s⁻¹ to 10⁶ s⁻¹. On average, when the strain rate was raised by one order of magnitude, the yield stress increased by 74% and the steady-state flow stress increased by 35%. Zbib et al.¹⁸ used DD simulation to study the evolution of dislocation structure under constant strain rate of 1 s⁻¹ in single crystal molybdenum. Even at smaller strain rates, DD simulations are used to predict the temperature and strain rate dependency of the strain rate exponent during Andrade's creep. To consider the strain rate dependency of the dissipation on the dislocation interactions, Sullivan et al.¹⁹ developed a model to study the behavior of thermally activated dissipative systems at different strain rates.

In this work, phase-field dislocation dynamics (PFDD) simulations are performed in nc Ni at strain rates in the range of 10⁶–10⁹ s⁻¹ at 0 K. The PFDD approach has the capability of investigating the collective behavior of dislocations with a wide spectrum of grain sizes, strain levels, and strain rates.^20,21 Moreover, the PFDD model can rigorously capture partial dislocations and stacking faults besides full dislocations by incorporating the material dependent gamma surface obtained from atomistic simulations.^22–24 The effect of strain rate, varying three orders of magnitude, on the deformation process is studied in a nc Ni sample with an average grain size of 15 nm. In particular, its impact on the stress-strain curve and the evolution of partials and full dislocations will be investigated.

This paper is organized as follows. In Sec. II, the formulation of the PFDD model and the geometry of the nc sample used in the simulations are described. The results of the simulations are reported in Sec. III, including dislocation structures, stress-strain curves, and density of stacking faults. Section IV contains the discussion and conclusions.

The plastic distortion $βijp(r)$ due to crystallographic slip can be expressed as²⁵ 

where d(x₁, x₂) is the displacement jump in the slip plane, x₁ and x₂ are the coordinates in the slip plane S_α, m^α is the slip plane normal, $δSα$ is a Dirac delta function confining the slip to the plane S_α, and the sum is over all slip systems α. The phase fields ξ^α that represent a slip density are defined as

where s^α is the direction of the Burgers vector and $|b|$ its magnitude, and l is the distance between slip planes. Therefore, the plastic distortion can be written in terms of the phase fields as

and the total distortion field can be obtained directly from the plastic part using the elastic Green's function, G(r),²⁶ as

where C_klmn are the components of the stiffness tensor, * represents the convolution operator, and the isotropic elastic Green's function is

The evolution of the phase fields, ξ^α(r), follows from the minimization of the energy through a system of time-dependent Ginzburg-Landau equations of the form

where N is the number of slip systems and L is a kinetic coefficient determined by the dislocation velocity.²² The energy contains three terms: the strain energy, Eⁱⁿ, the energy from an external driving force, E^ex, and the misfit energy, E^mf, associated with the gamma surface.²⁷ The strain energy is

where V is the simulation domain. By replacing Eqs. (3) and (4) in Eq. (7) and using Parseval's identity, the strain energy becomes²⁵ 

where ^∧ denotes the Fourier transform, * is the complex conjugate, k is the wave number vector in Fourier space, and

The external energy under strain controlled conditions is²⁴ 

where $ϵ¯ij$ are the components of the applied strain. The gamma surface is incorporated in the misfit energy by integrating the planar energy density functional, described in  Appendix A, over each slip plane

Replacing Eqs. (8), (10), and (11) in Eq. (6) and taking the variational derivatives (see  Appendix B) produces a system of coupled equations

As shown in Eq. (4), the elastic part is solved using the static Green's function, which enforces equilibrium of the elastic part of the deformation at each applied strain increment. The time-dependent Ginzburg-Landau equation, Eq. (12), is used to evolve the plastic strain. To calculate the strain rate, the distance traveled by one dislocation within one strain increment is divided by the dislocation velocity, which is assumed to be 1000 m/s.²² This is an approximation given that the velocity of the dislocation is related to the applied stress. By varying the strain increment in different simulations, the strain rates realized in the simulations range from 1 × 10⁶ s⁻¹ to 1 × 10⁹ s⁻¹.

The nc sample used in the simulations has a size of 75 × 75 × 40 nm³ and an average grain size of 15 nm as shown in Fig. 1. This polycrystalline sample is created from the atomic positions from MD simulations.¹² The coordination number of each atom is calculated and the non-12-coordinated atoms are considered as GB atoms, indicated in blue in Fig. 1. These GB atoms are mapped to the PFDD grid.

Random arrays of immobile and mobile dislocations are assigned to represent misorientations in the GB.^28,29 Upon loading the energetics of the system drives the emission of GB dislocations and their interaction with dislocations gliding in the grain interior. The PFDD simulations are conducted under strain-controlled loading condition in which the sample is subjected to an increasing applied strain ϵ₁₃ at a constant strain increment. Periodic boundary conditions are applied in all three directions. The Burgers vector is used as the grid size to resolve the dislocation core structure. The material properties of Ni used in this work are shown in Table I.

The phase fields, ξ^α(r), represent the slip in units of the Burgers vector in the $〈110〉$ directions. Accordingly, slip in any direction on the slip plane can be represented as³⁰ 

Unslipped regions have zero displacement; regions slipped by an extended full dislocation have a displacement of the Burgers vector of a perfect dislocation, $1/2〈110〉$⁠; and regions slipped by a partial dislocation have a displacement of $1/6〈211〉$⁠. The total slip, $Δ$⁠, is projected in the $[11¯2]$ direction to identify partial dislocations and stacking faults

The contour plots of Δ_p show partial dislocations and stacking faults. This visualization scheme is similar to the common neighbor analysis,³¹ which is commonly used to visualize stacking faults, dislocations, and grain boundaries in MD simulations. In addition, the contribution to strain by extended full dislocations can be quantified by projecting $Δ$ in the [110] direction

The contours of Δ_f show the areas slipped by partial and extended full dislocations.

The detailed evolution of individual dislocations observed in the PFDD simulations is illustrated in this section. Figure 2 shows the interaction between two extended full dislocations with opposite signs at a strain rate of 1 × 10⁹ s⁻¹. Each extended full dislocation is formed by a leading partial and a trailing partial (in green), separated by a stacking fault (in yellow). As the applied strain is increased, the extended full dislocations attract each other and annihilate.

In Fig. 3, a sequence of snapshots of PFDD simulations (at a strain rate of 1 × 10⁷ s⁻¹) is shown. A leading partial dislocation glides, recedes, and resumes gliding as the applied strain is increased. Changes in the direction of dislocation glide are caused by local variations of the stress distribution due to plastic activity on neighboring grains. It can be observed also that the dislocation is arrested by the GBs.

The impact of the strain rate on the stress-strain curves is analyzed in this section. For each simulation, a fixed strain increment in ϵ₁₃ (see Fig. 1 for axis directions) is applied to the domain. Fig. 4 shows the simulated stress-strain curves, with strain rates in the range of 1 × 10⁶ s⁻¹–1 × 10⁹ s⁻¹. The figure exhibits a significant growth of the stress with increasing strain rate and an overshoot in agreement with experiments³² and simulations.^{5,11,12,17,33} The stress overshoot is reduced with decreasing strain rate, as the dislocations have enough time to glide and reduce the local stress between strain increments.

Another prominent feature of Fig. 4 is the jerky character of the stress-strain curves at lower strain rates, 1 × 10⁶ s⁻¹ and 1 × 10⁷ s⁻¹. This behavior indicates microscopically heterogeneous deformation processes.^34–36 It is important to note that, in agreement with avalanche models, the jerky behavior is only observed at lower strain rates. At slow loading rates, dislocation activity is more spread in time and space, so each glide event is manifested as a distinct stress drop. At fast loading rates, in contrast, different dislocation glide events overlap in time. Therefore, the stress-strain curve shows the average of many of these single events, resulting in a smoother curve.

The grains in the PFDD simulations are oriented such that the stress σ₁₃ is the resolved shear stress in the plane $(1¯11)$⁠. Therefore, the critical stress calculated using the 0.2% strain offset shown as a dashed straight line in Fig. 4, is the critical resolved shear stress (CRSS). The predicted CRSS stress as a function of the strain rate is compared to experimental results in nc Ni with an average grain size of 40 nm (Ref. 32) in Fig. 5. The experimental CRSS is calculated with a Taylor's factor, m = 3, as CRSS = σ_f/m, from the uniaxial tensile stress, σ_f, reported by Schwaiger et al.³² Figure 5 shows a strong upturn in the CRSS in agreement with experimental observation³⁷ and constitutive models.^14,16,38

The significant differences among the stress-strain curves in Fig. 4 suggest that the activation of dislocations exhibit different features at different strain rates. The microscopic details of the PFDD simulations are used in this section to quantify these differences. Figures 6–9 show the area glided by leading partial dislocations (stacking faults) in yellow and extended full dislocations in gray at different strain rates.

At a strain rate of 1 × 10⁶ s⁻¹, Fig. 6 shows a glide event in a large region involving multiple grains. After being emitted from a GB, the first dislocation glides and activates dislocations at neighboring GBs. As the applied strain is increased, trailing partials start gliding forming extended full dislocations.

At higher strain rates, 1 × 10⁷ s⁻¹ to 1 × 10⁹ s⁻¹, the onset of dislocations gliding occurs at several GBs simultaneously, as shown in Figs. 7–9. The number of events increases with the strain rate. This can be explained as follows. Dislocation glide from a GB starts when the dislocation can overcome the energy barrier. At lower strain rates, the dislocation can glide over a large area and relax the local stress, reducing the number of events in other grains. This scenario is different at high strain rates in which the local stress rises before the dislocations can glide and, therefore, more leading partial dislocations start gliding from GB to accommodate the applied strain. This process is also responsible for the initial limited activity of trailing partials to form extended full dislocations at high strain rates.

To examine the relative contribution of partial and extended full dislocations to plastic strain at different strain rates, the area slipped by partial and extended full dislocations is measured. Figure 10 shows A_p, the stacking fault area, and A_f, the area swept by extended full dislocations, as a function of the applied strain at different strain rates. At strain rates 10⁸ s⁻¹ and 10⁹ s⁻¹, the activity of partial dislocations reaches a maximum and then decreases when trailing dislocations start to glide to form an extended full dislocations. In contrast, A_f increases monotonically with the applied strain in all the cases. The plastic strain is proportional to the sum A = A_f + 1/2A_p, where the 1/2 arises from the fact that for the same slipped area the strain accommodated by a partial dislocation is half of a full dislocation. Figure 10(c) displays A as a function of the applied strain. It can be noted that the slope of the curve A versus the applied strain is the same for all strain rates after plasticity is driven by gliding of full dislocations (after the stress drop).

Figure 10(c) shows a delay on the dislocation activity at higher strain rates. This delay is also observed by Brandl et al.¹¹ in MD simulations of nc Al when the stain rate was raised from 10⁷ s⁻¹ to 10⁹ s⁻¹ and it is explained by thermal activation. In the current simulations at 0 K, this effective delay in dislocation activity can be explained as follows. The onset of gliding of a dislocation is determined by the local resolved shear stress. At lower strain rates, partial dislocations are able to glide larger distances between load increments such that trailing partial dislocations can start gliding to form extended full dislocations during one strain increment step in the simulation. These events increase the plastic strain and, therefore, reduce the resulting stress. On the other hand, at larger strain rates the distance glided by leading partial dislocations is limited during each strain increment and no trailing partials are observed. The stress is increased rapidly in the next loading step, inducing the gliding of other leading partial dislocations that have a small contribution to the plastic strain. This leads to the effective delay in plastic strain observed in Fig. 10(c) and to the large stress observed in Fig. 4. It is important to note that in the simulations the dislocation velocity is assumed constant. Including a linearly growing dependency of the velocity on the resolved shear stress may change the predicted yield stress, but an increase of the yield stress with strain rate will still be observed.

To illustrate the relation between the delay in dislocation glide and the yield stress, Figure 11 shows the stress-strain curves together with the evolution of A as a function of the applied strain. The vertical dashed lines indicate the strain at which the slipped area reaches the value A = 75 × 75 nm², which is the area of the cross section of the simulation domain. This event coincides with the drop of the stress in the stress-strain curves for all strain rates. It is also interesting to note that there are multiple plateaus and jumps during the evolution of A at low strain rate (1 × 10⁶ s⁻¹). The plateaus in A indicate an elastic regime (no dislocation gliding), while the jumps denote dislocation gliding that correspond to a drop of the stress in the stress-strain curve.

PFDD simulations are performed in nc Ni with a grain size of 15 nm under strain rates from 10⁶ to 10⁹ s⁻¹ at 0 K. In the stress-strain curve, an increase in the yield stress is observed when the strain rate is increased, in agreement with experiments³² and simulations.^{5,11,12,17,33} While this high strain rate sensitivity is attributed to thermal activation in MD simulations, the same effect is observed in current simulations at 0 K. At lower strain rates, partial dislocations can glide larger distances and trailing dislocations become active to form extended full dislocations during one strain increment in the simulation. At larger strain rates, the area glided by partial dislocations is limited during each strain increment bounding its contribution to the total slip. This leads to an effective delay in plastic strain illustrated in Fig. 11, which is also related to the increase of the yield stress. The sharp stress upturn is universally observed in several materials³⁷ for strain rates beyond a critical range. A large number of constitutive models^14,16,38 have been proposed to explain this phenomenon. For example, the transition from thermal activation of dislocation glide at low strain rates to phonon drag at high strain rates,³⁸ dislocation nucleation at the shock front,¹⁶ and strain-rate-induced non-Arrhenius behavior.¹⁴ Nevertheless, the mechanisms for the sharp upturn are still under debate.

The dislocation microstructures that form at different strain rates also have striking differences. At low strain rate, the initial dislocation activity is limited to a few GBs and glide events with large areas slipped by extended full dislocations are observed. In contrast, dislocation activity is observed simultaneously at several GBs across the entire sample at high strain rates but the activity is limited to only leading partial dislocations with small areas slipped reducing the plastic strain.

In addition, in agreement with avalanche theories,^34–36 the simulations capture a jerky stress-strain response at low strain rates. This jerky behavior indicates a heterogeneous deformation process in which fast dislocation glide is separated by slow elastic processes.

The authors are grateful for the support from the U.S. Department of Energy, Office of Basic Energy Sciences under Contract No. DEFG-02-07ER47398.

The model takes into account the crystal disregistry by including the full gamma surface,^22,23,27,39 in the total energy. The planar energy density function, $ϕ(x)$⁠, describes the parameterization of the gamma surface for one slip plane, and it is written in terms of the 3 phase fields on each slip plane as^23,39,40

The material-dependent coefficients, c₀ – c₄, a₁, a₃, are determined by fitting the gamma surface obtained with atomistic simulations. The values of the coefficients for Ni used in the PFDD simulations are listed in Table II.^23,39

The system of time-dependent Ginzburg-Landau equations in Eq. (12) requires the calculation of the variation of the total energy with respect to all the phase fields. The steps to obtain Eq. (12) are summarized in this section. The chain rule is used to calculate the variational derivative of the strain energy

Equation (8) is used then, to calculate the variation of the strain energy with respect to the plastic distortion

A more compact form is obtained through using the symmetry of A_mnuv and change of variable

By using the definition of the plastic distortion in Eq. (1), the last term in Eq. (B1) can be expressed as

By replacing Eqs. (B4) and (B3) in Eq. (B1), the first term on the right hand side of Eq. (12) can be obtained

The second term on the right hand side of Eq. (12) is the variational derivative of the external energy with respect to the phase field

The surface integral is replaced by a volume integral to calculate the variational derivative of the misfit energy

and, therefore, the variational derivative is