Length scale effects and multiscale modeling of thermally induced phase transformation kinetics in NiTi SMA

MD simulations are performed employing the Finnis–Sinclair many body interatomic potential [4, 8, 11, 14, 15, 43], developed based on the embedded atom method (EAM) and considering a tight binding second moments approximation (TBSMA) theory. Quasi-static MD simulations are performed in LAMMPS [44], for temperature driven phase transformation process. Brief description of the simulation process is provided subsequently, while more details are available in recent studies by the authors [8, 15]. This EAM potential was initially proposed by Lai and Liu [43], based on the second-moment approximation of tight-binding theory (TBSMA). According to this potential, the total energy of the system over all atoms is expressed as

$$\begin{eqnarray}&&E=\sum _{i}\left\{\displaystyle \sum _{j\ne i}{A}_{\alpha \beta }\exp \left[-{p}_{\alpha \beta }\left(\displaystyle \frac{{r}_{ij}}{{d}_{\alpha \beta }}-1\right)\right]-\sqrt{\sum _{j\ne i}F({r}_{ij})}\right\}.\end{eqnarray} \tag{ 1 }$$

The first term in (1) indicates the pairwise interaction energy and the second term represents many-body effects. Recently the many-body interaction part (second term in (1)) of Lai and Liu potential has been modified by Zhong et al [4, 45], where the second term in (1) is expressed through

$$\begin{eqnarray}F({r}_{ij})=\,\left\{\begin{array}{ll}{\xi }_{\alpha \beta }^{2}\exp \left[-2{q}_{\alpha \beta }\left(\displaystyle \frac{{r}_{ij}}{{d}_{\alpha \beta }}-1\right)\right] & {r}_{ij}\leqslant {r}_{1}\\ {c}_{0,\alpha \beta }+{c}_{1,\alpha \beta }({r}_{ij}-{r}_{1})+{c}_{2,\alpha \beta }{({r}_{ij}-{r}_{1})}^{2}+{c}_{3,\alpha \beta }{({r}_{ij}-{r}_{1})}^{3} & {r}_{1}\leqslant {r}_{ij}\leqslant {r}_{c}\end{array}\right..\end{eqnarray} \tag{ 2 }$$

In (1) and (2), $i$ and $j$ denote atoms, $\alpha ,\beta$ denote the type of atoms (i.e. Ni or Ti), and ${r}_{1}$ and cut-off radius $({r}_{c})$ are adopted as 4.0 Å and 4.2 Å, respectively. Different parameters of this potential are determined through fitting the properties (i.e. lattice parameters, cohesive energy, elastic constants and vacancy formation energy) of the B2 phase at 0 K temperature, from first principles calculations [4]. More details about this EAM potential parameters are provided in Zhong et al [4].

To study length scale effects, free boundary conditions in all three directions of the simulation cell are imposed. In doing so, periodic boundary condition is used in all sides, with vacuums (empty space) all around the simulation cell are considered. In every direction, a minimum 20 Å empty space is provided, while the potential cutoff radius $({r}_{c})$ is 4.2 Å. This eliminates the possibility of interaction between surface atoms, due to the periodic boundaries and thus enables the capturing of length scale effects. Based on such boundary conditions, at constant temperature, at first NPT simulation, followed by NVT equilibration is performed. The size and thickness of the examined cells are described in the results and discussion section. Initially, pristine B2 crystal structures, i.e. 3D simulation cells with the crystallographic orientation [110], [$\bar{1}$10] and [001] are considered. For the horizontal plane, directions along thickness [110] and along one side length [$\bar{1}$10] are considered. For the vertical direction, the [001] crystal orientation is considered. The simulation process is as follows. First, a pristine austenite sample is relaxed using an energy minimization method. Then the sample is subjected to a velocity distribution (at a constant temperature, controlled via proper thermostat and barostat) and held for 360 ps, under an NPT stabilization where the systems pressure is kept constant at 0 bars. The NPT stabilized sample is then subjected to NVT conditions for an additional 240 ps. The atomic coordinates were averaged over 20 ps in order to avoid effects from thermal fluctuations. In each case of NPT or NVT stabilization, the adopted time step is ${\rm{\Delta }}t=0.0005$ ps. During simulation, the temperature and pressure are controlled by a Nose–Hoover thermostat and Berendsen barostat. It has been verified that the simulated cells are well-equilibrated and the fluctuations in the total energy and temperature are smaller than the 0.001 eV/atom and 0.50 K, respectively. Finally, the well equilibrated NVT simulations (i.e. averaged over last 20 ps) is used to report the results.

To distinguish austenite B2 phase from the different martensite B19' phase variants, the bond length based order parameter proposed by Mutter and Nielaba [8, 11, 12, 15, 46] was adopted. This order parameter has a value of −1.0 and +1.0 for pristine B2 and B19' structures, respectively. It is expressed as

$$\begin{eqnarray}&&\chi ({d}_{0},{d}_{1})=\displaystyle \frac{[{d}_{0}\{{d}^{({\rm{B}}2)}+{d}_{1}^{({\rm{B}}19^{\prime} )}\}]-[{d}_{1}\{{d}^{({\rm{B}}2)}+{d}_{0}^{({\rm{B}}19^{\prime} )}\}]}{[{d}^{({\rm{B}}2)}\{{d}_{0}^{({\rm{B}}19^{\prime} )}-{d}_{1}^{({\rm{B}}19^{\prime} )}\}]},\end{eqnarray} \tag{ 3 }$$

where bond length parameters are ${d}^{({\rm{B}}2)}=2.62\,{\mathring{\rm{A}}},\,{d}_{0}^{({\rm{B}}19^{\prime} )}=2.46\,{\mathring{\rm{A}}}\,{\rm{and}}\,{d}_{1}^{({\rm{B}}19^{\prime} )}=2.65\,{\mathring{\rm{A}}}.$ The order parameter becomes +1.0 for the martensite B19' phase when ${d}_{0}=\,{d}_{0}^{({\rm{B}}19^{\prime} )}$and ${d}_{1}=\,{d}_{1}^{({\rm{B}}19^{\prime} )}.$ For the austenite B2 phase ${d}_{0}=\,{d}_{1}=\,{d}^{({\rm{B}}2)}$ and thus the value of the order parameter becomes −1.0. However, the lattice parameters of B2 and B19' phases change due to the linear thermal expansion or contraction $(\alpha {\rm{\Delta }}T),$ and due to the internal strains $({\varepsilon }^{* })$ that develop during the phase transformation process. Thus the order parameter shows a narrow Gaussian distribution around −1.0 and +1.0 for the B2 and B19' phases, respectively [8, 15, 24]. More details about this order parameters can be found in [8, 11, 12, 15, 24], and in the supplementary material.

PhF simulations, a widely used and validated methodology [7, 16, 24], is employed to simulate the MPF evolution, based on the Ginzburg–Landau type polynomial energy function [7, 8]. In PhF simulations, austenite B2 transforms into 12 different variants of monoclinic martensite B19', with all field variables evolving stochastically over the simulation domain. The growth or evolution rate of the ith field variable ${\psi }_{i}$ can be expressed through the time dependent Ginzburg–Landau (TDGL) kinetic equation as [7, 8]

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\psi }_{i}}{\partial t}=\displaystyle \sum _{j=1}^{12}[{\bar{\vartheta }}_{ij}({G}_{j}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}+{G}_{j}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}+{G}_{j}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}})]+{\eta }_{i}({\bf{x}},t),\end{eqnarray} \tag{ 4 }$$

where ${\bar{\vartheta }}_{ij}$ is the kinetic coefficients matrix, ${G}_{j}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}};$ ${G}_{j}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}$ and ${G}_{j}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}$ represent the driving force acting on the jth field variable due to the local free energy density $({g}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}),$ gradient energy density $({{g}}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}})$ and elastic energy density $({g}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}),$ respectively, and ${\eta }_{i}({\bf{x}},t)$ is the Langevin noise term to incorporate the stochasticity in the evolution of the B19'variants.

For $i\,\ne j,$ the driving force on the ith PhF variable is considered that it does not affect the evolution of jth PhF variable, thus ${\overline{\vartheta }}_{ij}=\,\vartheta {\delta }_{ij},$ where ${\delta }_{ij}$ is the Kronecker delta. Further, from the fluctuation–dissipation theorem [7], the correlation function of the Langevin noise term in the spatiotemporal domain can be expressed as a function of temperature $\langle {\eta }_{i}({\bf{x}},t){\eta }_{i}({\bf{x}}+{\bf{s}},t+\tau )\rangle =2{k}_{{\rm{B}}}T\,{\overline{\vartheta }}_{ij}{\delta }_{ij}\delta ({\bf{s}})\,\delta (t),$ where ${k}_{{\rm{B}}}\,$is the Boltzmann constant, $T$ is simulation temperature and $\delta (\cdot )$ is the Dirac delta function. With this, (4) results into

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\psi }_{i}}{\partial t}=\,\vartheta \displaystyle \sum _{i=1}^{12}[({G}_{i}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}+{G}_{i}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}+{G}_{i}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}})]+{\eta }_{i}({\bf{x}},t).\end{eqnarray} \tag{ 5 }$$

