Thermodynamically consistent non-isothermal phase-field modelling of elastocaloric effect: indirect vs direct method

TL;DR

This paper develops a thermodynamically consistent non-isothermal phase-field model to simulate elastocaloric effects, comparing indirect and direct methods for calculating temperature changes, and demonstrates its effectiveness for designing elastocaloric cooling devices.

Contribution

It introduces a novel thermodynamic phase-field model for elastocaloric effect simulation, incorporating both indirect and direct computational methods, and validates their accuracy and efficiency.

Findings

Maximum adiabatic temperature change of 9.5 K in single crystal alloy.

Discrepancy between methods is within 10% at low stress.

Both methods are feasible for low-stress elastocaloric evaluation.

Abstract

Modelling elastocaloric effect (eCE) is crucial for the design of environmentally friendly and energy-efficient eCE based solid-state cooling devices. Here, a thermodynamically consistent non-isothermal phase-field model (PFM) coupling martensitic transformation with mechanics and heat transfer is developed and applied for simulating eCE. The model is derived from a thermodynamic framework which invokes the microforce theory and Coleman--Noll procedure. To avoid the numerical issue related to the non-differentiable energy barrier function across the transition point, the austenite-martensite transition energy barrier in PFM is constructed as a smooth function of temperature. Both the indirect method using isothermal PFM with Maxwell relations and the direct method using non-isothermal PFM are applied to calculate the elastocaloric properties. The former is capable of calculating both isothermal entropy change and adiabatic temperature change ($\Delta T_{\text{ad}}$), but induces high computation cost. The latter is computationally efficient, but only yields $\Delta T_{\text{ad}}$. In a model Mn-22Cu alloy, the maximum $\Delta T_{\text{ad}}$ ($\Delta T_{\text{ad}}^{\text{max}}$) under a compressive stress of 100 MPa is calculated as 9.5 and 8.5 K in single crystal (3.5 and 3.8 K in polycrystal) from the indirect and direct method, respectively. It is found that the discrepancy of $\Delta T_{\text{ad}}^{\text{max}}$ by indirect and direct method is within 10% at stress less than 150 MPa, confirming the feasibility of both methods in evaluating eCE at low stress. The results demonstrate the developed PFM herein, combined with both indirect and direct method for eCE calculations, as a practicable toolkit for the computational design of elastocaloric devices.

Tables

Table 1. Table 1: Material parameters and simulation parameters of Mn-22Cu for improper and proper MT PFM

parameter	name	value
$c_{11}$	elastic constant	76.588 GPa
$c_{12}$	elastic constant	14.588 GPa
$c_{44}$	elastic constant	31 GPa
$T_0$	chemical equilibrium temperature	245 K
$Q$	latent heat for improper MT	$4.84\times {10}^7$ J/m3
$L^{\eta }$	kinetic coefficient for improper MT	50 m3/s/J
$L^e$	kinetic coefficient for proper MT	1 m3/s/J
${\beta }^{\eta }$	gradient energy coefficient for improper MT	2.5$\times {10}^{-9}$ J/m
${\beta }^e$	gradient energy coefficient for proper MT	1$\times {10}^{-9}$ J/m
$K$	thermal conductivity	40 J/m/s/K
$\rho$	density	7500 kg/m3
$c_{\text{v}}$	specific heat per unit volume	2.64$\times$106 J/m3/K
$\alpha$	thermal expansion coefficient	10-5 K-1

Equations

$$\rho {\left(\frac{\partial S}{\partial \sigma }\right)}_T={\left(\frac{\partial \epsilon }{\partial T}\right)}_{\sigma }\text{or}\rho {\left(\frac{\partial S}{\partial \epsilon }\right)}_T=-{\left(\frac{\partial \sigma }{\partial T}\right)}_{\epsilon }.$$

$$\Delta S_{\text{iso}}(0\to \sigma )={\int }_0^{\sigma }\frac{1}{\rho }{\left(\frac{\partial \epsilon }{\partial T}\right)}_{\sigma }\text{d}\sigma ,$$

$$\Delta T_{\text{ad}}(0\to \sigma )=-{\int }_0^{\sigma }\frac{T}{\rho c}{\left(\frac{\partial \epsilon }{\partial T}\right)}_{\sigma }\text{d}\sigma$$

$${\sigma }_{ji,j}+b_i=0\text{in}\Omega ,$$

$$u_i={\hat{u}}_i\text{on}\partial {\Omega }_u,{\sigma }_{ji}{n}_j={\hat{t}}_i\text{on}\partial {\Omega }_{\sigma }$$

$${\sigma }_{ij}={\sigma }_{ji}.$$

$${\int }_{\partial \Omega }{\xi }^{\chi }{n}_i\text{d}A+{\int }_{\Omega }{\zeta }^{\chi }\text{d}V+{\int }_{\Omega }{\zeta }^{\chi \text{ex}}\text{d}V=0.$$

$${\xi }_{,i}^{\chi }+{\zeta }^{\chi }+{\zeta }^{\chi \text{ex}}=0.$$

