Molecular Dynamics Simulation Study of the Mechanical Properties of Nanocrystalline Body-Centered Cubic Iron

: In the present work, the mechanical properties of nanocrystalline body-centered cubic (BCC) iron with an average grain size of 10 Å were investigated using molecular dynamics (MD) simulations. The structure has one layer of crystal grains, which means such a model could represent a structure with directional crystallization. A series of uniaxial tensile tests with different strain rates and temperatures was performed until the full rupture of the model. Moreover, tensile tests of the models with a void at the center and shear tests were carried out. In the tensile test simulations, peak stress and average values of ﬂow stress increase with strain rate. However, the strain rate does not affect the elasticity modulus. Due to the presence of void, stress concentrations in structure have been observed, which leads to dislocation pile-up and grain boundary slips at lower strains. Furthermore, the model with the void reaches lower values of peak stresses as well as stress overshoot compared to the no void model. The study results provide a better understanding of the mechanical response of nanocrystalline BCC iron under various loadings.


Introduction
Nanocrystalline (NC) materials in which the grain size is in the range of nanometers can have interesting and compared to macrocrystalline materials unpredictable values of material properties, such as high modulus of elasticity, ultimate strength, hardness, flow stress and high ductility [1][2][3][4][5][6].One of the reasons for such properties is a high ratio of grain boundaries (GBs) which act as an obstacle for dislocation movements [1,3,[7][8][9].The reason for high elasticity moduli and ultimate strengths in molecular dynamics (MD) simulations can also be caused by high strain rates [2][3][4], which need to be applied due to limited computational resources when investigating large models.It is found in the literature that the strain rates in the tensile test should be held under 5 × 10 8 s −1 to get reasonable structure responses and stress values [4,10].The most common stress-strain curves of NC solids in MD simulations have an overshoot of stress values, which are nearly constant as loading continues (at high strain values) [1,4,11].
Usually, modeled BCC structures in MD simulations have a low ratio of BCC lattice in crystal grains compared to the volume of GBs [3,5,12,13].In this work, we increased the ratio of BCC lattice in grains.Therefore, the models of one layer of grains with periodic boundary conditions were constructed.With such a model, the effect of GBs to act as seeds for dislocation movements was lowered [8].This should lead to lower stress values during the tensile test.Moreover, mechanical properties of NC metals are highly dependent on grain size [1,[14][15][16].Jeon et al. [1] found that the highest values of strength in NC BCC iron are reached at the grain sizes of approximately 14.7 nm, which is the size where the Hall-Patch behavior becomes inverse.
In MD simulations of NC iron, high plastic strains can be reached.Generally, three mechanisms leading to plastic deformation occur GB slip, dislocation glide and diffusion connected processes [17].It is investigated that increasing grain size ratio of dislocation glide mechanisms gets higher than GB slip [7].Since the diffusion of atoms has a high effect on deformation behavior [18][19][20][21], it is presented that the diffusion coefficient is highly increased by prestraining, especially at lower temperatures [22].Furthermore, it is investigated that high strain rates cause GB shearing to be a more superior deformation mechanism than diffusion-based ones [10,23,24].Mechanical properties are also affected by the number of grains, grain rotation, grain orientation and nanoscale impurities [1,10].The initial structures which are usually modeled in MD simulations can be described as perfect with no vacancies or impurities, which also contributes to the difference in values of derived mechanical properties compared to the ones of macrocrystalline iron.Despite the differences in derived mechanical properties, the beauty of MD simulations is that one can still obtain the influence of different parameters on the mechanical properties of the nanostructure, which also have the same effect on the macrostructure.Moreover, in the recent literature, there have been first-principles MD methods proposed to eliminate the intrinsic limitations of classical MD simulations [25,26].
This work aims to study the influence of different temperatures and strain rates on nanocrystalline BCC iron in the uniaxial tensile test and shear test.In the tensile test, the dependency of stress to strain is derived.In the shear test, we observed shear stress and shear strain relationship.Compared to the tensile tests in [1] and following the Hall-Petch rule [24], our model has the average grain size in the range at which we can expect the highest values of stresses for the chosen strain rates.In the tensile test it is expected the stresses increase with the strain rate [2][3][4]27] and decrease with the temperature [22,27].In the second part, the influence of a hole in the NC BCC iron during the tensile test is investigated.It is expected void leads to stress concentrations, which initiates dislocation glides at lower external loads.Moreover, in the shear tests, the peak volumetric stress and flow stresses increase with the strain rate [28] and decrease with temperature [29].

Model Generation
The initial atom coordinates were programmed in Wolfram Mathematica software [30].Firstly, x and y coordinates of crystal grain corners were chosen.Hence, the grain shapes and sizes were determined, as shown in Figure 1a.A small gap up to 2 Å was created between the grains.Therefore, there was no instability and excessive heating at the start of the MD simulations.If the gap between the grains kept smaller than 2 Å, then there were possibilities of overlap of electron clouds, i.e., nuclei are less well shielded.Therefore, the atoms from one grain would have repelled strongly with the neighboring atoms from other grains, which resulted in high kinetic energy.The thermodynamic temperature of the system arises from the kinetic energy of the vibrational motion of its atoms.Thus, one may observe the excessive heat in the MD system.We constructed the model of nine crystal grains.For each grain, the lattice orientation was chosen.Secondly, x and y coordinates of atoms with corresponding orientation were put into each grain (see Figure 1b).Therefore, we modeled a two-dimensional (2D) layer of each grain.Lastly, each grain's corner and central atoms were extruded, respectively, in the normal direction of the initial 2D structure, with the step of lattice length.Such grain orientation could represent directional crystallization since the material structure in the extruded direction does not change.OVITO software is used to visualize the structure and results [31].

Simulation Details
All the MD simulations were performed in LAMMPS [32], developed by Plimpton.In this work, the interatomic potential of BCC iron developed by Mendelev was used [33].In all three Cartesian directions, periodic boundary conditions (PBCs) were applied.Firstly, we performed the energy minimization by using the steepest descent algorithm.In order to get the best-conformation structure without residual stresses, we carried out the annealing heat cycle followed by equilibration.For the first equilibration, canonical (NVT) ensemble was used where Langevin thermostat was set to the temperature of 300 K. Annealing cycle was done in NVT ensemble where the structure was heated to 650 K in 25 ps, held at this temperature for 25 ps and cooled down to 300 K in 25 ps.The second equilibration was done in NVT, microcanonical (NVE) and isothermal-isobaric (NPT) ensemble, respectively, each lasting 25 ps.Lastly, the structure was heated to the desired temperature (300, 500 and 700 K) and deformed till rupture.The motion of each atom was determined by the forces produced by Newton's equation of motion.In the mechanical tests, the normal and shear deformation at different strain rates (0.004, 0.006, 0.008, 0.010, 0.012 ps −1 ) were applied.During deformation in the tensile test, the PBCs and zero pressures were assigned to the other two directions.In the shear test, pressure normal to the shearing plane is set to zero.Visualization and dislocation analysis (DXA) were done in OVITO visualization software [31].

Tensile Test
In Figures 2 and 3, the snapshots of the model at different strains during the tensile test for different strain rates and temperatures are shown.In both figures, the 1st snapshot was taken at a strain of 0.0, the 2nd in the elastic region of deformation at the strain of 0.03, and the 3rd at the strain where each structure reached peak stress.The 4th and 5th snapshots are taken in the plastic region of flow stresses at strains of 0.15 and 0.25.For analysis and visualization, OVITO's module dislocation analysis (DXA) [34] was used.The module allows us to observe ratios of different lattice configurations and identification of all dislocation line defects.The purple atoms are arranged in the BCC lattice while the white represents the dislocation types of <111>, <100> and <110> (see Figures 2 and 3).In all tensile simulations, the two most common deformation mechanisms were recognized, GB slip and dislocation glide [13].Moreover, grain rotation [17] and growth were noticed.The undeformed state at a strain of zero represents the structure right after the equilibrations, which was obtained after the annealing cycle.After those cycles, initial structure with 9 flawless BCC monocrystals changed to the polycrystalline structure of 9 grains containing a different density of dislocations (see numbered grains in Figure 2a).It was observed that the highest density of dislocation in the initial structure developed in grain number 4 and the lowest in the grain 2 (using DXA tool of the OVITO software [31]).Due to model construction, the ratio of BCC lattice before deformation in all the samples equals 84.3%, which is higher than most of the literature [3,12,13].In all samples, we noticed a decrease in the BCC ratio with deformation, where the ratio change in elastic compared to the plastic region is negligible.With increased strain rates, a significant decrease of BCC lattice was noticed at the same strain values.Even if only one layer of atoms, such effect is observed as a higher increase in a number of white atoms in the top layer in Figure 2b from strain 0 to 0.25, compared to Figure 2a.This means that the fraction of stacking faults in the grain lattice and fraction of GBs increases with increased strain rate in the flow stress region.From the tests done at different temperatures and the same values of strain rate, we found that the BCC ratio is also reduced with an increase in temperatures.Still, the effect of different strain rates was a more dominant factor in the BCC lattice ratio change (see Figure 3).
During the tensile test some of the initial dislocations (yellow arrows in Figures 2a and 3a,b) started to move to the boundaries and piled up there, presumably because of high energy GBs [13].Others, which ended on GBs of neighboring grains with a low difference in lattice angle, initiated dislocations in the neighboring grain.At high values of strain, some GBs started to lose their improper lattice structure.Such regions changed to the boundary at which neighboring grains share a plane of atoms.GB that went through such a process of generating twin planes from strain 0 to 0.25 are shown in Figure 2a,b (green ellipses).The main difference in GB configurations happened after the peak stresses.The reason for such processes could be more frequent GBs slides at the high strain values and increased diffusion coefficient at the higher prestrains [22].Furthermore, before some grains grew in a way that the GB between neighbor grains disappeared, and grains started to share a plane of atoms.In Figure 2b, grain growth was observed gradually from strain 0 to 0.15.Newly formed GB is marked with a black ellipse.It is seen that grain number 8 is dissolved and became part of neighboring grains 7 and 9. Therefore, both the grains sizes were increased.The grain rotation can be observed in Figures 2a  and 3b from strain 0 to 0.25 in grain number 2, which is marked with red lines.
Figure 4 shows the stress-strain curves at different temperatures and strain rates.All the stress-strain curves have a linear elastic region, which is followed by a stress overshoot.After the peak stresses, nearly constant values of flow stress were observed.The fluctuation of the curves in the flow stress region is dependent on the distribution, orientation and size of grains.The averaging of more stress-strain curves with a change in the parameters, as mentioned above, leads to more constant flow stress values.Tensile stress (GPa) 0.006 ps -1 at T = 500 K 0.008 ps -1 at T = 500 K 0.006 ps -1 at T = 700 K 0.008 ps -1 at T = 700 K

b)
Tensile strain (-) 0 0.04 0.08 0.12 0.16 0.2 0.24 In Figure 4a, tensile tests were performed at the temperature of 300 K, and strain rates from 0.004 to 0.012 ps −1 with a step size of 0.002 ps −1 were applied.It is noticed that the change in strain rate does not affect the elastic modulus.However, the peak stresses increase significantly with strain rates.Moreover, the flow stresses increased moderately with strain rates.The reason for the increase in the flow stress is that the density of partial dislocations in the region before the strain at peak stresses decreases with increased strain rate.This leads to increased values of delayed partial dislocations with increased strain rate, causing stress overshoots after the overshoot dislocation density starts to rise, which leads to stress reduction to the constant value of flow stress.Similar observations were made by Rida et al. [28].
Figure 4b demonstrates the stress-strain curves of tensile tests, which were carried out at temperatures 500 and 700 K with strain rates of 0.006 and 0.008 ps −1 .It is observed that the elastic modulus decreases with increased temperature.Moreover, a decrease of the elastic modulus with temperature leads to a lower ratio of peak to flow stresses, which is seen as lower overshoot.Furthermore, higher temperatures lead to lower peak and flow stresses.As mentioned before, not only higher strain rates but also higher temperatures decreased the BCC lattice ratio in the model when being deformed.In addition, we compared strains at a temperature of 300, 500 and 700 K at strain rates of 0.006 and 0.008 ps −1 .It is found that high temperatures decrease the energy needed for GBs sliding more than for dislocation glide, making GBs dominant deformation mechanism, which is also discussed by Ghaffarian et al. [35].The reduction of energy needed to initiate GBs slide makes plastic deformation easier, which leads to lower peak and flow stresses.

