Interplay of plasticity and phase transformation in shock wave propagation in nanocrystalline iron

Fe is a prototypical material showing phase transformation both induced by temperature (from the low-temperature bcc to the high-temperature fcc phase) and by pressure. In the latter case, Fe transforms from the low-pressure bcc phase (termed α) at a transition pressure of around 13 GPa to the high-pressure hcp (or ) phase [1, 2]. The physics underlying this transformation, its dynamical and kinetical aspects, as well as its connection to the magnetic nature of the Fe electronic system have been explored in depth [2, 3].

Shock waves in materials have long been known to feature the elastic → plastic splitting for sufficiently high shock strengths [4]; for much higher shock strengths, however, a single overdriven wave would result. In phase-transforming materials, in addition, the occurrence of the phase transformation leaves a novel fingerprint in the shock wave. The wave can split into three parts (so-called 'three-wave structure') because in addition to the elastic and plastic wave, a phase-transformed wave could also appear [5, 6]. Fe has developed into the prototypical material for studying these specific transformation shock waves [7–11].

Using atomistic computer simulation, shock waves in Fe were the subject of a series of detailed studies [12–14]. However, up until now the typical three-wave structure could not be found in the simulation [15]; more precisely, recent molecular-dynamics (MD) simulations [14] show a phase transition lacking previous dislocation activity in the Fe samples.

Recently, interest in the transient features of shocks and in shocks in thin films increased. Thus, Crowhurst et al [16] found experimentally that phase transformation begins within 100 ps after the shock started and needs a similar time span for completion. Deviatoric stress exceeds 3 GPa, while transition stresses above 25 GPa are needed to initiate the shock. Similar conclusions are found by Smith et al [11]; these authors stress that 'over-pressurization' of the phase boundary is a common feature in high-strain-rate compression.

In the present paper, we use MD simulations to study the propagation of shock waves through a nanocrystalline Fe sample with the aim of identifying the interplay of plasticity and phase transformation. Our work goes beyond previous modeling in several respects. We use a finite shock loading time and a large nanocrystalline sample; both features allow dislocations to nucleate more easily. In addition, we employ an interatomic interaction potential (the Ackland potential, see below), which provides a better description of both dislocation and phase transformation properties than previously employed potentials [17]. By comparing our results to those of a potential that does not feature the phase transformation (the Mendelev potential, see below), we characterize the decisive influence of phase transformation on the shock wave dynamics. In particular, here we study the influence of shock wave velocity on the transformation, and thus go beyond our previous results that were restricted to a fixed piston velocity of 0.7 km s⁻¹ [15].

Fe belongs to the elements most studied using atomistic simulations, and many interaction potentials are available to simulate various materials aspects. The pressure-induced bcc → hcp transition in Fe is closely connected to the magnetic nature of Fe [18]; empirical interatomic interaction potentials nevertheless are able to capture the essential mechanical properties of such phase transitions.

A recent review [19] analyzes the ability of several often-used Fe potentials to reproduce experimental and ab initio results of dislocation properties, including dislocation core structure, line energy, etc. It was shown that the Mendelev family of potentials [19–21] provides a reasonable description of dislocation properties. Here we use the potential described in [20]. We shall denote it as the Mendelev potential. Note, however, that this potential does not properly describe the pressure-induced bcc → hcp phase transition since its barrier is high in this potential and the transition pressure amounts to about 60 GPa.

In addition, we also employ here a new embedded-atom-method potential, similar to the Machová and Ackland potential [22], but fitted specifically to provide the phase transition. This potential has been described in detail recently [17] and gives a bcc $\to$ hcp transition at 13.75 GPa. We shall denote it as the Ackland potential.

Our sample contains approximately 8 million atoms. It has a length of 85.6 nm in the direction of the shock (z direction) and an edge length of 30 nm perpendicular to it (x and y directions). The sample was constructed with a Voronoi construction scheme [24] with a log-normal size distribution and a mean grain size of 7.5 nm. This translates into typically 4 grains in the lateral direction (30 mm), and 16 grains along the z direction. The complete sample contains 256 grains.