$$\{\begin{array}{ll}{\xi }^{\chi }={\xi }_i^{{\eta }_I}& \\ {\zeta }^{\chi }={\zeta }^{{\eta }_I}& \text{for}\chi ={\eta }_I\\ {\zeta }^{\chi \text{ex}}={\zeta }^{{\eta }_I\text{ex}}& \end{array},\{\begin{array}{ll}{\xi }^{\chi }={\xi }_{ijk}& \\ {\zeta }^{\chi }={\zeta }_{jk}& \text{for}\chi ={\epsilon }_{kj}^0\\ {\zeta }^{\chi \text{ex}}={\zeta }_{jk}^{\text{ex}}& \end{array}.$$

$${\int }_{\Omega }\dot{e}\text{d}V={\int }_{\partial \Omega }({\sigma }_{ji}{\dot{u}}_j+{\xi }^{\chi }\dot{\chi }-j_i^{\text{h}})n_i\text{d}A+{\int }_{\Omega }(b_i{\dot{u}}_i+{\zeta }^{\chi \text{ex}}\dot{\chi }+q^{\text{h}})\text{d}V,$$

$$q^{\text{h}}=\{\begin{array}{ll}Q{\sum }^I\dot{{\eta }_I}& \text{for}\chi ={\eta }_I\\ Q^e({\dot{\epsilon }}_{11}^0-{\dot{\epsilon }}_{22}^0)& \text{for}\chi ={\epsilon }_{kj}^0\end{array},$$

$${\int }_{\Omega }\dot{e}\text{d}V={\int }_{\Omega }({\sigma }_{ji}{\dot{\epsilon }}_{ij}+{\sigma }_{ji,i}{\dot{u}}_j+{\xi }_{,i}^{\chi }\dot{\chi }+{\xi }^{\chi }{\dot{\chi }}_{,i}-j_{i,i}^{\text{h}}+b_i{\dot{u}}_i+{\zeta }^{\chi \text{ex}}\dot{\chi }+q^{\text{h}})\text{d}V,$$

$$\dot{e}={\sigma }_{ij}{\dot{\epsilon }}_{ij}+{\xi }^{\chi }{\dot{\chi }}_{,i}-{\zeta }^{\chi }\dot{\chi }-j_{i,i}^{\text{h}}+q^{\text{h}}.$$

$${\int }_{\Omega }\dot{s}\text{d}V+{\int }_{\partial \Omega }\frac{j_i^{\text{h}}}{T}{n}_i\text{d}A-{\int }_{\Omega }\frac{q^{\text{h}}}{T}\text{d}V\ge 0$$

$$\dot{s}+\Big{(}\frac{j_{i}^{\text{h}}}{T}\Big{)}_{,i}-\frac{q^{\text{h}}}{T}\geq 0.$$

$$f=e-Ts$$

$$\frac{1}{T}\frac{\partial f}{\partial t}=\frac{1}{T}\frac{\partial e}{\partial t}-\frac{s}{T}\frac{\partial T}{\partial t}-\frac{\partial s}{\partial t}.$$

$$\frac{1}{T}\frac{\partial f}{\partial t}\leq\frac{1}{T}\sigma_{ij}\dot{\varepsilon}_{ij}+\frac{1}{T}\Big{(}\xi^{\chi}\dot{\chi}_{,i}-\zeta^{\chi}\dot{\chi}\Big{)}-\frac{s}{T}\frac{\partial T}{\partial t}+j_{i}^{\text{h}}\Big{(}\frac{1}{T}\Big{)}_{,i}$$

$$f=f({\epsilon }_{ij},\chi ,\dot{\chi },{\chi }_{,i},T,T_{,i})$$

$$\frac{\partial f}{\partial t}=\frac{\partial f}{\partial {\epsilon }_{ij}}\frac{\partial {\epsilon }_{ij}}{\partial t}+\frac{\partial f}{\partial \chi }\frac{\partial \chi }{\partial t}+\frac{\partial f}{\partial \dot{\chi }}\frac{\partial \dot{\chi }}{\partial t}+\frac{\partial f}{\partial {\chi }_{,i}}\frac{\partial {\chi }_{,i}}{\partial t}+\frac{\partial f}{\partial T}\frac{\partial T}{\partial t}+\frac{\partial f}{\partial T_{,i}}\frac{\partial T_{,i}}{\partial t}.$$

$$\frac{1}{T}\Big{(}\frac{\partial f}{\partial\varepsilon_{ij}}-\sigma_{ij}\Big{)}\dot{\varepsilon}_{ij}+\frac{1}{T}\Big{(}\frac{\partial f}{\partial\chi}+\zeta^{\chi}\Big{)}\dot{\chi}+\frac{1}{T}\frac{\partial f}{\partial\dot{\chi}}\frac{\partial\dot{\chi}}{\partial t}+\frac{1}{T}\Big{(}\frac{\partial f}{\partial\chi_{,i}}-\xi^{\chi}\Big{)}\dot{\chi}_{,i}+\frac{1}{T}\Big{(}\frac{\partial f}{\partial T}+s\Big{)}\dot{T}+\frac{1}{T}\frac{\partial f}{\partial T_{,i}}\dot{T}_{,i}-j_{i}^{\text{h}}\Big{(}\frac{1}{T}\Big{)}_{,i}\leq 0$$

$${\sigma }_{ij}=\frac{\partial f}{\partial {\epsilon }_{ij}},{\xi }^{\chi }=\frac{\partial f}{\partial {\chi }_{,i}},\frac{\partial f}{\partial \dot{\chi }}=0,s=-\frac{\partial f}{\partial T},\frac{\partial f}{\partial T_{,i}}=0.$$

$$\frac{1}{T}\Big{(}\frac{\partial f}{\partial\chi}+\zeta^{\chi}\Big{)}\dot{\chi}-j_{i}^{\text{h}}\Big{(}\frac{1}{T}\Big{)}_{,i}\leq 0.$$

$$\dot{\chi }=-L\left(\frac{\partial f}{\partial \chi }+{\zeta }^{\chi }\right)$$

$$j_i^{\text{h}}={\kappa }_{ij}{\left(\frac{1}{T}\right)}_{,j}$$

$$\dot{\chi }=-L\left[\frac{\partial f}{\partial \chi }-{\left(\frac{\partial f}{\partial {\chi }_{,i}}\right)}_{,i}-{\zeta }^{\chi \text{ex}}\right].$$

$$j_i^{\text{h}}=-K_{ij}{T}_{,j}$$

$$f({\epsilon }_{ij},\chi ,{\chi }_{,i},T)=f^{\text{chem1}}+f^{\text{chem2}}+f^{\text{grad}}+f^{\text{ela}}.$$

$$f^{\text{chem1}}=A(T)\sum^{I}\eta_{I}^{2}-B(T)\sum^{I}\eta_{I}^{3}+C(T)\Big{(}\sum^{I}\eta_{I}^{2}\Big{)}^{2}$$

$$f^{\text{grad}}=\frac{1}{2}{\beta }^{\eta }(T)\sum ^I{\eta }_{I,i}{\eta }_{I,i}$$

Full text

††journal: ******

1 Introduction

The solid-state cooling technology features itself as an environmentally friendly and efficient alternative to the traditional vapor compression refrigeration, thus helping move towards carbon neutrality. Typical solid-state cooling by magnetocaloric [1, 2, 3], electrocaloric [4, 5], barocaloric  [6, 7] and elastocaloric [8, 9, 10, 11] effects has been widely investigated, among which elastocaloric cooling stands out with the large temperature change [9, 12], large working temperature window [10], and excellent coefficient of performance [13, 14]. The elastocaloric cooling is realized by means of elastocaloric effect (eCE), which originates from the latent heat associated with martensitic transformation (MT) in shape memory alloys (SMAs). Upon loading, the exothermic austenite-martensite transformation (also called the conventional MT) would cause a temperature increase in the adiabatic process. Upon unloading, the endothermic martensite-austenite transformation (also called the inverse MT) occurs, and a rapid drop in temperature arises [8]. The eCE can be quantified by the adiabatic temperature change ($\Delta T_{\text{ad}}$) or the isothermal entropy change ($\Delta S_{\text{iso}}$).

Modelling eCE plays an important role in the computational design of elastocaloric devices and can be an essential complement to experiments. In general, phenomenological constitutive models and phase-field model [15, 16] (PFM) are utilized to calculate the eCE in SMAs. The phenomenological Tanaka-type model [17, 18, 19, 20, 11, 21] has advantages in the description of temperature change and stress-strain behavior. Very good agreement between the measured and calculated values of $\Delta T_{\text{ad}}$ is obtained. Analogously, the phase transformation kinetic model [22, 23, 24] and crystal plasticity-based constitutive model [25, 26] are also proposed to predict the eCE. These phenomenological constitutive models can be used to obtain the macroscopic elastocaloric properties, but have difficulties in simulating the spatial and temporal evolution of microstructure details during the MT and thus are hardly applicable to the optimization of eCE by microstructure engineering.

eCE in SMAs is intrinsically ascribed to MT, which is a first-order phase transformation and can be induced by external stress/strain field or by temperature field. The PFM developed from Landau’s theory of phase transformation is widely used to simulate the stress- and temperature-induced MT [27]. Wang et al. [28] built a three-dimensional PFM of MT, taking into account the transformation-induced strain, which comprehensively described a generic cubic-to-tetragonal MT. In the PFM, each of the martensitic variants is described by an order parameter $\eta_{I}$ ($I=1,2,...,n$ with $n$ as the total number of various crystallographically equivalent martensitic variants), and the evolution of each of variants is governed by the Ginzburg-Landau equations. Several studies also considered hexagonal-to-orthorhombic [29], cubic-to-tetragonal [30], and tetragonal-to-monoclinic [31] transitions. Over the last decades, the PFM for MT has been constantly improved to include/study more physical phenomena in SMAs. For example, the martensitic reorientation, the temperature-induced transformation, and the stress-induced MT have been of high interests. Levitas et al. [32, 33, 34] constructed a thermodynamically consistent PFM for transformations between austenite and martensitic variants and martensitic reorientation. Since then the phase-field (PF) theory is extended for the cases of surface stresses [35, 36], large strain [37, 38], and eCE [39, 40, 41, 42]. The research that explores eCE by PFM has recently emerged, owing to the rise of efficient solid-state cooling technology by eCE as well as the PFM’s advantage in simulating the evolution of microstructure and temperature during the MT process.

In 2015, Levitas et al. [43] proposed a multiphase phase-field theory for temperature- and stress-induced MT, which allows for a presence of the third phase at the interface between the two other phases. Then, Cui et al. [44] employed a non-isothermal PFM to study the stress- and temperature-induced MT, as well as the corresponding latent heat in SMAs. However, the discontinuous piecewise energy barrier between austenite and martensite is adverse to the calculation of $\Delta S_{\text{iso}}$ or $\Delta T_{\text{ad}}$, which possibly makes the Maxwell relations based eCE analysis problematic around the transition point. In real materials, the austenite-martensite transition energy barrier varies continuously and increases with temperature in the whole temperature range. Furthermore, Sun et al. [45, 46] modified this energy barrier as a continuous piecewise function of temperature in PFM and discussed the effect of grain size, crystal orientation, and loading rate on $\Delta T_{\text{ad}}$. Cissé et al. [40] introduced an exponential function to express the energy barrier and computed qualitatively the eCE in CuAlBe, including $\Delta T_{\text{ad}}$ calculated by the direct method, coefficient of performance and cyclic deformation. Their energy barrier equations are continuous but not differentiable at the transition point (see Fig. 1 for more details). In addition, Xu et al. [47, 48] constructed a similar PFM by introducing an extra grain boundary energy to investigate the grain size dependent super-elasticity and eCE using the direct method in nanocrystalline NiTi SMAs. A mesoscale model [49] based on the Müller-Achenbach-Seelecke theory without PF order parameters was also developed to describe the evolution of local temperature and strain in SMA films during elastocaloric cycling. Meanwhile, the latent heat effect [44, 46], size effect [45, 50, 47], plasticity [40, 41] , and microstructure design [48, 42] are considered to regulate eCE. Nevertheless, all the above mentioned works are focused on the application of PF simulations to directly evaluate $\Delta T_{\text{ad}}$ during loading and unloading by solving the heat-transfer equation.

In general, the indirect and direct methods are two basic approaches to calculate the key parameter $\Delta T_{\text{ad}}$ [51]. Directly measuring the (adiabatic) temperature change distribution by thermography without thermodynamic calculations is called the direct method. The difficulty of ensuring adiabatic conditions usually leads to the inevitable heat loss and thus the minor $\Delta T_{\text{ad}}$ calculated by the direct method. Hence, researchers also adopt the indirect method [52, 53, 54] to calculate $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ by means of the Maxwell relations (Eq. 1) or Clausius-Clapeyron equation. $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ can be calculated from the average stress-strain responses at each temperature or average strain-temperature responses at each stress. The Maxwell relations for eCE evaluation read as

[TABLE]

Using Eq. 1, $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ between the initial and final stress are expressed as [8]

[TABLE]

[TABLE]

where $S$ is entropy, $T$ temperature, $\varepsilon$ strain, $\sigma$ stress, $\rho$ material’s density, and $c$ specific heat.

The indirect method is initially extensively utilized to calculate the magnetocaloric effect [2, 3]. Bonnot et al. [8] first introduced the Maxwell relations (Eq. 1) to evaluate $\Delta S_{\text{iso}}$ associated with the MT in Cu-Zn-Al single crystal and found that the indirectly calculated $\Delta S_{\text{iso}}$ agrees well with the directly measured value. This work has brought about the widespread use of the indirect method for evaluating eCE experimentally. Chen et al. [53] found the directly measured $\Delta T_{\text{ad}}$ in nanocrystalline Ti-Ni-Cu SMA is consistent with $\Delta T_{\text{ad}}$ from the indirect method. However, Pataky et al. [55] demonstrated that the directly measured $\Delta T_{\text{ad}}$ is much lower than the temperature change calculated by the Maxwell or Clausius-Clapeyron relations. Similarly, Qian et al. [56] measured the eCE in CuAlZn and CuAlMn SMAs under compression, and showed that $\Delta T_{\text{ad}}$ estimated by the indirect method is nearly three times as large as the directly measured one. Experimental results indicate the discrepancy in $\Delta T_{\text{ad}}$ from the indirect and direct method, resulting in some misunderstandings and disputes about the eCE in the same material, similar observations have been reported for electrocaloric effect of ferroelectric materials  [57]. It implies that if the stress-strain curves across the transition point is not measured with a sufficiently small incremental step of temperature, the indirect method may lead to unrealistic $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ owing to the possibly incorrect estimation of $\partial\varepsilon/\partial T$ [20, 51]. In short, the literature review indicates the calculation method and its effectiveness as a key issue for evaluating eCE by PF simulations. The indirect method has the merit of calculating both $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$, whereas the direct method only attains $\Delta T_{\text{ad}}$. However, the indirect method based on the Maxwell relations, which is widely adopted in experimental studies of eCE, is relatively unexplored in terms of PF simulations. Whether the eCE evaluated by the indirect and direct method using PFM is consistent with each other remains yet unknown.

In this work, we aim to develop a thermodynamically consistent non-isothermal PFM coupling MT with mechanics and heat transfer to evaluate eCE using indirect and direct method. The model is derived from a thermodynamic framework which invokes the microforce theory and the Coleman–Noll procedure. In order to account for the possibly non-sharp transition in real materials and avoid the problematic calculation due to the non-differentiable energy barrier function across the transformation temperature, the austenite-martensite transition energy barrier in PFM is introduced as a smooth function of temperature. The PFM details and their derivation accomplished by a systematic use of thermodynamic principles are presented in Sect. 2 for both improper and proper MT. Sect. 3 gives an application of non-isothermal PFM to Mn-22Cu alloy. In Sect. 4, the PFM for both single crystal and polycrystal is numerically implemented by the finite element method. In Sect. 5, benchmark simulations for reproducing the ferroelastic behavior and elastocaloric cycle are performed to verify the PFM. In addition, the elastocaloric properties in Mn-22Cu SMAs are obtained through the indirect and direct method, and the comparison of $\Delta T_{\text{ad}}$ calculated from these two methods are discussed. Sect. 6 gives the conclusive summary.

2 Thermodynamically consistent phase-field model

The phase-field method emerges as a powerful tool for modelling microstructure evolution and predicting eCE in elastocaloric materials. It describes a microstructure using a set of conserved or non-conserved order parameters (OPs) that are continuous across the smooth interfacial regions [58], which could implicitly track the positions of interfaces. For deriving the partial differential equations for the non-isothermal PFM, we consider a closed system of volume $\Omega$ in which a pure material undergoes a first-order phase transformation between austenite and martensite. Noticeably, the MT could be divided into improper and proper MT. The improper MT is characterized by the displacements of atoms within a unit cell and is described by soft optical displacement modes (e.g., ceramic materials), where strain is generated as a secondary effect. The proper MT is a homogeneous stress-free (eigen) strain that characterizes a change in crystal lattice parameters [59, 60], e.g., fcc to bcc MT in Fe alloys. Therefore, the proper and improper MT in SMAs should be distinguished and modelled by PFM using different OPs. In 1995, Kartha et al [61] proposed a phase field theory for the proper MT by using the strain tensor, which is based on a straightforward approach that uses the Helmholtz free energy defined by the OPs of eigen strain $\varepsilon_{ij}^{0}$. In addition, the non-conserved OPs $\eta_{I}$ are also chosen to distinguish various phases in SMAs regardless of the improper or proper MT [44].

In this work, non-conserved OPs $\chi$ is chosen to describe various phases during MT in SMAs. For the improper MT, $\chi$ is $\eta_{I}$ which is chosen to distinguish various phases: $\eta_{I}=0$ represents austenite and $\eta_{I}=1$ represents the $I^{\text{th}}$ martensitic variant. For the proper MT, the OPs $\chi$ is the eigen strain $\varepsilon_{ij}^{0}$, which can directly represent the martensitic variants and austenite. The following subsections present the complete thermodynamic consistent derivation that involves balance law, Coleman–Noll analysis, constitutive relations and evolution equations, and Helmholtz free energy.

2.1 Balance law

Balance of linear momentum

For the body $\Omega$ with a boundary $\partial\Omega$, the quasi-static mechanical equilibrium equation is described by

[TABLE]

[TABLE]

where $\sigma_{ji}$ is the Cauchy stress tensor and $b_{i}$ is the body force. Here, we assume $b_{i}=0$. The Latin indices ($i$, $j$) run over the range of 1–3. $\hat{u}_{i}$ is the displacement prescribed on the boundary $\partial\Omega_{u}$, $n_{j}$ is the outward surface unit vector, and $\hat{t}_{i}$ is the surface traction on the boundary $\partial\Omega_{\sigma}$ [62].

• 2.

Balance of angular momentum

The stress tensor is symmetric, i.e.,

[TABLE]

• 3.

Balance of microforce associated with OPs $\chi$

Based on the Gurtin’s microforce theory [63], we assume that there exist a set of forces that accounts for the phase transition dynamics. These forces are called microforces because they are involved with the local transformation of the material, rather than the macroscopic movements. In this work, the microforce is associated with the non-conserved OPs $\chi$, describing the phase transformation between austenite and martensite. The relevant microforce definitions are as follows.

$\xi^{\chi}$: the microstress,

$\xi^{\chi}n_{i}$: the surface microforce with $n_{i}$ as the unit outer normal to $\partial\Omega$,

$\zeta^{\chi}$: the internal microforce,

$\zeta^{\chi\text{ex}}$: the external microforce.

The balance of microforce associated with the phase transformation is written as

[TABLE]

By means of the Gauss law, we can obtain the equivalent local microforce balance as

[TABLE]

Specifically, for the improper MT ($\chi=\eta_{I}$) and proper MT ($\chi=\varepsilon_{kj}^{0}$), the microforces can be written as [64]

[TABLE]

• 4.

Balance of energy

The principle of the balance of energy can be stated as: the rate of change specific internal energy density $e$ is equal to the sum of power due to external forces and the heat input to the system. In solid with negligible inertia, the thermal heat flux $j_{i}^{\text{h}}$ and the heat source per unit volume $q^{\text{h}}$ are considered, and the internal microforce does not contribute to the energy change.

The first law of thermodynamics can be presented in the form of energy balance equation as

[TABLE]

where $u_{j}$ is the displacement, and $n_{i}$ is the unit outer normal to $\partial\Omega$. During loading or unloading, the elastocaloric material will release or absorb the latent heat associated with the stress-induced MT. Herein, $q^{\text{h}}$ is the release/absorb rate of internal heat source from MT, which is only assumed to be related to the evolution rate ($\dot{\chi}$) of phase transformation [41, 65]. For example, the $q^{\text{h}}$ for both $\eta_{I}$ and $\varepsilon_{kj}^{0}$ in cubic to tetragonal MT can be written as

[TABLE]

where $Q$ and $Q^{e}$ are the latent heat and play a critical role in the temperature change during MT. By Gauss law which converts the surface integration into volume one, Eq. 10 can be rewritten as

[TABLE]

where $\dot{\varepsilon}_{ij}$ is the total strain rate with $\dot{\varepsilon}_{ij}=\frac{1}{2}(\dot{u}_{i,j}+\dot{u}_{j,i})$. Then, combining Eqs. 4, 6 and 8, as well as considering that the equation holds for any arbitrary volume, the energy balance equation in the local form reads

[TABLE]

2.2 Constitutive relations and evolution equations

Second law of thermodynamics

For the non-isothermal system, the second law of thermodynamics or entropy inequality, combining the global entropy balance with the Clausius–Duhem inequality [66] for the volume $\Omega$, is expressed as

[TABLE]

where $s$ is the specific entropy per unit volume and $T$ is the temperature. Converting the surface integration in Eq. 14 into volume integration and considering its validity in any volume, we can obtain

[TABLE]

• 2.

Free energy imbalance

Herein, the Helmholtz free energy is chosen as a proper thermodynamic potential in the non-isothermal system. The Helmholtz free energy density $f$ per unit volume is defined by

[TABLE]

Taking time derivatives at both sides, we get the relation

[TABLE]

Substituting the internal energy balance equation Eq. 13 and the second law of thermodynamics Eq. 15 into the above relation Eq. 17, we can get an inequality

[TABLE]

the inequality in Eq. 18 is referred as the free energy imbalance. It plays an analogous role to Eq. 15 in placing restrictions on the constitutive relations.

• 3.

Coleman–Noll type analysis

In order to close the model, the constitutive relations for the Cauchy stress, the internal energy density, the entropy density, the heat flux, and the microforces should be provided. In this part, we derive the explicit form of the constitutive relations and the kinetic equations in terms of a thermodynamic potential. In this derivation, the Coleman–Noll argument [67] is applied so that the resulting constitutive relations will be thermodynamically consistent. Invoking Truesdell’s principle of equipresence [68], it is reasonable to assume that $f$, $s$, $e$, $\sigma_{ij}$, $j_{i}^{\text{h}}$, $\xi^{\chi}$ and $\zeta^{\chi}$ depend on $\varepsilon_{ij}$, $\chi$, $\dot{\chi}$, $\chi_{,i}$, $T$, $T_{,i}$. Specifically, the Helmholtz free energy density $f$ can be written as

[TABLE]

The time derivative of $f$ and the chain rule lead to

[TABLE]

Now substituting Eq. 20 into the free energy imbalance Eq. 18, and grouping terms together, we can get

[TABLE]

Here we provide an analysis of Eq. 21 by invoking the arguments made by Coleman and Noll [67], and notice that Eq. 21 is linear with respect to $\dot{\varepsilon}_{ij}$, $\partial\dot{\chi}/\partial t$, $\dot{\chi}_{,i}$, $\dot{T}$, $\dot{T}_{,i}$. Hence, to satisfy Eq. 21 in any admissible thermodynamics process, the coefficients of linear terms must vanish and thus these constitutive relations can be obtained

[TABLE]

With the relations in Eq. 22, $f$ is independent of $\dot{\chi}$ and $T_{,i}$, and so Eq. 19 can be rewritten as $f=f(\varepsilon_{ij},\chi,\chi_{,i},T)$. Then, the left nonlinear terms in Eq. 21 are reduced to

[TABLE]

• 4.

Evolution equations

The inequality in Eq. 23 can be satisfied by assuming the following relations

[TABLE]

[TABLE]

where positive $L$ and positive semi-definite $\kappa_{ij}$ are coefficients with respect to temperature. Further, combining Eqs. 8, 22 with 24, the evolution equations for non-conserved OPs $\chi$ can be written as

[TABLE]

The Eq. 26 coincides with the general Allen–Cahn equation [69]. In addition, the heat flux equation reads as

[TABLE]

where $K_{ij}=\kappa_{ij}/T^{2}$ is a tensor representing the thermal conductivity.

2.3 Helmholtz free energy

In the non-isothermal PFM, Helmholtz free energy density $f$ of the system with microstructure evolution consists of the chemical free energy density $f^{\text{chem1}}+f^{\text{chem2}}$, the gradient free energy density $f^{\text{grad}}$, and the elastic energy density $f^{\text{ela}}$, i.e.,

[TABLE]

2.3.1 Improper MT

The chemical free energy or Landau free energy $f^{\text{chem1}}$ of a closed system containing austenite and martensite is determined by the distribution of OPs $\eta_{I}$. The minimum energy states at $\eta_{I}=0$ and $\eta=1$ represent the austenite phase and the corresponding martensite variant $I$, respectively. The associated chemical energy can be constructed as a Landau polynomial [70], i.e.,

[TABLE]

where $A(T)$, $B(T)$, and $C(T)$ are temperature-dependent coefficients.

The gradient energy is often constructed to account for the interface between different phases, which can be expressed in terms of OPs’ gradient, i.e.,

[TABLE]

where $\beta^{\eta}(T)$ is the gradient energy coefficient depending on temperature.

The total strain tensor $\varepsilon_{ij}$ is given as

[TABLE]

where $\varepsilon_{ij}^{\text{ela}}$ is the elastic strain, $\varepsilon_{ij}^{\text{tr}}$ is the transformation strain generated by structural transformation, and $\varepsilon_{ij}^{\text{th}}$ is the thermal strain caused by temperature change. Therefore, the elastic energy can be expressed as

[TABLE]

in which C_{ijkl}=\sum^{I}\Phi(\eta_{I})C_{ijkl}^{\eta_{I}}+\big{(}1-\sum^{I}\Phi(\eta_{I})\big{)}C_{ijkl}^{A} with the interpolation function $\Phi(\eta_{I})=\eta_{I}^{3}(10-15\eta_{I}+6\eta_{I}^{2})$. $C_{ijkl}^{\eta_{I}}$ and $C_{ijkl}^{A}$ are the elastic tensors of the $I^{\text{th}}$ martensitic variant and austenite, respectively.

The chemical energy $f^{\text{chem2}}$ represents the main concave term of the free energy and is related to the heat conduction (e.g., the internal energy density from temperature change). The complete and general form of the $f^{\text{chem2}}$ can be formulated as

[TABLE]

where $c_{i}=\sum^{I}\Phi(\eta_{I})c_{i}^{\eta_{I}}+(1-\sum^{I}\Phi(\eta_{I}))c_{i}^{A}$ with $c_{i}^{\eta_{I}}$ and $c_{i}^{A}$ as the $i^{\text{th}}$ coefficients for the specific heat (Eq. 34) of $I^{\text{th}}$ martensitic variant and austenite, respectively. Meanwhile, the specific heat per unit volume is expressed as a polynomial function of temperature, i.e.,

[TABLE]

In addition, in order to derive the general kinetic equation for temperature $T$, combining Eqs. 8, 13, 22 and 27, the evolution equation of $e$ reads as

[TABLE]

Neglecting the first three terms on the right side of Eq. 35, one yields the common governing equation for temperature as

[TABLE]

2.3.2 Proper MT

As an example here, the OPs of eigen strain $\varepsilon_{ij}^{0}$ is considered to model the proper MT. The corresponding chemical free energy density $f^{\text{chem1}}$ describing the proper MT can be represented as [71]

[TABLE]

where $Q_{1}$ and $Q_{5}$ are bulk and shear modulus, respectively. The coefficients $Q_{2}$, $Q_{3}$, and $Q_{4}$ are Landau constants determining the transition temperature $T_{0}$ and the transformation strains in the product phase. $e_{i}$ are the symmetry-adapted strain defined in term of the transformation strains as

[TABLE]

In addition, the $f^{\text{chem2}}$ is similar to the case of the improper MT in Eq. 33.

The austenite and martensite can be described by the eigen strain $\varepsilon_{ij}^{0}$. For instance, $\varepsilon_{23}^{0}=\varepsilon_{13}^{0}=\varepsilon_{12}^{0}=0$ is set for a cubic to tetragonal MT. Further, the austenite (cubic) is represented as ($\varepsilon_{11}^{0}=0$, $\varepsilon_{22}^{0}=0$, $\varepsilon_{33}^{0}=0$), the three tetragonal martensitic variants are described by tet1 = $(-\varepsilon_{0},~{}1/2\varepsilon_{0},~{}1/2\varepsilon_{0})$, tet2 = $(1/2\varepsilon_{0},~{}-\varepsilon_{0},~{}1/2\varepsilon_{0})$, and tet3 = $(1/2\varepsilon_{0},~{}1/2\varepsilon_{0},~{}-\varepsilon_{0})$, where $\varepsilon_{0}$ is the magnitude of the spontaneous strain at a given temperature.

Similarly, the gradient energy density can be written as

[TABLE]

where $\beta^{e}(T)$ is the strain gradient coefficient. The elastic strain energy density can be written as

[TABLE]

3 Application of non-isothermal PFM to Mn-22Cu alloy

In this section, we apply the above non-isothermal phase-field framework to a model Mn-22Cu alloy and present the detailed formulations for the total free energy, constitutive relations, and governing or evolution equations. In Mn-22Cu SMA, there exist a face-centered cubic high-symmetry austenitic phase at high temperature and three variants of a face-center tetragonal low-symmetry martensitic phase at low temperature [72, 73].

3.1 Improper MT

The three martensitic variants are energetically equivalent. Herein we choose non-conserved OPs $\eta_{I}$ ($I=1,2,3$) to represent the improper MT in PFM, and the value of $\eta_{I}$ varies from 0 to 1.

3.1.1 Chemical free energy density

The chemical free energy $f^{\text{chem1}}$ represents the chemical driving force of the MT in a stress-free Mn-22Cu alloy, which can be expressed as a Landau 2-3-4 polynomial

[TABLE]

where $A(T)$, $B(T)$ and $C(T)$ are positive temperature-dependent coefficients, expressed as $A(T)=16\Delta G^{*}$, $B(T)=A(T)-4\Delta G_{\text{m}}$ and $C(T)=0.5A(T)-3\Delta G_{\text{m}}$. $\Delta G^{*}$ is the temperature-dependent energy barrier between austenite and martensite. $\Delta G_{\text{m}}$ is the driving force of MT. The temperature-dependent profile of $f^{\text{chem1}}$ in non-isothermal system and the different models for austenite-martensite energy barriers are shown in Fig. 1a and b, respectively. $\eta_{I}=0$ or 1 represents the system stable or metastable from the principle of minimization of the free energy. More explicitly, OPs $\eta_{I}=1$ indicates the $I^{\text{th}}$ martensitic variant and $\eta_{I}=0$ ($I=1,2,3$) represents austenite.

Since the MT is a first-order diffusionless structural transformation and is not sharp in real materials [74], there are several ways in the literature [40, 44, 46] to formulate the energy barrier function $\Delta G^{*}$. For instance, Cui et al. [44] proposed that

[TABLE]

Sun et al. [46] proposed that

[TABLE]

and Cissé et al. [40] proposed that

[TABLE]

where $Q$ is the specific latent heat and $T_{0}$ is the chemical equilibrium temperature. It can be seen from Fig. 1b that the first derivative of all these $\Delta G^{*}$ functions [40, 44, 46] are discontinuous at $T_{0}$, leading to the problematic calculation of driving force and eCE at $T_{0}$. Herein, we formulate the energy barrier as a smooth function of temperature. Specifically, we take advantage of the hyperbolic tangent function to modify the piecewise $\Delta G^{*}$ function, i.e.,

[TABLE]

$\Delta G_{\text{m}}$ is also a continuous function of temperature, i.e.,

[TABLE]

It should be noted that $\delta T$ is a new parameter associating with the energy barrier, which could be adjusted according to the transformation temperature window from experimental results. In this work, we assume a moderately sharp transition and set $\delta T=2$ K. As shown in Fig. 1b, compared to the functions proposed in literature [40, 44, 46], our modified energy barrier function is continuous and differentiate at $T_{0}$. This modification could resolve the difficulty of calculating eCE by indirect and direct method at the transition point.

In addition, the internal energy density from temperature change $f^{\text{chem2}}$ reads

[TABLE]

where $c_{\text{v}}=\sum_{I=1}^{3}\Phi(\eta_{I})c_{\text{v}}^{\eta_{I}}+(1-\sum_{I=1}^{3}\Phi(\eta_{I}))c_{\text{v}}^{A}$ is the specific heat per unit volume. $c_{\text{v}}^{\eta_{I}}$ and $c_{\text{v}}^{A}$ are the specific heat per unit volume of the $I^{\text{th}}$ martensitic variant and austenite, respectively. The difference between $c_{\text{v}}^{\eta_{I}}$ and $c_{\text{v}}^{A}$ is insignificant for Mn-22Cu alloy so that they approximately take the same values.

3.1.2 Gradient energy density

The gradient energy density can be expressed as a function of the OPs’ gradient [60], i.e.,

[TABLE]

where $\beta^{\eta}(T)$ is the gradient energy coefficient related to the interfacial energy and interface thickness [75].

3.1.3 Elastic strain energy density

The elastic strain energy density in Mn-22Cu alloys can be given as

[TABLE]

where $C_{ijkl}$ is the component of the fourth-order elastic tensor. In this work, $c_{ij}$ in Table LABEL:t1 denotes the Voigt notation of $C_{ijkl}$, i.e., $c_{11}=C_{1111},c_{12}=C_{1122},c_{44}=C_{1212}$.

We incorporate $\varepsilon_{ij}^{\text{tr}}$ and $\varepsilon_{ij}^{\text{th}}$ into the PFM. $\varepsilon^{\text{tr}}_{ij}$ can be expressed as the stress-free eigen strain, which is in general defined as

[TABLE]

$\varepsilon^{00}(i)$ ($i=1,2,3$) is determined by the orientation relationship and lattice distortion between martensite and austenite with regard to the FCC-FCT MT. It can be simplified to [60]

[TABLE]

where $\varepsilon_{1}=(a-a_{c})/a_{c}$ and $\varepsilon_{3}=(c-a_{c})/a_{c}$. We choose $\varepsilon_{1}=0.01$ and $\varepsilon_{3}=-0.02$ [44] in this work.

In addition, the thermal strain is computed as

[TABLE]

where \alpha_{ij}=\sum_{I=1}^{3}\Phi(\eta_{I})\alpha^{\eta_{I}}\delta_{ij}+\big{(}1-\sum_{I=1}^{3}\Phi(\eta_{I})\big{)}\alpha^{A}\delta_{ij} is a tensor representing thermal expansion, in which $\alpha^{\eta_{I}}$ and $\alpha^{A}$ are thermal expansion coefficients for the $I^{\text{th}}$ martensitic variant and austenite, respectively. $T^{\text{ref}}$ is the reference temperature at which there is zero thermal strain. In the PF simulations, $T^{\text{ref}}$ is the initial temperature.

3.1.4 Constitutive relations for PFM of Mn-22Cu alloy

Using Eq. 22, the constitutive relations for PFM of Mn-22Cu alloy is given by

[TABLE]

The internal energy density $e$ for Mn-22Cu alloy can be simplified as

[TABLE]

3.1.5 Governing equations for PFM of Mn-22Cu alloy

We adopt Eq. 26 to govern the spatial and temporal evolution of $\eta_{I}$, which is analogous to the time-dependent Ginzburg-Landau (TDGL) kinetic equation [76, 77]. Further, using Eqs. 35 and 54, and assuming quasi-static mechanics, the temperature evolution equation in non-isothermal PFM can be obtained. Finally, The governing equations are deduced and summarized as

[TABLE]

where $L^{\eta}$ is the kinetic coefficient characterizing the interfacial migration.

3.2 Proper MT

For the proper MT, the OPs $\varepsilon_{ii}^{0}$ are used to directly describe the austenite and martensite. The corresponding chemical free energy or Landau free energy density for 2D domain in cubic-tetragonal MT of Mn-22Cu can be represented as [78]

[TABLE]

where $e_{1}=(\varepsilon_{11}^{0}+\varepsilon_{22}^{0})/\sqrt{2}$, $e_{2}=(\varepsilon_{11}^{0}-\varepsilon_{22}^{0})/\sqrt{2}$, and $e_{3}=(\varepsilon_{12}^{0}+\varepsilon_{21}^{0})/2$. These coefficients ($Q_{1}$ to $Q_{5}$) in Landau free energy could be obtained by fitting the experimental results. For cubic to tetragonal MT, the eigen strains for martensitic variants are

[TABLE]

where $\varepsilon_{0}=0.02$ is the transformation strain in Mn-22Cu. Here, ($e_{1}=0,~{}e_{2}=0,~{}e_{3}=0$) represents for austenite, and ($e_{1}=0,~{}e_{2}=0.04,~{}e_{3}=0$) represents for martensitic variant 1 (V1), and ($e_{1}=0,~{}e_{2}=-0.04,~{}e_{3}=0$) represents for martensitic variant 2 (V2). For a first-order transition, the coefficients of $Q_{4}$ and $Q_{5}$ is positive and is needed for stability [79]. Here, we assume that the value of $Q_{2}(T),Q_{4}(T),Q_{5}(T)$ are similar to $A(T),B(T),C(T)$ in this work owing to the lack of experimental value for Mn-22Cu. Thus, $Q_{2}(T),Q_{4}(T),Q_{5}(T)$ are assumed as

[TABLE]

Similar to the case of the improper MT, the internal energy density, gradient energy density, and elastic strain energy density can be written as

[TABLE]

respectively.

The constitutive relations of proper MT PFM are analogous to Eq. 53, and the governing equations can be written as

[TABLE]

4 Finite-element implementation

Herein, we use finite element method to solve the governing equations in Eqs. 55 and 60, and convert the strong forms into weak forms by introducing a test function. Note that the degrees of freedom $\eta_{I}$ and $\varepsilon_{ii}^{0}$ are replaced by the generalized degree of freedom $\chi$ in the finite-element implementation. Therefore, the degrees of freedom are set as $u_{1},u_{2},u_{3},\chi,T$. Assuming $K_{ij}$ is a constant $K$ in Eqs. 55 and 60, the weak forms are formulated as

[TABLE]

where $\phi_{i}$, $\psi$ and $\vartheta$ are the test function for $u_{i}$, $\chi$ and $T$, respectively. Note that the surface terms ($\sigma_{ij}n_{j}$ and $T_{,i}n_{i}$) in Eq. 61 represent the surface traction and heat flux boundary conditions. $n_{i}$ is the normal vector of the boundary $\partial\Omega$.

By introducing the shape functions for independent variables and test functions, the discretized equations can be written as

[TABLE]

where $L$ denotes the node number. $N^{L}$ is the shape function. Here we assume quasi-static mechanics, and neglect $\zeta^{\chi\text{ex}}$, and do not consider the dynamic PFM. After the insertion of Eq. 62 into Eq. 61, the following elemental residuals can be obtained

[TABLE]

With regard to the time dependence of the residuals, we use the implicit backward Euler method to realize the time discretization [80]. The residual equation for the current time step $t_{n+1}$ is

[TABLE]

where $(\underline{\mathbf{d}}_{n+1}^{J}-\underline{\mathbf{d}}_{n}^{J})/\Delta t=\underline{\mathbf{\dot{d}}}_{n+1}^{J}$ and $\Delta t$ is time step. $\underline{\mathbf{d}}_{n+1}^{J}$ should be solved in this equation. For solving these non-linear equations, the Newton iteration scheme is performed at each time step. The corresponding iteration matrix is

[TABLE]

where $\underline{\mathbf{K}}^{LJ}$ is the stiffness matrix and $\underline{\mathbf{D}}^{LJ}$ is the damping matrix. We adopt the open source Multiphysics Object Oriented Simulation Environment (MOOSE) [81] to implement this PFM.

In order to evaluate the eCE of polycrystalline material, we make use of the rotation matrix to establish the polycrystalline PFM. The rotation matrices spinning around the $x$-axis, $y$-axis and $z$-axis are given by

[TABLE]

where $\alpha,\beta,\gamma$ are the Euler angles. The three-dimensional rotation matrix $\underline{\mathbf{Q}}$ is given by

[TABLE]

The relationship between the global and the local variables are

[TABLE]

in which $u$, $\varepsilon$ and $\sigma$ are the values in the global coordinate, while $u^{\prime}$, $\varepsilon^{\prime}$ and $\sigma^{\prime}$ are values in the grain local coordinate. $\underline{\mathbf{K}}_{\sigma}$ and $\underline{\mathbf{K}}_{\varepsilon}$ are the stress rotation matrix and strain rotation matrix, respectively. The eigen strain ($\varepsilon^{\text{tr}}_{ij}$) and elastic tensor ($C_{ijkl}$) should be rotated by $\underline{\mathbf{K}}_{\sigma}$ and $\underline{\mathbf{K}}_{\varepsilon}$. It is not difficult to deduce $\underline{\mathbf{K}}_{\varepsilon}$ and $\underline{\mathbf{C}}^{\prime}$ from Eq. 68 that

[TABLE]

where $\underline{\mathbf{F}}$ is a diagonal matrix with diagonal element values (1, 1, 1, 2, 2, 2). The stress rotation matrix ($\underline{\mathbf{K}}_{\sigma}$) is written in A.

5 Simulation results and discussions

In this work, we utilize a 2D domain to calculate the eCE for reducing computation cost. The Mn-22Cu material parameters used for PF simulations are taken from literature [44] and summarized in Table LABEL:t1. The material parameters of austenite and martensite in Mn-22Cu are regarded to be the same owing to their negligible differences. The finite element mesh size should be smaller than the minimum value of interface thickness ($\delta=\sqrt{\beta/2\Delta G^{*}}\approx 14.8$ nm) between austenite and martensite or between different austenite variants, thus the mesh size is chosen as $\Delta l=10$ nm. In the calculation of eCE by indirect and direct method, the uniaxial compressive stress is applied in the $y$-direction (favor the martensitic variant 2 (V2) formation), the left and right boundaries are set mechanically free (inserts in Fig. 6a and 7a), and adiabatic boundary conditions are specified by assuming zero heat flux.

5.1 Comparison of two models for improper and proper MT

The MT behavior and eCE in Mn-22Cu are modelled through two models, the ”OP-$\eta$” model using OPs $\eta_{I}$ for improper MT and ”OP-$\varepsilon$” model using OPs $\varepsilon_{ii}^{0}$ for proper MT. As shown in Fig. 2, under a 500 MPa uniaxial compressive stress, the stress-strain behavior and the temperature change during MT between these two models show good agreement. The consistency obviously prove the validity and generality of PF simulation in MT and eCE between the two models. The distinction in the temperature change during the beginning of loading is caused by the difference in PF kinetic coefficient $L$. In addition, the results of stress-strain behavior and microstructure evolution are consistent between these two models. In fact, Sun et al. [45, 46] have utilized the above two PFMs to simulate the cubic-tetragonal MT in SMA. Besides, the OPs $\eta_{I}$ are utilized in PFM built by Khachaturyan et al. [59, 82, 60] to model the proper and improper MT. This shows the correctness and validity of using OPs $\eta$ to model the proper or improper MT at least in term of phenomenological results.

The main discrepancy between the two models is the choice of OPs. The eigen strain $\varepsilon_{ij}^{0}$ for the two models is $\varepsilon_{ij}^{00}(I)\eta_{I}$ and $\varepsilon_{ii}^{0}$, respectively. Therefore, the transformation strains in the two models are all 0.02 ($\eta_{I}=1$ for the former, and $\varepsilon_{ii}^{0}=-0.02$ for the latter), when the austenite to martensite transformation completely occurs. For the cubic-tetragonal MT in Mn-22Cu, the two models are equivalent in PF simulation results. Note that the following results of MT behavior and eCE in PF simulation are based on the ”OP-$\eta$” model.

5.2 Benchmark simulation of ferroelastic behavior

Phase transformation, mainly including stress- and temperature-induced MT, would occur in Mn-22Cu SMAs under external fields, and leads to extraordinary macroscopic behaviors, e.g., the shape memory effect and superelasticity. In order to validate the model, benchmark simulations including that stress- and temperature-induced MT in isothermal (see Fig. 3) and non-isothermal PFM (see Fig. 4), and thermodynamic cycle (see Fig. 5) are carried out. For the single grain case in Fig. 3, austenite (A) with a small portion of V2 martensite embryo is set as the initial condition (inset of Fig. 3c). Figs. 3a and b show the stress- and temperature-induced MT in isothermal PF simulation, respectively. In Fig. 3a, a large enough compressive stress, along the easy axis of V2, would transform the initial A into V2. The mean value of $\eta_{2}$ ($\overline{\eta}_{2}$) increases as the stress is imposed and decreases as the stress is released, indicating the conventional MT during loading and the inverse MT during unloading. In Fig. 3b, upon cooling, high-temperature austenite would turn into low-temperature martensite, or vice versa. The superelasticity effect during stress-induced MT is also clearly presented in Fig. 3c. Upon stress loading, the small martensitic embryo (red stripes) grows to be a large martensitic domain owing to the stress-induced MT. Upon unloading, the austenitic embryo (blue stripes) grows and the original state restores.

Simulating, the superelasticity effect also can be used as a benchmark for validating the non-isothermal PFM for polycrystalline, as shown in Fig. 4. The polycrystalline model size is $1000~{}\text{nm}\times 1000~{}\text{nm}$ and has nine grains with random orientations. The maximum compressive stress (along the easy axis of V2) is 300 MPa and the initial temperature is 245 K. During the loading, the strain increases linearly with the stress at first, corresponding to the elastic response of the austenite. In this regime, there is no MT and temperature and $\overline{\eta}_{2}$ change a little. As the stress increases to 230 MPa, the stress-induced MT occurs, V2 bands start to grow, and temperatures apparently changes (inset of Fig. 4a). As shown in Fig. 4b, the temperature and $\overline{\eta}_{2}$ rapidly rise and then recover to the initial values upon unloading. The temperature is increased from 245 to 255 K, resulting in $\Delta T_{\text{ad}}=10$ K. The temperature distribution is almost homogeneous because of the high thermal conductivity and no thermal barrier in the domain wall. The temperature is found to be approximately proportional to order parameters, which is consistent with the results reported in [44].

The refrigeration cycle of eCE is schematically shown in the insert of Fig. 5. When a SMA in the austenitic (cubic) phase is axially stressed or strained, an exothermic austenitic-martensitic transformation occurs, which under adiabatic conditions makes the material heated up (① in Fig. 5). This heated material then releases heat to the surroundings and cools down to the ambient temperature (② in Fig. 5). When the stress is removed, the crystal structure transforms back to the austenitic phase (③ in Fig. 5). Finally, the material cools down and is now able to absorb heat from the surroundings (④ in Fig. 5). The simulated temperature vs time profile under loading and unloading is also shown in Fig. 5, showing good agreement with the above process from ① to ④. These simulation results of MT and thermodynamic cycle indicate that eCE can be soundly handled by the isothermal and non-isothermal PFM.

5.3 Indirect method

The indirect method is usually based on the data of superelastic response, such as the stress-strain curves at various temperatures or strain-temperature curves at various stresses. According to the Eqs. 2 and 3, indirect method is used to calculate the isothermal entropy change $\Delta S_{\text{iso}}$ and adiabatic temperature change $\Delta T_{\text{ad}}$ by means of the isothermal PFM. We investigate the eCE under 50, 100, 150, 200, 250, and 300 MPa compressive stress, which is lower than the maximum stress that Mn-22Cu alloy can endure.

5.3.1 eCE of single crystal

The single crystal model size is $300~{}\text{nm}\times 500~{}\text{nm}$. The typical stress-strain curves at different temperatures from 215 to 275 K are shown in Fig. 6a, which clearly indicates three stages involving elastic stage in austenite, transformation stage, and elastic stage in martensite. A minimum critical transformation stress is found at T = 245 K, because the energy barrier between austenite and martensite vanishes there. Fig. 6b shows the strain-temperature curves at different compressive stresses from 50 to 300 MPa. During the loading of 100 MPa stress, the strain firstly notably increases and then significantly decreases, indicating that both the inverse and conventional MT occur. By means of the Eqs. 2 and 3, $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ can be calculated, as summarized in Figs. 6c and d, respectively. Large $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ observably appear near the transformation temperature (245 K). This is caused by the first-order phase transformation which releases or absorbs lots of heat at the transformation temperature.

The values of $\Delta S_{\text{max}}$ (J/kg/K) and $\Delta T_{\text{max}}$ (K) at varied stresses are: $\Delta S_{\text{max}}=-3.3$ and $\Delta T_{\text{max}}=2.3$ for $\sigma=50$ MPa, $\Delta S_{\text{max}}=-13.5$ and $\Delta T_{\text{max}}=9.5$ for $\sigma=100$ MPa, $\Delta S_{\text{max}}=-27.8$ and $\Delta T_{\text{max}}=19.5$ for $\sigma=150$ MPa. The maximum $\Delta T_{\text{ad}}$ 9.5 K at 100 MPa is consistent with the experimentally measured value of 11.6 K [56]. It is found that $\Delta S_{\text{max}}$ and $\Delta T_{\text{max}}$ do not increase further for the stress beyond 150 MPa, because of the already complete MT under 150 MPa stress. However, the operating temperature window is effectively broadened if the stress increases. In addition, there exists an inverse/negative eCE, i.e., positive $\Delta S_{\text{iso}}$ or negative $\Delta T_{\text{ad}}$ under an external stress. This is a consequence of the inverse MT occurring in the case of co-existed austenite and martensite when the ambient temperature is below $T_{0}$ [8, 83, 84]. The further explanation in terms of microstructure for the negative eCE will also be discussed in Fig. 7e.

5.3.2 eCE of polycrystal

Fig. 7 shows the simulation results for eCE of polycrystal. $\Delta T_{\text{ad}}=3.5$ K for $\sigma=100$ MPa is close to the measured 3.9 K under 4.0$\%$ strain  [56]. In general, $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ in polycrystalline are lower than that in single crystal. The internal interaction caused by the grain boundary, grains with adverse orientations, and the large negative eCE together contribute to the low $\Delta S_{\text{iso}}$ and $\Delta T_{\text{ad}}$ in polycrystal. Besides, the large negative $\Delta T_{\text{ad}}$ of polycrystal (-13 K for 200 MPa) in Fig. 7d results from the inverse MT and the local large compressive stress due to the grain boundary, as shown the microstructure evolution in Figs. 7e and f. Below 245 K, the initial phase is martensite, and the inverse MT would occur under a large enough compressive stress. In the case of 150 MPa and 241 K (Fig. 7e), V1 is changed into A at first, and then into V2. This inverse MT absorbs an amount of heat, resulting in the large inverse eCE. At 248 K (Fig. 7f), the initial austenite phase is transformed into V2 accompanied with a small portion of V1.

5.4 Direct method

In experiment or numerical simulation, the straightforward way to assess eCE is the direct method, which directly measures the $\Delta T_{\text{ad}}$ during loading and unloading. Cissé et al. [40] and Sun et al. [46] used the non-isothermal PFM to directly calculate $\Delta T_{\text{ad}}$ in CuAl11Be2 and FePd SMAs, respectively. Here we use the modified non-isothermal PFM to evaluate eCE in Mn-22Cu alloys based the direct method. The applied compressive stress is 100 MPa, the loading rate is 0.2 MPa/ns, and other conditions are consistent with the indirect method for a comparison.

5.4.1 eCE of single crystal

For single crystal, the strain-stress response and the adiabatic temperature change calculated by the direct method are shown in Fig. 8. It can be found in Fig. 8a that the sample almost stays at the linear elastic regime at 250 K greater than the equilibrium temperature. As the initial temperature decreases to 245 K, the critical transformation stress decreases and the area with MT increases, leading to in a larger $\Delta T_{\text{ad}}$ (Fig. 8b). The transformation strain is 0.8$\%$ at 245 K under 100 MPa compressive stress and the corresponding $\Delta T_{\text{ad}}$ is around 8.5 K. $\Delta T_{\text{ad}}$ in Fig. 8b shows a strong temperature dependence. In other words, controlling the ambient temperature to obtain high transformation strain is a feasible neat idea to improve the eCE of SMAs.

5.4.2 eCE of polycrystal

The settings of polycrystalline model are consistent with those in the indirect method. Fig. 9a and b show the stress-strain response and $\Delta T_{\text{ad}}$ at different temperatures, respectively. Under the same stress, the maximum transformation strain and $\Delta T_{\text{ad}}$ occurs at $T$= 245 K and the corresponding $\Delta T_{\text{ad}}$ is around 3.8 K. Polycrystal is more closer to the microstructure of experimental samples and thus its maximum $\Delta T_{\text{ad}}$ of 3.8 K by the direct method and 3.5 K by the indirect method is more closer to the experimental data (3.9 K) [56]. Compared to the maximum $\Delta T_{\text{ad}}$ of 8.5 K in single crystal, the smaller $\Delta T_{\text{ad}}$ in polycrystal is probably a consequence of domain wall and adverse orientation grain in the polycrystalline model. The grain misorientation may hinder the martensitic nucleation and growth, as shown in Fig. 4a. In addition, comparing Fig. 8a with Fig. 9a, it is obvious that a smaller hysteresis exists in polycrystal.

5.5 Indirect vs direct method

In order to compare the indirect and direct method which are based on PFM, in Fig. 10 we present the $\Delta T_{\text{ad}}$-$T$ curves calculated by both methods under 100 MPa compressive stress. It is clear that the overall trend of $\Delta T_{\text{ad}}$ varying with $T$ is consistent for both methods. For single crystal in Fig. 10a, $\Delta T_{\text{ad}}$ calculated by the two methods at different ambient temperatures is close. Under 100 MPa, the maximum $\Delta T_{\text{ad}}$ ($\Delta T_{\text{ad}}^{\text{max}}$) calculated by indirect and direct method is 9.5 and 8.5 K (3.5 and 3.8 K) for single crystal (polycrystal), respectively. This indicates the consistency between indirect and direct method and the reliability of calculating eCE by PF simulations. However, the peak temperature corresponding to $\Delta T_{\text{ad}}^{\text{max}}$ calculated by the two methods is slightly different. For the direct method, the peak temperature is around 245 K (the chemical equilibrium temperature), which is lower than the indirect method (247 K for single crystal and 248 K for polycrystal).

The insets in Fig. 10 present the $\Delta T_{\text{ad}}^{\text{max}}$ as a function of the applied compressive stress. It can be seen that eCE can be significantly tuned by the uniaxial compressive stress. Results from the direct method is essentially in agreement with those from the indirect method when the stress is less than 150 MPa, with a discrepancy of $\Delta T_{\text{ad}}^{\text{max}}$ within 10$\%$. At higher compressive stress (above 150 MPa), $\Delta T_{\text{ad}}^{\text{max}}$ calculated by the indirect method is apparently larger than that form the direct one. This could be attributed to the temperature-induced MT in the non-isothermal PFM simulation, i.e., the relatively large temperature increase at high stress could in turn remarkably hinder the austenite-martensite transition and thus lower the temperature [85]. Cheng et al. [86] also report the similar phenomena in electrocaloric effect. In addition, $\Delta T_{\text{ad}}/\sigma_{\text{max}}$ is usually utilized as a parameter to evaluate the field normalized caloric effect, also known as the specific adiabatic temperature [87]. It is found that $\Delta T_{\text{ad}}/\sigma_{\text{max}}$ is 0.131 and 0.091 K/MPa for single crystal (0.054 and 0.037 K/MPa for polycrystal) by the indirect and direct method, respectively. A large $\Delta T_{\text{ad}}/\sigma_{\text{max}}$ would be beneficial to enhance the overall performance and efficiency of the elastocaloric refrigerator.

The agreement between the calculated and experimental results is satisfactory, implying that the direct method based on the non-isothermal PFM and the direct method based on isothermal PFM are reliable for the calculation and evaluation of eCE in SMAs. Nevertheless, the computational cost and the available information on eCE are also critical factors to be considered. The great advantage of the indirect method is that all the $\Delta S_{\text{iso}}$, $\Delta T_{\text{ad}}$, and refrigerating capacity can be readily obtained. But the indirect method depends on the temperature dependent strain-stress curves with small temperature intervals for an accurate integration, leading to heavy computation. If one is only interested in $\Delta T_{\text{ad}}$, the direct method is the best choice since it doe not require the overall stress-strain-temperature data and thus is computationally efficient.

6 Conclusion

In summary, we have developed a thermodynamically consistent non-isothermal PFM for the simulation of eCE. The PFM couples MT with mechanics and heat transfer to evaluate eCE by using the indirect and direct method. The model considering both improper and proper MT is derived from a thermodynamic framework which invokes the microforce theory accommodating non-local effects, thermodynamic laws, and the Coleman–Noll procedure. In order to avoid the problematic calculation due to the non-differentiable energy barrier function across the transformation temperature and consider the possibly non-sharp transition in real materials, the austenite-martensite transition energy barrier in PFM is introduced as a smooth function of temperature by using the hyperbolic tangent function. On one hand, all constitutive relations are represented in terms of a thermodynamic potential. Therefore, the PF modelling work is reduced to the design of a proper form of the thermodynamic potential. On the the hand, the framework automatically satisfies the first, second, and third laws of thermodynamics.

After being numerically implemented by the finite element method, the developed PFM is demonstrated to be capable of recapturing the microstructure response and calculating eCE via both indirect and direct method in a model material (i.e., Mn-22Cu SMA) under an external loading. Under a compressive stress of 100 MPa, $\Delta T_{\text{ad}}^{\text{max}}$ calculated by the indirect and direct method is 9.5 and 8.5 K for single crystal (3.5 and 3.8 K for polycrystal), respectively. A large $\Delta T_{\text{ad}}$ exists in single crystal, but the working temperature window is narrow and can be improved by increasing the compressive stress. Besides, negative eCE caused by the inverse MT is found, especially for polycrystal, which reduces $\Delta T_{\text{ad}}^{\text{max}}$. Overall, the direct method based on the non-isothermal PFM and the direct method based on isothermal PFM are reliable for the calculation and evaluation of eCE in SMAs. But there are still some differences between these two methods. $\Delta T_{\text{ad}}$ calculated by the indirect and direct method shows tiny discrepancy (within 10$\%$) under a low stress ($\leq$150 MPa). At higher stress, $\Delta T_{\text{ad}}^{\text{max}}$ from the indirect method is apparently larger than that from the direct one, mainly owing to the fact that the relatively large temperature increase at high stress in the non-isothermal PFM could in turn remarkably hinder the austenite-martensite transition and thus lower the temperature. The indirect method is computational expensive, but can yield all the $\Delta S_{\text{iso}}$, $\Delta T_{\text{ad}}$, and refrigerating capacity. The direct method is computationally efficient, but only yields $\Delta T_{\text{ad}}$. The developed non-isothermal PFM along with the indirect and direct method could be a computational toolkit to unveil strategies for designing high-performance elastocaloric devices.

CRediT authorship contribution statement

Wei Tang: Formal analysis, Investigation, Visualization, Data curation, Writing – original draft. Qihua Gong: Conceptualization, Resources, Supervision, Investigation, Data curation, Writing – original draft & review & editing. Min Yi: Conceptualization, Resources, Supervision, Project administration, Funding acquisition, Writing – original draft & review & editing. Bai-Xiang Xu: Conceptualization, Supervision, Writing – original draft & review. Long-Qing Chen: Conceptualization, Supervision, Writing – original draft & review.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors acknowledge the support from National Natural Science Foundation of China (NSFC 12272173, 11902150), 15${}^{\text{th}}$ Thousand Youth Talents Program of China, Fundamental Research Funds for the Central Universities (1001-XAC21021), the Research Fund of State Key Laboratory of Mechanics and Control of Mechanical Structures (MCMS-I-0419G01), and a project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.

Appendix A Rotation matrix

[TABLE]