Specimen with Void
At the nanoscale, continuum approximations based on the finite element method are no longer valid since the technique can not consider dislocation movement and grain slides.From the continuum mechanics, a problem of the plate with a hole is known, where a stress concentration is three times higher on the hole edge than the pressure applied on the plate edges [36].In this section, the influence of a hole on the mechanical properties of NC BCC iron is discussed.In the [37] deformation process of NC monocrystal with a hole is analyzed, where the NC model response can be described as follows 1. Dislocations nucleate near the void; 2. With additional deformation, dislocations initiated near void propagate to the specimen edges; 3. Reduction of specimen stress resistance; 4. Necking of specimen seen at higher strain values.
In our simulations, the initial model (see Figure 1) was taken, and a through-hole with a diameter of 30 Å was made in the center of the model, as it is shown in the most left snapshot in Figure 5.We were interested in the influence such a hole would make on the stress-strain curve in a uniaxial tensile test.We need to be aware of periodic boundary conditions used, which means the necking will not form the same way as in the work of Potirniche et al. [37].Here, MD simulations were performed at a temperature of 300 K and strain rates of 0.006 and 0.008 ps −1 .Snapshots of the deformed structure at different strains are shown in Figure 5.In Figure 5, at the strain value of zero, the model after minimizations, equilibration and annealing is shown.We can comprehend the structure of 9 grains in the model.The snapshot at a strain of 0.015 shows von-Misses stress of each atom in linearly elastic structure response.One can notice that the highest stresses are formed in the GBs, however, a small stress increase formed in the BCC lattice next to the hole edge (yellow ellipses) was observed.At a strain value of 0.08, we recognized stress concentration in BCC lattice near the edge on both sides of a hole increases (red ellipses).This does not only correspond to stress concentration in NC monocrystal [37] but also to the continuum mechanics model.The next snapshot at a strain of 0.12, the structure is loaded with the highest stress.Here, we can observe that grain with number 6 starts to slide into the void.Furthermore, crack initiation is seen (black ellipse).With further increase in strain, the grain 6 kept sliding into the hole while the hole in pulling axis direction kept increasing its size.
Figure 6 depicts the stress-strain curves of the models with and without holes at the center.The model with a hole reached lower peak stress than the model without a hole.However, the flow stresses of both models are nearly the same values for the considered range of strain rates.The significant difference was also observed in the low strain region.The model without a hole has a linearly elastic response over an extended range of strain 0.04, while the hole model curve starts to divert from the linearity at a strain of 0.018.With further investigation of change in neighbor atom distances, we concluded that curve flattening at low strains happens due to small GBs slips, which became more frequent with further deformation, leading to additional stress-strain curve flattening.
Tensile strain (-) 0 0.04 0.08 0.12 0.16 0.2 0.24 0 We observed that the hole first initiated stress concentrations, leading to dislocation nucleation at smaller strains.The hole also represented a volume to where the grains could be pushed more easily, making GB slips possible at lower strains.When such a slip happened, structure relaxation occurred.The described process leads to lower overshoot compared to the one without the hole model.

