On the deformation-induced grain rotations in gradient nano-grained copper based on molecular dynamics simulations

In this paper, the mechanical behavior of gradient nano-grained copper under uniaxial deformation was investigated using molecular dynamics simulations. The stress response was found to be different in the regions with different grain sizes, which was attributed to the different dislocation activities due to the dislocation-grain boundary synergies. The phenomenon of grain rotation was observed and a program was developed to accurately evaluate the grain rotation and explore its dependence on the grain size and the initial crystal orientation. It is found that all grains tend to rotate to the 30° orientation, consistent with the activation theory of the slip systems under the uniaxial deformation. The rotation magnitude is larger for larger grains, but the rotation rate is more diversely distributed for smaller grains, indicating more disturbance from grain boundary mechanisms such as the grain boundary sliding and the grain boundary diffusion for smaller grains. The effect of temperature on the grain rotation is also investigated, showing an increase of the dispersion of grain rotation distribution with the increase of temperature. This paper aims at providing insights into the synergistic deformation mechanisms from dislocations and grain boundaries accounting for the exceptional ductility of the gradient nano-grained metals.


Introduction
Grain rotation, referred to as the change of lattice orientation of grains, has received increasing attention since the ignition of research interests in nanocrystalline metals in 1990s. Although grain rotation is not treated as an important mechanism in coarse-grained metals, it is believed to play an important role in the plastic deformation of nanocrystalline metals in which the volume percentage of grain boundaries upsurges as the grain size decreases to nanometers [1][2][3]. Indeed, the phenomenon of grain rotation has been reported from a significant number of experiments and simulations in the nanograined films and bulk materials. From the in situ tensile test of Pt thin films under a high-resolution TEM, Wang et al. [4] directly showed that when the grain size is below 6 nm, the plastic deformation mechanism transits from dislocation glide to grain rotation coordinated in multiple grains mediated by the climb of grain boundary (GB) dislocations. Li et al. [5] investigated the softening of nanocrystalline materials considering the grain rotationbased nonhomogeneous plastic deformation using a phase mixture model, in which the nanocrystalline materials were treated as composites consisting of grain interior and GB phases. Experiments and simulations [6][7][8] have indicated that as the grain size decreases from micrometer to a few nanometers, the dominating mechanism of plastic deformation evolves from full dislocation activity in large grains, to partial dislocation activity in smaller grains, and finally to GB-mediated mechanisms including grain creep, grain sliding, and grain rotation. It indicates that the Hall-Petch relation [1,9] that prevails for the conventional metal materials can no longer be applied for predicting the strength of the nanocrystalline metals [10]. The size dependence has been found for grain rotation. Some authors showed that grain rotation is much more evident in finer nanocrystals [11,12], while others reported a reversal size dependence of grain rotation [13]. Theoretical models were also established to calculate the grain rotation rate considering the mass diffusion and GB sliding, e.g., the efforts by Kim et al. [14,15], Harris et al. [16], and Wheeler [17].
Grain rotation is usually accompanied with other mechanical or thermal processes such as cracking, annealing, and grain growth. Cracking is a special mechanical scenario in which the grain rotation mechanism may take effect. It has been found that the deformation at the crack tip is accommodated by grain rotation in nanocrystalline metals. A number of in situ experiments [18] and theoretical modeling efforts [19,20] have been devoted to investigate the effect of grain rotation on the crack growth and blunting. For example, some authors (Li et al. [21,22], Zhang and Zhou [23]) showed that the nanograin rotation can significantly enhance the dislocation emission from the crack tip, leading to evident crack blunting in nanomaterials. Liu et al. [19] demonstrated that grain rotationinduced deformation mechanism can lead to an increase of the critical crack intensity factor. Grain rotation is also involved in the process of grain growth. The synergies among the grain rotation, the GB migration, the dislocation activities, and the grain growth have also been explored in a number of researches through the theoretical approach [20,24,25], the in situ TEM testing [26,27], the phase field simulations [28,29], etc. Several interesting phenomena and mechanisms have been reported. For example, Liu et al. [30] demonstrated that the cooperation between the nano-grain rotation and the GB migration could lead to the enhancement of the dislocation emission. Moldovan et al. [31,32] demonstrated that the grain coalescence induced by grain rotation leads to the power-law rate of grain growth with a universal scaling exponent. Farkas et al. [33] reported a linear kinetics for grain growth accompanied by grain rotation in the nanocrystalline Ni, which is different from the square root kinetics observed in the coarsegrained counterparts. Vuppuluri and Vedantam [34] deduced the average grain growth rate as R ∼ t 1/3 under the coupled mechanisms of GB migration and grain rotation coalescence. Besides cracking and grain growth, grain rotation can also be mediated by other mechanical or thermal processes. For example, Danilenko et al. [35] observed annealing-induced grain rotation in ultrafinegrained aluminum alloy from in situ TEM studies. Wang et al. [36] investigated the synergistic effects of rolling and annealing on the grain rotation and the resulted texture in nickel.
Most of the existing researches relating to grain rotation were carried out for the conventional metals with uniform grain size distribution. It is widely known that the conventional metals have long faced a dilemma of the strength-ductility tradeoff, i.e., an increase in strength leading to a decrease in ductility and vice versa [1,37,38]. Recently, a novel approach for overcoming this dilemma was proposed by Lu [39] and other researchers [40,41] who introduced gradient into the grain size distribution, referred to as gradient nano-grained (GNG) metals. The GNG metals have nonuniformly distributed grain size that increases gradually from a few nanometers near the surface to several micrometers in the interior. They possess several exceptional properties such as high strength, good strain hardening ability, satisfactory ductility, and excellent erosion resistance and can be economically produced through mechanical surface treatments such as mechanical surface attrition, frictional sliding, and high energy shot peening [42][43][44]. These advantages endow the GNG metals with a broad application potential and have attracted sustained research interests in investigating the underlined mechanisms behind their intriguing mechanical behaviors [45,46].
As reviewed by Chen et al. [47], grain rotation has been extensively studied through theoretical models, molecular dynamics simulations, and experimental investigations. However, few efforts have been devoted to grain rotation in the GNG metals. In our previous work [48], the grain rotation behavior during the uniaxial tension of nanocrystalline copper was investigated using the crystal plasticity finite element method (CPFEM). It was found that the grain rotation rate does not show any dependence on grain size. Noting that the CPFEM framework is based on the dislocation slip systems and does not incorporate the GB mechanisms, the CPFEM results may deviate from the reality. On the other hand, atomistic-based method such as molecular dynamics (MD) is a capable tool that can capture various mechanisms in the microstructure evolution of nanocrystalline metals [49][50][51]. In this paper, we carried out a thorough MD study to investigate the grain rotation behavior in gradient nanocrystalline copper, which is a widely investigated material in the literature [52][53][54]. This work is aimed at promoting the understanding of the synergistic deformation mechanisms of the nonuniform nano-grained materials.