Next, the energy density functions can be expressed as [7, 8]

$$\begin{eqnarray}&&{g}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}={g}_{0}+{\rm{\Delta }}g(T)\left\{\tfrac{A}{2}\left(\displaystyle \sum _{i=1}^{12}{\psi }_{i}^{2}\right)-\tfrac{B}{3}\left(\displaystyle \sum _{i=1}^{12}{\psi }_{i}^{3}\right)+\tfrac{C}{4}{\left(\displaystyle \sum _{i=1}^{12}{\psi }_{i}^{2}\right)}^{2}+\tfrac{D}{4}\left(\displaystyle \sum _{i=1}^{12}{\psi }_{i}^{4}\right)\right\},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{g}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}=\tfrac{1}{2}\displaystyle \sum _{i=1}^{12}{\beta }_{mn}({\psi }_{i})\left(\displaystyle \frac{\partial {\psi }_{i}}{\partial {{\bf{x}}}_{m}}\right)\left(\displaystyle \frac{\partial {\psi }_{i}}{\partial {{\bf{x}}}_{n}}\right),\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&{g}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}=\tfrac{1}{2}{({\boldsymbol{\varepsilon }}-{{\boldsymbol{\varepsilon }}}^{* })}^{T}{\boldsymbol{C}}({\boldsymbol{\varepsilon }}-{{\boldsymbol{\varepsilon }}}^{* }).\end{eqnarray} \tag{ 6c }$$

In (6a), ${g}_{0}$ represents the free energy density of the austenite phase, ${\rm{\Delta }}g(T)$ the difference in the free energy between austenite and martensite phase at temperature $T,$ and $A,\,B,\,C\,{\rm{and}}\,D$ are constants that control the shape of the free energy function. In (5b), for isotropic gradient energy coefficient ${\beta }_{mn}({\psi }_{i})$ is assumed the same for all the field variables, i.e. ${\beta }_{mn}({\psi }_{i})=\beta {\delta }_{mn}$ and thus (6b) will can be rewritten as [7, 8]

$$\begin{eqnarray}&&{g}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}=\tfrac{1}{2}\displaystyle \sum _{i=1}^{12}\beta {\delta }_{mn}\left(\displaystyle \frac{\partial {\psi }_{i}}{\partial {{\bf{x}}}_{m}}\right)\left(\displaystyle \frac{\partial {\psi }_{i}}{\partial {{\bf{x}}}_{n}}\right)=\tfrac{\beta }{2}\displaystyle \sum _{i=1}^{12}{\left(\displaystyle \frac{\partial {\psi }_{i}}{\partial {{\bf{x}}}_{m}}\right)}^{2}=\tfrac{\beta }{2}\displaystyle \sum _{i=1}^{12}{({\rm{\nabla }}{\psi }_{i})}^{2}.\end{eqnarray} \tag{ 7 }$$

In (6c), ${\boldsymbol{\varepsilon }}$ is the linearized strain tensor, i.e. ${\{{\varepsilon }_{11},{\varepsilon }_{22},{\varepsilon }_{33},2{\varepsilon }_{23},2{\varepsilon }_{31},2{\varepsilon }_{12}\}}^{T};$ ${{\boldsymbol{\varepsilon }}}^{* }$ is the nonlinear stress-free strain due local phase transformation, i.e. ${\{{\varepsilon }_{11}^{* },{\varepsilon }_{22}^{* },{\varepsilon }_{33}^{* },2{\varepsilon }_{23}^{* },2{\varepsilon }_{31}^{* },2{\varepsilon }_{12}^{* }\}}^{T},$ and ${\boldsymbol{C}}$ is the 6 × 6 isotropic elasticity tensor in the Voigt notation system [7, 8]. Elastic interactions play a critical role in the evolution of MPF, thus, the accurate estimation of the elastic energy density is important. Since in the present study no external strain is applied (i.e. ${\boldsymbol{\varepsilon }}=0$), (6c) can be written as [7, 8]

$$\begin{eqnarray}&&{g}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}=\tfrac{1}{2}{({{\boldsymbol{\varepsilon }}}^{* })}^{T}{\boldsymbol{C}}({{\boldsymbol{\varepsilon }}}^{* })=\tfrac{1}{2}{\left(\displaystyle \sum _{{i}=1}^{12}{\psi }_{{i}}\cdot {{\boldsymbol{\varepsilon }}}_{{i}}^{0}\right)}^{T}{\boldsymbol{C}}\left(\displaystyle \sum _{{i}=1}^{12}{\psi }_{{i}}\cdot {{\boldsymbol{\varepsilon }}}_{{i}}^{0}\right),\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{\varepsilon }}}_{i}^{0}$ represents the transformation strain tensor associated with the ith PhF variable [7, 8]. More details about the estimation of ${{\boldsymbol{\varepsilon }}}_{i}^{0}$ is provided in the supplementary material. Here, ${{\boldsymbol{\varepsilon }}}^{* }$ represents the stress-free strain (i.e. internal strain tensor in B19'), which arises due to local MPT from B2 to B19', expressed as

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}^{* }=\displaystyle \sum _{i=1}^{12}{\psi }_{i}\cdot {{\boldsymbol{\varepsilon }}}_{i}^{0}.\end{eqnarray} \tag{ 9 }$$

Therefore, the associated driving forces for the evolution of jth field variable can be described through the following set of equations

$$\begin{eqnarray}&&{G}_{i}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}=-\displaystyle \frac{\partial {g}^{{\rm{l}}{\rm{o}}{\rm{c}}{\rm{a}}{\rm{l}}}}{\partial {\psi }_{i}}={\rm{\Delta }}g(T)\left\{-A{\psi }_{i}+B{\psi }_{i}^{2}-C{\psi }_{i}\left(\displaystyle \sum _{j=1}^{12}{\psi }_{j}^{2}\right)-D{\psi }_{i}^{3}\right\},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{G}_{i}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}=\displaystyle \frac{\partial {g}^{{\rm{g}}{\rm{r}}{\rm{a}}{\rm{d}}{\rm{i}}{\rm{e}}{\rm{n}}{\rm{t}}}}{\partial {\psi }_{i}}=\beta ({{\rm{\nabla }}}^{2}{\psi }_{i}),\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&{G}_{i}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}=-\displaystyle \frac{\partial {g}^{{\rm{e}}{\rm{l}}{\rm{a}}{\rm{s}}{\rm{t}}{\rm{i}}{\rm{c}}}}{\partial {\psi }_{i}}=-\left[\displaystyle \sum _{i=1}^{12}{\psi }_{i}\cdot {({{\boldsymbol{\varepsilon }}}_{i}^{0})}^{T}\right]\cdot ({{\boldsymbol{C}}}^{T}{{\boldsymbol{\varepsilon }}}_{i}^{0})={F}_{i}^{{\rm{e}}{\rm{l}}}.\end{eqnarray} \tag{ 10c }$$

With all the terms in TDGL kinetic equation in (5) obtained, the stochastic spatiotemporal evolution of any field variable (e.g. the jth field variable) can be expressed as [7, 8]

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\psi }_{i}}{\partial t}=\vartheta \left[{\rm{\Delta }}g(T)\left\{-A{\psi }_{i}+B{\psi }_{i}^{2}-C{\psi }_{i}\left(\displaystyle \sum _{j=1}^{12}{\psi }_{j}^{2}\right)-D{\psi }_{i}^{3}\right\}+\beta ({{\rm{\nabla }}}^{2}{\psi }_{i})+{F}_{i}^{{\rm{e}}{\rm{l}}}\right]+{\eta }_{i}({\bf{x}},t).\end{eqnarray} \tag{ 11 }$$

Equation (11) is easily normalized in terms of length and time scale. By employing evenly spaced grid size $({\rm{\Delta }}{\bf{x}})$ and time steps $({\rm{\Delta }}t)$ in (11), it follows that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\psi }_{i}}{\partial \tilde{t}}=\left[\left\{-A{\psi }_{i}+B{\psi }_{i}^{2}-C{\psi }_{i}\left(\displaystyle \sum _{j=1}^{12}{\psi }_{j}^{2}\right)-D{\psi }_{i}^{3}\right\}+\tilde{\beta }({\tilde{{\rm{\nabla }}}}^{2}{\psi }_{i})+\displaystyle \frac{{F}_{i}^{{\rm{e}}{\rm{l}}}}{{\rm{\Delta }}g(T)}\right]+{\tilde{\eta }}_{i}(\tilde{{\bf{x}}},\tilde{t}),\end{eqnarray} \tag{ 12 }$$

