Molecular dynamics study of phase transformations in NiTi shape memory alloy embedded with precipitates

The present study utilizes molecular dynamics simulations to study the athermal and stress-induced martensitic transformation of B2 to B19′ phase in a Ni-Ti alloy and the effect of precipitates on the phase transformation. The simulations demonstrate the existence of an intermediate B19 phase between martensite and austenite phases. The Nickel-Titanium shape memory alloy with precipitates is studied by introducing Ni3Ti, NiTi2 and Ni4Ti3 precipitates individually using atomistic simulations. The results show that the phase transition temperature decreases in the presence of a large volume fraction of precipitates. A blended precipitates model with Ni3Ti and NiTi2 is simulated to study the phase transformation in equiatomic NiTi alloy. The results indicate that the precipitates would initiate the emergence of the B19 phase and reduce the transition temperature. In addition, the variation of Nickel content by embedded precipitates would lead to a change in the microstructural phenomena.


Introduction
Ni-Ti shape memory alloy (SMA) is extensively used in biomedical treatments (Duerig et al 2003), automation, manufacturing and aerospace (Hartl and Lagoudas 2007) applications due to their diverse mechanical properties. The mechanism of recoverable strain in SMA is controlled by the martensitic transformation with applied stress or change in temperature (Wei et al 1998). The early study by Miyazaki and Otsuka (1986) shows that SMA properties could be affected by heat treatment, alloy composition, precipitation and substitution of other elements, etc. The present work focuses on the factors influencing Ni-Ti SMA properties during athermal and stress induced martensitic transformation.
The primary phase transformation which leads to the shape memory effect is between the austenite and martensite phases. The formation of various phases is a manifestation of different atomic structures. Therefore, to understand the mechanism of the phase transformation, it is necessary to generate a configuration of phases at the nanoscale. Molecular dynamic (MD) simulations have been successfully used in the investigation of phase transformations at the nanoscale (Allen 2004). Ishida and Hiwatari (2007) demonstrated the formation of B19 phase during martensitic transformation and changes in the transformation temperature with Ni concentration. Zhong et al (2011) and Chen et al (2018) studied the stress-induced martensitic transformation and phase recovery phenomenon of B2 and B19′ phase in Ni-Ti SMA at the nanoscale. MD simulations using a high accuracy interatomic potential, derived from Lai and Liu (2000), showed the formation of B19 phase in the middle of martensitic transformation (Ko et al 2015). Additionally, Ling and Roy (1981) used electrical resistance (ER) value based on the temperature to define the middle phase as the R phase, followed by other research works (Shamimi et al 2018, Zheng et al 2008, Zhang and Sehitoglu 2004, Pelosin and Riviere 1998) that observed R phase in both cooling and heating process by using differential scanning calorimetry (DSC). However, Sewak and Dey (2019) had recognized the intermediate phase as the B19 phase instead of R phase, according to the crystal structure parameter and the perturbed angular correlation (PAC).