Molecular dynamics modelling of GNG copper under uniaxial deformation
In real practice, the grain size of the GNG metals usually changes from several nanometers on the surface to hundreds of nanometers at a depth of 100 μm. Considering the forbidden computational cost in MD that a cube of just 1 μm size contains 80 billion atoms, it is not feasible to simulate a real GNG sample whose length is around 100 μm. Therefore, we adopted an economic scheme to save the computation cost by accelerating the gradient so that the grain size changes from several nanometers to tens of nanometers across a length of less than one micrometer. Besides, the sample possesses a quasi-threedimensional configuration in which the thickness along the out-of-plane direction is very small and the shape of grains is actually columnar instead of equiaxed with the axes along the out-of-plane direction. Figure 1(a) shows one of the computational samples of GNG copper. The computational samples were created using the homemade codes. Firstly, a MATLAB program was developed to generate the gradient polycrystal topology using the Voronoi method. To incorporate the gradience in the grain size, the positions of the grain center seeds were adjusted according to a linear varied distribution along the depth direction. After the generation of the GNG topology, the polycrystals were filled with copper atoms arranged in face-centered cubic (FCC) lattices whose orientations vary from one grain to another. A FORTRAN program was developed to fill in the atoms. In this work, all the grains were orientated with their [111] crystal direction along the out-of-plane direction. The purpose of choosing the [111] orientation is that in that orientation, the three slip systems that could be possibly activated are ( ) , whose slip directions are all within the plane (111) which is also the specimen plane. The resulted grain rotation due the activation of these slip systems is always around the [111] axis, causing the [111] crystal direction kept along the out-of-plane direction. In this way, the rotation angle can be easily determined from the relative positions of atoms on the close-packed (111) plane. The initial orientation of the grains followed a uniform random distribution, resulting in GBs with randomly distributed misorientation angles. The dimension of the MD sample is about 1,000 × 2,000 × 18.75 Å 3 . In this work, the gradient of the model is represented by the combination of the minimum grain size and the maximum grain size, e.g., 3-25 representing the smallest grain size of 3 nm and the largest of 25 nm. Three MD samples with gradients of 3-50, 3-25, and 3-10 were created. The smallest grain size was 3 nm and the largest grain size was 50 nm. The reason for choosing this grain size range is that the smallest grain size of 3 nm is believed to be the lower bound for nanocrystalline metals and the largest grain size of 50 nm corresponds to the most computational cost that we can afford.
The simulations of uniaxial tensile deformation were carried out using the open source code LAMMPS. The embedded atom potential (EAM) was used to model the interaction between atoms. The period boundary condition was applied to all the sample boundaries except for the boundaries along the gradient direction. The dynamic equations were explicitly integrated using the velocity Verlet algorithm with a time step of 1 fs. The pressure along the out-of-plane direction was kept at 1 atm using the isothermal-isobaric (NPT) ensemble, while the pressure along the gradient direction was left uncontrolled at self-equilibrated state. Before loading application, the samples were relaxed at the temperature of 300 K for about 0.5 ns to iron out any possible energy-unfavorable defects. Then a uniaxial affine deformation was applied to the direction perpendicular to the grain size gradient at a constant strain rate of 0.5 × 10 9 s −1 with the temperature controlled at 300 K. This strain rate was chosen as a result of the balance of the computational cost and the computation fidelity. As a common limitation, molecular dynamic simulations suffer from the tiny integration time step and the resulting extremely high strain rate. Decreasing the strain rate would lead to a significant increase of the computational cost. A case simulation for the 50 nm grain-sized sample straining to 0.4 under the rate of 0.5 × 10 9 s −1 took 43.88 h using 256 nodes of Shanghai Supercomputer Center clusters (Intel Skylake Xeon, 2.6 GHz), which was close to the limit of computational cost we can afford. Although the adopted strain rate is still several orders of magnitude higher than the experimental strain rate of usually within 10 −3 to 10 2 , the simulations are believed to be able to reflect the dislocation-dominated grain rotation mechanisms. Besides, this strain rate was also adopted in the literature [10]. The maximum strain is unanimously set as 0.4 for all the samples so that adequate significant plastic deformation could take place.
To investigate the effect of grain size on the deformation behavior of the nanocrystalline copper, a series of uniform grained samples with different grain sizes were created (see Figure 2) and deformed under the same loading conditions. The results of the uniform grained samples were compared with those of the GNG sample to provide insights into the synergistic deformation mechanisms in GNG metals. Figure 3 compares the stress versus strain curves for the uniform grained samples with different grain sizes. It shows that the Young's modulus is smaller for the smaller grain size of 3 nm. This can be attributed to the larger fraction of GB region for the 3 nm sample. As reported by the literature [55,56], the GB region is much less dense and thus has a smaller modulus compared with the grain interior region. Therefore, a smaller grain size corresponds to a smaller Young's modulus. From Figure 3, the sample with smaller grain size also has a lower peak stress, while the tendency is opposite in terms of the plateau flow stress which is larger for the smaller grain size. In general, the stress-strain curve is less undulant for the smaller grain size. On the other hand, the strain hardening effect is more obvious for the larger grain size as can be seen by the steeper slope of the fitted lines of the plastic section of the stress-strain curve (see the dash lines fitted for the strain range 0.08-0.4). This stronger strain hardening can be attributed to the more evident dislocation activity in the larger grain size, as shown later in Section 3.2. The material parameters obtained from the stress-strain curves are listed in Table 1. The strain hardening modulus of the 50 nm grain size is three times as large as that of the 3 nm grain size.   Figure 4 shows the snapshots of the samples during the uniaxial deformation at strains 0.1, 0.2, and 0.4. In these graphs, the green lines indicate the Shockley partial dislocations, the red lines indicate the stair-rod dislocations, the atoms-colored blue are at the GBs, and all other atoms are not displayed for clarity. To make the dislocation lines more clearly seen, only part of the samples are shown. Evident dislocation activities can be seen from all these snapshots. However, some differences still exist between the small grain (3 nm) sample and the large grain (50 nm) sample in terms of the dislocation-GB synergy and the grain coalescence/partition. In the small grain (3 nm) sample, the dislocations are relatively nonuniformly distributed, i.e., the dislocations are abundant in some grains, while scarce in others. On the other hand, abundant dislocation activities can be seen in nearly all the grains in the large grain (50 nm) sample. In the small grain sample, the complete extended dislocations are seldom displayed. The dislocations displayed are usually leading partials (see solid arrows) which are nucleated from the GBs, slip across the whole grains, and disappear in the opposite GBs. On the other hand, the complete extended dislocations (see hollow arrows) can be clearly seen in the large grains, with the leading partials, the trailing partials, and the stacking faults in between. In addition, dislocations in the large grains can also be nucleated from the grain interior besides the GBs. These phenomena indicate that the impediment effects of the GBs on the dislocation activities are more evident in the small grain sample. Another interesting phenomenon is that in the small grain sample the grain coalescence process is evident, while in the large grain sample the grain  partition is the dominant process. This can be attributed to the abundant dislocation activities in the large grains that the dislocations have large probability to intersect with each other, leading to the generation of sub-GB structures.