The sample is relaxed using high-temperature annealing at 80% of the melting temperature for 100 ps. As in our previous work [17], we found this essential for equilibrating the grain boundaries (GBs). The shock wave simulations are performed at a temperature of 10 K. This low temperature was chosen in order to minimize thermal noise and thus to facilitate the analysis of dislocation formation.

The public-domain MD code LAMMPS [23] was used in this study to perform the simulations. Shock waves are initiated by the so-called piston-driven algorithm [25, 26]. Here, a thin slab (thickness of a few atomic planes of Fe, here, one lattice constant) of the simulation volume is treated as the piston in the simulation; atoms are assigned the pre-determined piston velocity U_p in z direction, and zero velocity transverse to it. We study piston velocities varying from U_p = 0.5 to U_p = 0.9 km s⁻¹. In this way, a piston is created that drives into the simulation cell. Note that these atoms are not subjected to the forces of the surrounding atoms. MD simulates the remainder of the system volume under NVE conditions.

The piston velocity is increased linearly from 0 to the shock velocity U_p during a ramp time of 5 ps. We found that a finite ramp time is essential to ease dislocation nucleation [15, 26], since lower strain rates give dislocations more time to nucleate. After the described shock-loading time, the simulation is continued until a total simulation time of 15 ps. We apply periodic boundary conditions in the x and y directions and free boundaries along the shock direction.

To analyze the results, we determined various material properties as a function of the z position in the sample. These profiles were calculated by dividing the sample into small regions. Each of these equally sized regions is a small slab with a thickness of 0.8565 nm along the z-direction. We monitor the velocity profile in z direction, v_z, the stress in z direction, p_zz, the temperature T, and the shear stress ${{p}_{{\rm shear}}}$. The latter is defined as

$$\begin{eqnarray}&&{{p}_{{\rm shear}}}=\frac{1}{2}\left( {{p}_{zz}}-{{p}_{{\rm trans}}} \right),\end{eqnarray} \tag{ 1 }$$

where p_ij denote the components of the stress tensor, and the transverse pressure is defined as

$$\begin{eqnarray}&&{{p}_{{\rm trans}}}=\frac{1}{2}\left( {{p}_{xx}}+{{p}_{yy}} \right).\end{eqnarray} \tag{ 2 }$$

Finally, we analyze the local crystal structure using adaptive common-neighbor analysis [27] in order to quantify the fraction of phase-transformed material. This gives us the fraction of atoms in the sample that are locally in a bcc, hcp or fcc environment. We shall denote the fraction of non-bcc atoms as ${{f}_{{\rm non}-{\rm bcc}}}$; it serves as a quantitative measure for phase-transformed material. For the Ackland potential, we also determine the fraction of close-packed (cp) material, ${{f}_{{\rm cp}}}$.

The samples were visualized using common-neighbor analysis [28, 29] and OVITO [30].

For a qualitative overview, we display in figure 1 the changes induced in our sample by the passage of the shock wave. Figure 1(a) shows the results for the Mendelev potential, which do not allow for phase transition, while the generation of plasticity in the bcc material and changes in the grain structure of the nanocrystalline material are described realistically. The structure after passage of the 0.5 km s⁻¹ shock (top frame) corresponds well to the initial, virgin structure of the material. It is characterized by the GB network, while only a few point defects are visible in the grains. At 0.7 km s⁻¹ (middle frame) only a few changes can be discerned; an example is given by the big grain on the lower side at z = 20 nm. Such features are due to dislocation formation. The GB network is unchanged. Figure 2(a) presents a snapshot of the plasticity developing in this system. Besides straight-screw dislocations crossing grains, shear loops and vacancies are seen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At 0.9 km s⁻¹ (bottom frame of figure 1(a)) the changes induced by plasticity are considerable. Between the piston and around 70 nm, many GBs have broadened and acquired complex shapes. Dislocation lines spanning an entire grain have appeared (see the grain at 40 nm, the middle of the specimen, or at 60 nm in the lower part of the specimen). We conclude that widespread dislocation activity has occurred for this strong shock.

