Molecular Dynamics Study of Bulk Properties of Polycrystalline NiTi

: Molecular dynamics (MD) simulations were carried out to study the bulk polycrystalline properties of NiTi. Thermally driven phase transitions of NiTi between martensite and austenitewere simulated using single crystalline simulation domains. With external stress boundary conditions, MD simulation successfully predicted experimentally observed phase transformation temperatures of bulk polycrystalline. Elastic characteristics of NiTi martensite were simulated using polycrystalline simulation domains that consist of realistic disorientations and grain boundary structures. The existence of grain disorientation and grain boundary lowered the potential energy of the simulation domain, which led to more realistic elastic modulus prediction. Analysis of simulation domains that predicted realistic bulk polycrystalline properties showed that the major difference between single crystalline and polycrystalline structures is atomic stress distribution. simulation domain. The results indicate that the polycrystalline structure properties can be predicted


Introduction
Shape-memory alloys (SMAs) have been utilized in many fields due to their attendant shape-memory effect and superelasticity. Applications of SMAs have been widely developed in the biomedical, aerospace, construction, and automotive industries [1]. Specifically, nickeltitanium alloy (i.e., Nitinol) has been one of the most promising shape-memory alloys for commercial applications. Due to this reason, many molecular dynamics (MD) studies focused on understanding and predicting its atomic structure and related phenomena.
Although many aspects of NiTi were revealed, most of the studies focused on the characteristics of single crystalline NiTi. However, some of the predicted properties from MD simulations of single crystalline NiTi do not accurately represent the characteristics of bulk polycrystalline NiTi, including thermally induced phase transformation and elastic characteristics. For thermally induced phase transformation, Ko et al. [2] used the secondnearest-neighbor (2NN) MEAM potential to show phase transitions between twinned B19 structure and B2 structure, but the prediction of too high austenite finish temperature (A_f above 480 K) made the result unrealistic. Zhong et al. [3] used the Finnis-Sinclair (FS) potential in heating and cooling simulation of NiTi, but only showed abrupt changes of a sum of all the shear components without structural information of martensite structure. Mutter and Nielaba [4] showed abrupt changes in both bonding angle and edge lengths indicating phase transitions during heating/cooling, but the martensite structure with different edge lengths and bonding angles was formed after one heating/cooling cycle. Moreover, phase transformation temperatures were only calculated in the forms of T M = 0.5(M s + M f ) and T A = 0.5(A s + A f ), which are difficult to compare with experimental data. For elastic characteristics of NiTi, the prediction of elastic properties of martensite NiTi has been an issue. Many simulation studies from single crystalline NiTi calculated elastic constants and associated Young's modulus using Hill methods [5]. While experimentally observed Young's modulus was 20-40 GPa [6], predictions from MD and first-principle studies showed much higher values of 133-526 GPa [7][8][9]. NiTi has been widely used since its Metals 2021, 11, 1237 2 of 15 mechanical properties can be tailored by composition and post heat treatment. Therefore, precise prediction of phase transformation temperatures and elastic properties is important for designing mechanical properties of specific purposes.
Recently, MD studies using polycrystalline NiTi simulation domains were conducted to achieve a more accurate prediction of bulk properties. Stress-induced martensitic transformations of NiTi were studied using polycrystalline simulation domains that consist of nano-sized grains [10][11][12][13]. Although very high strain rates (10 7 -10 9 s −1 ) were used, the above studies successfully showed the elastic modulus of austenite as 60-80 GPa, which agreed well with experimental data [14]. It should be noted that elastic modulus is independent of strain rate [15]. However, no attempts have been made to predict realistic elastic properties of martensite NiTi using polycrystalline microstructure. In addition, the predicted A f by MD simulation was 200 K and 300 K from 10 nm and 15 nm polycrystalline simulation domains, respectively [10]. Considering that the experimentally observed A f is 380 K [16], it can be concluded that either a different approach or a more realistic polycrystalline simulation domain is necessary.
Applying proper stress boundary conditions needs to be carefully considered since studies showed that the role of hydrostatic pressure and shear stress is important when predicting the properties of martensite NiTi. Bakhtiari et al. [17] showed that stable martensite (B19 ) can be achieved with the hydrostatic pressure above 2 GPa. Huang et al. [18] reported that the B19 structure is unstable and can be stabilized only with shear stresses between 0.77 and 1.4 GPa. Wagner and Windl [19] stated that shear stress above 1 GPa is necessary to stabilize B19 structure. While the above studies showed the importance of applied stress, there remain a number of issues to be considered. Such large stress states from both hydrostatic pressure and shear stress cases for stabilizing B19 structure can only exist during loading (normal or shear directions) situations or after cold deformation processes such as cold rolling or extrusion. Considering that martensite NiTi can also be formed in stress-free conditions [20], it is doubtful that such large stress is required in the prediction of bulk NiTi properties.
The effects of grain structure should also be considered since MD simulation results showed that both stress-induced and thermally induced phase transformations are influenced by grain structures. Zhang et al. [12] showed that simulation domains of different polycrystalline structures experienced different stress-induced martensitic transformation behaviors, which emphasized the importance of grain size, grain orientation, and grain boundary. Wang et al. [21] compared single crystalline and polycrystalline structures in MD simulations to show the effects of grain boundaries on stress-induced phase transformations. Li et al. [10] showed that MD simulation domains with different grain structures exhibited different thermally induced phase transformation temperatures. The above studies showed the importance of proper grain structure in the prediction of bulk polycrystalline properties.
This study is focused on predictions of bulk polycrystalline NiTi properties using MD simulations by considering the effects of stress boundary conditions and grain structures. Thermally induced phase transitions of NiTi under different external stresses were simulated in order to properly predict precise phase transformation characteristics. External stresses were applied in both normal and shear directions on a single crystalline simulation domain to represent the bulk material state. In this way, bulk polycrystalline properties can be predicted without spending extremely large computational resources. Elastic characteristics of polycrystalline NiTi martensite were calculated using different grain structures including experimentally observed grain disorientations and grain boundaries. Simulation domains that successfully predicted bulk polycrystalline properties were analyzed to understand bulk polycrystalline structure. While properties of NiTi are well known, but it is also well known they are sensitive to the Ni percentage and grain size. Being able to predict such affects would be beneficial to the future design and applications of NiTi.

