Deformation and reverse phase transformation mechanism of high-pressure HCP iron during unloading process

Iron is the most abundant metallic element on earth, the most important metal, and the most stable nucleus on the periodic table. The phase transformation and plasticity of iron under high pressure have been paid much attention. In 1956, Bancroft et al.¹ first observed the phase transition of shock-compressed polycrystalline iron under shock compression at approximately 13 GPa through wave structure measurements. Subsequent static high-pressure x-ray diffraction experiments confirmed the phase transformation and identified the structure of the high-pressure HCP structure.² Laser-driven loading experiments of polycrystalline iron revealed that the α–ɛ phase transition has a strong impact on the structural transition of iron.³ After the α–ɛ phase transformation, the spall surface is smoother and has higher density twinning. The samples without phase transformation showed brittle fracture characteristics and lack twinning in the recovery. In another study of the recovered samples after shock loading, there were rod-like structures in the iron.^4,5 Due to the rapid phase transformation and the structural transition process, it is experimentally difficult to directly observe the above process at the atomic level in ultra-short spatiotemporal scales.

In recent years, benefiting from the rapid development of computer technology, molecular dynamics (MD) has been widely used to study the structural transformation of iron. This method offers an atomic-scale perspective to understand the microscopic mechanism and macroscopic dynamic behavior. It can offer insight where experiments are unavailable. In 2005, Kadau et al. reported the structural transformation mechanism of iron from α to ɛ under compression by MD simulation and proposed the shuffling mechanism for compression along the direction [001].⁶ In addition, a series of factors affecting the structural transformation, such as phase transformation,⁷ dislocations,⁸ loading orientations,^9–11 loading methods,^11,12 twinning,^13–16 temperature,^17–19 defects,^18,20–24 polycrystals,^25–28 and strain-rate effect,^11,29,30 have also been well researched. Despite these advances, the α–ɛ–α reversibility of iron remains underexplored, and the dynamic pathway and the microscopic mechanism of its reverse transformation are yet to be clarified. In 2015, Gunkelmann et al. confirmed that polycrystalline iron would leave higher density twinning after the α–ɛ–α phase transition through the different potentials.²⁵ In 2020, Shao et al. revealed the microstructure differences caused by different loading orientations.^9,31 In 2021, Guo et al. studied the compression of single crystal iron under uniaxial and triaxial conditions, and the results showed that the former showed good reversibility, while the latter showed poor reversibility with a pronounced strain-rate dependence.³² Further investigations showed the difference in microstructure evolution under different compression paths.³³ In addition, the research on nanowires also showed that there are different deformation mechanisms under different loading paths.³⁴ However, how to effectively establish the link between the unloading path and the microstructure transition is still an open question.

All previous research on the reversibility of structural transformation of iron starts from the initial BCC phase. However, loading along different crystal directions results in the formation of complex composite phase structures, which will bring great challenges to the subsequent research on the reverse phase transition process. Thus, a more effective and simplified model could facilitate a deeper understanding of the above mechanism. To achieve this goal, our simulations start from the high-pressure HCP single crystal of iron for the first time, unloading along three different orientations to simulate the possible unloading paths in experiment and to research the physical phenomena and mechanisms involved in this process, which will provide a physical reference for further understanding the structural transformation of iron.

Section II provides details of the model construction. Section III shows the simulation results, primarily including the waveform, the microstructure view, and the phase fraction. Section IV, the coupling between phase transformation and plasticity, the formation and evolution of twinning, and the polycrystal and damage fracture process of the HCP single crystal iron are discussed.

The modified version of the Ackland potential^35,36 is used in our simulations, which specifically fit to the BCC–HCP phase transition and the reverse phase transition of iron. Due to the significant anisotropy of HCP structures under different unloading paths, we have established three models with different unloading orientations. As shown in Fig. 1, the unloading orientations are along the normal directions of the $( 1 ¯ 2 1 ¯ 0)$⁠, $( 10 1 ¯ 0)$⁠, and $(0001)$ planes, corresponding to unloading along the x, y, and z directions in the Cartesian coordinate system.

