Strain-Rate Dependence of Plasticity and Phase Transition in [001]-Oriented Single-Crystal Iron

Abstract

Non-equilibrium molecular dynamics simulations have been used to investigate strain-rate dependence of plasticity and phase transition in [001]-oriented single-crystal iron under ramp compression. Here, plasticity is governed by deformation twinning, in which kinetics is tightly correlated with the loading rate. Over the investigated range of strain rates, a hardening-like effect is found to shift the onset of the structural bcc-to-hcp phase transformation to a high, almost constant stress during the ramp compression regime. However, when the ramp evolves into a shock wave, the bcc–hcp transition is triggered whenever the strain rate associated with the plastic deformation reaches some critical value, which depends on the loading rate, leading to a constitutive functional dependence of the transition onset stress on the plastic deformation rate, which is in overall consistence with the experimental data under laser compression.

1. Introduction

Over the past few decades, there has been an increasing interest in using dynamic compression at extremely high strain rates to investigate both plasticity and pressure-induced phase transition [1,2,3]. Usually, such studies are performed under the conditions of uniaxial strain, which involve large deviatoric stresses. High strain rates affect material behavior in both the mechanisms of plastic deformation and the kinetics of polymorphic phase transitions [4,5]. They can even affect melting temperatures at high pressures [6]. Shock-loading experiments were classically realized under planar plate impacts where the strain rates are typically about 10⁴ to 10⁶ s⁻¹ [7,8,9]. Then, iron was found to plastically yield at a stress of 0.92 to 1.3 GPa, while the ground-state body-centered cubic (bcc) structure ($\alpha$ phase) was found to transform into the high-pressure hexagonal close-packed (hcp) structure ($\epsilon$ phase) at a stress ${\sigma }_T$ of about 12.88 to 14.26 GPa, which is higher for single-crystal than polycrystalline iron [7,9,10]. Higher strain rates can be achieved under laser ramp compression. Indeed, strain rates from 3 to 9 × 10⁷ s⁻¹ were reported by Amadou et al. [11] under nanosecond laser compression, where the phase transformation onset stress, ${\sigma }_T$, was found to vary from 11 to 25 GPa with a constant completion time of 1 ns, suggesting the existence of an isokinetic regime over the explored range of strain rates. Furthermore, an extended analysis of strain-rate effects, up to 10⁸ s⁻¹, on both plastic flow and phase transition kinetics in iron was reported by Smith et al. [12,13], where the elastic limit (${\sigma }_E$) in polycrystalline iron was found to vary from $~$1 to $~$5.5 GPa, while ${\sigma }_T$ ranged from $~$14 to $~$40 GPa. Moreover, a constitutive functional dependence on the plastic deformation strain rate, ${\dot{\epsilon }}_P$, was evidenced for both ${\sigma }_E$ and ${\sigma }_T$. At strain rates higher than about 5 × 10⁶ s⁻¹, a sharp increase in ${\sigma }_E$ with increasing ${\dot{\epsilon }}_P$ was interpreted as a transition in the plastic flow regime from thermally activated to phonon drag dislocation dynamics, which affects the structural phase transition kinetics by limiting the new phase growth through energy dissipation [13]. In this regime, ${\sigma }_T$ was found to scale as ${\dot{\epsilon }}_P^{0.18}$. A strain rate on the order of about 10⁸ s⁻¹ was also reported by Hawreliak et al. [14]. However, ${\sigma }_T$ was found to be much lower, at about 15 GPa, with no obvious correlation with the strain rate. Even higher strain rates up to about 10⁹ s⁻¹ are currently available under picosecond laser compression [15,16]. In this ultra-fast deformation regime, Crowhurst et al. [17] and Hwang et al. [18] found ${\sigma }_T$ values of 25 and 34 GPa, respectively.