Simulation Method
The open-source MD program LAMMPS [22] and the FS potential developed by Zhong et al. [3,9] were used. Two MD simulations were conducted: i) thermally induced phase transformation under external stress conditions, and ii) elastic constant calculations of polycrystalline with misorientations and grain boundary.
Thermally induced phase transformation of NiTi was simulated by changing the temperature of the system. Atom positions were generated and velocities were sampled from the isothermal-isobaric ensemble with constant N, p, and T (NPT) controlling temperature and external stress tensors, where N is the number of particles, p is pressure, and T is temperature. B2 (austenite) structure in an orthogonal shape box was constructed and first stabilized at 450 K for 500 ps (250,000 time steps). The temperature was decreased from 450 K to 10 K and then increased again to 450 K with a rate of 0.1 K/ps. External stresses on the domain were applied in both normal and shear directions (x, y, z, xy, xz, and yz). Each stress component was controlled independently. In this study, seven different external stress settings were applied on simulation domains: 0 GPa, 10 −4 GPa, 10 −3 GPa, 0.01 GPa, 0.05 GPa, 0.1 GPa, 0.5 GPa, and 1.1 GPa. The time step of ∆t = 0.002 ps was used under periodic boundary conditions. Additional simulation cases were conducted only with hydrostatic pressure applied (normal stresses in a coupled manner) in order to provide comparison data.
The adiabatic elastic constants at a finite temperature were calculated from the stress tensor fluctuation during elastic deformation of the simulation domain. The compliance tensor (S ijkl ) and the elastic constants (C ijkl ) are calculated using the Equations (1) and (2).
x abi x abj x abk x abl av (2) where ε ij is the strain tensor, V 0 is the equilibrium volume of the N particle system, P ij is the internal stress tensor, x ab is the vector joining a and b of length r ab , and f is the function of f = r −2 u − r −1 u with the pair potential denoted by u(r). The ensemble average is denoted by the bracket, av . The detailed methodology is referred to in the work of Kimizuka et al. [23]. Bulk mechanical properties of polycrystalline structures were calculated from elastic constants. Macroscopic elastic modulus E, shear modulus G, and bulk modulus B were calculated by using Voigt [24], Reuss [25], and Hill averaged values [5], as described by the following Equations (3)-(9): 15 G R = 4(S 11 + S 22 + S 33 ) − 4(S 12 + S 23 + S 31 ) + 3(S 44 + S 55 + S 66 ) where B and G are Voigt-Reuss-Hill average of bulk modulus and shear modulus, respectively. Subscripts R and V represent the Reuss and the Voigt averages. E V and E R are the elastic modulus calculated by the Voigt and Reuss methods. Prediction of bulk polycrystalline properties was conducted using simple, but realistic simulation domains. As shown in Figure 1, the simulation domain was constructed with different grain distributions and grain boundaries, which were reported as important factors affecting the elastic modulus calculation [26,27]. Two adjacent grains were arranged following the experimentally observed sigma disorientations (Σ3, Σ5, Σ9, and Σ15) according to the coincident site lattice (CSL) theory [28]. While the left-side grain 1 had a constant orientation angle of 0 • , orientation angles for the right-side grain 2 were changed to match sigma disorientation values of Σ3, Σ5, Σ9, and Σ15. Multiple simulation domains were constructed with different initial grain boundary (GB) thickness values of 0-15 Å. It was set in the range to represent realistic grain structure since the grain boundary thickness of NiTi was reported as approximately 1 nm [29,30]. The maximum grain boundary thickness was set as 15 Å to consider widely accepted grain boundary thickness, which is equivalent to 2-3 atomic layers [31]. The total numbers of atoms inside each grain for the simulations of thermally induced phase transformation and elastic constant calculation were selected to be 60,016 and 221,050 on the basis of domain size sensitivity test results. Prediction of bulk polycrystalline properties was conducted using simple, but realistic simulation domains. As shown in Figure 1, the simulation domain was constructed with different grain distributions and grain boundaries, which were reported as important factors affecting the elastic modulus calculation [26,27]. Two adjacent grains were arranged following the experimentally observed sigma disorientations (Σ3, Σ5, Σ9, and Σ15) according to the coincident site lattice (CSL) theory [28]. While the left-side grain 1 had a constant orientation angle of 0°, orientation angles for the right-side grain 2 were changed to match sigma disorientation values of Σ3, Σ5, Σ9, and Σ15. Multiple simulation domains were constructed with different initial grain boundary (GB) thickness values of 0-15 Å. It was set in the range to represent realistic grain structure since the grain boundary thickness of NiTi was reported as approximately 1 nm [29,30]. The maximum grain boundary thickness was set as 15 Å to consider widely accepted grain boundary thickness, which is equivalent to 2-3 atomic layers [31]. The total numbers of atoms inside each grain for the simulations of thermally induced phase transformation and elastic constant calculation were selected to be 60,016 and 221,050 on the basis of domain size sensitivity test results.