The model uses two unloading methods: one is unloading from both sides of the boundary (free unloading model) and the other is unloading at a certain strain rate (strain-rate unloading model). The latter can well reflect the evolution of phase transformation during the unloading process and help to identify lattice orientation.

For the free unloading model, the length of the unloading direction is set to 120 nm (individually for x, y, and z), while the length of the other two directions is set to 25 nm, containing a total of about 7.5 × 10⁶ atoms. We set the initial pressure of the model to 15, 25, and 40 GPa, corresponding to lower than, slightly higher than, and far higher than the HCP–BCC equilibrium pressure, respectively.^35,37 Before unloading, the system is first equilibrated under the conditions of 60 K and 15, 25, and 40 GPa pressure in the isothermal–isobaric ensemble (NPT) for 40 ps. Then, the boundary of unloading orientation is set to a free-boundary condition for 60 ps so that structural transformation is fully completed. For the strain-rate unloading model, the model size is $30×30×30 n m 3$⁠, containing about 2.7 × 10⁶ atoms. The initial pressure of this model is set to 40 GPa, the temperature is 60 K, and the relaxation reaches equilibration in the isothermal–isobaric ensemble (NPT) for 40 ps. Afterward, unloading is performed at a strain rate of $4× 10 9 s − 1$ along the x, y, and z directions, respectively.

In the analysis, the OVITO program is used to obtain visualization models. The “common neighbor analysis” (CNA)³⁸ is used to identify the phase atom types, the “Polyhedral template matching” (PTM) and “Grain segmentation” are utilized to count the number of grains, and the “dislocation extraction algorithm” (DXA) is utilized to analyze dislocation. Different colors are used to distinguish phase atoms: HCP (red), BCC (blue), FCC (green), and “other” atoms (gray, represent disordered atoms).

As shown in Fig. 2, for the free unloading model along the x direction, the final products at different initial pressures are consistent. Except for the BCC phase, there are rod-like structures composed of disordered atoms, which result from twinning. These residual rod-like structures were reported in previous experiments,^4,5 and a detailed explanation of the reason for the formation of the rod-like structure will be given in Sec. IV A.

For the free unloading model, the unloading area can be divided into two parts. Taking Fig. 2(b) as an example, before 10 ps, the unloading front starts from the free surfaces at both ends and gradually propagates inward. This is Part-I. The second part starts at about 11 ps. At this moment, although the phase transformation wave has not yet propagated to the middle region, the stress value in this region has already decreased to a certain value. Therefore, the unloading process from HCP to BCC occurs in this region, also leaving the rod-like structure. This is Part-II. We count the reverse transition threshold separately for the two parts, as shown in Table I.

Interestingly, for the middle unloading region (Part-II), as shown by the waveform diagrams and the corresponding microstructure in Fig. 2(a), the stress had already decreased to 0.5 GPa at 11 ps; yet, the HCP structure was still maintained without a transition. Subsequently, the unloading process from HCP to BCC occurred at 12 ps, with the stress increasing to about 2 GPa. This phenomenon is called “overunloading,” which means that the stress has already dropped below the threshold of a structural transition, but the transformation does not occur until the stress further decreases. This is because the inward propagating unloading waves from both sides of the boundary will first meet in the middle region of the model to form a low-pressure zone, but at this time, the phase transformation wave has not yet propagated to this region. The perfect HCP single crystal structure used in the model lacks nucleation sites, which ultimately leads to the overunloading in the middle region. Adding initial defects into the model can effectively avoid the occurrence of overunloading in the middle region because it provides favorable nucleation sites. After that, a negative stress zone is formed in the middle region, as shown in Figs. 2(a) (29 ps), 2(b) (22 ps), and 2(c) (21 ps). As shown in Table I, as the initial stress of the model increases, the phase transition threshold of Part-II also increases, specifically at 0.5, 4.2, and 9.0 GPa, respectively.