These experimental ultra-high strain rates are comparable to those currently available in molecular dynamic simulations. Thus, Gunkelmann et al. [19], using the modified version of the Ackland EAM potential [20] to study the behavior of polycrystalline iron at a strain rate of about 10⁹ s⁻¹, showed that yielding occurs through dislocation activities around the grain boundary at 10 GPa with the onset stress of a structural bcc-to-hcp phase transformation ${\sigma }_T$ of about 23 GPa. A strain-rate regime ranging from 10⁹ s⁻¹ to 10¹¹ s⁻¹ was explored by Wang et al. [21] using a modified analytic embedded atom model [22], and they found that ${\sigma }_T$ varied from 25 to 38 GPa with an exponential dependence in ${\dot{\epsilon }}_P$,  ${\sigma }_T\propto {\dot{\epsilon }}_P^{0.196}$. Finally, using Voter–Chen potential [23], singularized by a lack of plasticity before the structural phase transformation [24], Shao et al. reported a linear dependence of ${\sigma }_T$ on ${\dot{\epsilon }}_P$ [25] over a strain rate range from 10¹⁰ s⁻¹ to 10¹¹ s⁻¹ with ${\sigma }_T$ ranging from 15 to 25 GPa. Thus, despite extensive research, the dynamic response of iron, including the bcc–hcp phase transition, its kinetics, and its strain-rate dependence under high strain rates still remains an open issue. In this context, using the modified version of the Ackland iron potential [20], we have reported in previous papers that defect-free single-crystal iron at a 50 K initial temperature subjected to ramp compression along the [001] direction exhibits a hardening-like effect, which has been shown to inhibit the nucleation of the hcp phase so that the onset of the phase transformation is shifted to very high pressures (on the order of 100 GPa) [5,26,27,28]. This remarkable effect, never reported before, deserves extensive study because it may contribute to the strong variations in the transition onset pressure with the loading rate classically observed experimentally. Here, we go further to explore the strain-rate dependence of both plasticity and phase transition and to investigate the influence of the hardening-like effect on the ${\sigma }_T-{\dot{\epsilon }}_P$ relationship. Because we reported that this hardening-like effect is more important when the ramp wave propagates along the [001] direction than along other low-index crystallographic directions [26], we choose to simulate the response of a [001]-oriented single-crystal iron.

2. Method and Computational Details

