Dynamics of Shape Memory Alloys Patches with Mechanically Induced Transformations

A mathematical model is constructed for the modelling of two dimensional thermo-mechanical behavior of shape memory alloy patches. The model is constructed on the basis of a modified Landau-Ginzburg theory and includes the coupling effect between thermal and mechanical fields. The free energy functional for the model is exemplified for the square to rectangular transformations. The model, based on nonlinear coupled partial differential equations, is reduced to a system of differential-algebraic equations and the backward differentiation methodology is used for its numerical analysis. Computational experiments with representative distributed mechanical loadings are carried out for patches of different sizes to analyze thermo-mechanical waves, coupling effects, and 2D phase transformations.

1. Introduction. The Shape Memory Effect (SME) in certain metallic systems, known as Shape Memory Alloys (SMA), has been the subject of considerable theoretical and experimental research efforts. Upon the action of thermal, magnetic, electrical, hydrostatic and other fields, the SMAs can recover their original shape after being permanently deformed [19]. These materials are an intrinsic part of the smart material and structure technology. They can directly transduce thermal energy into mechanical and vice versa, making them very attractive in micro-sensor and actuator applications. Their application areas include also biomedicine, communication industries, robotics to name just a few.
New promising engineering applications of SMAs and their ongoing proliferation out paced, at least initially, the development of adequate mathematical models and physical theories. At the same time, it is clear that a better understanding of the dynamic behavior of SMA structures is paramount for further advances in engineering applications of SMAs. Such an understanding can be gained with mathematical modelling based on physically plausible theories and experimental observations, and increasing research efforts are directed to this area [2]. At present, the unique properties of the SMA are viewed to be a consequence of its unique ability to undergo reversible phase transformations and reversible martensitic reorientations when it is subjected to appropriate loadings. Hence, a series of concerted efforts have been directed to the analysis of such phase transformations and reorientations in order to bridge the gap between engineering applications and our ability to predict the response of the material under various loadings.
However, even for the one dimensional case, the analysis of the dynamics of SMA structures is quite involved due to multivariant martensitic phase transformations, strong nonlinearity in the mechanical field, and a strongly coupled nonlinear pattern of interaction between mechanical and thermal fields (e.g., [2,18] and references therein). To analyze thermo-mechanical interactions numerically is a highly nontrivial task even in the case where the phase transformations are not included into the model due to the fact that strong nonlinearity and thermo-mechanical coupling have still to be dealt with. For many practical applications, one dimensional models are able to satisfactory approximate the dynamics of the SMA rods and wires. However, for many other applications a better understanding of the dynamics of SMA structures in dimensions higher than one is required. For instance, a number of practical applications lead to the usage of SMA membranes embedded in a microvalve, SMAs thin film, and so on. In such cases, it is essential to develop new models that are applicable for the description of nonlinear thermo-mechanical waves and phase transformations in higher dimensions.
Up to date, many instructive investigations have been carried out to understand the dynamics during the process of martensitic phase transitions. They provide a firm background for the application development, in particular in the one-dimensional cases where the modelling of shape memory alloys can be based on the Landau-Ginzburg free energy functional (e.g., [3,7,17] and references therein). Although considerable efforts have been devoted to the development of a general theory for the dynamics of SMAs, most of the results are applicable to the one dimensional case only ( [7,21] and references therein), while the results concerning two or three-dimensional cases are rarely available in the literature in the context of modelling thermo-mechanical waves and phase transformations in SMAs. Several theoretical contributions to the development of free energy functionals that play an essential role in the modelling of SMA dynamics have been reported in the literature (e.g., [5,21]. Such results were reported for the three dimensional case and were developed on the basis of a modified Landau-Ginzburg theory, but their numerical implementation was not reported. Recently, for the simulation of 2D microstructures in ferro-elastic materials, several models were proposed by reduction procedures applied to the general 3D model [10,11]. No numerical analysis of thermo-mechanical coupling was offered and only static mesoscale simulations of microstructures with fixed temperature were presented.
In this paper, we propose a two-dimensional macroscopic dynamical model for the description of thermo-mechanical behavior of SMA structures. Our free energy function, constructed for 2D dynamical models that are applicable at the mesoscopic level, is similar to those previously discussed in [1,13,14]. We apply this new model to the numerical analysis of thermo-mechanical behavior of 2D SMA structures under various mechanical loadings. We demonstrate that phase transformations in SMAs can be successfully predicted with the proposed model.

2.
Square to Rectangular Transformations. The main purpose of the numerical analysis of SMAs is to gain a better understanding of the dynamics of SMA structures. Since phase transformations are responsible for SMA unique properties, this implies that the developed mathematical models should be capable to capture three-dimensional cases. Examples of the alloys exhibiting such transformations include N b 3 Sn, InT l, F eP d, some copper based SMAs and many others (e.g., [11,22,14] and reference therein). Hence, studies of the specific transformation chosen here are often regarded as a prerequisite for a better understanding of the cubic to tetragonal and tetragonal to orthorhombic transformations. Up to this point, a number of numerical studies of the dynamics of phase transformations have been carried out, including those concerned with the formation and growth of microstructures. They explained various phenomena such as the formation of tweed microstructures observed by a variety of techniques ( [12,22,11,10,14] and references therein). These studies concentrated primarily on mesoscale pattern formation and microstructures evolution.
According to Landau's theory of phase transformations, the basis of any nonlinear continuum-thermodynamics model for phase transformations is a non-convex free energy functional. This idea has been used by many authors (e.g., [5,1,14] among others). The global and local minima of the free energy potential with respect to the strain tensor (or deformation gradients) correspond to the stable and mesostable states at a given temperature. In the context of SMAs, in particular when we are concerned with 1D simulations, we can use the Helmholtz free energy Ψ and its related strain energy W, as it has been demonstrated by a number of authors ( [6,7,18,16,3] and references therein): where ψ 0 (θ) models thermal field contributions, ψ 1 (θ)ψ 2 (ǫ) models shape memory contributions, and ψ 3 (ǫ) models mechanical field contributions. The last two contributions can be simulated by one free elastic energy function which is dependent on the temperature. Recall that ǫ = ∂u/∂x (with u being the displacement)is the strain in this case and it is chosen as the only order parameter for the phase transformations. The functions ψ 0 could take the following form: where c v is the specific heat constant, θ 0 is the reference temperature for the transformation.
To characterize the austenite at high temperature and the martensite variants at low temperature in the SMAs, a free elastic energy functional F was established earlier for the square to rectangular transformations, in which the Landau free energy function F l [10,11] was modified. Note that for the square to rectangular transformation, there is only one order parameter needed [7,7,14] to characterize the martensite variants in a 2D space, as shown by Fig. 1(a). As it has been discussed in detail before (e.g., [14,10,11,22,12] and references therein), a simple free elastic energy functional could be chosen as a function of the order parameter in the following form: . . , 6, d 2 , and d 3 are the materialspecific coefficients, and e 1 , e 2 , e 3 are dilatational, deviatoric, and shear components of the strains, respectively, which are defined as follows: The Cauchy-Lagrangian strain tensor η is given by its components in the standard manner with the repeated-index convention used: where u i is the displacement in the i th direction in the coordinate system, and x is the coordinates of a material point in the domain of interest. In this free energy functional, the deviatoric strain e 2 is chosen as the order parameter.
In the above elastic free energy function F , the Ginzburg term F g is the term proportional to the square of the strain gradients. It produces an energy cost for deviations from spatial uniformity such as in the presence of domain walls. This term is necessary for the simulation of microstructures and phase growth ( [11,14]and references there in), but here we are interested in the macroscopic behavior of the SMAs. Since only the dynamical properties of the material before and after the phase transformation are of interest, the simulation scale does need to be as fine as being able to reproduce the domain wall movement, and therefore the Ginzburg term can be ignored in this case.
At the same time, the Landau free energy function F l should be able to be converted into a triple well, double well, or convex potential, depending on the specific temperature and the material-specific constants. To include the temperature dependency of the free energy functional, the material parameter A 2 is assumed to be linearly dependent on the material temperature A 2 = a 2 (θ − θ 0 ). Now, if we assume that the shape memory contributions and mechanical field contributions ψ 1 (θ)ψ 2 (ǫ) + ψ 3 (e i ) are represented by the above mentioned Landau free elastic energy functional, and take the thermal contribution ψ 0 the same as in the 1D case, the final form of the Helmholtz free energy function for the square to rectangular transformations will take the following form: This representation can be easily interpreted in a standard way. In Fig.1 (b), we plotted the Landau free energy function for a SMA sample based on the material Au 23 Cu 30 Zn 47 as a function of the deviatoric strain e 2 for different values of temperatures. It is easy to see that the function has two local minima at lower temperatures (210 o K), which correspond to two rectangular martensite variants, while only one minimum at the center corresponds to the square austenite phase when the temperature is high (270 o K). When the temperature is in between (e.g., around 245 o K in Fig. 1b), there are three local minima, which indicates the existence of metastable phases.
In the next section, the above free energy functional will be used to derive governing equations for the dynamics of SMA patches with square to rectangular transformations.

A 2D Model for SMA Dynamics.
To obtain the governing equations describing the dynamical behavior of SMA patches, we start our discussion from physical conservation laws. Due to coupling between mechanical and thermal fields, we have to take into account the mass, momentum, and energy balances simultaneously. In its general form, the system describing coupled thermo-mechanical interactions for the first order phase transformations in a 2D SMA structure can be written as follows [21,19]: where ρ is the density of the material, u = {u i }| i=1,2 is the displacement vector, v = ∂u/∂t is the velocity vector, σ = {σ ij } is the stress tensor, q is the heat flux, e is the internal energy, f = (f 1 , f 2 ) T and g are mechanical and thermal loadings, respectively.
Having constructed the free energy functional Ψ as in Eq.2.6 in section 2, we note that the internal energy e is connected with the free energy function Ψ by the following relationship: This gives the internal energy as function of e i and θ: Finally, we have to specify the constitutive relationships that couple stresses, strains (deformation gradients), temperature and heat fluxes. We follow [19] in doing so: where it is implicitly assumed that these relations may involve spatial and temporal derivatives of the functions. For the constitutive relationship Φ 2 , the Cattaneo-Vernotte model can be used to capture a finite speed of the thermal wave propagation: An approximation to this general case is provided by the generalized form of the Fourier law, but all simulations reported in this paper has been performed for the classical case of q = −k∇θ (3.12) with constant heat conductivity k. Further we assume that the constitutive relationship Φ 1 between stresses and strains are given by the following equation [7,21]: By substituting the free energy functional defined by Eq.(2.6) into the above equation, all the entries in the stress tensor can be obtained as follows: ρ(a 1 e 1 + a 2 (θ − θ 0 )e 2 − a 4 e 3 2 + a 6 e 5 2 ), . (3.14) By assembling all the equations discussed above together, the governing equations for the dynamics of the shape memory alloys patches with square to rectangular transformation can be written as: It is expected that this 2D mathematical model should be capable to capture thermo-mechanical interactions and phase transformations in the 2D SMA structures. This is verified in Section 5 by a series of computational experiments. Model (3.15) is supplemented by appropriate boundary and initial conditions which are problem-specific and will be discussed in the next sections.
4. Numerical Methodology. Since the SMA can be in a high temperature phase (austenite) as well as in a low temperature phase (martensite), at the computational level one faces a fairly complex task of dealing with different equilibrium configurations of the metallic lattice simultaneously. Therefore, it is important to preserve intrinsic properties of model (3.15). As it is shown in [16] for the 1D model, a conservative scheme can be constructed and rigorous a priori estimates for the solution of such a scheme can be derived. A similar numerical approximation can be derived by applying formally the Finite Volume Method (FVM) which is easier generalizable to a higher dimensional case. For convenience, two velocity components are introduced into the model, and then the displacement components are replaced by the dilatational and deviatoric strains as follows: (4. 16) In the 2D model, there are three strain components e 1 , e 2 , e 3 , which are dependent on the two displacements u 1 , u 2 . So, after the replacement, we need one more equation to close the model. This is given by the compatibility relation in terms of the strains [11,14]:

LINXIANG WANG AND RODERICK MELNIK
After simple transformations, the model (3.15) can be written as follows:  The 2D model given by Eq.(4.18) is a differential-algebraic system. It is obtained by keeping the constitutive relations as algebraic equations with the stress components kept as independent variables to be solved for. The idea of simulating the thermo-mechanical waves by the differential-algebraic approach is stimulated by the idea of [17] where the same approach was adopted for the simulation of phase transformations in SMA rod.
The DAE system (4.18) is solved numerically by the method of lines. First, the model is spatially discretized, and then a time integrator is employed to the resultant system. Velocity components v 1 and v 2 are discretized at (x i+1/2 , y j+1/2 ), x i+1/2 = (i + 1/2)h x , y j+1/2 = (j + 1/2)h y where h x and h y are the grid sizes along the x and y directions, respectively, i = 1, 2, ..., m x and j = 1, 2, ..., m y with m x and m y being the numbers of discretization points in the x and y directions. Variables e 1 , e 2 , θ, σ 11 , σ 12 , σ 22 are discretized at (x i , y j ). After the discretization in space, the system can be re-cast in the following form: de 1 (i, j) dt = (I y D x v 1 (i + 1/2, j + 1/2) + I x D y v 2 (i + 1/2, j + 1/2))/ √ 2, de 2 (i, j) dt = (I y D x v 1 (i + 1/2, j + 1/2) − I x D y v 2 (i + 1/2, j + 1/2))/ √ 2, dv 1 (i + 1/2, j + 1/2) dt = I y D x σ 11 (i + 1, j + 1) + I x D y σ 12 (i + 1, j + 1) + f 1 , dv 2 (i + 1/2, j + 1/2) dt = I y D x σ 12 (i + 1, j + 1) + I x D y σ 22 (i + 1, j + 1) + f 2 , (4.19) where D x and D y are the first order difference operators in the x and y directions, respectively, while I x and I y are discrete interpolation operators in the x and y directions, △ is the discrete Laplace operator. The above DAE system (4.19) can be written in the following operator-matrix form: with matrix A = diag(a 1 , a 2 , ..., a N ) having entries "one" for differential and "zero" for algebraic equations for stress-strain relationships, and vector-function H defined by the right hand side parts of (4.19) based on the spatial discretisation of the original PDE system. This (stiff) system is solved with respect to the vector of unknowns U having 6 × m x × m y + 2 × (m x + 1) × (m y + 1) components by using the second order backward differentiation formula (e.g., [8]): where n denotes the current time layer. We note that to deal with strong (quintic) nonlinearities in the order parameter, a smoothing procedure similar to that proposed in [20,16] has been employed during the integration process. In particular, we have used the following expansions: where y n is the unknown variable on the current time layer n, while y n−1 is the unknown variable on the previous time layer n−1 (for the 2D case, y = e 2 , while for the 1D case y = ǫ = ∂u/∂x). We summarize this smoothing procedure as follows.
Nonlinear terms are averaged here in the Steklov sense, so that for nonlinear function f (y) (in particular, y 3 and y 5 ), averaged in the interval [ǫ n−1 , ǫ n ], we have g(y n−1 , y n ) = 1 y n − y n−1 yn yn−1 f (η)dη, y n−1 = y(t n−1 ), y n = y(t n ). (4.23) Finally, we note that the analysis of mathematical models for the description of SMAs and numerical methodologies for their solution are the subject of a separate discussion that goes beyond the present paper. For one dimensional models we refer the interested reader to recent papers [9,16] while mathematical models for the general 3D case have been recently discussed in [15].

Computational Experiments on
Copper-based SMAs. The above constructed mathematical model and the numerical methodology have been applied to modelling the dynamics of SMAs patches, including the thermo-mechanical interactions, as well as the nonlinear thermo-elastic phase transformations. The capability of the mathematical model has been demonstrated by the numerical experiments reported below. First, we report the results of simulations for Au 23 Cu 30 Zn 47 . For this specific material, its physical parameters for 1D Falk model are available from the literature [6,18,20]: k 1 = 480 g/ms 2 cmK, k 2 = 6 × 10 6 g/ms 2 cmK, k 3 = 4.5 × 10 8 g/ms 2 cmK, θ 1 = 208K, ρ = 11.1g/cm 3 , C v = 3.1274g/ms 2 cmK, k = 1.9 × 10 −2 cmg/ms 3 K.
For the 2D model given by Eq. (4.18), there are no 2D experimental values for the physical parameters. Here, we take all the parameters in the Landau free energy function in the 2D model identical to those reported previously in [16,24], which gives us: a 2 = k 2 , a 4 = k 2 , a 6 = k 3 . Then, we take the values a 1 = 2a 2 and a 3 = a 2 , as suggested in [11,4]. All thermal parameters are taken the same as in the 1D case.
The first numerical experiment aims at the numerical analysis of the dynamical thermo-mechanical response of the SMA patch under varying distributed mechanical loadings, small enough not to induce phase transformations. The purpose of this experiment is to verify the capability of the model to capture the nonlinear thermomechanical behavior. The SMA patch covers an area of 1 × 1cm 2 .
The initial conditions for this simulation are taken as: The second numerical experiment aims at the simulation of square to rectangular phase transformations under a periodic mechanical loading in the same SMA patch. The one-period profile of the loading is taken as follows: 4 ≤ t ≤ 6, sin(π(t − 2))/3), 6 ≤ t ≤ 10, 0, 10 ≤ t ≤ 12. (5.26) The computational results in the 1D case have previously shown that the mechanical loading (5.26) is strong enough to induce martensitic transformations (e.g., [18]). The initial and boundary conditions for this experiment were taken identical to (5.24) and (5.25).
Similarly to the analysis of phase transformations in the 1D case, we can calculate the deviatoric strain that corresponds to austenite and martensite variants in a SMA patch. In doing so, we have to minimize the Landau free energy functional, and the condition ∂F l /∂e 2 = 0 yields:    where dθ is the difference between the current material temperature and the transformation temperature. We observe that e 2 = 0 corresponds to the austenite. Let us denote (a 4 + a 2 4 − 4a 2 dθa 6 )/2a 6 by e m , then e 2+ = + √ e m or e 2− = − √ e m are the strains that correspond to the two martensite variants. Let us call them martensite plus and martensite minus, respectively. For the material in hand, with the given initial temperature, we estimate that dθ = 42 o , and this, in its turn, gives e 2+ = 0.12 and e 2− = −0.12.
The numerical results of this simulation are presented in Fig.3. We present the strain e 2 , as well as u 1 , u 2 , and θ on the horizontal line as functions of time. We provide also two representative snapshots of e 2 , and θ. From the time variation of the order parameter e 2 , we observe a periodic pattern in the phase transformations induced by periodic loading. Analyzing the e 2 plot at t = 2 we note that the SMA patch is divided into two sub-domains: in the upper-left triangle-like area, the simulated deviatoric strain is close to e 2+ which corresponds to the martensite plus, while on the opposite side, the deviatoric strain is close to e 2− which corresponds to the martensite minus. Between these two areas there is an interface. At t = 8, the martensitic transformation occurs again, but the martensite combination is switched due to reverse pattern of loading.
Our next example concerns a SMA sample of rectangular size 1 × 0.5 with all conditions identical to the previous experiments. We took only 19 nodes in the x direction and 9 nodes in the y direction. We expect that after the transformation, the martensite combination should be similar to its counterpart in the experiment 2, but scaled along the y direction. This is demonstrated in Fig.4. If we apply the mechanical loading in the x direction only, the symmetry along the diagonal will be broken, but the symmetry along the middle vertical line x = 0.5 will be preserved. This has been also confirmed in our computational experiments.
Finally, we apply our methodology to modelling the AlCuZn patches. The physical parameters for this material are available in the literature for the 1D case (e.g., [3]). In the numerical experiment reported here, we enlarged the specific heat by factor of 10. This material has a small specific heat and the temperature increases quickly upon mechanical loadings, if the boundary are insulated, which prevents us from observing a stable martensite otherwise. Our results relate to a 1 × 1cm 2 AlCuZn patch with the following initial conditions: and the boundary conditions identical to (5.25). The loading follows the pattern of (5.26) with factor 6000 replaced by 1200. The time span for this simulation is [0, 24] which corresponds to two periods of loading. There are 15 nodes in each direction, and the time stepsize is set to 1 × 10 −4 . The deviatoric strain e 2 , temperature θ, and displacements u 1 , u 2 on the central horizontal line are presented in the upper four plots of Fig.6. Two snapshots of the distribution of e 2 over the whole patch are also presented in the bottom of Fig.6. Although the material does not exhibit cubic-to-tetragonal or tetragonal-to-orthorhombic transformations in the 3D case, in the 2D case conclusions similar to those derived for AuCuZn patches can be drawn.