The unloading wave exhibits an obvious double wave. The first wave (close to the free boundary) is the phase transition and the plastic coupling wave, while the second wave is an elastic unloading wave. The region between these two waves represents the process of the transition from elastic to plastic. As shown in Figs. 2(b) at 6 ps and 2(c) at 5 ps, when the unloading process does not produce rod-like structures, the interface between HCP and BCC maintains relatively smooth, and the waveform of the first wave maintains a smooth state. In this case, the first wave only represents the phase transition wave from HCP to BCC. Conversely, when rod-like structures are produced during the unloading process, both phase transformation and twinning of plastic deformation occur simultaneously in the first wave, resulting in waveform fluctuations, as shown in Figs. 2(b) at 9 ps and 2(c) at 7 ps. As shown in Table I, the reverse phase transition threshold of Part-I increases with the initial pressure, specifically at 4.1, 6.8, and 11.0 GPa. Similarly, the reverse phase transition threshold for both Part-I and Part-II increases with increasing initial stress, and this phenomenon also exists for another two unloading directions. This indicates that for the high-pressure HCP single crystal iron, the threshold of the unloading reverse phase transition is dependent on the initial stress.

For the model freely unloading along the y direction, the unloading volume can be divided into two parts, exhibiting the “overunloading” phenomenon mentioned above. The overunloading stress can reach −2.3 GPa, as shown in Fig. 3(a) at 10 ps. Interestingly, the twinning structure of the unloading product exhibits excellent directionality, with consistent orientation and forming lamellar twinning grains. As the unloading wave meets to form the negative pressure zone, the original twinning structure undergoes detwinning and reconstruction. Then, a new twinning structure parallel to the unloading direction will be formed, leading the grain orientation to change. The specific mechanism of this process is analyzed in detail in Sec. IV A. It should be noted that during the formation of the twinning structure, there is an obvious fluctuation in the stress waveform, corresponding to the structural transition, as shown in Figs. 3(a) (11 ps), 3(b) (10 ps), and 3(c) (9 ps). As shown in Table I, reverse phase transition thresholds for Part-I are 1.0 GPa, 4.9 GPa, and 9.4 GPa, while those for Part-II are −2.3 GPa, −1.2 GPa, and 2.0 GPa. Both sets of thresholds increase with the increasing of initial stress. Additionally, the unloading wave exhibits an obvious double wave structure, comprising the phase transition and the plastic coupling wave followed by the elastic unloading wave.

For the model freely unloading along the z direction, the phenomenon of “overunloading” is more pronounced, as shown in Fig. 4(a) at 10 ps. The overunloading stress even reaches −10.5 GPa without undergoing any structural transition. As shown in Table I, the unloading reverse phase transition threshold along the z direction is lower than that for the other two directions, indicating that it is more difficult to unload along the z direction. Meanwhile, as shown in Fig. 4(a), for the model with an initial stress of 15 GPa, unloading starts from the middle part rather than from the free surface, indicating that the unloading process in this direction is influenced by the initial stress. In addition, the unloading products along the z direction differ from those in the other two directions, and the unloading process transitions through HCP to FCC and finally to the BCC phase structure. This process leaves many grain boundaries, resulting in polycrystalline formation. The reason for this phenomenon is explained in detail in Sec. IV B.

When the structural transition occurs, an obvious double wave structure is observed, comprising the phase transition and the plastic coupling wave followed by the elastic unloading wave. Interestingly, for Figs. 4(b) at 6 ps and 4(c) at 4 ps, there is a significant time gap between the elastic unloading wave and the subsequent phase transition and the plastic coupling wave. This indicates that the structural transformation along this unloading direction lags the stress unloading process.