Dislocation activities
The effect of the grain size on the dislocation activity can be more clearly seen from the variations of Shockley dislocation density and the fraction of hexagonal closepacked (HCP) atoms during the deformation as shown in Figure 5a and b, respectively. From Figure 5a, the density of Shockley dislocation is larger for a smaller grain size. This can be explained by the larger fraction of GBs acting as the dislocation sources in the small grain sample. Nevertheless, the fraction of HCP atoms is less for the sample with smaller grain size, as shown in Figure 5b. Noting that HCP atoms are usually the atoms at the stacking faults, the proportion of HCP can be taken as a measure of the intensity of dislocation activity. It indicates that the larger grains are more favorable medium for the motion of dislocations. From Figure 5b, it also shows that the HCP proportion reaches a saturated value at a certain strain and will decrease as the strain further increases. This decrease of HCP proportion can be attributed to the absorption of Shockley partial dislocations by the GBs when the extended dislocations slide across the whole grain.

Grain rotations
In our previous work, the phenomenon of grain rotations was observed in the deformation simulations of GNG metals using the CPFEM. It was found that the magnitude of grain rotation strongly depends on the initial lattice orientation of the grain, but does not depend on the grain size. This conclusion may not hold true for real situations since the CPFEM framework only incorporates the dislocation slip mechanisms and cannot reflect the GB mechanisms such as GB mass diffusion and GB sliding. In this work, the phenomenon of grain rotation is investigated from our MD simulations. For this purpose, a program was developed to track the change of orientation for each grain based on the statistics on the relative positions of all the atoms in that grain. The variation of the orientation angle with the strain is shown in Figure 6. Here, the orientation angle is defined as the angle between the[ ] 110 direction and the horizontal direction which is the tensile axis. Each curve in Figure 6 corresponds to one grain in  the sample. The orientation angle is wrapped to the range between 0 and 60 o since the (111) crystal plane is triple symmetric. The breakings on the curves are caused by the annihilation of that grain due to its coalescence with other grains, or the generation of a new grain due to the dislocation activities. Figure 6 clearly shows that the orientation angle tends to converge to 30°, indicating the occurrence of texture. This phenomenon is consistent with the theoretical model of the grain rotation under the uniaxial deformation due to the activation of different slip systems, as demonstrated by Hosford [57]. The slip system will be activated once the resolved shear stress on it, measured by the Schmid factor, exceeds that on other slip systems. As the plastic deformation goes on, the grain tends to rotate to the orientation that the slip direction of the activated slip system is parallel to the tensile axis. The Schmid factor μ is calculated as in equation (1) = | | μ λ ϕ cos cos (1) where, λ is the angle between the slip direction and the tensile axis, and ϕ is the angle between the normal of slip plane and the tensile axis. For the configuration adopted in this work that the [111] crystal direction is always along the out-of-plane direction, the possible slip systems that can be activated are the 9 slip systems, i.e., ( (1) for the involved slip directions and slip planes, respectively. From Tables 2 and 3, the Schmid factor μ for any of the 9 slip systems can be easily obtained. For example, from Tables 2 and 3 Figure 7. As can be seen, the Schmid factor is highly dependent on the initial grain orientation.
When the initial grain orientation is between 0°and 30°, the slip system ( ) 111 [ ] 011 has the largest Schmid factor among all slip systems and thus will be activated, leading to the anti-clockwise grain rotation. While for the

[ ] 011
( / + )/ π θ cos 6 3 Table 3: ϕ cos of the involved slip planes initial grain orientation between 30°and 60°, the slip system ( ) 111 [ ] 110 has the largest Schmid factor and is the activated system, leading to the clockwise grain rotation. From Figure 7, the Schmid factor takes its peak values at the orientations 15°and 45°, respectively, for the slip systems ( ) 111 [ ] 011 and ( ) 111 [ ] 110 , indicating the largest driving forces of the anti-clockwise and clockwise grain rotation at these two orientations. It also indicates that the rotation driving force becomes nullified at the orientations 30°, 0°, and/or 60°because of the simultaneous activation of two of the three slip systems ( ) , and ( ) 111 [ ] 101 whose effects on the rotation are opposite to each other at these orientations. Indeed, from our MD simulations the grain rotation rate is positive in the orientation range (0-30°) with its maximum at 15°, while is negative in the orientation range (30-60°) with its minimum at 45°, as shown in Figure 8, which is from the 20 nm sample at strain 0.2. Here, statistics on the grain orientation are made at every 0.02 strain increment, and the grain rotation rate is measured as the change of grain orientation during each strain increment. The variation tendency in Figure 8 shows an example of the dependence of grain rotation rate on the initial grain orientation, which is representative for all the investigated samples. The resultant effect of the grain rotation is that all grains tend to the orientation of 30°.
The results of the grain rotation Δθ versus the initial grain orientation θ 0 can be fitted by a sinusoidal function as equation (2).
The amplitude of the fitted curve A represents the rate of the grain rotation, while the standard deviation of the fitting D represents the diversity of the grain rotation distribution. For the specific sample shown in Figure 8, A and D are 2.16°and 1.35°, respectively. Figure 9 shows the fitted amplitude varying with the strain. It shows that larger grains rotate faster compared with the smaller grains, owing to the more evident dislocation activities in the larger grains. This can lead to faster convergence of the grain orientation of larger grains. Indeed, if we plot the standard deviation of the grain orientation from its stable value 30°versus the strain in Figure 10, the slope of the curve is larger for the larger grain.   On the other hand, Figure 11 shows that the standard deviation of the grain rotation rate is larger for the smaller grains, especially for the elastic deformation range, indicating the larger diversity of the grain rotation distribution for the smaller grains. This phenomenon may be attributed to the effects of GB mechanisms such as the GB sliding and the GB diffusion on the grain rotation. As the grain size decreases, the GB mechanisms play a more important role, leading to the larger proportion of the GB-induced grain rotation driven by the GB energy and the external work [2]. As demonstrated in the literature [2], the GB-induced grain rotation mechanism is not a cause of texture, but can be taken as a disturbance to the sinusoidal type grain rotation distribution shown in Figure 8. It is noted from Figure 11 that as the strain reaches the plastic stage, the dependence of the standard deviation of grain rotation on the grain size becomes opaque, which is different from the elastic stage. This could be attributed to the dominance of the dislocation slip-induced grain rotation mechanism which favors the large grains in the plastic stage. After reaching the plastic stage, the abundant dislocation activities give rise to the significant slip-induced grain rotation, whose fluctuation surpasses the contribution of the GB-induced grain rotation. Therefore, the smaller grains can no longer exhibit a larger fluctuation of grain rotation than the larger grains.
As a summary of the effects of grain size on the grain rotation, Table 4 lists the amplitude A and the standard deviation D for fitting the grain rotation distribution at 0.06 strain, as well as the slope k of the grain orientation deviation from 30°versus the strain (see Figure 10) for different grain sizes. It is noted that the investigated grain size range is from 3 to 50 nm. Within that range, the results show a monotonic increase of grain rotation amplitude and a monotonic decrease of the dispersion of the grain rotation distribution as the grain size increases. Whether these trends are suitable for grain sizes outside of this range is remained for future study.

Effect of temperature
We further investigated the effect of temperature on the stress-strain curve and the grain rotation behavior. For this purpose, additional simulations were carried out for straining the sample of 12 nm grain size at the temperature of 10 K. The stress-strain curve is shown in Figure 12. It is seen that the yield stress and the plateau stress at 10 K temperature are higher than those at 300 K temperature, indicating a harder mechanical response at a lower temperature. Similar to the grain rotation behavior at 300 K, the dependence of grain rotation on the initial grain orientation at 10 K can also be fitted by a sinusoidal curve. The fitting results are shown in Figure 13. Figure 13(a) shows that the fitted amplitude is not sensitive to the temperature, while Figure 13(b) shows that the standard deviation of the fitting is larger   at 300 K than that at 10 K. This can be possibly attributed to the stronger contribution of GB-induced grain rotation at the increased temperature, leading to a larger dispersion in the grain rotation distribution at 300 K compared to that at 10 K. On the other hand, the magnitude of grain rotation induced by the dislocation slip is not much sensitive to the temperature change.

Conclusions and discussions
In this paper, grain rotations induced by the deformation in the GNG copper was explored using MD simulations. First, the computational samples of GNG copper were created using Voronoi method. A systematic study was carried out to investigate the stress-strain curves, the dislocation activities, and the grain rotation through simulations of different grainsized samples. The synergistic effects of the dislocation mechanisms and the grain boundary mechanisms on the grain rotation were evaluated. The following conclusions can be drawn.
(1) Among the investigated grain size from 3 to 50 nm, the 3 nm grain size corresponds to the smallest Young's modulus, but the largest plateau flow stress.
On the other side, the 50 nm grain size corresponds to the strongest strain hardening capability. (2) The larger grains are more favorable for the motion of dislocations, leading to larger fraction of HCP atoms during the deformation.
(3) Grain rotation was observed and its rate was found to depend on the initial crystal orientation. It makes the orientations of all grains tend to 30°, leading to the texture during the uniaxial tensile deformation. (4) The grain rotation magnitude is larger for larger grains, but the rotation rate is more diversely distributed for smaller grains, indicating more disturbance from the GB mechanisms such as the GB sliding and the GB diffusion for smaller grains. (5) The grain rotation magnitude at 10 K temperature is comparable to that at 300 K, while the rotation rate is less diversely distributed at the lower temperature.
Although we did not investigate other nanocrystalline metals, we believe that some of the conclusions obtained from the simulations of copper in this paper are applicable to other FCC metals. For example, dominated by the dislocation slip mechanism, the grain rotation magnitude is larger for larger grains, while the rotation rate is more diversely distributed for smaller grains due to the influence of the GB mechanisms. However, whether these conclusions are applicable for other types of metals are difficult to be judged. Besides, the current study is based on the conventional nanocrystalline metals, and the conclusions may not be suitable for nano-twinned crystals. The complex interactions between the twin boundaries and the dislocations [58,59] will possibly make the grain rotation behavior different from that without twin boundaries. These topics are remained for future study. Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Research funding: The authors acknowledge the support from National Natural Science Foundation of China (11772231), Shanghai Science and Technology Innovation Plan (20ZR1462700), and Shanghai Supercomputer Center.