Phase Transformation Temperature
Effects of external stress on thermally induced phase transition were studied from two simulation groups: normal stress (hydrostatic pressure) only and both normal and shear stresses. Potential energy (PE) profile results from the two groups during cooling and heating simulations are shown in Figure 2 and Figure 3, respectively. Abrupt changes in PE profiles indicate clear transitions between austenite and martensite structures at phase transformation temperatures. Phase transformation temperatures (M , M , A , and A ) from two simulation groups were acquired from tangential lines drawn from PE profiles where abrupt changes were observed and compared, as shown in Figure 4. It was observed that phase transformation temperatures were not strongly affected by increased normal stress within the 0-1.1 GPa range. Phase transformation temperatures decreased by ~4 K and ~11 K as normal stress increased to 0.5 GPa and 1.1 GPa, as shown in Figure 4b. Wan et al. [32] reported similar results such that phase transformation

Phase Transformation Temperature
Effects of external stress on thermally induced phase transition were studied from two simulation groups: normal stress (hydrostatic pressure) only and both normal and shear stresses. Potential energy (PE) profile results from the two groups during cooling and heating simulations are shown in Figures 2 and 3, respectively. Abrupt changes in PE profiles indicate clear transitions between austenite and martensite structures at phase transformation temperatures. Phase transformation temperatures (M s , M f , A s , and A f ) from two simulation groups were acquired from tangential lines drawn from PE profiles where abrupt changes were observed and compared, as shown in Figure 4. It was observed that phase transformation temperatures were not strongly affected by increased normal stress within the 0-1.1 GPa range. Phase transformation temperatures decreased by~4 K and~11 K as normal stress increased to 0.5 GPa and 1.1 GPa, as shown in Figure 4b. Wan et al. [32] reported similar results such that phase transformation temperatures decreased by 5-10 K as hydrostatic pressure increased from 0 GPa to 1 GPa. On the other hand, results from normal + shear stresses showed a different trend. As shown in Figure 4c, phase transformation temperatures increased as external stress increased. In addition, phase transformation temperatures were more strongly affected by increased external stress when both stresses were applied. Especially, A f temperatures increased by 43 K and 56 K respectively as external stress was increased to 0.05 GPa and 0.1 GPa. No abrupt changes in PE profiles were observed, indicating that phase transformation was suppressed when external stresses in both normal and shear directions exceeded 0.1 GPa. The comparison of phase transformation temperatures from the two groups clearly showed that the presence of shear stress affected phase transformation temperatures. temperatures decreased by 5-10 K as hydrostatic pressure increased from 0 GPa to 1 GPa. On the other hand, results from normal + shear stresses showed a different trend. As shown in Figure 4c, phase transformation temperatures increased as external stress increased. In addition, phase transformation temperatures were more strongly affected by increased external stress when both stresses were applied. Especially, A temperatures increased by 43 K and 56 K respectively as external stress was increased to 0.05 GPa and 0.1 GPa. No abrupt changes in PE profiles were observed, indicating that phase transformation was suppressed when external stresses in both normal and shear directions exceeded 0.1 GPa. The comparison of phase transformation temperatures from the two groups clearly showed that the presence of shear stress affected phase transformation temperatures.   temperatures decreased by 5-10 K as hydrostatic pressure increased from 0 GPa to 1 GPa. On the other hand, results from normal + shear stresses showed a different trend. As shown in Figure 4c, phase transformation temperatures increased as external stress increased. In addition, phase transformation temperatures were more strongly affected by increased external stress when both stresses were applied. Especially, A temperatures increased by 43 K and 56 K respectively as external stress was increased to 0.05 GPa and 0.1 GPa. No abrupt changes in PE profiles were observed, indicating that phase transformation was suppressed when external stresses in both normal and shear directions exceeded 0.1 GPa. The comparison of phase transformation temperatures from the two groups clearly showed that the presence of shear stress affected phase transformation temperatures.