Figure 5 shows the variation of the phase fraction along different unloading directions. The results show that the BCC phase fractions in the final product for unloading along the x, y, and z directions are about 90%, 80%, and 70%, respectively, with the remaining portions primarily consisting of disordered atoms. These disordered atoms form rod-like structures, lamellar twinning, and polycrystalline boundaries in the three different unloading directions. As shown in Fig. 5(d), during the unloading process (before about 20 ps), only unloading along the z direction will produce a significant fraction of the FCC phase. This is because unloading in this direction involves a transition from HCP to FCC, followed by FCC to BCC, which is confined to regions near the unloading surface. Therefore, the overall proportion of FCC in the model does not exceed 9%. The strain-rate unloading model could better reflect this process. When the unloading waves meet to form a tensile zone, a very small amount of FCC structure will be formed at the intersection of the structures under the action of tension.

As the initial stress increases, with 60 ps as the final moment, the fraction of the BCC phase in the final product decreases slightly, while the fraction of disordered atoms increases slightly, with both changes not exceeding 5%. The FCC phase generated during the unloading process increases with the increase of initial stress. After unloading, the remaining FCC phase is less than 0.5% for unloading along the x–y direction and less than 2.1% for unloading along the z direction.

Taking the 25 GPa initial stress as an example, as shown by the arrows in Figs. 5(a2)–5(d2), an inflection point appears in the curve, indicating an acceleration in the unloading rate. This is because that the structural transformation in the middle region occurs. As shown in Figs. 5(a1)–5(a3), when the HCP fraction drops to about 0%, it signifies complete unloading. It can be found that under the same initial stress conditions, the unloading rates along the x and y direction are faster compared to that along the z direction.

Figure 6 shows the simulation results of strain-rate unloading along the x, y, and z directions. All the stress curves have a clear inflection point, where the stress (the red line) rapidly increases under a very small strain. This inflection point serves as a signal of structural transformation. It should be noted that we take the values at both ends of the inflection point as the value of the phase transition threshold. For the x, y, and z directions, the respective phase transition thresholds are (10.2, 12.7), (6.3, 16.0), and (−7.6, 8.2) GPa. These correspond to the phenomenon of “overunloading” mentioned above. At these points, the microstructural transformation views are shown in Figs. 6(a2), 6(b2), and 6(c3). For unloading along the x and y directions, the HCP to BCC phase transition occurs first, followed by the formation of rod-like or lamellar twinning structures. In contrast, unloading along the z direction first occurs with the formation of disordered atomic clusters in the HCP structure, which act as phase nucleation sites for the HCP to FCC phase transition. (This requires a large tensile force to destroy the lattice, leading to form disordered atoms, so the stable HCP structure will still exist when stress is negative.) Finally, it will form the polycrystalline BCC structure. This indicates that when unloading along the z direction, the phase transition needs to pass through the FCC intermediate state.

For the maximum shear stress (maximum absolute value, the blue line) in different unloading directions, the x direction is about −3.9 GPa and the y direction is about −7.3 GPa. Interestingly, the z direction is about −16.6 GPa, much higher than other directions, which is the direct reason for the rapid increase in disordered atoms observed in Fig. 7(d). In addition, during the unloading along the y and z directions, the trends in shear stress and stress exhibit a strong correlation, with inflection points coinciding for these. However, for the unloading along the x direction, no structural transformation occurs when the shear stress reaches its minimum value.

As shown in Fig. 7, for the fraction of the BCC phase in the final product, the x direction is the highest, followed by the y direction, and the lowest is the z direction. This is due to the formation of regular rod-like and twinning structures during unloading along the x and y directions. However, the polycrystalline structure is produced during unloading along the z direction, having a large number of grain boundaries. In addition, the phase transition from HCP to FCC during the unloading along the z direction leads to the maximum fraction of the FCC reaching about 29%.