where dimensionless time and time steps are expressed in terms of real time $t$ and time steps ${\rm{\Delta }}t,$ i.e. $\,\tilde{t}=t\vartheta {\rm{\Delta }}g(T)$ and ${\rm{\Delta }}\tilde{t}={\rm{\Delta }}t\vartheta {\rm{\Delta }}g(T),$ and the dimensionless spatial coordinate system is expressed in terms of the grid size ${\rm{\Delta }}x$ i.e. $\tilde{{\bf{x}}}={\bf{x}}/{\rm{\Delta }}x,$ and $\tilde{\beta }=\beta /{\rm{\Delta }}g(T){({\rm{\Delta }}x)}^{2}$ denotes the normalized isotropic gradient energy coefficient. Finally, the normalized Langevin noise term ${\tilde{\eta }}_{i}(\tilde{{\bf{x}}},\tilde{t})$ is generated as a random number from the spatiotemporally uncorrelated normal distribution with zero mean and variance of $2{k}_{{\rm{B}}}T/{\rm{\Delta }}g(T){({\rm{\Delta }}x)}^{3}.$

For the simulation process, initially 2D B2 simulation cells of different sizes (from 32 nm to 2048 nm) are considered. The crystallographic orientation of the B2 simulation cells are adopted as [110] (along x direction) and [001] along (y direction) and then the simulation cells are stabilized at constant temperature, as detailed subsequently. The spatial and temporal resolutions are adopted as ${\rm{\Delta }}{\bf{x}}=1\,{\rm{nm}}$ and ${\rm{\Delta }}t=0.001\,{\rm{ns}},$ respectively, and simulations are performed for a total time of $20\,{\rm{ns}},$ out of which, for the first $5\,{\rm{ns}},$ the simulation is perturbed by uncorrelated Langevin noise. Convergence of the simulation process is ensured by comparing different statistical measures of the MPF between two consecutive iteration steps, i.e. iteration is stopped once ${{\rm{\Delta }}}_{Y}\,\leqslant {\delta }_{Y}.$ Here, for any variable $Y,$ the difference in value between two consecutive iteration steps is expressed as ${{\rm{\Delta }}}_{Y}={\rm{abs}}\{({Y}_{n+1}-{Y}_{n})/{Y}_{n}\}$ and the tolerance limit is adopted as ${\delta }_{Y}=0.001.$ Further details about the PhF simulation process is provided in [8].

Heat transfer in the simulation domain is modeled with the LBM that has received significant attention in simulation of transport processes [47–53]. Considering an isothermal process for heat transfer, Peng's implementation of thermal LBM [48] is a simplified form of He et al [47] thermal LBM. It considers that heat diffuses through the NiTi medium, and heat is dissipated during phase transformation. For the transient case, the heat diffusion equation with energy dissipation can be expressed as

$$\begin{eqnarray}&&\tfrac{\partial T}{\partial t}=\left(\tfrac{1}{\rho {c}_{p}}\right)[\{{\rm{\nabla }}\cdot (k{\rm{\nabla }}\cdot T)\}-\dot{Q}],\end{eqnarray} \tag{ 13 }$$

where $\rho ,$ ${c}_{p},$ $k,$ and $\dot{Q}$ are the density, thermal conductivity, heat capacity and heat energy dissipation, respectively.

Here, the LBM solves the generalized heat diffusion equation via propagation and collision steps, employing specific so-called elements [8, 51–53]. In the present study, the D2Q9 lattice configuration is used for the transport process [8, 51–53]. The temperature dependent thermal conductivity, heat capacity and heat energy dissipation are obtained from different experimental studies [54–56] and further details are available in [8]. In the LBM framework for heat transfer through a solid, all velocity components are zero, i.e. ${\bf{u}}=0.$ It is considered that during the heat diffusion process, heat energy dissipates due to phase transformation, thus, for the transient condition, the generalized heat diffusion equation solved by LBM contains a heat source term. In two-dimensions (2D), both mass and heat transport can be described by the lattice Boltzmann equation, expressed as [8, 47–53]

$$\begin{eqnarray}&&{f}_{\alpha }({\bf{x}}+{{\bf{e}}}_{\alpha }{\delta }_{t},t+{\delta }_{t})=\left(1-\tfrac{1}{{\tau }_{f}}\right){f}_{\alpha }({\bf{x}},t)+\left(\tfrac{1}{{\tau }_{f}}\right){f}_{\alpha }^{{\rm{e}}{\rm{q}}}(\rho ,{\bf{u}}),\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{g}_{\alpha }({\bf{x}}+{{\bf{e}}}_{\alpha }{\delta }_{t},t+{\delta }_{t})=\left(1-\tfrac{1}{{\tau }_{g}}\right){g}_{\alpha }({\bf{x}},t)+\left(\tfrac{1}{{\tau }_{g}}\right){g}_{\alpha }^{{\rm{e}}{\rm{q}}}(T,{\bf{u}})-{\varpi }_{\alpha }\left(1-\tfrac{1}{2{\tau }_{g}}\right)\dot{Q}\left(\tfrac{1}{\rho {c}_{p}}\right),\end{eqnarray} \tag{ 14b }$$

where ${\bf{x}}$ the space coordinates or position vector; ${{\bf{e}}}_{\alpha }$ the discrete velocity components, ${\bf{u}}$ the velocity vector; ${\delta }_{x}={\delta }_{y}$ the lattice constant (grid size)/spatial resolution; ${\delta }_{t}$ the time step/temporal resolution; ${\tau }_{f}=\left(\tfrac{3\upsilon }{c{\delta }_{x}}+\tfrac{1}{2}\right)$ and ${\tau }_{g}=\left(\tfrac{3k}{2c{\delta }_{x}\rho {c}_{p}}+\tfrac{1}{2}\right)\,$are the dimensionless relaxation time for mass transport and heat energy transfer, respectively, $\upsilon$ is the kinematic viscosity of the medium, $c$ the speed in the lattice that depends on the spatial and temporal resolution, expressed as $c={\delta }_{x}/{\delta }_{t}$ and is related to the speed of sound in the lattice through ${c}_{s}^{2}=RT=\tfrac{{c}^{2}}{3},$ where $R$ is the universal gas constant and $T$ is absolute temperature of the system. Further, ${f}_{\alpha }$ and ${g}_{\alpha }$ denote the density and heat energy distribution function, respectively, ${f}_{\alpha }^{{\rm{e}}{\rm{q}}}$ and ${g}_{\alpha }^{{\rm{e}}{\rm{q}}}$ denote the equilibrium density and the heat energy distribution function, respectively, and ${{\bf{e}}}_{\alpha }\,$is the discrete velocities components for the D2Q9 element (as shown in figure S4 of supplementary material available online at stacks.iop.org/MSMS/25/045002/mmedia). More details about those quantities is provided in the supplementary material.

Finally, the density, velocity and temperature at any node are calculated by [8]

$$\begin{eqnarray}&&\rho ({\bf{x}},t)=\sum _{\alpha }{f}_{\alpha }({\bf{x}},t),\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}&&{\bf{u}}=\tfrac{\displaystyle {\sum }_{\alpha }{{\bf{e}}}_{\alpha }{f}_{\alpha }({\bf{x}},t)}{\displaystyle {\sum }_{\alpha }{f}_{\alpha }({\bf{x}},t)},\end{eqnarray} \tag{ 15b }$$

$$\begin{eqnarray}&&T({\bf{x}},t)=\sum _{\alpha }{g}_{\alpha }({\bf{x}},t)-\dot{Q}\left(\displaystyle \frac{{\delta }_{t}}{2\rho {c}_{p}}\right).\end{eqnarray} \tag{ 15c }$$

In this study an important predictive feature is added to the CWM method, mainly based on the recent study by Gur et al [8, 41] on the wavelet based surrogate time-series generation for heterogeneous catalysis. The following 1D example illustrates the methodology [8, 41].

Considering a 1D discrete data vector $y(x)$ of length $n={2}^{N},$ the WT of $y(x)$ hierarchically decomposes $y(x)$ into $N(=\,{\mathrm{log}}_{2}(n)-{2}^{2})$ number of scales and returns ${2}^{N}$ wavelet transform (WT) coefficients. The first ${2}^{N-1}$ WT coefficients represent the finest characteristics of the original signal, the next ${2}^{N-2}$ WT coefficients represent the next level (coarser) information in the data and so on. This process is well documented in the literature [8, 37, 39–41]. Thus, the WT of $y(x)$ will yield ${W}_{y}(s,x),$ i.e. a set of wavelet coefficients, where $s$ represents scale and $x$ stands for spatial coordinate in the fine scale simulation. Now, ${W}_{y}(s,x)$ is decomposed as

$$\begin{eqnarray}&&{W}_{y}(s,x)={\phi }_{y}({s}_{0},{\rm{\Delta }}x)\oplus {\phi }_{y}({s}_{1},{\rm{\Delta }}x)\oplus {\phi }_{y}({s}_{2},{2}^{1}{\rm{\Delta }}x)\oplus \cdots \cdots \oplus \,{\phi }_{y}({s}_{N},{2}^{N-1}{\rm{\Delta }}x),\end{eqnarray} \tag{ 16 }$$