Martensite Twin Formation
Different trends of phase transformation temperatures should be explained by the atomic structure of martensite. Figures 5 and 6 show the atomic structures colored by local von Mises shear strain invariant [33] at 10 K from two simulation groups with applied stresses of 0.1 GPa for both. Martensite twins were observed from both cases, but their orientations were different. The orientation of the twin was expressed in the form of (h k l) [a b c], indicating the (h k l) shear plane along with the [a b c] direction. Orientations were (0 1 1) (1 0 0) and (1 0 0) [011] for normal stress (hydrostatic pressure) case and normal + shear stress case, respectively. Martensite twin of (0 1 1) was reported as type II

Martensite Twin Formation
Different trends of phase transformation temperatures should be explained by the atomic structure of martensite. Figures 5 and 6 show the atomic structures colored by local von Mises shear strain invariant [33] at 10 K from two simulation groups with applied stresses of 0.1 GPa for both. Martensite twins were observed from both cases, but their orientations were different. The orientation of the twin was expressed in the form of (h k l) [a b c], indicating the (h k l) shear plane along with the [a b c] direction. Orientations were (0 1 1) (1 0 0) and (1 0 0) [011] for normal stress (hydrostatic pressure) case and normal + shear stress case, respectively. Martensite twin of (0 1 1) was reported as type II  1 1) (1 0 0) and (1 0 0) [11] for normal stress (hydrostatic pressure) case and normal + shear stress case, respectively. Martensite twin of (0 1 1) was reported as type II twinning mode [34], which was dominantly observed in experiments [35,36]. These results show the effect of shear stress on the simulation domain, which not only increased phase transformation temperatures but also achieved a realistic form of the martensite twin (1 0 0) (0 1 1). twinning mode [34], which was dominantly observed in experiments [35,36]. These results show the effect of shear stress on the simulation domain, which not only increased phase transformation temperatures but also achieved a realistic form of the martensite twin (1 0 0) (0 1 1).  Simulation domains at 10 K under different normal + shear stresses are shown in Figure 7. Additional stress cases of 0.2 GPa and 0.3 GPa are also shown in Figure 7. It should be noted that martensite twins were not formed when large external stresses (equal to or greater than 0.2 GPa) were applied. The formation of the martensite twin was suppressed, and the slip was observed instead. Similarly, phase transformation between martensite and austenite was suppressed with large external stresses as shown in Figure  3. It is the known characteristics of NiTi that detwinning and slip occur beyond the critical stress [37]. Therefore, it is important to compare simulation results with experimental cases to validate whether this simulation study predicted a realistic critical stress value. Stress-induced martensite starts forming under a certain amount of stress (critical stress), twinning mode [34], which was dominantly observed in experiments [35,36]. These results show the effect of shear stress on the simulation domain, which not only increased phase transformation temperatures but also achieved a realistic form of the martensite twin (1 0 0) (0 1 1).  Simulation domains at 10 K under different normal + shear stresses are shown in Figure 7. Additional stress cases of 0.2 GPa and 0.3 GPa are also shown in Figure 7. It should be noted that martensite twins were not formed when large external stresses (equal to or greater than 0.2 GPa) were applied. The formation of the martensite twin was suppressed, and the slip was observed instead. Similarly, phase transformation between martensite and austenite was suppressed with large external stresses as shown in Figure  3. It is the known characteristics of NiTi that detwinning and slip occur beyond the critical stress [37]. Therefore, it is important to compare simulation results with experimental cases to validate whether this simulation study predicted a realistic critical stress value. Stress-induced martensite starts forming under a certain amount of stress (critical stress), Simulation domains at 10 K under different normal + shear stresses are shown in Figure 7. Additional stress cases of 0.2 GPa and 0.3 GPa are also shown in Figure 7. It should be noted that martensite twins were not formed when large external stresses (equal to or greater than 0.2 GPa) were applied. The formation of the martensite twin was suppressed, and the slip was observed instead. Similarly, phase transformation between martensite and austenite was suppressed with large external stresses as shown in Figure 3. It is the known characteristics of NiTi that detwinning and slip occur beyond the critical stress [37]. Therefore, it is important to compare simulation results with experimental cases to validate whether this simulation study predicted a realistic critical stress value. Stress-induced martensite starts forming under a certain amount of stress (critical stress), as indicated by plateau lines in stress-strain plots. Andani and Elahinia [38] experimentally showed critical stress values of 0.13-0.25 GPa when both tensional and torsional stresses were applied at the same time with axial and shear strain rates proportional to each other. This experiment result can be compared with simulation cases when external stresses were applied in both normal and shear directions. According to observed slips in simulation domains as shown in Figure 7, the critical stress predicted was 0.1-0.2 GPa, which indicates simulations captured the realistic phenomena of bulk NiTi structure. as indicated by plateau lines in stress-strain plots. Andani and Elahinia [38] experimentally showed critical stress values of 0.13-0.25 GPa when both tensional and torsional stresses were applied at the same time with axial and shear strain rates proportional to each other. This experiment result can be compared with simulation cases when external stresses were applied in both normal and shear directions. According to observed slips in simulation domains as shown in Figure 7, the critical stress predicted was 0.1-0.2 GPa, which indicates simulations captured the realistic phenomena of bulk NiTi structure.  Table 1 shows phase transformation temperatures from simulation cases with different external stresses in both normal and shear directions. Experimental data [16] and simulation results from other MD studies [39,40] are also included for comparison. It should be noticed that calculated phase transformation temperatures from MD simulations were lower (~40 K) with external stress of 0-0.01 GPa. As external stress increased to 0.05-0.1 GPa, simulation results agreed well with experimental data showing the A discrepancy decreased from 45.4 K to 4.6 K and 11.5 K for 0.05 GPa and 0.1 GPa cases, respectively. These results indicate two things: a) phase transformation temperatures were affected by external stress in normal and shear directions, and b) more realistic properties, bulk material properties in this case, required external stresses in MD simulations for precise prediction. Bulk material always contains internal stress caused by microstructural features such as point defects, concentration fluctuations, line defects, phase boundaries, and grain boundaries that are known to play important roles in transformation characteristics [19]. Then, the question is what is the typical amount of internal stress in bulk Nitinol that exhibits shape-memory ability. Residual stresses obtained from X-ray diffraction studies showed 0.06-0.13 GPa [41] and 0.1-0.2 GPa [42] from bulk NiTi alloys in the form of 1mm plates. The amount of stress observed in experiments agreed well with applied stresses in simulations when predicted phase transformation temperatures matched well with experimental data. Therefore, it can be concluded that applying the realistic stress leads to the correct prediction of bulk characteristics of Nitinol: (a) formation of martensite twins and (b) more accurate prediction of phase transformation temperatures.  Table 1 shows phase transformation temperatures from simulation cases with different external stresses in both normal and shear directions. Experimental data [16] and simulation results from other MD studies [39,40] are also included for comparison. It should be noticed that calculated phase transformation temperatures from MD simulations were lower (~40 K) with external stress of 0-0.01 GPa. As external stress increased to 0.05-0.1 GPa, simulation results agreed well with experimental data showing the A f discrepancy decreased from 45.4 K to 4.6 K and 11.5 K for 0.05 GPa and 0.1 GPa cases, respectively. These results indicate two things: a) phase transformation temperatures were affected by external stress in normal and shear directions, and b) more realistic properties, bulk material properties in this case, required external stresses in MD simulations for precise prediction. Bulk material always contains internal stress caused by microstructural features such as point defects, concentration fluctuations, line defects, phase boundaries, and grain boundaries that are known to play important roles in transformation characteristics [19]. Then, the question is what is the typical amount of internal stress in bulk Nitinol that exhibits shape-memory ability. Residual stresses obtained from X-ray diffraction studies showed 0.06-0.13 GPa [41] and 0.1-0.2 GPa [42] from bulk NiTi alloys in the form of 1mm plates. The amount of stress observed in experiments agreed well with applied stresses in simulations when predicted phase transformation temperatures matched well with experimental data. Therefore, it can be concluded that applying the realistic stress leads to the correct prediction of bulk characteristics of Nitinol: (a) formation of martensite twins and (b) more accurate prediction of phase transformation temperatures. Phase transformation temperatures of Ni-rich NiTi with different Ni contents were simulated and compared with experimental data [16], as shown in Figure 8. Simulation domains were under normal + shear stresses of 0 GPa, 0.05 GPa, and 0.1 GPa. For bulk Nirich NiTi alloys, phase transformation temperatures are not only affected by Ni contents but also by a post-aging heat treatment that produces precipitates within the NiTi matrix [43]. Therefore, experimental data of NiTi alloys without post-aging heat treatment were selected for comparison to exclude any effects from precipitates. Phase transformation temperatures from Ni-rich NiTi simulation cases agreed well with experimental data with applied stresses of 0.05 GPa and 0.1 GPa. Results again verified that external stress on simulation is necessary to predict the realistic physical properties of bulk NiTi material.  Phase transformation temperatures of Ni-rich NiTi with different Ni contents we simulated and compared with experimental data [16], as shown in Figure 8. Simulati domains were under normal + shear stresses of 0 GPa, 0.05 GPa, and 0.1 GPa. For bu Ni-rich NiTi alloys, phase transformation temperatures are not only affected by Ni co tents but also by a post-aging heat treatment that produces precipitates within the N matrix [43]. Therefore, experimental data of NiTi alloys without post-aging heat tre ment were selected for comparison to exclude any effects from precipitates. Pha transformation temperatures from Ni-rich NiTi simulation cases agreed well with e perimental data with applied stresses of 0.05 GPa and 0.1 GPa. Results again verified th external stress on simulation is necessary to predict the realistic physical properties bulk NiTi material.