As shown in Figs. 6(d)–6(f), during unloading along the x direction, the number of grains gradually increases with increasing strain at the beginning. This is because the continuously formed rod-like structures divide the original single crystal into polycrystals. Then, the number of grains has a clear increase in the highlighted area in Fig. 6(d). This is due to the occurrence of void nucleation in this area. During unloading along the z direction, the number of grains decreases with increasing strain in the highlighted area, which is due to grain merging. The specific mechanism will be analyzed in detail in Sec. IV B. The subsequent increase in the number of grains is also due to void nucleation and fracture. It is interesting that in the case of unloading along the y direction, the number of grains in the highlighted area is significantly higher compared to other strain areas. This is because the initially formed lamellar twinning structures reach a peak in this area. Then, the system undergoes detwinning and reconstruction processes, resulting in the merging of different grains, leading the number of grains to greatly decrease.

Previous research has shown that single crystal iron exhibits excellent reversibility during the loading–unloading process along the [001] crystal direction. After the initial BCC structure undergoes a phase transition to HCP under high pressure, it can almost fully revert to its original BCC structure after unloading without leaving any other structure.^35,37 However, the unloading behavior in other directions is not clear. This section analyzes and discusses the unloading process of HCP single crystal iron along three different crystal directions. The results show that it has highly anisotropic characteristics and cannot return to the perfect BCC phase structure. This process involves coupling between plasticity and phase transitions. For unloading along the x and y directions, the HCP phase directly transforms to the BCC phase, followed by forming the different types of twinning due to differences in the phase transition mechanisms at the lattice level, with subsequent deformation dominated by twinning. For unloading along the z direction, defects are initially formed due to significant overunloading and shear stress, resulting in many disordered atoms. Subsequently, the phase transition occurs starting from the defects, accompanied by rapid formation of dislocations and interactions between grain boundaries.

Figures 8(a) and 8(b) show the three-dimensional lattice changes of the HCP to BCC transition, while (a2), (a3), and (b2) provide top views corresponding to (a1) and (b1), respectively. The simulation reveals two different types of HCP to BCC phase transitions: type I and type II, with the former having two variants and the latter having only one variant. For type I, an initial force is applied in the y direction. First, atoms on the $( 10 1 ¯ 0)$ plane of the HCP undergo opposite shear (the black arrows), as in Fig. 8(a2). Subsequently, due to stretching in the y direction, the atoms on the $( 10 1 ¯ 0)$ plane will move along the y direction. Simultaneously, the second layer of atoms (the red triangles) under the action of the shear stress will undergo shear displacement (the green arrow) along the XY direction, forming two vertices of the BCC lattice. The above steps transform the (0001) plane of the HCP to the (011) plane of the BCC (the blue line). When the initial shearing step occurs in the opposite direction, another variant of the BCC structure is produced [Fig. 8(a3)], with an intersection angle of about 112° between the two variants, as shown in Fig. 9C.

For type II, the initial force is applied in the x direction. First, atoms on the $( 1 ¯ 2 1 ¯ 0)$ plane of the HCP undergo opposite shear forces (the black arrows), as shown in Fig. 8(b2). The second layer of atoms (the red triangle) under the action of the shear stress will undergo shear displacement (the green arrow) along the Y direction, forming the two vertices of the BCC lattice. The above steps change the (0001) plane of HCP to the (011) plane of the BCC (the blue line).

There are two types of twinning formed by the HCP–BCC transformation. The first type, as shown in Figs. 9(a)–9(d) and B, involves the formation of a lamellar BCC phase along the (010) plane during unloading, which also serves as the twinning plane. Enlarged views A and C show that all BCC structures generated through this phase transition belong to type I, characterized by two variants. Due to the different lattice orientations, these two variants have an intersection angle, as shown in C, which is about 112°, thereby forming the twinning structure. Figures 9(a2)–9(d2) show the strain tensor-XY, indicating that regions where the HCP–BCC transformation occurs are accompanied by the change of the strain tensor-XY value. Notably, adjacent twinning regions exhibit opposite strain tensor-XY values, corresponding to shear displacements along different directions for the second layer of atoms in type I [the red triangle in Fig. 8(a)]. In brief, for unloading along the y direction, there is only one shuffling mechanism for the HCP–BCC phase transformation, with two variants, forming lamellar twinning with different lattice orientations.