where $\oplus$ denotes scale wise association in the wavelet formalism, ${\phi }_{y}({s}_{i},{2}^{i-1}{\rm{\Delta }}x),$ $i=\mathrm{1,}\,2,\cdots ,\,N$ denotes the WT at scale ${s}_{i}$ with the scaling interval ${2}^{i-1}{\rm{\Delta }}x$ and ${\phi }_{y}({s}_{0},{\rm{\Delta }}x)$ represents the WT at the coarsest scale in fine scale simulation for a particular scaling function.

Based on the WT, the multiscale coupling of MPF evolution is performed via the CWM method, as discussed subsequently. If, for thermally induced martensite phase transformation the MPF, $\xi ,$ is only available at specific and slightly overlapping wavelet space scales, e.g. $\xi (x)$ at fine MD territory scales, ${a}_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}},$ and $\xi (X)$ at coarse PhF territory scales, ${a}_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}},$ given the WT of $\xi$ at fine scales, ${W}_{\xi }({a}_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}},x),$ and that at coarse scales, ${W}_{\xi }({a}_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}},X),$ with the $x$ and $X$ denoting the location in MD and PhF domain, then the wavelet synthesis or fusion process, implies that

$$\begin{eqnarray}&&{W}_{{\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}}(a,X)={W}_{\xi }({a}_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}},x)\cup {W}_{\xi }({a}_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}},X)\end{eqnarray} \tag{ 17 }$$

i.e. the union (with proper weights for overlapping scales) of the WT at fine and coarse scales. Next, the inverse wavelet transform $({\rm{IWT}})$ is performed on the compounded wavelet coefficients, to obtain the MPF vector after compounding, thus

$$\begin{eqnarray}&&{\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}(X)={\rm{IWT}}[{W}_{{\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}}(a,X)]={\rm{IWT}}[{W}_{\xi }({a}_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}},x)\cup {W}_{\xi }({a}_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}},X)].\end{eqnarray} \tag{ 18 }$$

The CWM multiscaling process was described above for simplicity in 1D, yet its extension to 2D is straightforward. In the 2D case, both coarse and fine scale simulations are performed in 2D, and so is the forward and inverse WT.

The relevant flow chart and more details about the CWM method with the component for predicting the fine scale wavelet coefficients is provided in [8]. However, the algorithm is provided here in table 1, showing the main relevant steps of the multiscale coupling method, where ${\mathscr{W}}\,$ and ${{\mathscr{W}}}^{-1}$ denote the WTs and the IWT, respectively, and ${ {\mathcal H} }_{({s}_{1},{s}_{2})}$ and refers to the Heaviside function defined as

$$\begin{eqnarray}{ {\mathcal H} }_{({s}_{1},{s}_{2})}(s)=\left\{\begin{array}{ll}1 & {\rm{if}}\,{s}_{1}\lt s\lt {s}_{2}\,\\ 0 & {\rm{otherwise}}\,\end{array}.\right.\end{eqnarray} \tag{ 19 }$$

Also, a represents the coarsest scale resolved by the PhF simulation method, c the smallest scale resolved by the MD simulation method and b the cutoff scale, chosen based on the dominant scales resolved by the coarse and fine (by comparing the energy at the different wavelet scales). For more details about the CWM based multiscale coupling method, we refer to [8, 35–41].

Table 1.  Algorithm of multiscale coupling framework via CWM method.

1.   Given the coarse and fine martensite phase fraction (MPF): ${\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}({\boldsymbol{X}})$ with ${\boldsymbol{X}}\in [{{\boldsymbol{X}}}_{n},{{\boldsymbol{X}}}_{n}+{\rm{\Delta }}{\boldsymbol{X}}]$ and ${\xi }_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}({\boldsymbol{x}})\,$with $x\in [{{\boldsymbol{X}}}_{n},{{\boldsymbol{X}}}_{n}+{\boldsymbol{x}}];$ where ${\boldsymbol{x}}\ll {\rm{\Delta }}{\boldsymbol{X}}$

2.   Compute wavelet transforms (WT) of MPF: ${W}_{{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}}({\rm{s}},{\boldsymbol{X}})={\mathscr{W}}[{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}({\boldsymbol{X}})]$ and ${W}_{{\xi }_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}}(s,{\boldsymbol{x}})={\mathscr{W}}[{\xi }_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}({\boldsymbol{x}})]$

3.   Apply window filter on the wavelet coefficients of MPF: ${W}_{{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}}^{{ {\mathcal H} }_{({a},b)}}(s,{\boldsymbol{X}})={ {\mathcal H} }_{({a},b)}\cdot {W}_{{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}se}}(s,{\boldsymbol{X}})$ and ${W}_{{\xi }_{fine}}^{{ {\mathcal H} }_{(b,c)}}(s,{\boldsymbol{x}})={ {\mathcal H} }_{(b,c)}\cdot {W}_{{\xi }_{fine}}(s,{\boldsymbol{x}})$

4.   Compute compounding i.e. the coefficient of compound wavelet matrix (CWM): ${W}_{{\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}}(s,{\boldsymbol{X}})={W}_{{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}\oplus {\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}}^{{ {\mathcal H} }_{(a,c)}}(s,{\boldsymbol{X}})={W}_{{\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}}^{{ {\mathcal H} }_{(a,b)}}(s,{\boldsymbol{X}})\oplus {W}_{{\xi }_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}}^{{ {\mathcal H} }_{(b,c)}}(s,{\boldsymbol{x}})$

5.   Compute inverse wavelet transform (IWT) to obtain compounded MPF: ${\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}(s,X)={\rm{CWM}}({\xi }_{{\rm{c}}{\rm{o}}{\rm{a}}{\rm{r}}{\rm{s}}{\rm{e}}}({\boldsymbol{X}}),{\xi }_{{\rm{f}}{\rm{i}}{\rm{n}}{\rm{e}}}({\boldsymbol{x}}))={{\mathscr{W}}}^{-1}[{W}_{{\xi }_{{\rm{c}}{\rm{w}}{\rm{m}}}}(s,{\boldsymbol{X}})]$

Phenomenological models based on extensive continuum level experimental data on thermally induced martensitic phase fraction have been proposed [8, 15, 42, 57–59]. Experiments are reported on annealed (at a temperature of 750 °C) polycrystalline austenite NiTi SMA (mainly [110] B2 crystal with weak [$\bar{1}$11] B2 crystal), with cooling and heating rate of 20 °C min⁻¹, under various concentrations of Ni and Ti [42]. The model considered herein, essentially captures the kinetics of phase transformation rate as a function of phase fraction. It is assumed that, (a) the austenite transformation rate at any given temperature is proportional to the nucleation rate of austenite per unit volume and (b) the size of the austenite nucleus increases linearly with the austenite phase fraction. According to the model, the temperature dependent austenite phase fraction evolution is described through the Richards equation [42], i.e.

$$\begin{eqnarray}&&\xi (T)=1-{[1+\exp \{-g\cdot \upsilon (T-{T}_{m})\}]}^{(-1/\upsilon )},\end{eqnarray} \tag{ 20 }$$

where $\xi (T)$ is the temperature dependent austenite phase fraction, $T$ the temperature of interest in K, ${T}_{m}({\rm{in}}\,{\rm{K}})$ the temperature corresponding to the maximum transformation rate, while $g({\rm{in}}\,{{\rm{K}}}^{-1})$ and $\upsilon$ (unitless) are the fit parameters which dictate the transformation or growth rate of the austenite phase. Equation (20) is adopted here to represent the martensitic phase fraction evolution as a function of temperature at coarse scales or continuum level. More details about this phenomenological model are available in [8, 15, 42, 57–59].

The shape memory effect and super-elasticity are two paramount properties of SMAs, controlled by the thermo-mechanically induced austenite to martensite (and vice versa) phase transformation process. Above the austenite finish temperature ${A}_{f}$ the SMA is in the austenite phase and upon loading it transforms into martensite and during unloading the martensite back transforms into austenite. Below ${A}_{f},$ depending on the temperature, the SMA leaves residual strain upon unloading which can be recovered by heating above ${A}_{f}.$ In the process of loading–unloading, the SMA dissipates a large amount of energy through its flag-shaped hysteresis loop which again depends on the temperature and strain rate. Here the SMA material model of Helm and Haupt [60] is adopted, which has been verified experimentally and captures all essential aspects of the SMA behavior. The transformation temperatures considered here are those of NiTi SMA: (i) for the martensite start ${M}_{s}$ and finish ${M}_{f}$ temperatures are 295 and 275 K, respectively; (ii) for the austenite start ${A}_{s}$ and finish ${A}_{f}$ temperatures are 310 and 335 K, respectively. For uniaxial stress, it is expressed as

$$\begin{eqnarray}&&\sigma ={ {\mathcal F} }_{1}(\varepsilon ,{\varepsilon }^{* },T),\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&{\dot{\varepsilon }}^{* }={ {\mathcal F} }_{2}(\sigma ,{\sigma }^{* },\xi ,T),\end{eqnarray} \tag{ 21b }$$

$$\begin{eqnarray}&&{\sigma }^{* }={ {\mathcal F} }_{3}(\sigma ,{\sigma }^{* },\xi ,\dot{\xi },T),\end{eqnarray} \tag{ 21c }$$

