Atomistic studies of shock-induced phase transformations in single crystal iron with cylindrical nanopores

Computational Materials Science 122 (2016) 1–10

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Atomistic studies of shock-induced phase transformations in single crystal iron with cylindrical nanopores Li Wu a, Kun Wang b, Shifang Xiao a,⇑, Huiqiu Deng a, Wenjun Zhu b, Wangyu Hu c a

Department of Applied Physics, School of Physics and Electronics, Hunan University, Changsha 410082, China National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, Mianyang 621900, China c College of Materials Science and Engineering, Hunan University, Changsha 410082, China b

a r t i c l e

i n f o

Article history: Received 19 January 2016 Received in revised form 9 May 2016 Accepted 11 May 2016

Keywords: Phase transformation Variants Cylindrical nanopores Transformation kinetics NEMD

a b s t r a c t Non-equilibrium molecular-dynamic simulations with our modified analytic embedded-atom model potential were carried out on single crystalline iron of idealized cylindrical nanopores to study the effect of defect on shock response. The results showed that the shock response and phase transformation were influenced by cylindrical nanopores. Reflection wave, ‘‘hot spot” can be clearly observed under [0 0 1], [1 1 0] and [1 1 1] loading directions, and the nucleation and growth of dislocations appeared for the shock along [1 1 0] and [1 1 1] directions. The critical stress of phase transformation, the nucleation sites and martensitic variants were obviously influenced by the cylindrical nanopores. In addition, stress assisted transformation mechanisms and strain induced transformation mechanisms were discussed, and the kinetics of the variants which caused by the presence of cylindrical pores were studied. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction The study of the mechanical and structural response behavior of solids under extreme high-pressure is important to understand geophysical, astronomical processes and technological application. Since iron is the main component of inner core of Earth and has been widely used in modern industry and defense, its highpressure phase transition has been investigated extensively. The martensitic phase transition a ? e of iron was first reported by Bancroft et al. in 1956 [1]. Shock pressure of the phase transformation of iron varies from 13 GPa to more than 30 GPa, depending on the stress states [2], strain rates [3,4] and the shear deformation of samples [5]. The mechanism of phase transition under high pressure has been a focus of interest for a long time. To date, three possible mechanisms of the a ? e phase transition of iron have been proposed. The first mechanism is that the uniaxial compression along [0 0 1]BCC crystal direction make the atoms on (1 1 0)BCC crystal plane arrange into hexagon, then the shuffling along alternate (1 1 0)BCC planes transform the crystal into HCP structure with c axis lying along the [1 1 0]BCC direction [6,7]. The second mecha 1 1  BCC direc 1 2ÞBCC planes along the ½1 nism is that, the shear of ð1 tion forces atoms in the (1 1 0)BCC planes into hexagonal patterns, then the relative shuffle between the (1 1 0)BCC crystal planes trans⇑ Corresponding author. E-mail address: [email protected] (S. Xiao). http://dx.doi.org/10.1016/j.commatsci.2016.05.010 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

form the crystal into HCP structure [8]. The third mechanism is similar to the second one in the first step which forces atoms in the (1 1 0)BCC planes into hexagonal patterns, but then a metastable FCC phase appears before the final HCP phase forms [8]. Recently, the relationship between plasticity and a ? e phase transition of iron has received extensive attention. The coupling between plasticity and martensitic transformation was widely discussed in the study of austenitic stainless steel, and two coupling modes, stress assisted coupling mode (SACM) and strain induced coupling mode (SICM), were proposed [9,10]. The two coupling modes correspond to two different phase transition mechanisms, namely, stress assisted transformation (SAT) and strain induced transformation (SIT) mechanisms in our work, respectively [11]. SIT contributes to the phase transition through the nucleation of martensite at the generated dislocation and shear bands [9], and the selection of martensitic variants of SIT meets Schmid factor criterion [11]. While SAT does not need the assist of new defects [12], and the selection of martensitic variants of SAT meets the strain work criterion [11]. It is noteworthy that, in real engineering materials, there are different kinds of defect inevitably or artificially, such as impurities, voids, dislocations and grain boundaries, and so forth. These preexisting defects are generally known as the original positions of phase transition nucleation or damage development, so the pre-existing defects may also affect phase transformation property of material under shock loading. Various types of defects have been

2

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

extensively studied experimentally and theoretically [13–31] to discuss their influence on the properties of materials. And most of these works were focused on the shock response, such as shock-induced void collapse and melting, dislocation emission from nanovoid, plastic deformation, the nucleation, growth and coalescence of void and vacancy emission. Certainly, the influence of various types of defects of iron also has been reported [32–38]. For example, interplay of plasticity and phase transformation in shock wave propagation in polycrystalline iron was studied [34], the results demonstrated that the phase transformation is preceded by dislocation generation at the grain boundaries (GBs), plasticity is mostly given by the formation of dislocation loops. The work of Cui et al. [35] presented the microview as well as the wave profile analyses on the phase transformation of iron with a nanovoid to understand the role of defect under shock compression. In addition, Shao et al. [37] also studied the model of single crystalline iron with a nanovoid to detect the dynamics of shockinduced structural transition. What’s more, the model of single crystalline iron with cylindrical nanopores under shock loading was studied [38]. However, up until now the micro-mechanism that how the cylindrical nanopores influences the phase transformation of iron is yet to be systematically clarified. Meanwhile, the study of such defect is meaningful because the material containing extended cylindrical pores can be used for the development of filters, detectors, cooling elements in nanoelectronics, etc. With the rapid development of computer science and technology, molecular dynamics (MD) simulation allows us to explore the influences of defects on phase transformation at the atomistic level. In this paper, we introduced a nanoscale cylindrical pore with different pore diameter into a single crystalline iron. Using the molecular dynamics method with our modified analytic embedded-atom model (MAEAM) potential, the influences of the cylindrical nanopores on the structural transformations as well as kinetics of the phase transition under shock loading were investigated.