Elastic Modulus Results
Elastic constants of martensite structure were calculated from two different simu tion cases: single crystalline (SC) and bi-crystalline (BC) with an initial grain bounda

Elastic Modulus Results
Elastic constants of martensite structure were calculated from two different simulation cases: single crystalline (SC) and bi-crystalline (BC) with an initial grain boundary thickness of 0 Å. Equi-atomic NiTi structures were used in simulations at 200 K, when entire structures exist as martensite. The elastic constants obtained from MD, experiments, and ab initio simulations are shown in Table 2. Two ab initio calculations, Abinit and VASP, were conducted on the super-cell containing (1 0 0) twins using the linear response function (RF) method and ground state (GS) calculations, respectively [44]. All constants from BC cases were smaller than those of SC cases such that BC simulations predicted C44, C55, and C66 values as small as 16-28 GPa while SC and ab initio results showed larger values of 35-55 GPa and 21-90 GPa, respectively. The importance of C55 on the elastic modulus of B19 was discussed by Šittner et al. [14], indicating that C55 has a large impact on the martensite elastic modulus of B19 in tension. Similar to C55, C44 and C66 also affect the shear characteristics of the structure. Therefore, it can be stated that BC structure with grain boundary is overall softer and especially softer in shear deformation. Elastic properties B, G, and E were calculated from elastic constants as shown in Table 3. BC simulation cases yielded lower elastic modulus (49-59 GPa) compared with the SC case (86 GPa) and ab initio (119-137 GPa), and are closer to experimentally observed values of 20-40 GPa [6]. However, it is still debatable whether BC simulation results are truly predicting bulk polycrystalline properties. For martensite NiTi, elastic modulus can be measured differently by different techniques: static loading/unloading [6], dynamic impulse excitation of vibration [45], and in situ synchrotron X-ray diffraction [14]. While the static method yielded low measured elastic modulus in the range of 20-40 GPa [6], the measured values by the dynamic method and the X-ray diffraction method were much higher at 55-70 GPa [45] and 120 GPa [14], respectively. BC simulation results agreed well with results from the dynamic method; however, this study aims to predict static bulk properties of NiTi. Moreover, it was reported that the dynamic method cannot depict stress history [46], which cannot be ignored easily in shape-memory applications. In addition to the disorientation of adjacent grains, grain boundary also affects the elastic modulus of the simulation domain. The BC simulation domains with different initial grain boundary thicknesses are shown in Figure 9. Atoms are colored by their potential energy, and it can be observed that the grain boundary region is composed of atoms with lower potential energy. The averaged potential energy distributions of Ni and Ti atoms across grain boundaries for different simulation cases are shown in Figure 10. The potential energy distribution of the SC case was also plotted for comparison. Grain boundaries consisting of atoms with lower potential energy were observed between two grains. It should be noted that Ni atoms in non-boundary regions of the BC domains showed lower potential energy compared with the SC case, while Ti atoms in both BC and SC domains showed similar potential energy. In addition to the disorientation of adjacent grains, grain boundary also affects the elastic modulus of the simulation domain. The BC simulation domains with different initial grain boundary thicknesses are shown in Figure 9. Atoms are colored by their potential energy, and it can be observed that the grain boundary region is composed of atoms with lower potential energy. The averaged potential energy distributions of Ni and Ti atoms across grain boundaries for different simulation cases are shown in Figure 10. The potential energy distribution of the SC case was also plotted for comparison. Grain boundaries consisting of atoms with lower potential energy were observed between two grains. It should be noted that Ni atoms in non-boundary regions of the BC domains showed lower potential energy compared with the SC case, while Ti atoms in both BC and SC domains showed similar potential energy. The calculated elastic modulus from simulation domains is shown in Figure 11. The average potential energy of atoms was calculated to represent each simulation domain. Due to atoms with lower potential energies in grain boundary regions, as shown in Fig-Figure 9. Simulation domains of BC with different initial grain boundary thicknesses. Atoms were colored by their potential energy.
The calculated elastic modulus from simulation domains is shown in Figure 11. The average potential energy of atoms was calculated to represent each simulation domain. Due to atoms with lower potential energies in grain boundary regions, as shown in Figures 9 and 10, BC simulation domains showed lower potential energy compared with the SC simulation domain. Moreover, the increased grain boundary thickness in BC cases led to even lower potential energy and elastic modulus. The trend line in Figure 11 shows the relationship between potential energy and elastic modulus. Elastic modulus decreased to 42.8 GPa as the potential energy of the simulation domain increased to−4.854 eV/atom, following the second-order polynomial trend curve. According to the result, experimentally observed elastic modulus of 20-40 GPa [6] is expected to be achieved with a large simulation domain consisting of many grains and grain boundaries. While it is difficult to simulate such a large simulation domain due to the required prohibitive computational cost at present, this study clearly shows the way of predicting the characteristics of polycrystalline NiTi. 42.8 GPa as the potential energy of the simulation domain increased to−4.854 eV/atom, following the second-order polynomial trend curve. According to the result, experimentally observed elastic modulus of 20-40 GPa [6] is expected to be achieved with a large simulation domain consisting of many grains and grain boundaries. While it is difficult to simulate such a large simulation domain due to the required prohibitive computational cost at present, this study clearly shows the way of predicting the characteristics of polycrystalline NiTi.