For the second type, as shown in Fig. 10, unloading along the x direction reveals two distinct HCP–BCC phase transition types, type I and type II, from the enlarged views of B, C, and D diagrams. These transitions lead to the formation of a twinning structure in D. The twinning planes are along $( 1 ¯ 3 3 0 )$ and $( 1 3 3 0 )$⁠. From Figs. 10(a)–10(d) and A, it is shown that the HCP to BCC phase transformation is first carried out with type II before ɛ = −0.089. Because of the same type of phase transition, the point-line BCC structure generated at ɛ = −0.097 will be connected into a long strip structure at ɛ = −0.089. Subsequently, the HCP to BCC (c1)–(d1) transition primarily involves the two variants of type I. As shown in Figs. 10(b1) and (b2) and 8(b2), the type II phase transition mechanism is accompanied by the strain in the y direction. However, type II has only one variant type, which cannot effectively balance the strain in the y direction, leading to the concentration of strain tensor-YY in the HCP region that has not undergone a phase transition, as shown in Figs. 10(a2)–10(d2). It will lead to the occurrence of type I of HCP–BCC in this region. As shown in Figs. 10(c3)–10(d3), when type I starts to occur, the strain tensor-XY value changes. This process forms twinning, which is the rod-like structure mentioned above. The reason for the formation of the rod-like structure had also been reported across the prior MD study.⁴ However, the mechanisms differ from those observed in our simulation. The supplementary material provides a detailed comparison between our results and prior studies. In brief, for unloading along the x direction, there are two different shuffling mechanisms for the HCP–BCC, with type II occurring first, followed by type I, and finally forming a rod-like twinning structure.

In previous experiments, it has been confirmed that the shock recovered polycrystalline iron will leave rod-like^4,5 structures (corresponding to unloading along the x direction) and high density twinning^3,39 (corresponding to unloading along the y direction) after loading and unloading, and this simulation result provides a microscopic explanation for this phenomenon. Here, we propose a possible mechanism for the formation of the rod-like structure and twinning. Of course, there may be other mechanisms that require further experimental and theoretical research.

As the simulation evolves further, the original positive pressure gradually transforms into negative stress (tensile), and the model unloading along the x and y directions will undergo detwinning and reconstruction, leading to a change in the number of grains. For unloading along the y direction, the formed (010) plane of twinning is perpendicular to the unloading direction. During the subsequent stretching stage, due to the resistance effect of the twinning structure, it is difficult to further stretch, resulting in detwinning and reconstruction processes. This ultimately results in the formation of (100) twinning planes parallel to the stretching y direction, as shown in Figs. 11(f) and B. The detwinning process begins with grain rotation, as shown in the yellow elliptical area of Fig. 11A. The BCC atoms belonging to different grains undergo rotation, which is the conversion between the two variants of type I, resulting in the disappearance of the original twinning boundary and the connection of the upper and lower regions, forming new (100) planes of twinning. This rotation eventually reaches 90°, destroying the original grain structure and causing the merging of grains, which explains the observed reduction in the grain count in the highlighted area of Fig. 6(e).

For unloading along the x direction, as shown in Figs. 12(a) and 12(b), part of the twinning structure disappears, showing that detwinning occurs under tensile stress. With further stretching, a new twinning structure along the (101) plane is formed. Interestingly, the newly generated twinning structure intersects with the pre-existing twinning, leaving a large number of disordered atoms at the intersection region. This results in an unstable structure, which provides potential sites for subsequent void nucleation. As shown in Figs. 12(A) and 12(B), the (011) plane of the BCC structure, marked by the black square, shows that the newly generated twinning structure is different from the initial twinning structure in terms of lattice orientation. Importantly, this transformation occurs without any phase transition, indicating that the formation of new twinning is due to the rotation of the BCC lattice.