Methodology
The simulation domain used to simulate phase transformation is shown in figure 1(a) with 30 lattices on each side. Periodic boundary conditions are considered and time step, Δt = 1fs. The simulation is initiated with a B2 (BCC) structure, where green and grey atoms denote Ni and Ti, respectively. In the simulation, temperature and pressure are controlled by the settings of individual thermostat and barostat in LAMMPS (Plimpton 1995), which is suitable for the study of phase transformation in MD method (Chen et al 2018). In this work, the NiTi alloy model is simulated using an optimised second nearest-neighbour modified embedded-atom method (2NN MEAM) potential proposed by Ko et al (2015). In the MEA method (Baskes and Johnson 1994), the basic total energy equation is described as: Where, r ī is the background electron density, F i is embedded energy, and f R ij ij ( ) is the pair interaction between atoms i and j with range R ij which is a value computed from the total energy and the embedding function. The 2NN-MEAM uses a modified background electron density equation as introduced by Lee and Baskes (2000).
The 2NN MEAM potential parameters for the pure Ni and Ti model optimized by Ko et al (2015) are shown in table 1 which contains the sublimation energy E c (eV/atom), the equilibrium nearest-neighbour distance r e (Å), the bulk modulus B (10 12 dyne/cm 2 ), the scaling factor for the embedding energy A, the exponential decay factors for the atomic densities β ( i) , the weighting factors for the atomic densities t ( i) and the smooth cut-off function C min and C max . These parameters are related to experimental values (Ishida and Hiwatari 2007). The potential parameters of the interaction potential of Ni-Ti originally published by Ko et al (2015) are archived in OpenKIM , Karls 2019, which are summarized in tables 1 and 2.
To investigate the athermal and stress-induced phase transformations, an equiatomic NiTi model is used in an isothermal-isobaric (NpT) ensemble command in LAMMPS (Plimpton 1995) and the Nose-Hoover thermostat and barostat are used to control the temperature and pressure by adjusting the atom velocity and position. The model settings of the equiatomic model by controlling the temperature and pressure are used as a control group, comparing the models with precipitates and different proportions of elements. A model of 30 × 30 × 30 of NiTi lattices (BCC, a = 3.015 Å) is used to simulate the athermal martensitic transformation   between 100 K and 500 K with heating and cooling rates at 1 K/ps. Before quenching, the model goes through a relaxation process for 100 ps to approach an equilibrium state. In order to observe the stress-induced phase transformation, the simulation temperature is set over the austenite start temperature (As) when the microstructure is all austenite and reach a dynamic balance which would be suitable to observe the phase transformation. Then a uniaxial compression is set along Z direction to apply loading and unloading between 0 MPa and 2 GPa at a rate of 2 MPa/ps. The stress-induced transformation model is used to investigate the interim phase behaviour in the no-precipitate model to compare with the precipitate embedded model and athermal transformation model.
In the precipitate embedded model, Ni 4 Ti 3 , NiTi 2 and Ni 3 Ti structure is considered as precipitates with an oblate spheroidal shape in the middle of a NiTi simulation cell in a treatment group, as shown in figures 2(d)-(f). The precipitate structure data is obtained from Materials Project (Jain et al 2013) and embedded in the simulation cell. The models embedded with different sizes and shapes of precipitates are for the study of athermal phase transformation with different Ni proportions. In addition, 50 at% Nickel models as shown in figure 2(c) with two precipitates are for the investigation of precipitate behaviour and exclude the effect of variation in alloy composition. Models with one kind of precipitate are set as shown in figure 2(a) including  Ni 4 Ti 3 , NiTi 2 and Ni 3 Ti with precipitate fraction from 0 at% to 25 at% and the double precipitate model with the precipitate fraction from 0 at% to 30 at%.
In this paper, a Ni 0.5 Ti 0.5 model is simulated as a control group. The treatment group of the precipitate model is studied for the effect of precipitates on the phase transformation in NiTi alloy. By involving precipitates (Ni 3 Ti, NiTi 2 and Ni 4 Ti 3 ), the atomic percentage of Ni can be modified to simulate the phenomenon with a different proportion of Nickel in a cell. An equiatomic model of Ni-Ti alloy includes Ni 3 Ti and NiTi 2 which are the structural decomposition of NiTi (Sewak and Dey 2019). Ideally, with a large-scale model involving precipitates, the NiTi alloy simulation would be more accurate to predict the transformation temperature, transformation stress, volume change and total energy data.

Results and discussion
3.1. Observation of martensite phase As the SME is related to martensite and austenite phase transformation, the material data can be gathered from a simulation where the system is cooled and heated which can reproduce the martensite and austenite phases. The transformation temperature can be obtained from the simulations. As shown in figure 3, the simulations predict the martensite start temperature (M s ), austenite start temperature (A s ), martensite finish temperature (M f ) and  austenite finish temperature (A f ) as 224 K, 425 K, 218 K and 435 K, respectively. The sudden changes shown in figure 3 indicate the occurrence of phase transformation. The change in atomic structure is shown in figures 4, (a)-(d) that show different phases in the range of 500 K and 100 K in cooling and heating respectively, where green and grey are identified as Ni and Ti atoms. The twinned structure of martensite shows consistency with the experimental data studied by Waitz et al (2005). Meanwhile, figures 4(e)-(h) are coloured based on radial distribution function (RDF) g(r) that represent different phases corresponding to figures 4(a)-(d), respectively. Based on the RDF results, the structure in red and blue is marked as austenite and martensite, respectively.
The data in the RDF plot is normalized as shown in figure 5(a) which represent B2 austenite and B19′ martensite, which are similar to the RDF data reported by Mutter and Nielaba (2010) as shown in figure 5(b). The peaks in g(r) at r = 2.6 Å, r = 3.0 Å, r = 4.25 Å and r = 5.0 Å in figure 5(b) show the average volume of distribution atoms at certain distances which represent B2 phase atomic distribution pattern. Identically, the g(r) for B19′ phase in figure 5(a) has a similar trend with the perfect structure, r = 3.0 Å and r = 4.0 Å in figure 5(b) that show a slight difference, which may be due to the thermal vibrations. Between the r = 3.0 Å and r = 4.5 Å, the martensite curve shows four peaks which are different from the perfect curve (Mutter and Nielaba 2010).  These differences could have resulted from a different potential data, which determines the atom behaviour in simulations. In spite of this, the structure identification cannot fully depend on RDF (Mutter and Nielaba 2010).
To further identify the phase structure, the lattice structure, before and after the phase transformation, is detailed in figure 6. Figure 6(a) is a lattice from 500 K and the martensite lattice at 100 K is shown in figure 6(b). A tetragonal lattice of B2 is set to correspond to a monoclinic B19′ phase (Ko et al 2015). The phase structure in figure 6(b) shows a good agreement with the B19′ phase derived by Ko et al (2015). Furthermore, this structure corresponds to one of the 12 variants of the B19′ phase (Xu et al 2019).