Atomic Stress Distribution in the Simulation Domain
The above results showed that realistic predictions of bulk polycrystalline properties (phase transformation temperatures and elastic modulus) were achieved by either applying increased external stress (up to 0.1 GPa) on the simulation domain or by constructing grain disorientations and grain boundaries. Although the two methods seem following the second-order polynomial trend curve. According to the result, experimentally observed elastic modulus of 20-40 GPa [6] is expected to be achieved with a large simulation domain consisting of many grains and grain boundaries. While it is difficult to simulate such a large simulation domain due to the required prohibitive computational cost at present, this study clearly shows the way of predicting the characteristics of polycrystalline NiTi.

Atomic Stress Distribution in the Simulation Domain
The above results showed that realistic predictions of bulk polycrystalline properties (phase transformation temperatures and elastic modulus) were achieved by either applying increased external stress (up to 0.1 GPa) on the simulation domain or by constructing grain disorientations and grain boundaries. Although the two methods seem

Atomic Stress Distribution in the Simulation Domain
The above results showed that realistic predictions of bulk polycrystalline properties (phase transformation temperatures and elastic modulus) were achieved by either applying increased external stress (up to 0.1 GPa) on the simulation domain or by constructing grain disorientations and grain boundaries. Although the two methods seem completely separate from each other, both methods changed atomic stress distributions in simulation domains, especially in the shear direction. Ti atoms with different stress XY values from both simulation domains (elastic modulus and phase transformation temperature simulation cases) were plotted in the form of a histogram, as shown in Figure 12. It should be noted that more Ti atoms with positive stress XY were observed from both simulation domains that predicted more realistic bulk polycrystalline properties: initial GB above 10 Å in Figure 12a and applied external stress of 0.1 GPa in Figure 12b. In addition, those cases showed very similar XY stress distribution of Ti atoms. The similarity can also be checked by comparing the average shear stress XY applied on simulation domains. Stress XY was 0.26 GPa and 0.32 GPa from elastic modulus simulation domains with an initial grain boundary thickness of 10 Å and 15 Å, respectively. The values are close to 0.1 GPa, which was applied as external stress on the phase transformation temperature simulation domain. The results indicate that the polycrystalline structure properties can be predicted by applying proper stress boundary conditions, especially in the XY direction. Therefore, it can be concluded that such different atomic stress distributions are the major difference between single crystalline and bulk polycrystalline structures.
GB above 10 Å in Figure 12a and applied external stress of 0.1 GPa in Figure 12b. In addition, those cases showed very similar XY stress distribution of Ti atoms. The similarity can also be checked by comparing the average shear stress XY applied on simulation domains. Stress XY was 0.26 GPa and 0.32 GPa from elastic modulus simulation domains with an initial grain boundary thickness of 10 Å and 15 Å, respectively. The values are close to 0.1 GPa, which was applied as external stress on the phase transformation temperature simulation domain. The results indicate that the polycrystalline structure properties can be predicted by applying proper stress boundary conditions, especially in the XY direction. Therefore, it can be concluded that such different atomic stress distributions are the major difference between single crystalline and bulk polycrystalline structures.