Shearing Test
In this section, shearing tests on NC BCC iron are discussed, focusing on the influence of different strain rates on the mechanical properties.The MD simulations were carried out with two different strain rates of 0.004 and 0.012 ps −1 , which defined the upper face of the simulation box's translatory motion.The bottom part of the box was fixed.We need to be aware that the chosen simulation boundary conditions do not allow the simulation box to be deformed orthogonally to the translatory direction of the upper cell face.This means the shearing component of stress is the main load for structure deformation at lower shear strains.However, as loading continues, the influence of normal stresses increases.
Figure 7 shows the snapshots of the model at different shear strains during shearing test simulation at a strain rate of 0.012 ps −1 .As mentioned earlier, the purple atoms represent the ones arranged in the BCC lattice, and the white represents the dislocation types of <111>, <100> and <110>.The first snapshot is taken at the shear strain equal to zero, which shows the structure after the minimizations, equilibration and annealing.As in the tensile test, 9 grains were formed with different density of dislocations.The shear strain of 0.06 represents the structure in its maximum shear strain in the linear elastic response, where no dislocation activity was found.Atom configuration at a shear strain of 0.24 shows the structure at its maximum stress load.As we compared it to the undeformed state, grain boundary slips were noticed by observing a change in neighbor atom distances.We found out that the ratio of BCC lattice stays constant-it fluctuates around the initial value of 82%, which is different compared to the tensile test, where the ratio decreases.The last two snapshots represent the structure response in the shear stress flow region where grain boundary slips were the main deforming mechanism recognized.Figure 8 depicts the shear stress-strain curves.The structure response is linearly elastic for small shear strains-for both strain rates up to shear strain of approximately 0.05.At the shear strain of 0.05, the curve for strain rate 0.012 ps −1 starts to flatten and makes a small stress overshoot over a long shear strain.In the case of a slow strain rate of 0.004 ps −1 , at a shear strain of 0.075, the stress suddenly drops.The reason for such a response was sudden GB slip, which happened in the model.If simulations results are averaged for different grain configurations, the curves do not have such fluctuations before the peak values of shear stress and also in the flow stress region.Both overshoots are followed by constant shear flow stresses.We concluded that higher shear strain rates lead to higher peak shear stresses and higher shear flow stresses.Moreover, the strain rate has no significant influence on the shear modulus.The cross-section of the specimen is the same as in the tensile test.It is seen that for the same cross-section and strain rate, the structure can be plastically deformed with lower flow stresses when shear stress is applied instead of a normal one.This work can be effectively used to investigate the high-speed forming of nanocrystalline metals, e.g., iron and aluminum, using MD simulations [38][39][40][41].