For unloading along the z direction, as shown in Fig. 6(c), structural transformation initiates at the stress of about −7.6 GPa and the shear stress of about −16.6 GPa. At this time, randomly distributed disordered atoms will be formed in the HCP single crystal [Fig. 13(a)]. These disordered atoms act as nucleation sites for the HCP to FCC phase transition. The stress then rapidly rises from −7.6 to 8.2 GPa with the stress reversal phenomenon, which is due to drastic structural transformation. At the same time, as shown in Fig. 13(d), under the action of high shear stress, this process is accompanied by significant plastic deformation, forming Shockley dislocations and grain boundaries. Previous research has shown that there are multiple phase transition paths from HCP to FCC.^40–43 Due to the randomness of the distribution of disordered atoms offering phase transition sites, multiple types of HCP to FCC paths occur simultaneously. Also, this process occurs extremely rapidly, with disordered atoms reaching the peak within the strain of 0.008 [blue line in Fig. 7(b)]. The random distribution of phase transition sites leads to the newly generated FCC lattice having different orientations, and this process is accompanied by rapid formation of Shockley dislocations, leading to the formation of grains. The alternating plane structure formed with the FCC and HCP planes in Fig. 13(h) is due to the phase transition mechanism of HCP to FCC. The slip planes formed are $( 1 ¯ 3 3 0 )$⁠, $( 1 3 3 0 )$⁠, and $(100)$ planes, all parallel to the z-axis. The simulation results show that the BCC phase first forms mainly at the grain boundaries and then propagates into the interior of the grain, as shown in Figs. 13(c) and 13(d). It indicates that the grain boundaries can be regarded as the phase nucleation sites, which favors the reverse phase transformation from HCP and FCC to BCC.

As shown in the black circle of Figs. 14(a) and 14(b), for unloading along the z direction, grain boundaries formed during the unloading process merge with the phase transition process. As shown in the black circle of Figs. 14(a1) and 14(a2), the merging of grains is attributed to the similar orientations of the alternating FCC and HCP plane structures in adjacent grains. The small grain shown by the black circles in Figs. 14(b1) and 14(b2) also exhibits the merging with the two large grains (regions A and B). However, unlike Fig. 14(a), the upper A and lower B grains have similar lattice orientations, while the small grains have different orientations from A and B. As further unloading occurs, the large grains A and B tend to merge, while the phase structure in the middle small grain gradually returns to disordered atoms and becomes assimilated by the A and B grains. The yellow highlighted region of Figs. 14(c) and 14(d) shows the grain merging process during the stretch stage. With tensile strain increasing, Fig. 14(c) represents the merging between two large grains, which also requires similar lattice orientations between them. Figure 14(d) shows that small grains are assimilated by their surrounding grains. The principle of grain merging in Fig. 14(a) is consistent with that in Fig. 14(c), and the principle of small grain merging in Fig. 14(b) is consistent with that in Fig. 14(d). In addition, as shown by the yellow arrow in Fig. 14(d), the grain boundary migrates during the stretching process. All of these lead to an increase in the grain size, which explains the decrease in the number of grains in the highlight of Fig. 6(f).

When a negative stress zone is generated inside the model, a tensile force is formed. When the threshold of void nucleation is reached, voids will be generated, which eventually leads to spallation. The fracture process is accompanied by obvious plastic deformation. As mentioned in Secs. IV A and IV B, for unloading along the x and y directions, the plastic evolution in the stretch stage is mainly twinning, detwinning, and reconstruction. As shown in Figs. 15(a) and 15(b), void nucleation mainly occurs at the zone of twinning intersection, which is consistent with our previous study.³⁷ This phenomenon occurs because the twinning intersection zone represents the intersection of lattices with different orientations, which leads to the destruction of lattice stability and provides preferential sites for void nucleation. After the spallation, the model unloading along the x and y directions will leave rod-like structures and higher density twinning, respectively, which is consistent with the previous experimental results.^3–5,25 For the model of unloading along the z direction during the stretch stage, the plasticity is mainly based on the behavior of grains, including grain growth, merging, and boundary migration. A large number of grain boundaries provide sites for void nucleation.