2. Modeling and simulation The molecular dynamics code LAMMPS [39] and our MAEAM potential of iron were used in this work to perform the simulations. The potential function of iron, which provides a better description of both dislocation and phase transformation properties than previously employed potentials, has been detailed in Wang’s paper [7] and well validated in high pressure applications on the study of single crystals [7,40] as well as polycrystals [11]. To investigate the effect of crystal orientation on shock response, we first constructed three perfect single crystalline iron samples, with the Cartesian  1 0, (X, Y, Z) axes corresponding to the ([1 0 0], [0 1 0], [0 0 1]), (½1  0; ½1 1 2;  ½1 1 1) crystal directions, and their [0 0 1], [1 1 0]) and (½1 1 sizes are 17.16  17.16  167.61 nm, 16.99  16.88  192.04 nm and 16.99  16.82  195.39 nm, respectively. And then we set a cylindrical pore at the center of a perfect sample with the cylinder axis along with X axis by removing the atoms within a certain cylinder radius. In addition, in order to explore the effect of defect size on the transformation rate during phase transition of single crystalline iron with a cylindrical nanopore, four samples with different sized nanoscale cylindrical pore are studied for each shock loading direction, and their cylinder radius are about 1/10, 1/8, 1/6 and 1/5 of the length of the simulation box along Y axis, namely 3.43 nm, 4.58 nm, 5.72 nm and 6.86 nm in diameter for [0 0 1] loading direction, and 3.2 nm, 4.2 nm, 5.6 nm and 6.4 nm in diameter for both [1 1 0] and [1 1 1] loading directions. Periodical boundary conditions were applied to the X and Y directions. To equilibrate our model systems, we first applied an energy minimization using the conjugate gradient method, and then the samples were equilibrated in NPT

ensemble with temperature from 800 K to 10 K for 68 ps. Shock wave loading was applied along the Z direction by taking a few surface atomic layers at one end of the sample as a piston and pushing the piston inward at a certain velocity (vp). Different shocking pressures were generated by setting different piston velocity of 0.2, 0.3, 0.5, 0.6 and 0.8 km/s. The shock simulations were performed in an NVE ensemble. The motion equations of the atoms are integrated by the velocity Verlet algorithm with a time step of 0.4 fs and the run durations are up to about 40 ps, thus allowing for sufficient time for plasticity and transformed phase. The system was initially 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. Samples were visualized using adaptive common-neighbor analysis (ACNA) [41] within OVITO [42]. The local atomic stress is defined as [43]

" # ab ab a a 1 X pi pj a 1 X @V r i r j ab rij ð~rÞ ¼ 0 K ð~ rÞ  I ð~ rÞ 2 a–b @rab ~ r ab X a m

ð1Þ

where Ka ð~ rÞ is 1 if atom a lies within average volume X0 centered on ~ rÞ is the fraction length of bond ab r otherwise is 0. Function Iab ð~ which lies within average volume and satisfies Iab (~ r) = Iba (~ r). The wave profiles are analyzed by uniformly dividing our simulation domain into many small bins (with a thickness of about 6 angstrom each) along Z direction, and the average values of atomic stress as well as velocity calculated statistically over these bins. 3. Results and discussion 3.1. Shock wave profile Using large-scale molecular dynamics simulation, we investigated shock response of a model single crystalline iron with a cylindrical nanopore. Fig. 1 has shown the stress profiles and particle velocity profiles at five certain moments along three low index crystallographic orientations with a piston velocity of 0.5 km/s, revealing wave propagation during shock compression and subsequent release from the surface of cylindrical nanopores. Our results showed a few common features among the three directions, for example, multiple wave structure can be observed under the three loading directions. As one can see from Fig. 1(a), a reflected wave appeared after the right travelling wave arrive at the cylindrical nanopore at about 16 ps when shock along [0 0 1] direction, and the reflected wave of shock compression is a tensile wave. The flow stress suddenly fell to a lower stress states behind shock front but rose up to the full peak stress states at the cylindrical nanopore, and the particle velocity suddenly increases behind shock front due to reflected wave, which was similar to the situation of [1 1 0] and [1 1 1] loading directions. In all these three loading directions, multiple-wave-structure can be clearly observed, but is different from traditional pictures, which propose that the first wave is associated with an elastic deformation, the second is attributed to plastic deformations, and the last results from a certain type of phase transformation according to the time evolution. Taking the stress profile, particle velocity and atomic configurations at 20 ps (Fig. 1(a)) as an example, we can divide the wave profile into three parts, corresponding to three typical regimes as shown in Fig. 1(a). Region I located at the very beginning of the shock front where only elastic compression occurs. Region II was the phase transition regime. Region III is associated with elastic compression process. According to the theory of shock waves, when incident wave propagates to the free surface, the stress keeps zero. And in this work, the local stress in the stress profile gradually decreases as a result of the reflection wave which can

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

3

Fig. 1. Profiles of stress, particle velocity and the corresponding atomic configuration with a particle velocity of 0.5 km/s, where red indicates HCP sites, light blue indicates BCC sites, yellow indicates FCC sites, and dark blue denotes the other atom types commonly associated with some types of crystal defects, such as point defects and GBs, in the atomic view: (a) is the profile for the shock along [0 0 1] direction containing a cylindrical pore 6.86 nm in diameter; (b) is the profile for the shock along [1 1 0] direction containing a cylindrical pore 4.2 nm in diameter; (c) is the profile for the shock along [1 1 1] direction containing a cylindrical pore 6.4 nm in diameter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

release the local stress. Plasticity was not clearly observed under shock along [0 0 1] direction. In comparison with the shock along [0 0 1] direction, the stress profiles for the shock along [1 1 0] and [1 1 1] directions can be divided into four typical regimes. Beside the three characteristic regions as in the [0 0 1] direction, an extra region labeled as II in Fig. 1(b) and (c) appeared behind elastic precursor, and this region is associated with plastic deformations process, which behaves as the nucleation, propagation and multiplication of dislocations. In addition, the existence of the cylindrical nanopore induced ‘‘hot spot” [35] and changed the macroscopic response feature. The temperature of the region around the cylindrical pore was much higher than the other regions. As we can see from Fig. 2, for the sample shocked along [0 0 1] direction containing a cylindrical pore 6.86 nm in diameter, when the shock wave arrived at the cylindrical nanopore at about 16 ps, the temperature of the region around the edge of cylindrical nanopore began to rise, and the temperature continued to increase during the continuous compression because of the focusing of the momentum and energy toward the nanopore [14]. When the nanopore just collapsed completely at about 21.6 ps, the temperature reached the highest value of about 800 K, then the temperature decreased with increasing time. For the shock along [1 1 0] and [1 1 1] directions, the temperature of the region around the edge of the cylindrical nanopore showed similar changing tendency, what’s more, the highest temperature can reach 800 K and 1100 K for [1 1 0] and [1 1 1] loading directions, respectively. In addition, when the size of cylindrical nanopore increased, it was found that the attained highest temperature around the pore rose for the shock along different directions with the same particle velocity. Shock-induced plasticity was clearly observed in our results thanks to the superiority of our new embedded-atom-model potential, and the plasticity observed in the single crystalline iron with a cylindrical pore is different from that in a perfect single crystalline iron. We have observed that dislocations were emitted from the surfaces of cylindrical pore for [1 1 0] and [1 1 1] loading direction, however, shock-induced plasticity only can be observed

(a)

for the shock along [1 1 0] direction in the perfect single crystalline [7]. The nucleation, propagation and multiplication of dislocations under [1 1 0] loading direction was reported in perfect single crystalline iron [7], but the nucleation and growth of dislocations under [1 1 1] loading direction has not been reported in other simulation work. For the shock along [0 0 1] direction, no obvious dislocation can be observed before phase transformation. As showed in Fig. 3, after the shock front swept over the cylindrical pore, the FCC and HCP structure atoms expanded outward to form a ‘‘butterfly” phasetransition zone, which is similar to that observed by Cui et al. [35], with the time going, the butterfly shape grows up gradually. In Fig. 4, for the shock along [1 1 0] direction, a highly disordered region was created, the nucleation, propagation and multiplication  direction can be clearly of dislocations along the [1 1 1] and ½1 1 1 observed after the shock front swept over the cylindrical pore. In addition, the generated dislocations act as the nucleation sites for the continuing phase transition. For the [1 1 1] loading direction (Fig. 5), after the shock front swept over the cylindrical pore, dislocation lines were emitted from cylindrical pore surface and induced phase transformation, the dislocation along [1 1 1] direction expanded outward and act as the nucleation sites for the following phase transition. 3.2. Influence of cylindrical pore on phase transformation The effect of the cylindrical pore on the phase transformation of single crystalline iron will be discussed from the following three aspects: the critical stress of phase transformation, the nucleation sites and martensitic variants (e phases). When the piston velocity was 0.3 km/s, no phase transformation was observed in the perfect single crystals, but phase transformation starting from the surface of cylindrical nanopores appeared in all defective samples. So, the existence of nanosized cylindrical pore decreased the pressure threshold of the phase transformation, similar to the shock compression of single crystalline iron with spherical nanopore inclusion [35]. Meanwhile, the initial collapse of cylindrical pores was

800

18.4ps 20 ps 20.8ps 21.6ps 22.4ps 23.2ps 26.4ps

Temperature(K)

600

400

200

0

45

50

55

60

65

70

75

80

85

90

95

100

Z(nm)

(b)

20 ps

20.8 ps

21.6 ps

22.4 ps

Fig. 2. Profile of temperature for 0.5 km/s piston velocity simulation at different time and the corresponding atomic configurations, the sample shocked along [0 0 1] direction containing a cylindrical pore 6.86 nm in diameter. (a) The temperature of the region around the cylindrical pore at different time; (b) atomic configuration at four certain moment which corresponding to (a).

5

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

Fig. 3. Atomic configuration in the sample containing a cylindrical pore 6.86 nm in diameter, at the impact velocity of 0.5 km/s along [0 0 1] direction. Only atoms with FCC structure (yellow), HCP structure (red) and amorphous structure (dark blue) are shown, where BCC atoms have been removed before visualization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ċ

ĉ Č Ċ

Fig. 4. Perspective view of the evolution of dislocations and phase transition processes in the sample containing a cylindrical pore 6.4 nm in diameter at the impact velocity of 0.5 km/s along [1 1 0] direction. The color map is the same as Fig. 3. The nucleation, propagation and multiplication of dislocations along the  direction can be clearly observed after the shock front swept over ½1 1 1 and ½1 1 1 the cylindrical pore. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

captured at the impact velocity of 0.3 km/s for different loading directions and different sized defects. In a perfect crystal, phase transformation first appeared behind the shock front. When a strong enough shock wave swept across the perfect sample, many HCP structure atoms nucleated uniformly behind the shock front, and these nuclei grew up larger and linked together, then phase transition structure extended from the behind of the shock front. However, the phenomenon of nucleation in the defective crystal was obviously different from the perfect crystal. In our work, we have found the dislocation generating from the initial voids. The activities of dislocations can lower the barrier of nucleation and growth of the new phase, which eventually lead to a rapid increase in the quantities of the transition products. As we can see from Figs. 3–5, the nucleation sites first appeared around the surface of the cylindrical pore, and then phase transition region gradually expanded outward. For instance, the sample with a cylindrical pore 6.86 nm in diameter at the impact

Fig. 5. Perspective view of the evolution of dislocations and phase transition processes in the sample containing a cylindrical pore 6.4 nm in diameter at the impact velocity of 0.5 km/s along [1 1 1] direction. The nucleation, propagation and multiplication of dislocations along the [1 1 1] direction can be clearly observed after the shock front swept over the cylindrical pore.

velocity of 0.3 km/s along [0 0 1] direction, after 25 ps of simulation, the phase transition nucleation around the surface of the pore were observed, and the abundance of the HCP structure atoms can reach about 20%, while no phase transition nucleation was observed in the perfect crystal region. Besides, the impact strength directly affects the nucleation rate, the greater the impact strength, the faster the nucleation rate. To better reveal the effect of cylindrical pore on the martensitic variants, we mainly paid attention to the phase transformation region around the pore. Here the martensitic variants were distinguished by analyzing the orientation of the base plane of HCP phase. Furthermore, the base plane of the HCP phase can be accurately determined by the c axis analysis technique and the distribution of the c axis can be illustrated in a pole figure [7,11]. The c axis analysis technique and pole-figure technique were suitably applied to the analysis of variants [7,11], and the details about c axis analysis technique and pole-figure technique could refer to the work of Wang et al. [7]. By choosing several typical cross sections near the cylindrical nanopore with respect to the loading direction behind the shock front and applying c axis analysis technique, the orientation of c axis of the transition products along three loading directions at vp = 0.5 km/s are showed in Table 1 and Fig. 6. In conjunction with Table 1, Figs. 6 and 7, we can get  0, ½1 1 0, ½1 1  0, ½1  1 0 and that c axis are mainly located in: ½1 1  0 and ½1 0; 1 for the shock along [0 0 1] direction; ½1 1 0, ½1 1  1 0 for the shock along [1 1 0] direction; [1 0 1], [1 1 0], [011], ½1  1 0 and ½1 1  0 for the shock along [1 1 1] direction. In other ½1

Table 1 c axis orientations after phase transition for the shock along three directions with vp = 0.5 km/s. The angle in parentheses denotes the equivalent value. Shock direction

c axis’ projection orientation (h)

Deviation of shocking direction (u0 )

Original BCC phase orientation

[0 0 1]

180 45 135

135 90 90

[1 0 1] [1 1 0]  0 ½1 1

[1 1 0]

0 (180 ) –

90 0 (180 )

 1 0 ½1 [1 1 0]

[1 1 1]

0 (180 ) 90 30 155

90 35 145 145

 0 ½1 1 [1 1 0] [0 1 1] [1 0 1]

6

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

Fig. 6. Three groups of cross sections transverse to the three loading directions are cut from three samples at the impact velocity of 0.5 km/s to judge the orientation of c axis: (a) shows three cross sections which are corresponding to the same position (Z = 70 nm) of the sample shocked along [0 0 1] direction at 24 ps, the cylindrical pore diameter is 6.86 nm. The first picture shows the cross section of the atoms configuration, color coding as in Fig. 1, the rest of the two pictures are colored by the angle h and u0 , respectively; (b) shows two cross sections which are corresponding to the same position (Z = 88 nm) of the sample with a cylindrical pore 5.6 nm in diameter shocked along [1 1 0] direction at 24 ps. The first picture shows the cross section of the atoms configuration, the second picture is colored by the angle u0 , and the attached h are equal to 0° (180°); (c) shows three cross sections which are corresponding to the same position (Z = 91 nm) of the sample with a cylindrical pore 5.6 nm in diameter shocked along [1 1 1] direction at 24 ps. The first picture shows the cross section of the atoms configuration, the rest of the two pictures are colored by the angle h and u0 , respectively; more details on this method can refer to the work of Wang about c axis analysis technique [7]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

words, the main martensitic variants (shuffle planes) are (1 1 0),  0Þ (1 0 1) variants for [0 0 1] loading direction; (1 1 0) and ð1 1  0Þ variants for [1 1 0] loading direction; (1 0 1), (1 1 0), (0 1 1) ð1 1  0Þ variants for [1 1 1] loading direction. In consequence, and ð1 1 the final orientations of the c axis along all three low index directions are determined to be h1 1 0iBCC, which is consistent with the transition mechanisms of iron.

In our previous work [11], we have discussed SAT and SIT mechanisms for the shock upon polycrystalline iron in detail. For the SAT, the plastic deformation of compression along [0 0 1] direction will accumulate strain energy, thus the lattice of parent phase becomes instable, and then shuffle processes will take place with the sufficient accumulative strain energy. This process of phase transformation corresponds to the first kind of phase transition

7

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

Fig. 7. The pole figure of the c axis within the sample reference frame for the shock along the three directions. Equal area projection, contours in the number of HCP atoms with the corresponding c axis orientation. The positions which have brighter colors in each pole figure represent the main orientations of the c axis. Each point of the pole figure, corresponding to a certain c axis orientation, is colored by the number of HCP atoms with the same c axis orientation. As we can see from the first picture, for the shock  0, ½1 1 0, ½1 1  0, ½1  1 0 and ½1 0 1; the second picture shows that the c axis orientation are ½1 1 0, ½1 1  0 and ½1  1 0 for the shock along [0 0 1] direction, c axis orientation are ½1 1  1 0 and ½1 1  0 for the shock along [1 1 1] direction. More details on this along [1 1 0] direction; The third picture shows that the c axis orientation are [1 0 1], [1 1 0], [0 1 1], ½1 method can refer to the work of Wang et al. [7]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

mechanism of iron as mentioned in the introduction. The SAT mechanism could occur not only in polycrystalline iron but also in perfect single crystalline iron, where the plasticity only partially couples with the phase transition processes, the plasticity couples with the phase transition through the compression along [0 0 1]. While for the SIT, the large local stress concentration will result  1 0 slip system, then a nucleus of in the activation of ð1 0 1Þ½0 HCP-like defects will be formed. With continuing loading, the  1 0 slip sysHCP-like defects grow with the slippages of ð1 0 1Þ½0 tem. Because of the energetic preference, the HCP-like defects could relax to form a stable HCP phase and HCP phase will grow  plane with continuing slipping. The SIT slice by slice along ð1 0 1Þ could take place in polycrystalline iron but not be observed in perfect single crystalline iron, where the plasticity completely couples with the phase transition processes, the plasticity couples with the phase transition through the compression along [0 1 0] and the  1 0 slip system to form HCP-like defects. activation of ð1 0 1Þ½0 Wang et al. [11] verified that martensitic variants generated through the SAT could be well predicted by maximum strain work criterion, while the selection of martensitic variants generated through the SIT mechanisms is governed by minimum Schmid factor criterion. By comparing the martensitic variants of the perfect single crystalline iron [7] and the defective iron containing a cylindrical nanopore under shock loading, the existence of nanosized cylindrical pore increased a martensitic variants for the three loading direction, namely ð1 0 1Þ variants for the shock along [0 0 1] direction,  0Þ variants for the shock along the [1 1 0] and [1 1 1] direction. ð1 1 In addition, the variants, which caused by the presence of cylindrical pores, mainly distributed in the region around the cylindrical nanopore with time evolution, while for the region far away from the cylindrical nanopore, the emerged variants are the same as those emerged in single crystalline iron. In order to study the variants, which caused by the presence of cylindrical pore, the variants selection rules will be discussed in terms of strain works and Schmid factor criterions. The calculation about the strain works and Schmid factors are the same as our previous work [11], except that the stress is not uniaxial stress along z direction but actual local stress around cylindrical nanopore. Due to the presence of the cylindrical nanopore, the local stress around the cylindrical nanopore is complex and largely weakens the leading role of uniaxial stress on the variants selection. For the variants around the cylindrical nanopore, it is actual local stress rather than the uniaxial stress along z direction that play the leading role on the variants selection. Here we mainly focus on the region near the cylindrical

nanopore. Due to the specific location where variants located in, we have chosen local stress located in different location. For the  0Þ variants mainly disshock along [0 0 1] direction, (1 1 0) and ð1 1 tributed in the upper and lower of nanopores (region I and II in Fig. 2(b)), and (1 0 1) variants mainly distributed in the left and right of nanopores (region III and IV in Fig. 2(b)). So, we take the actual local stress which locate in the up and left of the nanopore to calculate strain work and Schmid factor for each {1 1 0}h1 1 1i slip system. Similarly, for the shock along [1 1 0] direction, (1 1 0) variants mainly distributed in the upper and lower of nanopores  0Þ variants mainly distributed (region I and II in Fig. 4), and ð1 1 in the left and right of nanopores (region III and IV in Fig. 4); For  0Þ the shock along [1 1 1] direction, (1 1 0), (1 0 1), (0 1 1) and ð1 1 variants are mixed together. The strain works of all possible transformations are listed in Table 2. Additionally, the corresponding Schmid factors of {1 1 0}h1 1 1i slip systems are calculated as shown in Table 3. According to Tables 2 and 3, for the shock along [0 0 1] direction, based on the data of [0 0 1]up column, we can get that (1 1 0) variants meets minimum Schmid factors criterion. Namely, (1 1 0) vari 0Þ variants are generated through the ants and its equivalent ð1 1 SIT mechanisms. Based on the data of [0 0 1]left column, we can get that (1 0 1) variants meets maximum Schmid factors criterion. Because the local stress in the left is tensile stress instead of compressive stress, the (1 0 1) variants does not meet minimum Schmid

Table 2 Strain work done by the first step of the corresponding phase transition mechanisms. The top line of the table denotes different single crystalline iron with cylindrical nanopore samples, whose shock directions lie along their [0 0 1], [1 1 0] and [1 1 1], respectively. [0 0 1]up and [0 0 1]left denotes that the local stress we choose locates in the upper and left of the nanopore along [0 0 1] loading direction, respectively, and the notes is the same with [1 1 0] and [1 1 1] loading direction. The sample shocked along [0 0 1] direction containing a cylindrical pore 6.86 nm in diameter; the sample shocked along [1 1 0] direction containing a cylindrical pore 5.6 nm in diameter; the sample shocked along [1 1 1] direction containing a cylindrical pore 5.6 nm in diameter. The leftmost column denotes different martensitic variants represented by the basal plane of e phase or the corresponding shuffle plane of a phase. If not specified, the same conventions would hold for Table 3. Variants

[0 0 1]up

[0 0 1]left

[1 1 0]up

[1 1 0]left

[1 1 1]up

[1 1 1]left

(1 1 0) (1 0 1) (0 1 1) (1 1 0) (1 0 1) (0 1 1)

11.3418 11.6993 11.5958 11.3079 11.6542 11.7993

0.0483 0.3040 0.3026 0.3129 0.1616 0.0066

6.5893 5.9587 5.1503 7.8783 5.1545 5.9550

1.5010 1.5949 1.9681 0.8827 2.0153 1.6404

5.3990 6.5935 6.2770 11.8578 8.1910 9.5397

3.8164 4.1888 2.8866 6.3352 4.1395 4.9645

8

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

Table 3 Schmid factor for each {1 1 0}h1 1 1i slip system under actual local stress around cylindrical nanopore. Variants

Slip direction

[0 0 1]up

[0 0 1]left

[1 1 0]up

[1 1 0]left

[1 1 1]up

[1 1 1]left

(1 1 0)

[1 1 1] [1 1 1]

0.0049 0.0037

0.3019 0.3088

0.0016 0.0017

0.1557 0.0734

0.2906 0.2857

0.2893 0.2015

(1 0 1)

[1 1 1] [1 1 1]

0.0131 0.0279

0.4155 0.4580

0.4662 0.2130

0.0088 0.4548

0.8048 0.2284

0.7846 0.2968

(0 1 1)

[1 1 1] [1 1 1]

0.0149 0.0263

0.3531 0.0397

0.2020 0.4758

0.4908 0.0774

0.2334 0.8052

0.3536 0.7536

(1 1 0)

[1 1 1] [1 1 1]

0.0142 0.0130

0.1570 0.1061

0.0160 0.0097

0.1003 0.0502

0.0060 0.0044

0.1384 0.2413

(1 0 1)

[1 1 1] [1 1 1]

0.0230 0.0198

0.3243 0.2140

0.4758 0.2033

0.0952 0.4943

0.8096 0.3279

0.7296 0.1093

(0 1 1)

[1 1 1] [1 1 1]

0.0086 0.0316

0.4533 0.3384

0.4662 0.2117

0.0517 0.4470

0.8099 0.3248

0.7229 0.0043

factors criterion. In a word, for the shock along [0 0 1] direction,  0Þ and (1 0 1) variants are generated through the SIT (1 1 0), ð1 1 mechanisms. For the shock along [1 1 0] direction, based on the data of [1 1 0]up column, we can get that (1 1 0) variants meets minimum Schmid factors criterion and (1 1 0) variants is generated through the SIT mechanisms. Based on the data of [1 1 0]left column,  0Þ variants meets minimum strain work critewe can get that ð1 1 rion, because the local stress is tensile stress instead of compres 0Þ variants does not meet maximum Schmid sive stress, ð1 1 factors criterion. However, for the shock along [1 1 1] direction, the situation is different from that for the shock along [0 0 1] and [1 1 0] direction. For the shock along [1 1 1] direction, the local stress at the upper and left of the nanopore are all compressive stress, and the highest temperature is much higher than that of the other two loading directions, and the emerged martensitic variants are mixed together. So, the phase transition mechanisms for the shock along [1 1 1] direction is more complex, and it is dif 0Þ variants is generated by SAT or SIT in ficult to determine the ð1 1 our present simulation.

Fig. 8. The number of variants’ atoms as a function of shock wave propagation time for the shock along [0 0 1] direction, where (a) is corresponding to the new variants which is (1 0 1) variants, while (b) is corresponding to the old variants which are  0Þ variants. (1 1 0) and ð1 1

3.3. Kinetics of the phase transition According to the above discussion, the introduction of cylindrical nanopores can influence martensitic variants for different loading direction. In order to study the effect of cylindrical pore on kinetics of the phase transition, we divided the martensitic variants which emerged under the shock into two kinds, one includes the same variants as those occurred in single crystalline iron and the other includes what only emerged in this work. To avoid ambiguity, the first and the later one are called ‘‘old variants” and ‘‘new variants”,  0Þ respectively. For the shock along [0 0 1] direction, (1 1 0) and ð1 1 variants are old variants, and (1 0 1) variants is a new variant. For the shock along [1 1 0] direction, (1 1 0) variants is old variants, and  0Þ variants is a new variants. For the shock along [1 1 1] direction, ð1 1  0Þ variants (1 1 0), (1 0 1) and (0 1 1) variants are old variants, and ð1 1 is new variants. The number of variants’ atoms under the shock along [0 0 1], [1 1 0] and [1 1 1] directions of single crystalline iron containing different size of cylindrical pores have been analyzed, respectively, and the results are shown in Figs. 8–10. The number of the variants caused by the cylindrical pores is far smaller than that located in the region far away from the cylindrical pores. So, the number of the variants, which caused by the cylindrical pore but is the same as that happened in single crystalline iron, have no significant impact on the number of old variants. For the shock along [0 0 1] direction, the number of old variants’ atoms is affected by the cylindrical pores from ‘B’ to ‘E’ as shown in Fig. 8(b). The

Fig. 9. The number of variants’ atoms as a function of shock wave propagation time for the shock along [1 1 0] direction, where (a) is corresponding to the new variants  0Þ variants, while (b) is corresponding to the old variants which is which is ð1 1 (1 1 0) variants.

existence of cylindrical pores could largely slow down the increasing rate of the old variants and promote the nucleation of the new variants. However, the new variants near the cylindrical pores

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

Fig. 10. The number of variants’ atoms as a function of shock wave propagation time for the shock along [1 1 1] direction, where (a) is corresponding to the new  0Þ variants, while (b) is corresponding to the old variants variants which is ð1 1 which are (1 1 0), (1 0 1) and (0 1 1) variants.

cannot grow bigger after shock waves passed by until the pore size increases to 6.86 nm. If the pore size is larger this critical size, the new variants would be able to grow up after shock waves passed by. The reason for the slow increasing rate of the old variants is due to reflections of shock waves from the surface of cylindrical pores as shown in Fig. 1, where high stress states, caused by the shock compressions, are released by the reflections of waves which could lead to an e ? a reverse transition. For the cases of [1 1 0] and [1 1 1] directions, similar phenomena could also be observed. As shown in Figs. 9 and 10, the critical pore size is about 5.6 nm and 4.2 nm for the shock along [1 1 0] and [1 1 1] directions, respectively. Below the critical size, the increasing rate of the old variants is decreasing with the growing pore size. While above the critical size, its increasing rate rises with the growing pore size because of the interactions between the old variants and dislocations, as well as the new variants whose growing rate increases with the promoting pore size. Therefore, cylindrical pores have two basic effects on the kinetics of the phase transition: One is slowing down the increase rate of the old variants while promoting the nucleation of the new variants. The other is that there is a critical pore size which determines whether the nucleus of the new variants could grow up, which will give rise to a distinctly different kinetics of the phase transition. 4. Conclusions Shock-induced phase transformation in single crystalline iron with cylindrical nanopores was investigated using molecular dynamics simulations with our modified analytic embeddedatom model potential, and a variety of shock strength with piston velocity ranging from 0.2 to 0.8 km/s has been considered on shock along three low index crystallographic directions. The main findings were listed as follows: (1) The cylindrical nanopore influences the shock response. Reflection wave can be observed when shock wave swept over the cylindrical nanopore. In addition, the temperature of the region around the cylindrical pore was much higher than the other regions, what’s more, when the size of cylindrical nanopore increased, the attained highest temperature around the pore rose. Meanwhile, the nucleation and growth of dislocations can be clearly observed under [1 1 0] and [1 1 1] loading directions, and the surfaces of the cylindrical nanopore act as effective sources for dislocations.

9

(2) The cylindrical nanopore influences the phase transformation. The existence of cylindrical pore decreases the critical stress of phase transformation, and the cylindrical pore is the favorite nucleation site of the phase transformation of iron for different loading directions. In addition, the introduction of the cylindrical nanopore increases martensitic variants of ð1 0 1Þ for the shock along [0 0 1] direction,  0Þ for the shock along the [1 1 0] and [1 1 1] direction. ð1 1 (3) SAT and SIT were discussed in terms of the variants near the cylindrical nanopores. For the shock along [0 0 1] direction,  0Þ and (1 0 1) variants were generated through (1 1 0), ð1 1 the SIT mechanisms. For the shock along [1 1 0] direction, (1 1 0) variants was generated through the SIT mechanisms  0Þ variants was generated through the SAT mechaand ð1 1 nisms. For the shock along [1 1 1] direction, it was difficult to determine whether SAT or SIT contribute to the genera 0Þ variants. tion of the ð1 1 (4) The number of variants’ atoms were investigated. Cylindrical pores can slow down the increase rate of the old variants while promote the nucleation of the new variants. In addition, there is a critical pore size which determines whether the nucleus of the new variants could grow up, which will give rise to a distinctly different kinetics of the phase transition.

Acknowledgements This work is supported by the National Natural Science Foundation of China (NSFC-NSAF U1530151, NSFC 51571088 and NSFCNSAF 11076012), the Fundamental Research Funds for the Central Universities, National Key Laboratory Project of Shock Wave and Detonation Physics, the Science and Technology Foundation of National Key Laboratory of Shock Wave and Detonation Physics and Chinese National Fusion Project for ITER with Grant No. 2013GB114001. The authors would like to acknowledge the support of the computation platform of National Supper-Computer Center in Changsha (NSCC). References [1] D. Bancroft, E.L. Peterson, S. Minshall, J. Appl. Phys. 27 (1956) 291. [2] N.V. Barge, R. Boehler, High Press. Res. 6 (1990) 133–140. [3] R.F. Smith, C.A. Bolme, D.J. Erskine, P.M. Celliers, S. Ali, J.H. Eggert, S.L. Brygoo, B.D. Hammel, J. Wang, G.W. Collins, J. Appl. Phys. 114 (2013) 133504. [4] J.C. Crowhurst, B.W. Reed, M.R. Armstrong, H.B. Radousky, J.A. Carter, D.C. Swift, J.M. Zaug, R.W. Minich, N.E. Teslich, M. Kumar, J. Appl. Phys. 115 (2014) 113506. [5] K. Caspersen, A. Lew, M. Ortiz, E. Carter, Phys. Rev. Lett. 93 (2004) 115501. [6] W.W. Pang, P. Zhang, G.C. Zhang, A.G. Xu, X.G. Zhao, Sci. Rep. 4 (2014) 5273. [7] K. Wang, S.F. Xiao, H.Q. Deng, W.J. Zhu, W.Y. Hu, Int. J. Plast. 59 (2014) 180– 198. [8] W.G. Burgers, Physica 1 (1934) 561–586. [9] J. Talonen, H. Hänninen, Acta Mater. 55 (2007) 6108–6118. [10] J. Venables, Philos. Mag. 7 (1962) 35–44. [11] K. Wang, W. Zhu, S. Xiao, K. Chen, H. Deng, W. Hu, Int. J. Plast. 71 (2015) 218– 236. [12] A. Das, P. Chakraborti, S. Tarafder, Mater. Sci. Technol. 27 (2011) 366–370. [13] R.E. Rudd, Comput. Mater. Sci. 24 (2002) 148–153. [14] T. Hatano, Phys. Rev. Lett. 92 (2004) 015503. [15] S. Traiviratana, E.M. Bringa, D.J. Benson, M.A. Meyers, Acta Mater. 56 (2008) 3874–3886. [16] K.J. Zhao, C.Q. Chen, Y.P. Shen, T.J. Lu, Comput. Mater. Sci. 46 (2009) 749–754. [17] S.Z. Xu, Z.M. Hao, Y.Q. Su, Y. Yu, Q. Wan, W.J. Hu, Comput. Mater. Sci. 50 (2011) 2411–2421. [18] Y. Tang, E.M. Bringa, B.A. Remington, M.A. Meyers, Acta Mater. 59 (2011) 1354–1372. [19] A.M. He, S. Duan, J.L. Shao, P. Wang, C. Qin, J. Appl. Phys. 112 (2012) 074116. [20] Y. Tang, E.M. Bringa, M.A. Meyers, Acta Mater. 60 (2012) 4856–4865. [21] H.J. Chang, J. Segurado, O. Rodríguez de la Fuente, B.M. Pabón, J. Llorca, Modell. Simul. Mater. Sci. Eng. 21 (2013) 075010. [22] S.J. Fensin, S.M. Valone, E.K. Cerreta, J.P. Escobedo-Diaz, G.T. Gray, K. Kang, J. Wang, Modell. Simul. Mater. Sci. Eng. 21 (2013) 015011.

10

L. Wu et al. / Computational Materials Science 122 (2016) 1–10

[23] C.J. Ruestes, E.M. Bringa, A. Stukowski, J.F. Rodríguez Nieva, Y. Tang, M.A. Meyers, Comput. Mater. Sci. 88 (2014) 92–102. [24] J.L. Shao, P. Wang, A.M. He, S.Q. Duan, C.S. Qin, Modell. Simul. Mater. Sci. Eng. 22 (2014) 025012. [25] Y. Zhao, Q. Fang, Y. Liu, X. Zeng, Int. J. Solids Struct. 51 (2014) 1617–1629. [26] D. Tramontina, C. Ruestes, Y. Tang, E. Bringa, Comput. Mater. Sci. 90 (2014) 82– 88. [27] F.P. Zhao, B. Li, W.R. Jian, L. Wang, S.N. Luo, J. Appl. Phys. 118 (2015) 035904. [28] J.F. Rodriguez-Nieva, C.J. Ruestes, Y. Tang, E.M. Bringa, Acta Mater. 80 (2014) 67–76. [29] C.J. Ruestes, E.M. Bringa, A. Stukowski, J.F. Rodríguez Nieva, G. Bertolino, Y. Tang, M.A. Meyers, Scr. Mater. 68 (2013) 817–820. [30] E.T. Seppälä, J. Belak, R.E. Rudd, Phys. Rev. B 69 (2004) 134101. [31] L.P. Dávila, P. Erhart, E.M. Bringa, M.A. Meyers, V.A. Lubarda, M.S. Schneider, R. Becker, M. Kumar, Appl. Phys. Lett. 86 (2005) 161902. [32] F. Gao, D.J. Bacon, P.E.J. Flewitt, T.A. Lewis, Lewis J. Nucl. Mater. 249 (1997) 77–86.

[33] M. Itakur, H. Kaburaki, K. Kusunoki, Nucl. Instrum. Methods Phys. Res. B 153 (1999) 122–125. [34] N. Gunkelmann, D.R. Tramontina, E.M. Bringa, H.M. Urbassek, New J. Phys. 16 (2014) 093032. [35] X.L. Cui, W.J. Zhu, H.L. He, X. Deng, Y.J. Li, Phys. Rev. B 78 (2008) 024115. [36] K. Kadau, T.C. Germann, P.S. Lomdahl, R.C. Albers, J.S. Wark, A. Higginbotham, B.L. Holian, M. Elert, M.D. Furnish, R. Chau, N. Holmes, J. Nguyen, Shock Compression Condens. Matter 955 (2007) 313–316. [37] J.L. Shao, S.Q. Duan, A.M. He, P. Wang, C.S. Qin, J. Phys.: Condens. Matter 22 (2010) 355403. [38] B. Ma, Appl. Mech. Mater. 543–547 (2014) 3910–3913. [39] S. Plimpton, Comput. Mater. Sci. 4 (1995) 361–364. [40] K. Wang, S.F. Xiao, M. Liu, H.Q. Deng, W.J. Zhu, W.Y. Hu, Procedia Eng. 61 (2013) 122–129. [41] A. Stukowski, Modell. Simul. Mater. Sci. Eng. 20 (2012) 045021. [42] A. Stukowski, Modell. Simul. Mater. Sci. Eng. 18 (2010) 015012. [43] J. Cormier, J.M. Rickman, T.J. Delph, J. Appl. Phys. 89 (2001) 99–104.