Conclusions
A detailed study has been performed to investigate the mechanical properties of nanocrystalline BCC iron through all-atom molecular dynamics simulations.In this work, the structure of the model has one layer of crystal grains, which means such a model could represent a structure with directional crystallization.The mechanical properties of the nanocrystalline BCC iron MD model in uniaxial tension, tension with a void at the center of the model and shear were studied.The following most notable findings are summarized.
1.In the tensile test MD simulations, peak values of stress and average values of flow stress increase with strain rate, however, the strain rate does not affect the elasticity modulus.Moreover, an increase in temperature leads to a decrease in elasticity modulus, peak stress value and flow stress.2. When the model with the void at the center is under tension, stress concentrations in structure occur, which leads to dislocation pile up and grain boundary slips at lower strains.Furthermore, the model with the void reaches lower values of peak stresses as well as stress overshoot compared to the no void model.Furthermore, for different strain rate values, the stress flows of both models have approximately the same values.3.For high strain rates, in shear test simulations, the peak shear stress values, as well as flow stresses, were higher.However, the shearing modulus was not affected.

Figure 1 .
Figure 1.Initial grain configuration of the MD simulation model.(a) The front view and (b) the three-dimensional view of the model.Each different color represents a different grain.

Figure 2 . 1 Figure 3 .
Figure 2. Deformation snapshots of the specimen show the progress in tensile loading at constant temperature and different strain rates.

Figure 4 .
Figure 4. Stress-strain curves at (a) different strain rates and (b) temperatures.

Figure 5 .
Figure 5. Sanpshots of specimen with a hole during the tesile test.

Figure 6 .
Figure 6.Stress-strain curves for model with and without hole at the center.

1 Figure 7 .
Figure 7. Snapshots of the MD simulation model during shearing test.

Figure 8 .
Figure 8. Shear stress-strain curves at different strain rate loading.