For the Ackland potential, which allows for both dislocation activity and phase transformation, the situation has qualitatively changed considerably; see figure 1(b). The phase transition proceeds with the nucleation of the hcp phase crisscrossed with stacking faults, which appear as fcc phase. For 0.5 km s⁻¹, the phase transition occurs only partially in a few grains, while for 0.9 km s⁻¹, the transformation is nearly complete approximately 10 ps after the shock crosses a given grain.

Plastic activity occurs for the Ackland potential only in a narrow spatial range. Dislocations appear inside the grains just before the phase transition is starting. Figure 2(b) presents a snapshot of the plasticity occurring for the 0.7 km s⁻¹ shock. While dislocations start nucleating in the narrow plastic region ahead of the transformation wave, here we show dislocations in the partially phase-transformed region. The dislocations shown are still inside the bcc phase. We observe only intragranular, but not intergranular, dislocations. Plasticity is similar to the one observed for the Mendelev potential; straight screw dislocations crossing grains, shear loops and vacancies are observed. For the special case of the Ackland potential and the intermediate shock, 0.7 km s⁻¹, a study showing screw dislocations crossing grains in the plastic region and simulated x-ray diffraction has been published elsewhere [15]. The material directly adjacent to the piston ($z\lt 13$ nm) does not phase transform; this is due to the modeling of the piston, where rigid Fe atoms stabilize the bcc state.

The transformed material does not appear in a uniform hcp phase, but rather in a mixture of the two (cp) phases, hcp and fcc, since the free-energy difference between these two phases is small. Nucleation often starts at a GB, and only rarely in the grain interior. The contribution of hcp is dominant; often, the fcc phase only occurs as a stacking fault in the hcp material. Note also that regularly layered structures of hcp and fcc, as they are here seen in some grains, may be associated with nR (rhombohedral cp) structures, which characterize a distorted hcp phase. In a previous study [17], we relaxed such phase-transformed material at a high pressure; the material turned hcp with fcc only appearing as stacking-fault planes within the hcp phase. We note that a recent study [15] finds that the diffraction pattern of the transformed phase features hcp peaks prominently but not the fcc peaks, in agreement with experiment [10].

The phase transformation is analyzed quantitatively in table 1, which presents the phase fractions of the initial bcc and the transformed hcp and fcc phases. The data refer to the entire sample at the end of the simulation, t = 15 ps, and thus comprises both shocked and unshocked material. The column denoted by 'other' represents disordered atoms in GBs, point and line defects, and regions which are not bcc/fcc/hcp, including distorted versions of those structures.

For the Mendelev potential, there is no phase transformation. However, there is a negligible fraction (less than 1%) of cp phase due to defective atoms with a similar topology. The bcc fraction decreases to the favor of 'other' atoms. As the discussion of figure 1(a) above proves, this is directly connected to the buildup of dislocations and other defects while plasticity increases.

For the Ackland potential, the bcc phase is substantially reduced, mainly by the buildup of cp phases. As discussed above, hcp is more abundantly produced, roughly double in size compared to the fcc phase. However, also 'other' atoms are produced, similar in size as for the Mendelev potential. This proves that dislocation generation and GB activity are as abundant in this potential as in the Mendelev potential.

We discuss the dynamics within the shock wave by displaying its spatially resolved characteristics for various times after launching the shock. Figure 3 gives these characteristics for the Mendelev potential. After the loading time of the shock, 5 ps, the particle velocity, v_z, and the longitudinal pressure, p_zz, show a stable profile. This is likely due to the very short time (10 ps) that elapses from the time when ramp loading ends to the end of the simulation, and the relative insensitivity of these coarse mechanical measures to the underlying lattice deformation dynamics. The ramp wave is relatively short, and the resulting width is unlikely to shock-up further, but some long-time relaxation of the microstructure might change the wave profiles for a much thicker sample. We can determine the wave speed from these plots. The shock front moves at 5.8–6.0 km s⁻¹. We note that this is larger than the velocity obtained in a longer sample, where a wave speed of 5.2 km s⁻¹ was observed [15]; this is due to the fact that the wave has not yet reached a steady-state propagation. However, this does not change our general findings and will be appropriate for experiments in extremely thin films [16].