Observation of intermediate phase
With stress applied along the Z direction (schematic in figure 1(a)), three phases are observed in the simulation, as shown in figure 7. The total energy variation indicates that an interim B19 phase exists between austenite and martensite. The discrete jumps observed in figure 7 indicate the occurrence of phase transformation. Most of the experiments (Miyazaki et al 1984, Miyazaki and Otsuka 1986, Zhang and Sehitoglu 2004 demonstrate that the R-phase exist in the middle of the phase transformation with certain conditions. However, none of the literature on MD simulations shows the existence of the R-phase with the developed potential. By comparing the x-ray diffraction (XRD) analysis measurement between B19 and R phase, Sewak and Dey (2019) found that the intermediate phase is the B19 phase instead of the R phase. As described in the investigation by (Ko et al 2015), they observed the occurrence of the B19 phase in the middle of B2 and B19′ phase transformation using MD simulations. Therefore, by using their potential data, as shown in figure 8, the blue quadrilaterals in the atomic structure indicate the shape and the construction conversion of the lattices. The structure in figures 8 (a), (b) and (c) demonstrate the atomic shape of B2, B19 and B19′, in good agreement with the results of Ko et al (2015). The phase defined as B19 in figure 8 has orthorhombic lattices, which are different from the experimentally observed R-phase (rhombohedral) by Wang et al (2014).

Effect of precipitates on phase transformation
In the MD model with precipitates, the Ni 4 Ti 3 , Ni 3 Ti and NiTi 2 structures are embedded as precipitates in simulation cells individually. By involving one precipitate, the Nickel atomic percentage changes with the volume fraction of precipitates. In order to eliminate the effect of varying alloying element concentration, the Nickel percentage in the model is controlled by involving two kinds of precipitates. Ni 3 Ti and NiTi 2 precipitate structures are combined and adjusted in a cell to balance the Nickel concentration in the Ni 0.5 Ti 0.5 model. The simulation data of various Ni concentrations with precipitates and the resultant phase transformation temperatures are shown in tables A.1 and A.2.
Generally, the presence of precipitates will influence the phase transformation start temperature. In the cooling process, as shown in figure 9(a), with the increasing precipitates proportion, the B19′ start temperatures could vary by 150 K in the 30 at% blended precipitate model. The models with a single type of precipitate, i.e. either of Ni 3 Ti, NiTi 2 and Ni 4 Ti 3 , show a similar trend where the phase transformation start temperature decreases from 250 K to around 125 K with increasing precipitate proportion. In the heating process ( figure 9(b)), the B2 phase start temperatures for the models with single type of precipitate decrease from 420 K to 300 K with increasing precipitate proportion. Similarly, the transformation temperature for the blended precipitates model decreases from 420 K (0 at% of precipitate) to 250 K (30 at% of precipitate). In the blended precipitate model, different volume fractions of Ni 3 Ti and NiTi 2 are used to control Nickel percentage to be 50 at%. This precipitate model is simulated to investigate the precipitate behaviour in equiatomic NiTi alloy. The B19′ and B2 phase transformation start temperatures decrease with increasing precipitate proportion. On the contrary, the B19 phase transformation start temperature shows a dissimilar trend in different kinds of precipitate models. As shown in figure 10, with Ni 4 Ti 3 , NiTi 2 and Ni 3 Ti precipitates individually embedded in the models, the B19 phase (intermediate phase) is observed in the cooling simulations with precipitates over 5.6 at%. However, the B19 phase is not observed during the heating in any models which agrees with the simulation by Li et al (2021b). Under specific conditions in the experiment, the interim phase could arise with increasing temperature (Shamimi et al 2018). In the presence of Ni 3 Ti precipitate (including double precipitate model), the B19 start temperature increases with increasing precipitate proportion. In the presence of Ni 4 Ti 3 and NiTi 2 precipitates,the B19 transformation temperature do not vary significantly with increasing precipitate fractions. In figure 10, the presence of Ni 3 Ti increases B19 transformation temperature with increasing Figure 10. Effect of (a) Ni 3 Ti (b) Ni 4 Ti 3 (c) NiTi 2 (d) NiTi 2 and Ni 3 Ti precipitates on start and finish temperatures of B19 phase as well as on start temperature of B19′ phase. Figure 11. The schematic of Ni 3 Ti (in red circle) and NiTi 2 (in blue circle) in (a) the initial state of the atomic model, (b) relaxed model at 500 K for 100 ps and (c) model at 100 K with martensite phase which is indicated by the herring-bone type structure (red lines). precipitate proportion. A high precipitate volume fraction causes a deviation from the coherent relationship with the matrix phase, which leads to the emergence of the interim phase (Li et al 2021a). As shown in figures 11(b) and (c), compared to the NiTi 2 precipitate, Ni 3 Ti is more disordered and assimilate some B2 structure into disordered phase and absorb the interfacial energy. The loss of interfacial energy promotes the nucleation of the interim phase and contribute to the increase of transformation temperature. Moreover, the increased stress also restricts the transformation from B2 to B19′ phase and triggers the emergence of the B19 phase.
The precipitates impede martensite formation by restricting atomic slide on the interface. As shown in figure 11, the martensitic interface prefers to form tangential to the elliptical precipitate due to the high  interfacial energy of austenite with precipitates. As shown in figure 12(a), the B19′ forms in the bulk of austenite at a distance away from the interface between precipitate and austenite which agrees with the MD simulation by Li et al (2021b). Moreover, the B19′ formation area is also the region of maximum strain in the bulk area. Relatively, in figure 12(b), austenite nucleates at the interface of precipitate and martensite in the heating process which means higher interfacial energy can contribute to the formation of austenite. Due to the presence of precipitate, the interfacial energy impedes the nucleation of the martensite phase and leads to the nucleation of interim phase B19 instead of B19′ phase (Jiang et al 2015). Therefore, the martensite start temperature decreases with the increasing proportion of precipitate and causes the formation of interim phase.
In addition to the influence of precipitates, the variation in alloying element concentration would also alter the start temperatures of B19′ and B2 phases. As shown in figure 13, the B19′ start temperature is the highest at around 50 at% of Nickel and has a decreasing trend with increasing or decreasing Nickel atomic percentage, i.e. 150 K at 55 at% and 100 K at 47 at%. The B2 transition temperature shows the same trend as the B19′ phase with the highest transition temperature at 420 K and lowest at both sides with 310 K (Nickel 47 at%) and 320 K (Nickel 55 at%). However, the variation of Nickel proportion is instigated by the change of Nickel and Titanium microstructure from NiTi to NiTi 2 , Ni 3 Ti and Ni 4 Ti 3 . On comparing the plots of Ni 4 Ti 3 with Ni 3 Ti in figure 13, it is evident that the type of precipitate also influences the transformation temperature. With the same precipitate proportion in the simulation, the higher percentage of Nickel in precipitate results in lower transformation temperature.
In general, the increased volume fraction of precipitate would decrease the phase transformation temperature in both cooling and heating, however, the fraction of precipitates in the models is much higher than a real specimen. The phases other than B2, B19 and B19′ observed by Otsuka and Ren (1999) can be Ni 3 Ti and NiTi 2 , transferred from NiTi. Sewak and Dey (2019) observed Ni 3 Ti and NiTi 2 phases in the annealed specimen but not in PAC measurements. The divergence between experiments in different conditions and equipment is inevitable, therefore the conclusions will lead to a diverse understanding regarding the microstructural phenomena. On the other hand, the models based on potential functions may not be perfect to simulate all the occurrences in phase evolution. Due to the limitation of size and simulation time in a MD model, the simulation is more idealized in studying the phase transformation compared to experimental data.

Conclusion
In this investigation, MD simulations of athermal phase transformation in Ni-Ti alloy with embedded precipitates and equiatomic compound are demonstrated. B2, B19 and B19′ phases are observed in the simulations in the temperature range of 10 K to 500 K. The transformation from austenite to martensite shows that the lattice transforms from tetragonal to monoclinic, corresponding to the rearranged B2 cubic lattice to B19′ monoclinic lattice. The occurrence of interim phase B19 between the B2 and B19′ phases was observed with precipitate amount above 9.9 at% in the model with hybrid precipitates and above 5.5 at% in the model with single precipitate. The precipitate influences the martensite nucleation by impeding atomic movements and creating high interfacial energy which also contributes to the formation of the interim phase. In the MD simulations, the model with precipitates is consistent with experimental data, however, using precipitates to interfere with the Nickel proportion has limitations in simulating the full microstructure phenomenon.