Conclusions
This work presented an MD study using an FS potential to predict bulk properties of polycrystalline NiTi: phase transformation temperatures and elastic modulus. Thermally induced phase transformation between martensite and austenite was more accurately simulated with external stress in both normal and shear directions. With applied stress of 0.05-0.1 GPa, which agreed well with experimentally observed internal stress, the type II martensite twins (1 0 0) (0 1 1) were formed as observed from experiments. Simulation domains with applied stress boundary conditions predicted experimentally reported phase transformation temperatures of bulk polycrystalline NiTi with various Ni contents. While other simulation studies predicted the maximum A discrepancy of 45 K or above, it was reduced to less than 12 K in this study. Such results are important since phase transformation temperatures are very sensitive to the components, and the difference of ~50 K is equivalent to 0.5% of compositional difference [47]. Considering that a very small difference in Ni content has a large effect on precipitate formation and mechanical properties [43], the precise prediction of phase transformation temperatures is critical.
Elastic characteristics of NiTi were studied with single crystalline and bi-crystalline simulation domains with different grain structures: experimentally reported grain disorientations and grain boundaries thicknesses. The existence of grain boundaries showed the major effect on softening C55 and C66. Moreover, the potential energy of simulation domains was lowered with the existence of grain disorientations and grain boundaries.