$$\begin{eqnarray}&&\xi ={ {\mathcal F} }_{4}({\varepsilon }^{* },{\dot{\varepsilon }}^{* },T),\end{eqnarray} \tag{ 21d }$$

$$\begin{eqnarray}&&\dot{\theta }={ {\mathcal F} }_{5}(\sigma ,\varepsilon ,{\dot{\sigma }}^{* },{\dot{\varepsilon }}^{* },\xi ,T),\end{eqnarray} \tag{ 21e }$$

where, $\sigma$ is the developed stress, which depend on applied strain $\varepsilon ,$ ambient temperature $T,$ developed internal strain is ${\varepsilon }^{* }.$ Furthermore, ${\dot{\varepsilon }}^{* }$ is the developed internal strain rate, ${\sigma }^{* }$ and ${\dot{\sigma }}^{* }$ are the developed internal stress and its rate, $\xi$ and $\dot{\xi }$ are the MPF and its time rate and finally $\dot{T}$ is the time rate of temperature evolution. Thus, the material model is expressed as a system of strongly nonlinear equations that need to be solved iteratively at every time step. For detailed description of the SMA material model and the iterative solution algorithm of those system of nonlinear equations i.e. equation (23). (a)–(e), the reader is referred elsewhere [8, 46, 60]. The hysteretic energy dissipation capacity of SMAs and its temperature modulation is addressed below.

Figure 1 shows MD simulations results on the effects of side length $(s)$ and thickness $(t)$ of NiTi thin samples on the thermally induced MPF. To demonstrate the effect of side length $s,$ the thickness of NiTi thin samples are kept constant i.e. for thickness $t=8\,{\rm{nm}}.$ To show the effect of thickness $t,$ the side of the NiTi thin samples is given a constant value i.e. $s=8\,{\rm{nm}}.$ For each case, both of the austenitic (${M}\to {A},$ i.e. heating cycle) and martensitic (${A}\to {M},$ i.e. cooling cycle) thermally induced phase transformation process is shown, where, M and A are used to abbreviate the martensite and austenite phase, respectively. Figure 1(a) shows the effect of side length of NiTi thin samples on the temperature driven MPF curves for $s=\mathrm{2,16}\,{\rm{and}}\,64\,{\rm{nm}}.$ Similarly, figure 1(b) shows the effect of NiTi sample thickness on MPF evolution, for thickness $t=2,\,16\,{\rm{and}}\,64\,{\rm{nm}}$ with constant side length. It can be observed that during austenitic phase transformation, the effect of the simulation cell side length is almost the same for both the austenite start $({A}_{s})$ and finish $({A}_{f})$ temperatures, while in the MPT process, this effect is more pronounced in the martensite finish temperature $({M}_{f})$ than in the start temperature $({M}_{s}).$ However, change in sample thickness causes a similar level of shift in the transformation temperatures (i.e. austenite and martensite start and finish).

Zoom In Zoom Out Reset image size

Size effects are shown in figure 2, where both the side length and thickness of NiTi thin samples are changed proportionally, i.e. $s=t$ and all three dimensions are increased proportionally from 2 to 64 nm; both heating and cooling cycles (i.e. austenitic and MPT process) are shown. Figure 2(a) shows the thermally induced MPF evolution curves for $s=t=2,\,4,\,8\,{\rm{nm}},$ and figure 2(b) shows the same curves for $s=t=16,\,32,\,64\,{\rm{nm}}.$ It can be observed, at any temperature, with increasing size the MPF value decreases. Furthermore, the austenite start $({A}_{s})$ and martensite finish $({M}_{f})$ temperatures show more sensitivity to size than the austenite finish $({A}_{f})$ and martensite start $({M}_{s})$ temperatures. A discussion about the shape or the size effects in terms of the thermodynamics and available free surfaces is provided subsequently.

Zoom In Zoom Out Reset image size

As discussed here, the length scale effects are attributed to the presence of the free surfaces, which are defects. The available free surface area per unit volume can be expressed as $\gamma =A/V=(2/t)+(4/s).$ Thus for any fixed thickness $t$or side length $s,$ the ratio of available free surface area per unit volume is $\gamma \propto 1/s$ or $\gamma \propto 1/t,$ with increasing $s$ or $t$ this ratio quantity decreases inversely proportional to $\sqrt{s}$ or $\sqrt{t},$ i.e. ${\rm{d}}\gamma /{\rm{d}}s\propto -1/{s}^{2}$ or ${\rm{d}}\gamma /{\rm{d}}s\propto -1/{t}^{2}.$ The presence of free surfaces changes the free energy density (energy per unit volume) of phase transformation. Assuming that energy reduces due to the presence of free surface area as ${e}_{s}$ (per unit surface), the drop in the free energy density is $\{(2/t)+(4/s)\}{e}_{s}$ (per unit volume). Thus, increase in $s$ or $t,$ decreases the free energy density, which finally increases the total energy of the system. Therefore, to obtain similar level of phase transformation, more driving force is required with increasing $s$ or $t.$ In terms of thermodynamics, with increasing $s$ or $t,$ more temperature drop is required to obtain similar level of phase fraction. This ultimately shifts the phase evolution kinetics curves more towards low temperature and reduces the transformation temperature with increasing $s$ or $t.$

Similarly to temperature dependent phase fraction evolution, length scales have a strong effect on the transformation temperatures (i.e. $\,{A}_{s},$ ${A}_{f},$ ${T}_{mA},$ ${M}_{s},$ ${M}_{f}$ and ${T}_{mM}$), as shown in figure 3. Figures 3(a1)–(c1) show the transformation temperatures associated with the austenitic (heating cycle) phase transformation process, while figures 3(a2)–(c2) shows the ones associated with the MPT process i.e. for cooling cycle. In every case, with increasing side length and thickness of the NiTi thin sample, the transformation temperature decreases. Figure 3(d) shows the effect of NiTi thin samples thickness and side length on the width of thermal hysteresis between the austenitic and MPT process, which represent the difference of temperature between the heating and cooling cycle for the same MPF. Interestingly, the width of the thermal cycle initially decreases with increasing thickness for samples of small side length. However, with increasing side length, this trend changes, i.e. the width of the thermal cycle initially decreases and after a critical thickness it increases. In every case, solid lines represent the shape effect, where only the side length changes while keeping the thickness constant. Whereas, the dotted line represents size effects, where both the side length and thickness of thin samples change. A nonlinear least-square minimization process on the data in figures 3(a1)–(c2) (through MATLAB 2014) shows a power law variation for the transformation temperatures with respect to side length. Thus, all of the phase transformation temperatures decrease with increasing side length, following a certain power law and the power law constant and exponent shows logarithmic variation with thickness, i.e. they decrease logarithmically with increasing thickness. In general, the transformation temperature model can be expressed as (i.e. as a function of $s$ and $t$)

$$\begin{eqnarray}&&\displaystyle \frac{{\theta }_{T}(s,t)}{{\theta }_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}}={\delta }_{c}\,\mathrm{ln}\left(\displaystyle \frac{{t}_{c}}{t}\right){[s]}^{-\left\{{\delta }_{e}\mathrm{ln}\left(\displaystyle \frac{{t}_{e}}{t}\right)\right\}},\end{eqnarray} \tag{ 22 }$$

where ${\delta }_{c}$ (in nm) and ${t}_{c}$ (in nm) are constants, respectively, and ${\delta }_{e}$ (unitless) and ${t}_{e}$ (in nm) are constant for the power law exponent. Here, ${\theta }_{T}(s,t)$ represents the shape and size dependency of the transformation temperatures i.e. ${A}_{s},$ ${A}_{f},$ ${T}_{mA},$ ${M}_{s},$ ${M}_{f}$ and ${T}_{mM};$ and ${\theta }_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}$ represent the transformation temperatures associated with the bulk.

Zoom In Zoom Out Reset image size