Zoom In Zoom Out Reset image size

Shear stress shows the expected evolution with piston velocity; it remains constant at a level of around 5 GPa behind the shock front for the weak, nearly elastic, shock, while it decreases behind the front for higher piston velocities, due to plasticity and the phase transformation. Shear stresses cannot sustain values larger than around 6 GPa. These are typical stresses for dislocation nucleation from nanocrystalline GBs at high strain rates [15, 17]. The fraction of non-bcc atoms, ${{f}_{{\rm non}-{\rm bcc}}}$, rises in the region behind the maximum of the shear stress, indicating again the abundant formation of defects in the form of plasticity. The temperature rises considerably above the initial level of 10 K, and reaches values of 45, 130 and 220 K with increasing shock strength due to plastic heating.

The corresponding data for the Ackland potential are displayed in figure 4. The profiles of v_z and p_zz show more details in their structure than for the Mendelev potential; they exhibit two changes of slope, which become particularly pronounced for higher piston velocities, ${{U}_{p}}\geqslant 0.7$ km s⁻¹. These points are known as 'knees'; they occur at particle velocities of around 0.15 and 0.45 km s⁻¹ and divide the shocked material into three regions, thus establishing the three-wave structure of the shock. These knees are also visible in the longitudinal stress profile, p_zz, at 10 and 25 GPa, respectively.

Zoom In Zoom Out Reset image size

The three-wave structure corresponds to the splitting of the shock wave into the elastic, the plastically deformed and the phase-transformed regions of the shock, which becomes clear by comparing with the snapshots shown in figure 1(b). Consider, for instance, the strong shock, U_p = 0.9 km s⁻¹, at 15 ps. Phase transformation starts at around 80 nm, as is clear both from figure 1(b), and from the strong increase of the fraction of cp material, ${{f}_{{\rm cp}}}$, at this depth. This position corresponds to the upper knee in the v_z and the p_zz plots. The first knee, however, occurs at around z = 90 nm; this marks the location where the shear stress starts being relieved from its high value of 4 GPa, and, hence, the onset of dislocation formation. Here, the elastic wave enters into the plastic regime. Note that temperature starts increasing at this position; it remains high throughout the shock wave. Dislocation motion, but also GB sliding and phase transformation contribute to generating thermal energy. This temperature increase is even more pronounced than for the Mendelev potential. The profiles of non-bcc atoms, ${{f}_{{\rm non}-{\rm bcc}}}$, show a distinctly different behavior than for the Mendelev potential. Due to the ongoing phase transformations, values close to 100% are now obtained, demonstrating the almost complete transformation of the material. The transformation is directly captured in the fraction of cp material, ${{f}_{{\rm cp}}}$, which can be seen to increase steadily with piston velocity in figure 4.

The shock front travels at 6.3–6.5 km s⁻¹. We note that this is similar to the value found in an increased simulation volume with a length of 400 nm simulated for 70 ps, where a wave speed of 6.3 km s⁻¹ was observed [15]. Due to the smaller duration of the present simulation, it was not possible to deduce reliable values for the two other regions of the shock. In the larger simulation, they moved at 5.9 and 5.5 km s⁻¹, respectively. We note that shock speeds show little dependence on U_p in our simulations for the range of 0.5 km s⁻¹ $\leqslant {{U}_{p}}\leqslant 0.9$ km s⁻¹ covered by our simulations. Kadau et al report a similar finding in their MD simulations for piston speeds ${{U}_{p}}\lt 1$ km s⁻¹ [14]. In this regime, shock speed is roughly independent of U_p and assumes a value of 5.3 km s⁻¹ for a single crystal and 5.6 km s⁻¹ for a polycrystal. This behavior is typical in the regime where plasticity is negligible [31].