Samples with up to 28 million atoms, 100a₀ × 100a₀ cross-section, and a length up 1400a₀ (with a₀ = 2.87 Angstrom, the lattice constant) were simulated using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics code [29]. The interactions between atoms were modeled through the embedded atom model (EAM) formalism [30,31], where the iron model was a modified version of the Ackland potential [20], which was successfully used to address both plasticity and phase transition in iron [5,19,26]. The sample was thermalized at a 50 K initial temperature, then dynamic ramp-wave compression was realized by driving an effective infinite-mass wall piston with an imposed velocity v(t) along the z-axis, oriented along the [001] crystallographic direction, while periodic boundary conditions were used for the transverse directions. In order to investigate strain-rate effects, the piston velocity was increased linearly from 0 to 1600 m·s⁻¹ with a rise time varying from 15 to up to 150 ps. The corresponding loading rates were about 1–10 × 10⁹ s⁻¹, which is comparable with that reported under picosecond laser dynamic compression [15,16,17,18]. The local thermodynamic and mechanical properties, such as longitudinal stress ${\sigma }_z$, shear stress, temperature, etc., were evaluated within a spatial planar bin (of 3-lattice constant width) perpendicular to the wave propagation direction in the same manner as our previous work [5,26,27]. Finally, local structural analysis was performed by adaptative common neighbor analysis (CNA), centrosymmetry, and DXA, as implemented in OVITO software (https://www.ovito.org, accessed on 18 January 2023) [32].

3. Results and Discussion

Figure 1 shows snapshots of atom configurations for simulations where the ramp rise time was 30, 45, and 60 ps. At each loading rate, the compressed material exhibits four distinct regimes: an elastic compression of bcc-iron (yellow zone); a regime where defects associated with plastic deformation can be observed (green zone with brown micro-features identified as twins); a regime where these defects associated with plastic deformation progressively disappear (blue zone); and a regime where iron is transformed from bcc to hcp phase (violet zone). More details on wave propagation and splitting of the compression fronts can be found in ref. [5]. At each loading rate, the sample is found to yield at a stress ${\sigma }_z$ of about 12 GPa via the generation of micro-twins with no sizeable influence on the loading rate. The fact that twinning deformation is independent of the loading rate is consistent with the well-known existence of the twinning threshold in bcc materials, such as tantalum, reported under both experiments [1,33] and molecular dynamics simulations [34]. Thus, regardless of the loading rate, twinning is found to be the main deformation mechanism for the defect-free single-crystal iron ramp compressed along the [001] crystallographic direction.

Figure 2 shows the evolution of the twin fraction during the compression for a ramp rise time of 15, 30, 45, and 60 ps. After the nucleation period, the twin fraction grows rapidly, reaching its maximum within 3–4 ps. Then, the twin fraction starts to decrease, i.e., twins formed under moderate pressure are removed upon further compression, which we refer to as a hardening-like effect (see Figure 1, blue zone). The longer the ramp rise time, the more important the twins’ fraction and the duration of the subsequent receding phase. Thus, upon increasing the ramp rise time, the maximum twins’ fraction increases from 28% to a value as high as $\approx$67%, in consistence with the 50% twins’ fraction usually assumed under dynamic loading [1,33]. Furthermore, this increase in the twins’ fraction upon increasing the ramp duration is consistent with both experimental observations [35] and thermo-mechanical simulations, including a model of twins’ nucleation and growth [36]. In the region affected by the hardening-like effect (blue zone), where twins are almost fully removed and no dislocation activities can be detected, the shear stress increases with compression, as can be seen in Figure 3. Such elastic stiffening of the bcc matrix leads to a confinement effect inhibiting the nucleation of the hcp phase. This confinement effect has been interpreted in the context of classical nucleation and growth theory, and it has been shown that in order to activate new phase nucleation two different scenarios are possible [5]: (i) ramp compressing must be kept at a much higher pressure by increasing the piston maximum velocity or (ii) the ramp wave, which steepens during its propagation due to the increase in sound velocity with pressure in the bcc phase, should be allowed to evolve into a shock wave. Thus, for the simulations shown in Figure 1, both piston maximum velocity and propagation distance are enough that two nucleation fronts can be observed. The first front, on the left part of the sample, corresponds to the top of the ramp, where pressure is sufficiently high for the hcp nuclei to grow and reach a critical size, despite a rigid, confining bcc matrix due to the hardening-like effect. The second nucleation front, in the middle of the sample, corresponds to the second scenario above (ii), where the so-called P1 wave evolves into a shock wave during its propagation.

Thus, the resulting instantaneous thermodynamic driving force is high enough for hcp embryos to grow and reach a critical size almost immediately behind the shock front [37]. If the ramp rise time is increased to a sufficiently high value, the P1 wave does not form a shock upon the simulation duration. Thus, the second nucleation front can no longer be observed for a ramp rise time above 150 ps.

Figure 4 shows the structural phase transition onset stress, ${\sigma }_T$, as a function of the ramp rise time for both scenarios (i) and (ii) mentioned above. In the first scenario, upon increasing the ramp rise time, i.e., decreasing the loading strain rate, the phase transition onset stress can be observed to remain almost constant, at about 100 GPa (red triangles), so that no loading rate dependence can be observed. This observation can be interpreted as follows: As the hardening-like effect confines the hcp embryos, only the material strength determines the maximum pressure supported by the bcc matrix. Yet, it has been shown that above a critical strain rate of 10⁵ s⁻¹, the material strength becomes strain-rate independent [38,39]. This result highlights the crucial role of the strength of the parent phase matrix in the dynamics of structural phase transformation due to the generation of the elastic strains necessary to accommodate the difference between the parent and daughter crystalline structures.

Although the transition onset stress in scenario (i) remains almost constant after the hardening-like effect, the pressure relaxation accompanying the structural phase transformation is found to be rate-dependent (Figure 5). Indeed, this relaxation is due to the coexistence of both bcc and hcp phases in the sample, which causes a drop in the sound velocity that evolves as a result of the balance between the hcp nucleation and loading rates [13]. Moreover, increasing the ramp compression duration seems to enhance both the size and the ovaloid shape of the hcp nuclei (see Figure 6).

Ramp compression can be seen as a succession of infinitesimal compression waves, which all propagate at the speed of sound in the upstream medium. Because this speed usually increases with pressure (above the elastic limit), these elementary waves catch up with each other so that the ramp becomes steeper with increasing propagation distance and evolves naturally into shock waves. Its rise time and the associated plastic deformation strain rate, ${\dot{\epsilon }}_P$, is thought to be associated with the microscopic processes of dislocation multiplication and motion, twinning, vacancy production, precipitate alteration, etc. [40,41]. Figure 5 clearly shows such gradual steepening of the compressive wave beyond the elastic limit (Pz > 12 GPa). This wave, usually referred to as the P1 wave, is associated to the plastic compression of the bcc phase, up to the onset of the transition to hcp. The average plastic deformation strain rate, ${\dot{\epsilon }}_P$, associated with the P1 wave is calculated as follows [41,42]:

$${\dot{\epsilon }}_P=\frac{1}{∆z_{P_1}}\underset{∆z_{P_1}}{\overset{ }{{\displaystyle \int }}}\dot{\epsilon }\left(z\right)dz$$

(1)

where $∆z_{P_1}$ is the width of the P₁ wave between the elastic front and the bcc–hcp phase transition front in the scenario (i), as can be seen in Figure 7. The evolution of ${\dot{\epsilon }}_P$ as a function of time for various ramp rise times is presented in Figure 8. It ranges from a few times 10⁹ s⁻¹ to about 10¹¹ s⁻¹, in consistence with the values reported under the picosecond laser compression experiment [16,17]. As expected, when the ramp steepens, ${\dot{\epsilon }}_P$ increases as a function of time more rapidly for lower ramp rise time than for higher ramp rise time [43]. For each rise time, ${\dot{\epsilon }}_P$ shows a vertical asymptotic upward variation corresponding to the shock formation during the wave propagation through the sample. The lower the loading ramp rise time, the earlier the shock formation. Due to the complexity of the plastic response involved here, i.e., twins’ formation, hardening-like effect, and phase transition, including hcp phase nucleation and growth, the time taken for steepening is much lower than the ratio of ~1.4 to 1.6 times the rise time of the ramp found from the analytical models and MD simulations of simple materials such as Al and Cu [43,44].

The bcc-to-hcp phase transition is triggered in scenario (ii) whenever the plastic strain rate, ${\dot{\epsilon }}_P$, reaches a certain critical value that depends on the loading conditions (see Figure 8). This critical ${\dot{\epsilon }}_P^c$ ranges from 6.71 × 10¹⁰ s⁻¹ to 0.88 × 10¹⁰ s⁻¹ when the ramp loading rise time increases from 15 to 60 ps. Indeed, when these critical values are reached, the dynamic loading switches from ramp to shock compression, where the material is carried almost instantaneously from the thermodynamic stability conditions of the bcc phase to those of the hcp phase. The resulting driving force is high enough for hcp embryos to reach a stable critical volume almost instantly behind the shock front. Thus, the hcp phase starts to nucleate behind the elastic–plastic transformation front leading to a second phase transition front in the sample, referenced above as scenario (ii), as can be seen in Figure 1. In this scenario, the transition onset stress, ${\sigma }_T$, can be observed to decrease from 60 to 42 GPa upon increasing the ramp rise time from 15 to 60 ps. Thus, in contrast with scenario (i), ${\sigma }_T$ is rate-dependent. It increases with the strain rate, which is due to the kinetics of the transformation, and is consistent with experimental observations [11,13,45].

Figure 9 shows the variation in ${\sigma }_T$, as a function of ${\dot{\epsilon }}_P^c$ (red stars), in comparison with various data reported in the literature under both experiments and MD simulations. The black line is the prediction of the Swegle–Grady (SG) law ${\sigma }_T=\alpha {\dot{\epsilon }}_P^{\beta }$ [40], using the parameters reported in ref. [14]. The red curve presents the fit of our data using a similar law where $\alpha =44.88\pm 0.47$ and $\beta =0.15\pm 0.01$. MD simulations are from Wang et al. [21], using a MEAM potential where the plasticity is dominated by twins; Shao et al. [25], using the Voter–Chen EAM potential [23], which does not predict plastic deformation before the structural phase transformation; and Gunkelmann et al. [19], using the same potential as we did.

Experimental data come from Smith et al. [13], where plastic deformation was thought to be dominated by phonon scattering from defects, nanosecond laser ramp compression iron by Amadou et al. [11] and Hawreliak et al. [14], or picosecond laser compression by Crowhurst et al. [17] and Hwang et al. [18]. Although our data deviate from the predictions by both the original SG model and Smith et al.’s constitutive relationship [13] at an ultra-high strain rate, their extrapolation to lower strain rates using this updated SG fit is consistent with the experimental data, including those by Amadou et al. [11], Smith et al. [13], and Hwang et al. [18]. On the other hand, the magnitude of the power law exponent in our simulations is slightly different from the 0.196 value reported by Wang et al. [21] in single-crystal iron for the same crystallographic orientation under MD simulations using the MEAM interatomic potential. This discrepancy is thought to be due to the difference in the underlying plastic deformation micro-process predicted by the different potentials. Indeed, in our simulations, twinning is followed by a hardening-like effect under further compression, while such a post-twinning effect is not predicted by the MEAM potential used by Wang et al.

4. Conclusions

NEMD simulations were used to investigate the strain-rate dependence of plasticity and phase transition in the [001]-oriented single-crystal iron at 50 K under dynamic ramp compression between 10⁹ s⁻¹ and 10¹⁰ s⁻¹. Iron was found to yield at 12 GPa, regardless of the loading rate, through the generation of micro-twins, which grew rapidly to a peak fraction before receding upon further compression. The longer the ramp rise time, the higher the maximum twins’ fraction and the longer the subsequent receding regime, in consistence with both experimental and theoretical observations in iron. As reported previously, twin recession in the absence of any dislocation slip induced a hardening-like effect, which shifted the phase transition onset stress to a very high value of $\approx$100 GPa scenario (i) independently of the loading rate. On the other hand, the phase transition could be triggered at lower stresses when the ramp evolved into shock wave scenario (ii). Then, the onset stress ${\sigma }_T$ was evidenced to be strain-rate-dependent. Indeed, the transition was triggered whenever the strain rate associated with the plastic deformation reached some critical value, ${\dot{\epsilon }}_P^c$. The higher the ramp compression time, the lower both ${\dot{\epsilon }}_P^c$ and ${\sigma }_T$. Thus, the onset stress, ${\sigma }_T$, in this scenario has been shown to follow a Swegle–Grady power law type [40] in ${\dot{\epsilon }}_P^c$ with an exponent of 0.15, i.e., ${\sigma }_T\propto {\dot{\epsilon }}_P^{c 0.15}$, which is in overall consistence with some experimental data under laser compression.