A similar power law model of phase transformation temperature as the function of size has been proposed previously for other cases of nanostructures (such as for nanoparticles) [11, 13]. Those models, and the one proposed here, are only applicable at certain range of length scales, because any size below a critical length scale does not show phase transformation, as reported in [9–11, 13, 17–19]. In this context, Fu et al [9] reported that the growth of the surface oxide and its interfacial diffusion exert dominant constraining effects that produce high residual stress and low recovery capability. Experimental studies by Waitz et al [10, 17–19] show the effect of grain size and interaction of grain boundaries on the martensite phase transformation process due to the presence of nanograins for polycrystalline samples. The presence of grain boundaries and the interaction between grains play an important role in the evolution (nucleation and growth) and stabilization of different phases (B2 and B19' variants). In contrast to those studies, the present study, mainly shows results for single crystals (with crystal orientation of [110], [$\bar{1}$10], [001]) and thus any effect due to the residual stress or due the interaction between different crystals is not considered. However, to address the effect of grain size, grain boundaries, interaction between different grains and with the simulation cell size, PhF simulations for polycrystalline NiTi are performed. Also, an indirect method (based on the multiscale approach) is proposed to mimic the phase fraction evolution kinetics in polycrystalline NiTi. Here, the simulation results from different size single crystals i.e. the wavelet coefficients of MPF, for crystal orientation of [100], [110] and [111] are used in a unique way to reproduce the effects of polycrystals. Details about this are provided in section 4.4, for the polycrystalline NiTi case.

Therefore, it is important to note that, in the present study, the reported size effects are due to the presence of free surfaces, which act as defects and trigger phase transformation. With decreasing simulation cell size, the surface area per unit volume increases, and this creates the dependency of the phase transformation process (transformation temperatures and phase fraction) on size. Mutter and Nielaba [11] and Zhang et al [13] demonstrate the phase transformation process for single crystal nano-scale NiTi simulation cells, by employing an EAM potential similar to the one used in this study. Importantly, in the Zhang et al [13], it is found that below size (diameter) of 1.15 nm, no phase transformation occurs. In this study, the minimum size or side length of any simulation cell is adopted as 2.00 nm, which is above 1.15 nm. Finally, in the previous studies mentioned above, experiments and/or simulations are performed at constant temperature rates, which can play an important role in the phase transformation process, whereas in the present study quasi-static simulations (equilibrations) are performed at constant temperature. Thus, as potential future work, it will be very interesting to verify the model by employing experimental or simulation data for polycrystalline NiTi samples and obtain the model parameters.

The evolution of martensite phase faction and the development of internal strain in B19' due to the evolving martensite phase are the two most important quantities that govern the NiTi material model. Thus, it is important to accurately study the evolution of these two internal variables and transfer all relevant information to the continuum scales via multiscale coupling.

Figure 4 shows the length scale effects at the mesoscale and their linkage to microscale i.e. atomistic simulation results. At the mesoscale, 2D stochastic PhF simulation are performed and are linked to 3D atomistic (MD) simulations. For all MD simulations, the thickness $t=2\,{\rm{nm}},$ and a range of side lengths, $s=2\,{\rm{to}}\,128\,{\rm{nm}},$ is considered, with free boundary conditions in all directions. At the mesoscale level, the range of simulation cell sizes is from $s=32\,{\rm{to}}\,2048\,{\rm{nm}},$ with free boundary conditions in all directions (x, y and z). Therefore, the side length cases of $s=32\,{\rm{to}}\,128\,{\rm{nm}}$ correspond to the MD and PhF overlap domain, where the two simulation methods are cross checked, showing very good agreement with each other. The error bars in the figures 4(a)–(c) are the results of uncertainty in the evolution of different phases (B2 and B19' variants). The uncertainty is quantified by performing a total of 25 Monte Carlo (MC) simulations at each temperature and for each size. In MD simulations, depending on the temperature, random velocities (based on Gaussian distribution) are provided to each atom and thus at the end of each simulation, the atomic coordinates are not the same as those obtained from other simulations at the same temperature and same size. Similarly at the beginning of stochastic PhF simulation, depending on the temperature, a Langevin noise term is added at every material point [7, 8]. This causes small but noticeable differences in the evolution of MPF, which ultimately produces uncertainties in the development of internal strain. Therefore, to account for those uncertainties, MC simulations are performed.

Zoom In Zoom Out Reset image size

Another important parameter is the internal strain in B19' that governs the stability of different phases (i.e. B2 and different variants of B19'). Since strong length scale effects are present in the phase evolution kinetics of NiTi (at both atomic and meso scales), it can be expected that similar effects should also be present for the transformation induced internal strain. However, this is not examined in the present work. In the MD simulations, the lattice parameters (bond lengths a, b and c, monoclinic bond angle β) of the B19' unit cell phase are measured. Here, it is important to mention that, in our MD simulations, with the adopted crystal orientation of B2 (i.e. [110], [$\bar{1}$10] and [001]) only two B19' variants are observed, which are favorable due to the requirement of low twining energy with respect to B2, as also reported in other MD studies [4, 8, 13–15] employing a similar interatomic potential. It is found from our MD simulations that, for the twin B19' unit cell, the twinning energy is ${E}_{{\rm{B}}19}=-5.057\,{\rm{eV}}$ with the difference in energy (i.e. the excess energy in reference to the single phase of B19') of ${\rm{\Delta }}E=-0.035\,{\rm{eV}},$ which shows very good agreement with the values reported in the Zhong et al [4] study. Similarly, the twin boundary energy (i.e. the excess energy per unit boundary area of twin, in reference to the single phase of B19') is also estimated and found that its value is very small compared to the other B19' variants for the twining case, i.e. ${\gamma }_{t}=-0.009\,{\rm{J}}\,{{\rm{m}}}^{-2}$ [4]. Finally, all of this indicates the possible reason for evolution of those two B19' variants, in our, as well as in other MD studies. Then, the deformation gradient matrix $({{\bf{F}}}_{i})$ for the martensitic (B2 to any variants of B19') phase transformation is estimated based on the lattice parameters of the B19' unit cell phase and of the pristine B2 phase (a₀). Since the lattice parameters of the B19' phase change with length scale $(s)$ and temperatures $(T),$ ${{\bf{F}}}_{i}$ becomes a function of $s$ and $T.$ More information about ${{\bf{F}}}_{i}$ for B2 → B19' phase transformation is available in previous studies [4, 25]. Next, to determine the transformation induced strain tensor of different variants of the B19' phase, the Green–Lagrange strain tensor $({{\boldsymbol{\varepsilon }}}_{i}^{0})$ for each of the variants is estimated as

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}_{i}^{0}(s,T)=\tfrac{1}{2}[{\{{{\bf{F}}}_{i}(s,T)\}}^{T}\{{{\bf{F}}}_{i}(s,T)\}-{\bf{I}}].\end{eqnarray} \tag{ 23 }$$

Finally, knowing the fraction of B19' variants in the MD simulation cell (determined from the monoclinic bond angle β), at every temperature and for every size, the internal strain tensor in the B19' phase is estimated through

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}^{* }=\displaystyle \sum _{i=1}^{2}{\psi }_{i}\cdot {{\boldsymbol{\varepsilon }}}_{i}^{0}=\tfrac{1}{2}\displaystyle \sum _{i=1}^{2}{\psi }_{i}\cdot [{\{{{\bf{F}}}_{i}(s,T)\}}^{T}\{{{\bf{F}}}_{i}(s,T)\}-{\bf{I}}].\end{eqnarray} \tag{ 24 }$$

Now, the mean value of the internal strain tensor of the B19' phase $\langle {{\boldsymbol{\varepsilon }}}^{* }\rangle$ and associated uncertainty (standard deviation) are estimated from 25 MC MD simulations. Length scale effects in the average internal strain in the B19' phase is due to the differential fraction evolution among the different B19' variants. More information about the internal strain is provided in supplementary material. In PhF simulations, at any material point the internal strain associated with any B19' variant (${{\boldsymbol{\varepsilon }}}_{i}^{0}$ for the ith field variable) [7, 26, 27] and its fraction (${\psi }_{i}$ for the ith field variable), is used to estimate the total internal strain in B19' phase, as expressed trough ${{\boldsymbol{\varepsilon }}}^{* }=\displaystyle {\sum }_{i=1}^{12}{\psi }_{i}{{\boldsymbol{\varepsilon }}}_{i}^{0}$ (for the 12 variants). Here it is important to mention that, although during the simulation process the possibility of evolution of all 12 B19' variants is considered, only four B19' variants (1, 2, 6 and 8) are observed. It is due to the consideration of 2D simulation cell, as discussed in [8]. Interestingly, it is found here that for both MD and PhF simulations, the mean internal strain tensor $({{\boldsymbol{\varepsilon }}}^{* })$ is orthorhombic in nature. The mean value $\langle {{\boldsymbol{\varepsilon }}}^{* }\rangle =\left\{\left(\displaystyle {\sum }_{k=1}^{M}\displaystyle {\sum }_{l=1}^{N}{\varepsilon }^{* }\right)/(M\,* \,N)\right\}$ and its associated uncertainty is measured based on the number of MC simulations (here $M$ 25) and number of material points $(N)$ in each simulation cell. Since the fractional evolution of B19' variants depends on the temperature and size of the simulation cell, the evolution of the average internal strain in the B19' phase is a function of temperature and simulation cell size. Next, from the mean internal strain tensors $\langle {{\boldsymbol{\varepsilon }}}^{* }\rangle$ in the B19' phase, strains in three orthonormal directions (i.e. $\langle {\varepsilon }_{11}^{* }\rangle ,$ $\langle {\varepsilon }_{22}^{* }\rangle$ and $\langle {\varepsilon }_{33}^{* }\rangle$) are shown in figure 5. This provides the direct estimation of temperature and size dependent lattice strains that B2 experiences during the MPT process. A recent relevant study by Ahadi and Sun [23] provides similar measures of lattice strains (in terms of stretches) and demonstrates the effect of grain size and applied strain on the evolution of lattice strain in NiTi SMA. As shown in figure 5, strong interplay between the average internal strain in B19' and the length scale of the simulation cell is present. With increasing length scale, a strong decay is observed in all components of the maximum value of the average internal strain in B19' i.e. ${\langle {\varepsilon }_{ii}^{* }\rangle }_{\max }\propto 1/f(s)\,.$ It is clear that during the martensitic transformation process, the $[110]$ and $[001]$ directions of B2 experience compressive strain, while the $[\bar{1}10]$ of B2 experiences tensile strain (as shown in figure 5(c)). Interestingly, for any length scale, the temperature dependency of the normalized $\langle {\varepsilon }_{33}^{* }(s,T)\rangle$ in B19' (normalized with respect to the maximum value i.e. ${\langle {\varepsilon }_{33}^{* }(s)\rangle }_{\max }$) follows the same Richard's equation as the one for the evolution kinetics of the MPF, i.e.