We studied the phase-change kinetics in more detail by considering the temporal history of a specific slab of material located between z = 40 and 50 nm; see table 2. The fraction of bcc material continuously decreases with time during the passage of the shock wave. The time for the phase transformation is seen to be of the order of 10 ps. Note the sharp increase of the bcc fraction for the 0.9 km s⁻¹ shock at around 10 ps after shock launch. Such short times are qualitatively consistent with the results described by kinetic models of the $\alpha \to \epsilon$ phase transformation in Fe, such as the overdrive model proposed by Smith et al [11]. A quantitative comparison is not straightforward since the continuum kinetic models consider transition rates under fixed thermodynamic conditions (temperature and pressure) which are not present in our MD simulation.

In our simulations, pressures, p_zz, of the order of 25 GPa are needed to initiate the phase-transformation process. Note that this is well above the equilibrium phase transition pressure of 13 GPa, which is also implemented in our potential. We argue that the 25 GPa threshold and $\sim 10$ ps kinetic timescale are specific for the very short time and length scales considered in our simulations; they can be compared with thin-film experiments [11, 16] which also hint at an increased critical pressure necessary for phase transformation. Conventional macroscopic experiments [5] report the transition to occur at the equilibrium pressure of 13 GPa and a transformation time scale of tens of ns.

As a summary, in this paper we present large-scale shock simulations in polycrystalline Fe with a new interatomic potential that describes both the phase transition and dislocation-based plasticity in the bcc phase, as it has been recently shown for homogeneous compression of Fe samples [17]. The shock profiles show a three-wave structure as it is observed in experiments, corresponding to elastically compressed bcc, compressed bcc with plasticity, and phase-transformed cp phase. Most of this transformed phase has an hcp structure, but we still find a fraction of the fcc phase, mostly as planar faults inside hcp grains.

In detail, we find the following.

• (i)  
  The successful modeling of shock wave interaction with a phase-transforming material needs both dislocation sources (in this case GBs), and a finite shock-loading time to allow dislocation nucleation before the phase transformation occurs. As a result, the experimentally observed three-wave structure can also be found in simulation, where plasticity precedes phase transformation.
• (ii)  
  For strong shocks, the high shear stress generated is first relieved by dislocation generation. Phase transformation occurs in the defective bcc grains, starting from GB junctions.
• (iii)  
  Particle velocities in excess of 0.6 km s⁻¹, leading to longitudinal stresses of 25 GPa, are needed for a nearly complete phase transformation. This is considerably above the equilibrium transformation pressure of around 13 GPa, which is faithfully modeled by the interatomic potential we used. The high transition pressure is caused by the small time and length scales available and can be compared to experiments performed in thin-film systems [11, 16].
• (iv)  
  The phase transformation is completed within around 10 ps. This short time is due to the high transition pressure; during macroscopic experiments, the transformation takes tens of ns.

We note that recently simulated x-ray diffraction patterns could be obtained [15], which are consistent with the atomistic analysis, showing a clear hcp signal, and lack of a clear fcc signal, together with polycrystal ring broadening due to significant plasticity in the phase-transformed grains. Thus our simulations, together with a new generation of high-pressure experiments that allow 'femtosecond visualization of lattice dynamics in shock compressed matter' [32] permit an in-depth, atomistic understanding of the processes occurring during shock wave propagation in phase-transforming materials.

The results presented in this paper go beyond Fe and are of general interest to the shock-loading, high-pressure community. For example, the specific topic of reproducing the experimentally observed shock-induced plasticity was recently studied for ceramics; see [33, 34]. In addition, high-pressure loading can lead to novel materials with improved properties [35].

This work has been supported by the Deutsche Forschungsgemeinschaft via the Sonderforschungsbereich 926. EMB acknowledges support from CONICET, SeCTyP (U. N. Cuyo), and PICT-2008–1325.