Conclusions
This work presented an MD study using an FS potential to predict bulk properties of polycrystalline NiTi: phase transformation temperatures and elastic modulus. Thermally induced phase transformation between martensite and austenite was more accurately simulated with external stress in both normal and shear directions. With applied stress of 0.05-0.1 GPa, which agreed well with experimentally observed internal stress, the type II martensite twins (1 0 0) (0 1 1) were formed as observed from experiments. Simulation domains with applied stress boundary conditions predicted experimentally reported phase transformation temperatures of bulk polycrystalline NiTi with various Ni contents. While other simulation studies predicted the maximum A f discrepancy of 45 K or above, it was reduced to less than 12 K in this study. Such results are important since phase transformation temperatures are very sensitive to the components, and the difference of~50 K is equivalent to 0.5% of compositional difference [47]. Considering that a very small difference in Ni content has a large effect on precipitate formation and mechanical properties [43], the precise prediction of phase transformation temperatures is critical.
Elastic characteristics of NiTi were studied with single crystalline and bi-crystalline simulation domains with different grain structures: experimentally reported grain disorientations and grain boundaries thicknesses. The existence of grain boundaries showed the major effect on softening C55 and C66. Moreover, the potential energy of simulation domains was lowered with the existence of grain disorientations and grain boundaries. While single crystalline simulation predicted martensite elastic modulus to be too high (87 GPa), lower elastic modulus (as low as 43 GPa) was observed from BC cases, which was close to the experimentally observed 20-40 GPa. [6] The precise prediction of elastic modulus is necessary in designing mechanical properties of NiTi. The relationship between potential energy and elastic modulus clearly showed that a realistic prediction of polycrystalline NiTi martensite properties can be achieved with proper simulation domains.
Similarities in atomic stress distribution were observed from two different simulation domains that predicted realistic bulk polycrystalline properties: single crystalline with pressure boundary conditions and bi-crystalline with grain disorientations and grain boundaries. Results showed that the major difference between single crystalline and polycrystalline structures is atomic stress distribution, which is necessary for prediction of realistic material properties.
In summary, this paper provides a way of predicting realistic polycrystalline NiTi properties by using bi-crystalline MD. The new knowledge this work provides includes the proper amount of stresses to be applied to the MD simulation to correctly predict phase transformation temperatures, and how considering the boundary conditions can improve the prediction of elastic modulus. While properties of NiTi are well known, it is also well known they are sensitive to the Ni content and grain structure. Being able to predict such effects would be beneficial to the future design and applications of NiTi. This work provides a way of correctly predicting both phase transformation temperatures and elastic properties with the different Ni and Ti compositions.
Author Contributions: J.L.: conceptualization, methodology, software, validation, formal analysis, data curation and writing-original draft preparation; Y.C.S.: conceptualization, resources, writingreview and editing, supervision, project administration. All authors have read and agreed to the published version of the manuscript.