$$\begin{eqnarray}&&\xi (s,T)=\displaystyle \frac{\langle {\varepsilon }_{33}^{* }(s,T)\rangle }{{\langle {\varepsilon }_{33}^{* }(s)\rangle }_{{\rm{\max }}}}=1-{[1+\exp \{-g\cdot \upsilon (T-{T}_{m})\}]}^{(-1/\upsilon )}.\end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

The other components of internal strain tensor i.e. $\langle \,{\varepsilon }_{11}^{* }\rangle$ and $\langle {\varepsilon }_{22}^{* }\rangle$ does not follow a relation similar to equation (25). The relevant direct relation between the MPF and the internal strain tensor component ${\varepsilon }_{33}^{* },$ implies that ${\varepsilon }_{33}^{* }$ is mainly responsible for the stabilization of the B19' variants and B2. Using a fitting process, it is found that the different parameters of the evolution equation (equation (25)) i.e. $\,g,$ $\upsilon$ and ${T}_{m}$ vary with the length and can be expressed as

$$\begin{eqnarray}&&g(s)={g}_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}{\left(\displaystyle \frac{{\delta }_{g}}{s}\right)}^{{\beta }_{g}},\end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray}&&\upsilon (s)={\upsilon }_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}{\left(\displaystyle \frac{s}{{\delta }_{\upsilon }}\right)}^{{\beta }_{\upsilon }},\end{eqnarray} \tag{ 26b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{m}(s)}{{T}_{m,{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}}={\gamma }_{m}-{\beta }_{m}\,\mathrm{ln}\left(\displaystyle \frac{s}{{\delta }_{m}}\right),\end{eqnarray} \tag{ 26c }$$

where ${g}_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}\,({\rm{in}}\,{{\rm{K}}}^{-1})$ and ${\upsilon }_{{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}$ (unitless) are the constants corresponding to the bulk, ${T}_{m,{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}$ is the temperature corresponding to the maximum martensite transformation rate of the bulk, ${\beta }_{g}$ (unitless) and ${\beta }_{\upsilon }$ (unitless) are constant for the power law exponent, ${\delta }_{g}$ (in nm), ${\delta }_{\upsilon }$ (in nm) and ${\delta }_{m}$ (in nm) are constants representing the length above which the length scale effect on $g,$ $\upsilon$ and ${T}_{m}$ is negligible, and ${\gamma }_{m}$ and ${\beta }_{m}$ represent the intercept and slop of the $\left(\tfrac{{T}_{m}(s)}{{T}_{m,{\rm{b}}{\rm{u}}{\rm{l}}{\rm{k}}}}\right)$ versus $\mathrm{ln}\left(\tfrac{s}{{\delta }_{m}}\right)$ line.

Since the off diagonal terms of the mean internal strain tensor are negligible, the diagonal terms furnish the principal components of the internal strain tensor. Interestingly, it is found that while the diagonal terms are only present in the mean value of the transformation induced total internal strain tensor of the B19' phase, the mean value of the total volumetric strain developed due to the MPT process (at any temperature and size) is almost zero, i.e. $\langle {\varepsilon }_{{\rm{v}}}^{* }\rangle ={\rm{tr}}(\langle {{\boldsymbol{\varepsilon }}}^{* }\rangle )=\langle {\varepsilon }_{11}^{* }\rangle +\langle {\varepsilon }_{22}^{* }\rangle +\langle {\varepsilon }_{33}^{* }\rangle \approx 0.$

Figure 6 shows the effects of length scale on the phase transformation temperatures (martensite start, maximum transformation rate, and finish in figure 6(a)) and absolute maximum value of different elements of average internal strain tensor (i.e. ${| \langle {\varepsilon }_{11}^{* }\rangle | }_{{\rm{\max }}},$ ${| \langle {\varepsilon }_{22}^{* }\rangle | }_{{\rm{\max }}}$ and ${| \langle {\varepsilon }_{33}^{* }\rangle | }_{{\rm{\max }}}$) in B19' (in figure 6(b)). As can be observed from figure 6(a), all transformation temperatures decrease with increasing length scale, following a linear trend with log of the length scale i.e. ${\theta }_{T}(s)\propto \,\mathrm{ln}(1/s).$ Thus, each of the transformation temperatures can be described through an equation similar to equation (26c), where the slope of the line is the same yet the intercepts are different. Similarly to the transformation temperatures, the absolute maximum value of each of the components of the internal strain tensor in B19' phase decreases with length scale and can be expressed as

$$\begin{eqnarray}&&{| \langle {\varepsilon }_{{i}{i}}^{* }(s)\rangle | }_{\max }={\gamma }_{{i}{i}}\exp \left\{-{\beta }_{{i}{i}}\,\mathrm{ln}\,{\left(\displaystyle \frac{s}{{\delta }_{{i}{i}}}\right)}^{2}\right\}\,{\rm{for}}\,{i}=1,\,2\,{\rm{and}}\,3\end{eqnarray} \tag{ 27 }$$

i.e. asymptotically converging to zero as size increases. Here, ${\gamma }_{ii}$ (in %), ${\beta }_{ii}$ (unitless) and ${\delta }_{ii}$ (in nm) are constants.

Zoom In Zoom Out Reset image size

This section describes the multiscale/multiphysics coupling results for the thermally induced phase transformation process. Multiscale coupling via the CWM methodology is performed on the MPF evolution. Here, the statistical information of length scale effects on the MPF evolution kinetics (obtained via MD and PhF simulations), is transferred/propagated from micro and meso scales to the macro scale. Furthermore, the CWM results are compared with continuum level experimental data and with macroscopic constitutive relations for NiTi SMA. Essentially, the continuum level experimental data/evolution law predicts an exponential dependency of MPF evolution on temperature, while the macroscopic constitutive relations are linear. Therefore, in this section, the experimental data is referred to as the 'exponential' model and MPF obtained from macroscopic constitutive relations is referred as the 'linear' model.

Figure 7 demonstrates the ability of the CWM method to transfer the trend and uncertainties of MPF to the continuum level. In this context, it is important to noted that, the experimental sample is the polycrystalline with mostly contain [110] B2 crystals, i.e. similar like the CWM simulations. Thus, it can be expected that, in the CWM method, the effect of the grain boundary energy and internal strain, on the MPF evolution is captured sufficiently, as is the case for the experiments. Figure 7(a) shows the evolution of MPF as a function of temperature predicted by the CWM method, experimental data, described above in section 3.1, and from the continuum level thermo-mechanical constitutive relation (mentioned in section 3.2). It can be observed that the evolution of MPF predicted by the CWM method matches very well the experimental data/exponential model i.e. the blue curve [8, 42]. Also, the CWM method transfers the statistics of length scale effects, as shown in figure 7(a), in terms of the uncertainty bars on the CWM curve. However, the continuum level linear material model shows significant deviation from the experimental and/or CWM results. To demonstrate the deviation, the differences between values from experimental data and from CWM and continuum material model are shown in figure 7(b). The difference between experimental data/exponential model and CWM predicted values is small or even negligible when uncertainties are considered, whereas this difference becomes very large for the case of experimental data and linear model obtained from continuum level thermo-mechanical constitutive relations.

Zoom In Zoom Out Reset image size

The interplay of cell/specimen size and grain size on the thermo-mechanical properties such as phase transformation process in SMAs has received little attention despite its importance [16], perhaps because of the difficulties in addressing it. In some recent studies [1, 21–23], attempts are made to model or experimentally demonstrate the effect of grain size on stress induced MPT (i.e. for superelasticity) in NiTi SMA. However, the interaction between sample size and grain size and their coupled effects on the MPT process was addressed in those studies. Therefore, differing from the previous studies [1, 21–23], in this study, some attempts are made to shed the light on the effect of such interplay. Here, 2D PhF simulations are used to address the interplay of simulation cell size and grain size on the thermally induced MPF evolution, while for the multiscale modeling, simulation results (of different size single crystal) from both MD and PhF are adopted, as discussed subsequently.

In order to perform PhF simulation of polycrystalline NiTi, simulation cells with the polycrystalline grain boundaries or interfaces are generated via Poisson–Voronoi tessellations, employing the MATLAB toolbox. To introduce the effect of grain boundary a gradient energy term is added, as mentioned in the Ahluwalia et al study [16]. This gradient energy term is expressed as the function of interfacial energy $(\gamma =0.9482\,{\rm{J}}\,{{\rm{m}}}^{-2}),$ width of the grain boundary $({\rm{\Delta }})$ and field variable $({\eta }_{i}).$ The value of the interfacial energy term $(\gamma )$ for all different type of interfaces (such as interface between all B2 and B19' crystals, as well as B2-B19' interface) is considered the same, as assumed in other studies [7, 16]. However, the effect of different B2 crystal orientation on the MPF evolution is incorporated in the PhF model through the elastic strain energy term [7, 8, 16]. This elastic strain energy depends on the transformation strain tensors for B2 to B19' variants $({\varepsilon }_{i}^{0}),$ which again depends on the crystal orientation (with respect to [100] B2 crystal) [7, 8, 16]. Therefore, before using those transformation strain tensors, they are rotated through the rotation matrix $R;$ i.e. ${\varepsilon }_{i}^{0}$ is converted into $R{\varepsilon }_{i}^{0}{R}^{T}$ [7, 8, 16]. Here, this rotation matrix depends on the crystal orientation of a particular B2 grain, with respect to the [100] B2 crystal. Thus, based on the aforementioned process, the effects of interface/grain boundary and crystal orientation are incorporated in the PhF model.

An indirect method (based on the multiscale approach) is used to mimic the results of polycrystalline NiTi SMA. For that, the wavelet coefficients of MPF from different size single crystal simulations (with the crystal orientation of [100], [110] and [111]) are used in a very unique way. Here, to construct the coefficients of the CWM of the polycrystalline NiTi, wavelet coefficients of different single crystal NiTi are fused from different length scales, i.e. from different size samples. In case of single crystals (with orientations [100], [110] and [111]), depending on the size of the cells, simulations are performed employing MD and PhF, as mentioned in section 2. It is important to mention that, before mixing the wavelet coefficients (obtained from different size single crystals of NiTi), they are multiplied by appropriate weight factors. These weight factors are adopted based on the fraction of the relevant crystal present in the simulation cell. Adopting these weight factors is important, since they act as a step up or step down operator and thus amplify or reduce the value of the wavelet coefficients of single crystal NiTi, such that similar (to a close approximation) wavelet coefficients of polycrystalline NiTi are obtained. Since, in this indirect method the CWM is constructed, it is referred to as the 'Hybrid CWM' and this term is used subsequently.

To demonstrate the coupled effect of grain size and length scale on MPF evolution and its associated multiscale modeling, a simulation cell of size ${s}=512$ nm is adopted and simulations are performed at four different temperatures (T = 300, 313, 325 and 350 K) with different grain sizes. Figure 8 shows the effect of grain size on the thermally induced MPF evolution, for simulation cell size ${s}=512$ nm. The average values of the grain size and corresponding MPF are directly plotted in figure 8, while their corresponding uncertainties are represented through the error bars. In order to directly compare the single crystal and polycrystal simulations, MPF values (mean and standard deviation) are provided for the four different temperatures. Figure 8(a) shows the coupled effects of grain size and simulation cell size, on the thermally induced MPF evolution, at T = 313 K. Also, figure 8(a) demonstrates the efficiency of the hybrid CWM method in mimicking the MPF evolution of polycrystalline NiTi SMA as a function of grain size, for T = 313 K. It is seen that the MPF obtained from the direct polycrystal 2D PhF simulations and the hybrid CWM method are in very good agreement with each other. Further, those results are compared with a reported computational study [16] and an experimental study [42]. Results from the computational study [16] agree very well with our direct polycrystalline NiTi simulation, as well as with the hybrid CWM predictions. However, the experimental value stays a little off from all other cases. Such deviation can be attributed to size effects for the NiTi samples/cells and due to the constant temperature simulation. The experiments [42] were conducted at the continuum level, with a fixed temperature rate (20 °C min⁻¹ cooling and heating rate).

Zoom In Zoom Out Reset image size

Next, to conform and emphasize the efficiency and prediction capability of the hybrid CWM method, results are generated at different temperatures (T = 300, 325 and 350 K) and compared to the direct simulation results of polycrystalline NiTi SMA. As shown in figures 8(b)–(d), the hybrid CWM results agree very well with the direct polycrystalline simulation results. Further, the deviation of the experimental MPF values from the simulation ones (from both the direct simulation and hybrid CWM) decreases with increasing temperature. For all examined temperatures, when the results from polycrystalline NiTi are compared to the results from single crystal NiTi, a significant difference in the MPF value is observed. Such noticeable difference in MPF value arises from the presence of grain boundary. Further, the simulations confirm that with the reduction in grain size, the evolution of the MPF decreases drastically and thus below a critical grain size MPT will not be observed in case of the polycrystalline NiTi SMA, as reported in many other studies [1, 10, 17–19, 21–23].

The above results show the coupled effect of grain size and simulation cell size on thermally induced MPF evolution, but do not clearly show the interplay between grain size and cell size and their effect on the MPF evolution. Figure 9 shows such interplay for three different temperatures (T = 300, 325 and 350 K), the interplay becoming prominent for smaller grain sizes. To clearly investigate the interplay, direct polycrystalline PhF simulations are performed. For each temperature (and mainly for T = 300 and 325 K), with the increase in simulation cell size (i.e. removal of length scale effects), the grain size effect becomes more pronounced, causing a drastic decrease in MPF evolution with decrease in grain size. This can be explained as follows. For larger simulation cell with smaller grain size, the effective amount (length) of grain boundary is high, which increases the total grain boundary energy. High grain boundary energy provides a driving force opposing the MPT process and thereby prevents the nucleation and growth of martensite (different variants of B19') from austenite (B2). Interestingly, it can be observed that with the increase in simulation cell size, as the length scale effect is removed, polycrystalline NiTi simulation results start agreeing very well with experimental data. To further demonstrate the effect of grain size in polycrystalline NiTi, results of single crystal NiTi simulations are also shown in figure 9. Finally, the relation between normalized MPF and normalized grain size can be express as

$$\begin{eqnarray}&&\displaystyle \frac{{\xi }_{p}(s,{s}_{g},T)}{\xi (s,T)}={\gamma }_{p}+{\beta }_{p}\,\mathrm{ln}\left(\displaystyle \frac{{s}_{g}}{s}\right),\end{eqnarray} \tag{ 28 }$$

where, ${\xi }_{p}(s,{s}_{g},T)$ and $\xi (s,T)$ are the MPFs corresponding to polycrystalline and single crystal NiTi SMA, respectively, ${s}_{g}$ and $s$ are the grain size and simulation cell size, respectively, ${\gamma }_{p}$ is a unitless constant which represent the normalized MPF value for a sample with very large grain size, and ${\beta }_{p}$ is a constant representing the normalized MPF and normalized grain size (in log scale) curve.

Zoom In Zoom Out Reset image size

This section describes the application of multiscale coupling to macroscopic continuum scale constitutive relations. In order to improve the continuum thermo-mechanical material model of NiTi SMA described previously, the temperature dependent MPF and its spatial uncertainties (obtained from multiscale coupling) are transferred into the continuum constitutive model. In the present version of the continuum model it is assumed that the MPF varies linearly with temperature.

Figure 10 shows the stress strain hysteretic loops and the associated energies obtained from different methods and from experimental data. The linear evolution law for the thermally induced MPF shows noticeable deviation from the experimental stress–strain data, while the exponential evolution law (expressed in equation (20)) and the CWM provided thermally induced MPF show very good agreement with the experimental data. The same holds for energy dissipation capacity of NiTi SMA (figure 10(b) for single cycle and figure 10(c) for multiple cycles). Thus, the improvement in the material modeling of NiTi SMA is rather significant, at least for those cases where NiTi SMA is used as an energy dissipative medium, such as in NiTi SMA based vibration dampers or actuators.

Zoom In Zoom Out Reset image size

The improvement in the continuum material model is demonstrated for other temperature cases in figure 11. It can be observed that within the transformation temperature zone, the stress–strain hysteresis loops obtained from the linear evolution law deviates from that obtained from CWM and from exponential evolution law. However, outside the transformation temperature regime, such differences in the stress–strain hysteresis loops diminish. Here at T = 200 K, the MPF obtained from different models is the same, thus in figure 11(a) stress–strain hysteresis loop obtained from different models are indistinguishable.

Zoom In Zoom Out Reset image size

The uncertainties associated with thermally induced MPF obtained from the multiscale coupling are used to demonstrate the hysteretic stress–strain response of NiTi SMA at different temperatures. In figure 12, for each temperature, five curves are plotted, i.e. at mean phase fraction value and at mean with two different levels of uncertainties $(\mu ,\,\mu \pm \sigma ,\,\mu \pm 2\sigma ).$ It can be observed that within the transformation temperature regime, uncertainty in the thermally induced MPF has a pronounced effect on the stress–strain response. However, this effect of uncertainties diminishes outside the transformation temperature zone.

Zoom In Zoom Out Reset image size