Twinning and grain boundaries both favor void nucleation, leading to formation of many voids at multiple locations near the spall surface, making it easier to leave a smooth spall surface from a macroscopic view. This is consistent with the smooth spall surface left by Fe undergoing the α–ɛ–α phase transition in experiments.^3,39

Figure 16 shows the spall strength for different unloading orientations, methods, and initial stress. Under the simulation conditions we set, the strain-rate unloading model has lower spall strength. For the free unloading method, the spall strength increases with increasing initial stress. Among them, the spall strength unloading along the z orientation is the lowest because the polycrystalline structure favors void nucleation over twinning.

In this work, the microstructure evolution and related mechanical characteristics for three different unloading orientations are studied by constructing perfect high-pressure HCP single crystal iron. The following are the main findings:

1. The overunloading phenomenon exists under three different unloading orientations. Specifically, the high-pressure HCP structure has reached the unloading threshold, but the structural transition has not occurred. For unloading along the normal direction of the (0001) plane, there will be an anomalous existence of the HCP structure under negative stress.

2. Before the structural transition, unloading along the normal direction of the (0001) plane has the largest shear stress, followed by $( 10 1 ¯ 0)$ and $( 1 ¯ 2 1 ¯ 0)$ planes. However, unloading along the normal direction of the (0001) plane has the lowest spall strength.

3. Two types of twinning formed by HCP–BCC transformation are revealed. For unloading along the normal direction of the $( 10 1 ¯ 0)$ plane, only one HCP–BCC shuffling mechanism occurs (type I in the paper), which has two variants, forming the lamellar twinning with different lattice orientations. For unloading along the normal direction of the $( 1 ¯ 2 1 ¯ 0)$ plane, two HCP–BCC shuffling mechanisms occur (types I + II in the paper), leading to rod-like twinning. The subsequent evolution process is dominated by twinning. When the unloading process gradually evolved into a stretch process, detwinning and reconstruction would occur, accompanied by grain rotation. Void nucleation mainly occurs at twinning intersections.

4. For unloading along the normal direction of the $(0001)$ plane, the unloading phase transition needs to go through the FCC phase, forming $( 1 ¯ 3 3 0 )$⁠, $( 1 3 3 0 )$⁠, and $(100)$ slip planes parallel to the unloading orientation, and then return to the BCC phase. Due to overunloading, the stability of the lattice is destroyed, and randomly distributed defects (disordered atoms) are formed, providing phase transition sites. This process is accompanied by rapid formation of Shockley dislocations. Meanwhile, multiple transition paths from HCP to FCC result in formation of grain boundaries and ultimately leaving polycrystalline structures. This process is accompanied by the coupling between dislocations, grain merging, and phase transitions. Void nucleation mainly occurs at grain boundaries.

5. There are two reasons for the decrease in the number of grains. The first is the merging of the lamellar grains due to the detwinning and reconstruction for unloading along the normal direction of the $( 10 1 ¯ 0)$ plane. The second is the merging between grains with similar lattice orientations, or small grains assimilated and merged by surrounding grains for unloading along the normal direction of the (0001) plane.

See the supplementary material that shows the comparison of our MD results about a rod-like structure with the study of Ref. 4.

This work was supported by the National Natural Science Foundation of China-NSAF (Grant No. U2030117).

The authors have no conflicts to disclose.

Jinmin Yu: Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Jianli Shao: Conceptualization (equal); Data curation (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Hua Shu: Conceptualization (equal); Data curation (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Xuyang Ma: Software (equal). Xichen Zhou: Software (equal). Xiuguang Huang: Project administration (equal). Sizu Fu: Project administration (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.