arXiv:cs/0702168v1 [] 28 Feb 2007

Simulation of Phase Combinations in Shape Memory
Alloys Patches by Hybrid Optimization Methods

Linxiang Wang a and Roderick V. N. Melnik b
a Mads

Clausen Institute, Faculty of Engineering,
University of Southern Denmark,
Sonderborg, DK-6400, Denmark

b Mathematical

Modelling and Computational Sciences
Wilfrid Laurier University, Waterloo,
75 University Ave W, Canada N2L 3C5

Abstract
In this paper, phase combinations among martensitic variants in shape memory alloys patches and bars are
simulated by a hybrid optimization methodology. The mathematical model is based on the Landau theory of
phase transformations. Each stable phase is associated with a local minimum of the free energy function, and
the phase combinations are simulated by minimizing the bulk energy. At low temperature, the free energy
function has double potential wells leading to non-convexity of the optimization problem. The methodology
proposed in the present paper is based on an initial estimate of the global solution by a genetic algorithm,
followed by a refined quasi-Newton procedure to locally refine the optimum. By combining the local and
global search algorithms, the phase combinations are successfully simulated. Numerical experiments are
presented for the phase combinations in a SMA patch under several typical mechanical loadings.

Key words: Phase combinations, shape memory alloys, variational problem, genetic algorithm,
quasi-Newton methods.

Email addresses: wanglinxiang@mci.sdu.dk (Linxiang Wang), rmelnik@wlu.ca (Roderick
V. N. Melnik).

Seventh IMACS International Conference on Iterative Methods in Scientific Computing, May 5-8, 2005

2

L. X. Wang, R. V. N. Melnik

1 Introduction

Shape Memory Alloys (SMA) are materials with increasing range of applications in engineering,
aerosapce, and biomedical industries. They possess unique properties of being able to recover their
original shape after permanent deformations. These materials can directly transduce thermal energy into mechanical and vice versa. The key to “pseudo-elastic” behaviour of these materials and
the “shape memory effect” [8,26] is held by the first-order martensitic phase transformations in
these materials. Indeed, drastic changes in their properties are originated from their microstructure or phase combinations. Via the phase transformation, the microstructure of the material can
be switched among various combinations between austenite, martensite variants, or their mixtures.
Austenitic phase is a more symmetric phase of the crystallic lattice, prevailing at high temperature,
while martensite is a less symmetric, low-temperature phase [1,8,26]. Between these two critical
situations, austenite and martensite might co-exist providing a typical example of phase combinations. Upon external loading, one phase combination can be switched to another. If the material is constrained at the boundary, a specific (“self-accommodating”) combination of different
phases will be established such that the bulk energy in the constrained domain will be minimized
[1,18,19].
The mathematical framework for modeling phase combinations in shape memory materials is
based on the solution of the variational problem with respect to a frame-indifferent non-convex
free energy function φ(Y, θ) ([5,6,11,19,18] and references therein) :
W (Y ) =

Z

φ(Y, θ)dV,

(1)

Ω

where Ω is the reference configuration associated with the considered material, Y is the deformation tensor, and θ is the temperature of the material. Hence, by minimizing W (Y ) from (1), we
minimize the bulk energy of the considered structure, as a functional of the deformation tensor, at
temperature θ. This procedure, performed particularly often at the mesoscale [7,10,11,19,18], has
several known difficulties. It is known that in the general case the variational problem given by
Eq.(1) may have infinitely many minimizers ([7,19,10] and references therein). On the other hand,
a more precise definition of the free energy that allow us to account for interfacial energy effects
by introducing gradients of the order parameters is a highly non-trivial task ([1,18] and references
therein) connected with additional difficulties. The problem can also be regularized by assigning
simplified, e.g. affine, boundary conditions. For example, for the simulation of simple laminated
microstructure, the affine boundary conditions are constructed by assuming that the deformation
gradient of the material on the domain boundary is a linear combination of its equilibrium deformation gradients:
Y (x) = (λF0 + (1 − λ)RF1 )x, x ∈ ∂Ω,

(2)

3

Phase Combinations by Hybrid Optimizations

where R is a rotation matrix satisfying a twinning equation (see [7,10,11,19]) and λ is the thickness
of layers in the laminated microstructure. F0 and F1 are the deformation gradients that minimize
the local energy function (for square to rectangular transformation) when the deformation gradient
∇Y takes any of the following following values:
√
 1 + em / 2


F0 = 

0

0



√
 1 − em / 2


√  , F1 = 
1 − em / 2


0

0



√ ,
1 + em / 2


(3)

where em is a local minimum of the free energy function φ(Y, θ). In the general case, however, the
above affine boundary conditions may not be appropriate.
In solving problem (1), we have to face also numerical challenges connected with several local
minima, resulted from phase mixtures under low and moderate temperature regimes, and nonconvexity of the problem ([4,10,19,30] and references therein). Minimization procedures based
on conventional local search methods with randomly chosen initial guesses [9] may not lead to
a satisfactory result in those cases where the solution space is discretized with a large number of
node points.
In the present paper, we construct a mathematical model for the simulation of phase combinations
in SMA materials on the basis of the Landau theory. We associate the phase combination with
the global minimizer of the bulk energy of the SMA structure with the prescribed mechanical
boundary conditions. We develop a hybrid optimization strategy consisting of two main steps.
Firstly, we apply the Genetic Algorithm (GA) and its global exploration capacity to obtain an
initial estimation of the global minimizer. Then, we apply the quasi-Newton method to refine
such an estimation locally. Finally, the developed procedure is demonstrated by several numerical
examples simulating phase combinations in SMA materials.

2 Landau Free Energy Function and Variational Formulation

For the modeling of phase combinations in SMAs structures, the first task is to characterize different phases, which is different from one material to another. Here our mathematical model will
be based on the square to rectangular transformations. In Fig.1 (a) we give a schematic representation of this case where the square lattice is the austenite, while two rectangles are the martensite
variants. The square to rectangular transformation could be regarded as a 2D analog of the cubic to tetragonal or tetragonal to orthorhombic transformations observed in general 3D cases in
Nb3 Sn, InT l, F eP d alloys and some copper based SMAs ([17,18] and reference therein). The
analysis of this 2D transformation is a first step in understanding more complex cubic to tetrago-

4

L. X. Wang, R. V. N. Melnik

nal and tetragonal to orthorhombic transformations. Most numerical studies of the dynamics of the
phase transitions up to date have been concentrated on the analysis of the formation and growth
of the microstructure ([1,17,18] and references therein). Such studies have been focused on the
mesoscale under either periodic or affine boundary conditions.
In what follows, we base our consideration on the Landau theory of phase transformation. According to this theory, the basis of any nonlinear continuum thermodynamical model for phase
transformation is a non-convex free energy function [23,17,18]. The local minima of the free energy function with respect to the strain tensor (or deformation gradients) correspond to the stable
and mesostable state at a given temperature, while the microstructure in the domain of interest can
be described by the minimizer of the bulk energy in the domain. One of the simplest realization
of this idea in the context of SMAs is based on the Helmholtz free energy Ψ ([13,23,18,22]and
references therein):
Ψ(θ, ε) = ψ0 (θ) + ψ1 (θ)ψ2 (ε) + ψ3 (ε),

(4)

where ψ0 (θ) models thermal field contributions, ψ1 (θ)ψ2 (ε) models shape memory contributions
and ψ3 (ε) models mechanical field contributions, ε = ∂u/∂x (u is the displacement) is the strain
which is chosen as the only order parameter in the 1D case. The thermal field contributions ψ0 can
often be modelled as follows [13,28] :
ψ0 = −cv θ ln θ,

(5)

where cv is the specific heat constant.
The numerical analysis of a system of conservation laws based on the above representation of
the free energy function (and a more general one) has been recently reported in detail in [22]. A
conservative numerical scheme was constructed for the solution of the problem. It was noted that a
standard energy inequality technique, applied to the convergence analysis of the scheme, can lead
to quite restrictive assumptions. In [22] it was shown how such assumptions can be removed.
This work focuses on the practical development of an algorithm suitable for the simulation of
phase combinations. In this context we note that a free elastic energy functional (denoted further
by F ), similar to the one discussed above, was established earlier to characterize the austenite at
high temperature and the martensite variants at low temperature in SMA patches, specifically for
the square to rectangular transformation where the Landau free energy function Fl was modified
[17,18,1,28]. Recall that for the square to rectangular transformation we have to deal only with two
martensite variants and only one order parameter [13,17,18] in order to characterize the martensite
variants and austenite in a 2D domain. Following previous works on the subject ([13,17,18,28] and

Phase Combinations by Hybrid Optimizations

5

references there in), we have:
F = Fs + Fg ,
a1 2 a3 2
e + e3 + Fl
2 1
2
A2 2 a4 4 a6 6
e − e2 + e2 ,
Fl =
2 2
4
6
Fs =

(6)

3
3
d2 X
d3 X
2
Fg =
(∇ei ) +
(∇2 ei )2 .
2 i=1
4 i=1

where ∇ is the gradient operator, A2 , ai i = 1, . . . , 6, d2 , and d3 are material-specific coefficients,
and e1 , e2 , e3 are dilatational, deviatoric, and shear components of the strains, respectively, defined
as follows:
√
e1 = (η11 + η22 ) / 2,
√
e2 = (η11 − η22 ) / 2,

(7)

e3 = (η12 + η21 ) /2.
The Cauchy-Lagrangian strain tensor η is given by its components
ηij (x, t) =

!

∂ui (x, t) ∂uj (x, t)
+
/2,
∂xj
∂xi

(8)

where ui is the displacement in the ith direction in the Cartesian system of coordinates, x=
(x1 , x2 , x3 ) are the coordinates of a material point in the domain of interest. It is known that in
this formulation the deviatoric strain e2 can be chosen as the order parameter. More precisely, e2
and e3 are two-component order parameter strains that have been discussed before in [1,2,3].
In the above formulation the Ginzburg term Fg is the term proportional to square of strain gradients. This term is often included to account for the presence of domain walls. It is essential in simulating phase growth and several other phenomena ([17,18]and references there in), but it can be
ignored if we are interested only in the macroscopic phase combinations of the SMAs patch under
mechanical loadings (indeed, the simulation scale in this case is too coarse to capture mesoscale
structures). Furthermore, the energy contribution of this term is typically small compared to other
terms.

6

L. X. Wang, R. V. N. Melnik

In order to be able to model the entire range of different phase combinations that correspond to
different temperatures, the material parameter A2 is assumed to be temperature dependent A2 =
a2 (θ − θ0 ), where θ0 is a critical temperature, responsible for the appearance of an addtional
minimum and corresponding to the austenitic phase when temperature increases. In Fig.1 (b) we
present the plots of the Landau free energy function, defined in this way, for the entire range of
temperatures of interest (the material Au23 Cu30 Zn47 ). We observe that the function has two local
minima at low temperatures (210o ), which correspond to two (rectangular in the interpretation of
Fig.1 (a)) martensite variants, while only one minimum at the center corresponds to the (square)
austenite phase when the temperature is high (270o ). When the temperature is in between two
critical values (e.g., at around 245o in the figure), we observe that there are three local minima
demonstrating co-existence of metastable and stable phases.
Now, if we 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 transformation will take the following
form:
Ψ(θ, ε) = −cv θ ln θ +

a1 2 a3 2
e + e3 + FL ,
2 1
2

FL =

a2
a4
a6
(θ − θ0 ) e22 − e42 + e62 .
2
4
6

(9)

By substituting the free energy function into the variational problem given by Eq.(1), the phase
combination problem can be written as the following variational problem. Given temperature θ,
find the displacements ux and uy (in the x and y direction) that minimize the bulk energy:
W (ux , uy ) =

a3
a2
a4
a6
a1
−cv θ ln θ + e21 + e23 + (θ − θ0 ) e22 − e42 + e62 dV.
2
2
2
4
6


Z 

Ω

(10)

Under a given temperature, the contribution of the thermal field will not change the profile of
the local free energy function, but rather only shift it upwards or downwards. Hence, if the applied
external force in 2D is given by its components fx , fy , the final problem to solve can be formulated
as follows:
W (ux , uy ) =

a1 2 a3 2 a2
a4
a6
e1 + e3 + (θ − θ0 ) e22 − e42 + e62 − fx ux − fy uy dV → min .
2
2
2
4
6

Z 

(11)

Ω

The above model is reduced to the well-known Falk model in the 1D case [13] that can be written
with respect the only order parameter ǫ = ∂u/∂x:
W (u) =

Z 

Ω

a2
a4
a6
(θ − θ0 ) ǫ2 − ǫ4 + ǫ6 − f u dV,
2
4
6


(12)

Phase Combinations by Hybrid Optimizations

7

As discussed in Section 1, we supplement the model (11) by appropriate boundary conditions. For
all examples discussed in Section 4, these are clamped boundary conditions:
ux = uy = 0, at x = xl , xr , or y = yt , yb ,

(13)

where xl and xr are the left and right boundaries along the x direction, yt and yb are the top
and bottom boundaries along the y direction, as sketched in Fig.3. External forces vary and are
specified in Section 4.

3 Numerical Implementation Based on Hybrid Optimization

The above variational problem is non-convex and its solution can only be obtained by numerical
methods. The procedure developed in this section consists of two main steps: firstly, the variational
problem is converted into a nonlinear minimization problem by spatial discretization, and then the
solution of the resulting problem is sought by a hybrid optimization strategy.

3.1 Spatial Discretization Procedure

For the spatial discretization, we employ the Chebyshev pseudospectral approximation (e.g., [4])
on a set of 2D Chebyshev points (xi , yj ) in the 2D domain of interest Ω = [−1, 1] × [−1, 1]:
xi = cos(

πi
πj
), yj = cos( ),
N
N

i, j = 0, 1, . . . , N,

(14)

where N + 1 is the total number of nodes in one direction. Other structures of interest can be
mapped onto the domain Ω by a linear transformation. Using the constructed grid, the displacements in the patch can be approximated as follows:
f (x, y) =

N X
N
X

fi,j ξi (x)ξj (y),

(15)

i=0 j=0

where f (x, y) represents either function ux or uy , fi,j is the function value at (xi , yj ). Functions
ξi (x) and ξj (y) are the ith and j th Lagrange interpolating polynomials along the x and y directions
respectively.

8

L. X. Wang, R. V. N. Melnik

Having obtained (15), the derivatives of function f , ∂f (x, y)/∂x and ∂f (x, y)/∂y, can be obtained
by calculating ∂ξi (x)/∂x and ∂ξj (y)/∂y. Following the standard technique found, e.g., in [27,4],
all the differentiation operators in Eq.(11) (or Eq.12) can be written in the matrix forms:
F x = DxF ,

F y = Dy F ,

(16)

where F x and F are vectors collecting all values of the derivative ∂f /∂x and the function f
at (xi , yj ), respectively, and similarly for Fy . The differentiation matrices D x and D y can be
calculated using the approximation given by Eq.(15). For instance, the differentiation matrix Dx
for the Falk model takes the following form:

Dij =


2N 2 + 1






6






2N 2 + 1


−


6


xj



−



2(1 − x2j )






ci (−1)i+j




cj (xi − xj )

i = j = 0,
i = j = N,
(17)
i = j = 1, 2, . . . , N − 1,
i 6= j,

i, j = 1, 2, . . . , N − 1,

where ci = 2 for i = 0, N and ci = 1 otherwise. Such matrices have dimensionality (N + 1) ×
(N + 1).
The bulk energy is given by an integral operator with the local free energy function as its integrand. We use the same set of points as chosen for the derivative approximation for constructing a
numerical integration formula for integral. In particular, we use Chebyshev collocation nodes and
the resulting quadrature formula is constructed by using the Chebyshev-Lobatto rule [4,27]. For
example, the formula for integration in the x-direction, generically represented by
Z1
0

f (x)dx ≈

N
X

wi f (xi ),

(18)

i=0

is exact for any polynomials with an order less than 2N − 1. In (18) weight coefficients wi are
defined in the standard manner [4,15].
By substituting the approximations of all the differentiation operators and integral operators into
the variational problem, we convert the original problem into the following minimization problem:
find u1i,j , u2i,j to minimize: W (u1i,j , u2i,j ),

(19)

Phase Combinations by Hybrid Optimizations

9

where u1i,j and u2i,j stand for the values of ux and uy at node (xi , yj ), respectively. Note that
W (u1i,j , u2i,j ) is just the discretized bulk energy and is a nonlinear algebraic function of u1i,j and
u2i,j . The resulting problem has 2 × (N + 1) × (N + 1) variables in total. Given the prescribed
boundary conditions, we have 2 × (N − 1) × (N − 1) variables in total.
The problem to be solved is non-convex minimization problem with strong nonlinearity. It is not
easily amenable to conventional gradient-based minimization methodologies due to multiple local
minima. On the other hand, the genetic algorithms (GA) can be helpful in locating an approximation to the global minimum, giving an initial approximation to gradient-based procedures. In
what follows, we combine these two ideas by employing a hybrid optimization strategy that takes
advantage of the global exploring capability of the GA and the local refinement accuracy of the
quasi-Newton method for non-convex problems.

3.2 Genetic Algorithm Locates Initial Approximations

The GA is a well established methodology for global optimization ([16,20,14] and references
therein). We highlight here only its main features in the context of our problem. The GA maintains
a population of individuals (chromosomes), say P (n), for generation n and each chromosome consists of a set of genes, where each gene stands for a parameter to be estimated. One chromosome
represents one potential solution to the minimization problem. Each chromosome is evaluated to
give some measure of its fitness according to the bulk energy defined in the previous section.
Some chromosomes undergo stochastic transformations by means of genetic operations to form
new chromosomes. Recall that there are two transformations in the GA: crossover, which creates
new chromosome by combining parts from two chromosomes, and mutation, which creates a new
chromosome by making changes in a single chromosome. New chromosomes, called offsprings
S(n), are then evaluated. A new population is formed by selecting fitter chromosomes from the
parent population and the offspring population. After some generations, the algorithm converges to
the fittest chromosome, which represents an estimated optimal solution to the problem [16,20,14].
Generation of initial chromosomes: In most of the GAs, the chromosomes are created by randomly
choosing the genes in a given range. For the current problem, this is not an effective way. Indeed,
if the strain in solid structures under consideration is assumed to be not very large, while displacements may vary slowly and be represented by smooth functions, there is no good reason to include
high frequency oscillations in displacements. At the same time, if the displacement values are randomly chosen on our discrete grid, high frequency oscillations may well be pronounced. Hence,
we smooth the randomly chosen chromosomes by the following filter operation:
U s = I d I sU r ,

(20)

10

L. X. Wang, R. V. N. Melnik

where U r is the vector collecting randomly chosen displacement values at all discretization nodes,
while U s is the resultant smoother profile of the displacements after the filter operation. The
matrix I s is an interpolation matrix which maps the given data at the discretization nodes onto
a sparser grid, using the least square approximation, while the matrix I d is another interpolation
matrix which maps the data on the sparse grid back to the given discretization nodes. These two
matrices can be constructed by using the Chebyshev collocation method, as discussed above. The
product of these two interpolation matrices could be regarded as a filter to remove high frequency
components.
Crossover: Crossover operator in the GA is for producing new children chromosomes from chosen
parent chromosomes. For a given pair of parent chromosomes x1 and x2 , the offspring is obtained
by the so called Intermediate recombination [20,15]:
o1 = x1 α + x2 (1 − α),

(21)

where α is a vector with the same size as x, and all its entries are randomly chosen independently.
This operation is capable of producing variables slightly larger than the hypercube in the solution
space defined by the parents, but confined by the parameter α. A typical range for the parameters
in α is [−0.25, 1.25], which is used in the current paper.
Mutation: The mutation operation is randomly applied with a low probability, typically in the
range 0.001 to 0.01, and modifies genes in the chosen chromosome. Its role in the GA is to make
sure that the probability of searching any potential solution (any points in the solution space) is
nonzero, and is often regarded as a measure to recover good genetic material that might be lost
through the operation of selection and crossover [20,16,14,15].
In practice, the mutation operation is carried out by replacing a randomly chosen gene by a randomly chosen new value in the given range. The mutated chromosome itself is randomly chosen
from the current population with a given probability. In the current problem, the mutation operation is applied to one chromosome in one generation, which means the probability for each
chromosome is one divided by the number of chromosomes.
Generation-alteration method: Generation alteration method determines how to evolve the current
generation to the next, which means how to select pairs of parents for producing children by
crossover and mutation operator, and how to select parents in the currents population that survive
in the next generation.
Taking into account their fitness, it is obvious that each chromosome should have a probability
determined by the associated function value. Here a simple linear map method is used as follows:
rank all the chromosomes in the current generation in order of decreasing function values, then the

Phase Combinations by Hybrid Optimizations

11

probability of being chosen for a specific chromosome i is calculated as:
ρi =

i
,
2(N + 1)N

(22)

where i is the position of the chromosome in the rank, while N is the number of chromosomes
in the rank. Then the parents for survival and crossover are chosen using the above calculated
probability. While all the chromosomes have nonzero probability to be chosen, fitter chromosomes
have a better chance.
The results of application of this procedure to the simulation of SMA phase combinations are
reported in the next section.

3.3 Quasi-Newton Method Refines the Solution

The output of the GA is used as the initial guess for local search methods. In what follows we apply
the quasi-Newton method to the minimization of the bulk energy given in Eq.11. Let’s denote any
potential minimizer by x. The bulk energy W (x) is minimized by x∗ if it satisfies ∇W (x∗ ) = 0.
To achieve this, we can organize the following iterative process. Let at the general k th step of the
iteration, the potential solution is xk . Then, the task is to estimate the next xk+1 = xk + dk such
that ∇W (xk+1 ) = 0 which means that
∇W (xk + dk ) = q(xk + dk ) = q(xk ) + ∇q(xk )dk = 0,

(23)

where q(xk ) = ∇W (xk ) is the gradient of the bulk energy at xk , ∇q(xk ) = ∇2 W (xk ) is the
Hessian matrix, while dk is the search direction. To make ∇W (xk+1) = 0, the search direction
should satisfy:
dk = −(∇q(xk ))−1 q(xk ).

(24)

This standard Newton-type procedure is computationally expensive and has other well-known
drawbacks [4,30] that preclude us from using it in the context of our problem. Instead, we apply
the quasi-Newton procedure by constructing a matrix B k (at the k th iteration), an approximation
to the Hessian matrix that satisfies the following condition:
q(xk ) + B k dk = 0.

(25)

12

L. X. Wang, R. V. N. Melnik

Provided with the initial guess from the GA, the local search method using the quasi-Newton
method is organized in a standard manner: at the general k th iteration (xk is the current estimated
solution):
(a) Find a descent direction dk using B k by the following formula:
k
dk = −B −1
k ∇W (x ),

(26)

(b) Compute the acceleration parameter αk by line search, find αk such that W (xk + αk dk ) is
minimized.
(c) Update the potential solution:
xk+1 = xk + αk dk .

(27)

Finally, the update of the approximation of the Hessian matrix B k+1 is organized by using the
BFGS update (Broyden-Fletcher-Goldfarb-Shanno update, e.g., [4,30,21]):
B k+1 = B k −

(Bkdk)(Bkdk)T
y k y Tk
,
+
y Tk dk
dTk B k dk

(28)

where y k = ∇W (xk+1 )−∇W (xk ). To initiate the iteration process, B0 for the first step is taken as
the identity matrix that corresponds to the deepest descent methodology. For the current problem,
we did not apply the line search for the acceleration parameter αk . Instead, a small value, regarded
as a relaxation factor, was assigned to αk .

4 Numerical Examples: Phase Combinations in SMA with Hybrid Optimization Procedure

By combining the GA and the quasi-Newton method, the bulk energy given in Eq.(11) and Eq.(12)
can be minimized with respect to displacements. To demonstrate the capability of this hybrid
optimization method, three numerical experiments are reported here. For all three experiments, the
GA first evolves a given number of generations, and the fittest chromosome in the last generation
is selected. This fittest chromosome is then used as the initial guess for the quasi-Newton method,
and is refined iteratively. The termination criterion for the quasi-Newton method is based on the
assumption that the norm of difference between two consecutive potential solutions is smaller than
the predefined value δ (chosen in experiments as δ = 1 × 10−6 ):
ksk+1k = kxk+1 − xk k ≤ δ.

(29)

Phase Combinations by Hybrid Optimizations

13

All simulations reported here have been carried out for Au23 Cu30 Zn47 . For this specific material,
its physical parameters are available for the 1D case: [13,25]:
a2 = 480 g/ms2cmK,
θ0 = 208K,

a4 = 6 × 106 g/ms2 cmK,

ρ = 11.1g/cm3 ,

cv = 3.1274g/ms2cmK,

a6 = 4.5 × 108 g/ms2 cmK,
k = 1.9 × 10−2 cmg/ms3 K.

Since 2D experimental values are not available to us, we take all the parameters in the Landau
free energy function the same as above and complete parameterization of the model by assuming
a1 = 2a2 and a3 = a2 as suggested in [17,12]. Earlier, we confirmed numerically that the essential features of the 2D problem can be captured with this parameterization, at least in the case
of square to rectangular transformations considered here [28,29]. Recent studies presented in [3]
provided encouraging results also for more general cubic to tetragonal transformations. However,
experimental results on multi-dimensional SMA samples are still lacking. The next step and a natural development of the study presented here would be accounting systematically for the dynamics
of the thermal field. Numerical experiments pertinent to this generalization would be much more
involved.
The first experiment is the simulation of phase combinations in a 1D SMA wire with the bulk
energy given by Eq.(12). The physical interpretation of the problem is sketched in Fig.2. The
length of the wire is 1cm , and the applied mechanical force is f = 500g/ (ms2 cm2 ) which is
evenly distributed along the whole length. The initial condition for this case, u = (rand − 0.5)/5,
provides a random distribution of displacement in the interval [−0.1, 0.1]. 15 nodes are used for
the spatial discretization. For the filter operation in the GA, 7 nodes are used for the sparser grid to
remove the higher frequency components. The GA evolves 800 generations with 60 chromosomes
in each generation, each chromosome consists of displacement values within the range from −0.1
to 0.1 on internal nodes.
To analyze the performance of the hybrid optimization method for phase combination simulation,
the bulk energy given by Eq.(12) has been monitored in the GA, and plotted in Fig.4 (left). The
quasi-Newton iteration process has also been monitored. The actual update step size in the quasiNewton method is plotted in Fig.4 (right). The trends in the two curves demonstrate how the
optimization process evolves. As expected, the GA starts to converge to the global minimizer
gradually, but has difficulty to locate it precisely. The quasi-Newton method starts with the output
of the GA as its initial guess, and refined the global minimizer effectively.
The final estimation for the phase combination in this case is given in the right column in Fig.5:
the order parameter ǫ (upper plot) and the displacement (lower plot). We observe that the entire
domain is divided into two parts, one with the strain values ǫ ≃ 0.11 and the other with ǫ ≃ −0.11.
This phase combination agrees well with the previously reported results (e.g., [13,19,25]). For

14

L. X. Wang, R. V. N. Melnik

comparison purpose, we have also provided an approximation obtained with the GA at the initial
stage of the procedure (see the left column of Fig.5). Although the quality of the final solution has
been refined with the quasi-Newton procedure, we note that the two parts in the distribution of ǫ
are captured by the GA. Quantitative values are easily identified as being between 0.1 to 0.15 for
one of the domains and between −0.15 to −0.1 for the other.
The second experiment aims at simulating the phase combination in a 2D SMA patch sketched
in Fig.3. The size of the patch is [−1, 1] × [−1, 1]cm2 and the applied mechanical forces fx =
fy = 3000g/ (ms2 cm2 ) are distributed evenly through the entire patch. In this case, we use
ux = (rand − 0.5)/5 and uy = (rand − 0.5)/5 as the initial conditions. In this experiment,
we use different node numbers for the GA and the quasi-Newton method. In the GA, there are 9
nodes in each direction used for the discretization along the x and y axes. This leads to the total
number of optimization parameters being 128. The filter operation is carried out in the same way as
explained for the 1D experiment, but for both x and y directions, and the number of nodes in each
direction on the sparse grid is 6. The GA is firstly run 1500 generations with 120 chromosomes
in each generation. The output of the GA is then interpolated onto a denser grid with 15 nodes
in each direction, and is used as the initial guess for the quasi-Newton method. The interpolation
methodology is based on the pseudospectral method with the Chebyshev collocation points, as
discussed in Section 3.
The final phase combination estimated in this case are presented in Fig.6. The distribution of the
order parameter e2 is given in the upper right subplot. For comparison purpose, the estimated
distribution of e2 from the GA is presented in the upper left subplot. The final displacement distributions along the x and y directions are presented by the two lower subplots. Once again, the entire
domain can be divided by two parts with values e2 ≃ 0.11 e2 ≃ −0.11, respectively, referred to as
martensite plus and martensite minus [24,29].
In the last experiment, we consider the same patch [0, 1] × [0, 1]cm2 with modified distributed
forces fy = 0, fx = 2000g/ (ms2 cm2 ) (with the same initial conditions as in the previous experiment). All other computational parameters are taken the same as those in the previous experiment.
The numerical results are presented in Fig.7 in a way similar to already reported. In this case, the
entire structure is divided into martensite plus and minus in a different way, the interface is the
central vertical line due to the horizontal symmetric loading. The simulated order parameter is still
close to either 0.11 or −0.11 and the estimated e2 distribution obtained with the GA can capture
the essence of its final profile.
In the analysis that follows we explain why the order parameter takes values close to either 0.11
or −0.11. From the Landau free energy function, it is easy to estimate the order parameter value
which minimizes the local free energy by setting:
δFl
= 0.
δe2

(30)

Phase Combinations by Hybrid Optimizations

15

For the 1D case, the order parameter should be replaced by ǫ. A simple operation gives the following equation:
a2 (∆θ)e2 − a4 e32 + a6 e52 = 0,

(31)

which means:
e2 =

a4 ±

q

a24 − 4a6 a2 ∆θ
.
2a6

(32)

Note that e2 = 0 is also a solution to (31) and can be associated with the austenite phase which
is unstable at the current temperature. Therefore, martensitic phases in this case are of greater
interest. The temperature difference from the transformation temperature here is given as ∆θ =
θ − θ0 = 2o , and the local minima can be estimated as: e2 = ±0.1146.
All three numerical experiments have demonstrated that the distributions of the order parameter
are represented by a combination of two different values, e2 = ±0.11, which agrees well with the
prediction of the above analysis.

5 Conclusions

In the present paper, phase combinations in the 1D and 2D SMA structures have been analyzed
in the case of square to rectangular transformations. The phase combinations have been obtained
by minimizing the bulk energy in the considered SMAs structures, subject to given temperature
distribution, mechanical loadings, and boundary conditions. By combining the global and local
search techniques, the developed hybrid optimization method provides a promising strategy for
the solution of the associated non-convex problem.

References

[1] R. A HLUWALIA , T. L OOKMAN , AND A. S AXENA, Elastic deformation of polycrystals, Phys. Rev. Lett.
91(2003), pp. 055501.
[2] R. A HLUWALIA , T. L OOKMAN , A. S AXENA , AND S. R. S HENOY, Pattern formation in ferroelastic
transitions, Phase Transitions, 77(5-7) (2004), pp. 457–467.

16

L. X. Wang, R. V. N. Melnik

[3] R. A HLUWALIA , T. L OOKMAN , AND A. S AXENA, Dynamic strain loading of cubic to tetragonal
martensites, Acta Materilia, 54(2006), pp. 2109–2120.
[4] Q. A LFIO , S. R ICCARDO ,

AND

S. FAUSTO, Numerical Mathematics, Springer-Verlag, 2000.

[5] J. M. BALL , R. D. JAMES, Fine phase mixtures as minimizers of energy, Archive. Rat. Mech. Anal.
100(1) (1987), pp. 13–52.
[6] J. M. BALL , AND R. D. JAMES, Proposed experimental tests of the theory of fine microstructure and
the two-well problem, Philos. Trans. R. Soc. Lond., Ser. A, 338 (1992), pp. 389–450.
[7] S. BARTELS , T. ROUBICEK, Linear-programming approach to non-convex variational problems,
Numerische Mathematik, 99(2)(2004), pp. 251–287.
[8] V.B IRMAN, Review of mechanics of shape memory alloys structures, Appl.Mech.Rev., 50(1997),
pp. 629–645.
[9] K. B HATTACHARYA , B. L I , AND M. L USKIN, The simply laminated microstructure in martensitic
crystals that undergo a cubic to orthorhombic phase transformation, Arch. Rat. Mech. Anal. 149(1999),
pp. 123–154.
[10] C. C ARSTENSEN, Ten remarks on nonconvex minimisation for phase transition simulations. Comput.
Methods Appl. Mech. engrg., 194(2005), pp. 169–193.
[11] S. C OLLINS , M. L USKIN , AND J. R IORDAN, Computational results for a two-dimensional model
of crystalline microstructure, In.D. Kinderlehrer, R. James, M. Luskin, and J. L. Ericksen (Eds),
Microstructure and Phase Transition, The IMA Volumes in Mathematics and its Applications, Vol.54,
Springer, (1993), pp. 51–56.
[12] S. H. C URNOE , AND A. E. JACOBS, Time evolution of tetragonal-orthorhombic ferroelastic, Phys.
Rev. B, 64(6)(2001), pp. 064101.
[13] F. FALK, Model free energy, mechanics, and thermomechanics of shape memory alloys. Acta
Metallurgica, 28(1980), pp. 1773–1780.
[14] S. F ORREST, Genetic algorithms: principles of atural selection applied to computation, Science,
261(1993), pp. 872–878.
[15] G. G AUTSCHI, Orthogonal polynomials and quadrature,Electronic Transaction on Numerical
Algorithm, 9(1999), pp. 65–76.
[16] D. E. G OLDBERG, Genetic Algorithm in Search, Optimization and Machine Learning, AddisonWesley, 1989.
[17] A. E. JACOBS, Landau theory of structures in tetragonal to orthorhombic ferroelastics, Phys. Rev. B,
61(10)(2000), pp. 6587–6595.
[18] T. L OOKMAN , S. S HENOY, D. R ASMUSSEN , A. S AXENA ,
and strain compatibility, Phys. Rev. B, 67(2003), pp. 024114.

AND

A. B ISHOP, Ferroelastic dynamics

Phase Combinations by Hybrid Optimizations

17

[19] M. L USKIN, On the computational of crystalline microstructure, Acta Numerica, 5(1996), pp. 191–
256.
[20] M. M ITCHELL, An Introduction to Genetic Algorithms, MIT Press, 1996.
[21] J. M. M ARTINEZ, Practical quasi-Newton methods for solving nonlinear systems, J. Comp. App.
Math., 124(2000), pp. 97–121.
[22] P. M ATUS , R. V. N. M ELNIK , L. X. WANG , AND I. RYBAK, Applications of fully conservative
schemes in nonlinear thermoelasticity: modelling shape memory materials, Mathematics and Computers
in Simulation, 65(4-5)(2004), pp. 489–509.
[23] R. V. N. M ELNIK , A. J. ROBERTS , AND K. A. T HOMAS, Computing dynamics of copper-based SMA
via centre manifold reduction of 3D models, Computational Materials Science, 18(3-4)(2000), pp. 255–
268.
[24] R. V. N. M ELNIK , A. J. ROBERTS , AND K. A. T HOMAS, Coupled thermomechanical dynamics
of phase transitions in shape memory alloys and related hysteresis phenomena, Mechanics Research
Communications, 28(6)(2001), pp. 637–651.
[25] M. N IEZGODKA , AND J. S PREKELS, Convergent numerical approximations of the thermomechanical
phase transitions in shape memory alloys, Numerische Mathematik, 58(1991), pp. 759–778.
[26] K. OTSUKA , C. M. WAYMAN(editors), Shape Memory Materials, Cambridge University Press, 1998.
[27] L. N. T REFETHEN, Spectral methods in Matlab, SIAM, 2000.
[28] L. X. WANG , AND R. V. N. M ELNIK, Dynamics of shape memory alloys patches, Materials Science
and Engineering A, 378(1-2) (2004), pp. 470–474.
[29] L. X. WANG , AND R. V. N. M ELNIK, Thermomechanical Waves in SMA Patches under Small
Mechanical Loadings, Lecture Notes in Computer Science, 3039(2004), pp. 645–653.
[30] Y. X. Y UAN ,
Chinese).

AND

W. Y. S UN, Theory and methods of optimization, Chinese Science Press, 1997 (in

18

L. X. Wang, R. V. N. Melnik

Variant 1

Austenite
Martensite

Variant 2
150

270o

Free energy Fl

100

50
o

245

230o

0

−50
210o
−100
−0.2

−0.15

−0.1

−0.05
0
Strain e2

0.05

0.1

0.15

Fig. 1. (a) Sketch of the square to rectangular transformation, (b) The temperature dependency of the free
energy function for the transformation.

19

Phase Combinations by Hybrid Optimizations

f

L
Fig. 2. Sketch of an one dimensional SMA wire under mechanical loading

fy
fx

Ly

Lx
Fig. 3. Sketch of a two dimensional SMA patch under mechanical loading

1.7

0
Best = 36.7267

−1

log10(f(x))

log10(  sk  )

1.65

−2

−3

−4

1.6

−5

0

200

400

600

generation

800

1000

−6
0

10

20

30

40

50

60

70

80

Iteration Number

Fig. 4. Convergence of the hybrid optimization method. Left: reduction of the bulk energy in the GA. Right:
reduction of the step size in the quasi-Newton method

20

L. X. Wang, R. V. N. Melnik

0.1

0.1

0.05

0.05

0

0
ε

0.15

ε

0.15

−0.05

−0.05

−0.1

−0.1

−0.15

−0.15

−0.2
0

0.2

0.4

0.6

0.8

−0.2
0

1

0.2

0.4

x

0.6

0.8

1

0.6

0.8

1

x
0.06

0.05

0.05

0.04

0.04

u

u

0.06

0.03

0.03

0.02

0.02

0.01

0.01

0
0

0.2

0.4

0.6

x

0.8

1

0
0

0.2

0.4

x

Fig. 5. Mechanically induced phase combination in a SMA wire.

21

Phase Combinations by Hybrid Optimizations

0.2

0.2
0.1

0

e2

e2

0.1

−0.1

−0.1
−0.2
1

−0.2
1

y

−1

0

−0.5

0.5

1

0.5

0

0

−0.5

y

−0.5
−1 −1

x

0.1

x

0.1

2

0.05

u

1

1

0.5
0

u

0

0.05

0

0

1

1

0.5

0.5
1
0

0.5

1
0

0

−0.5

−0.5

y

−1

0.5
0

−0.5

−0.5

y
x

−1

x

Fig. 6. Phase combination in a SMA patch induced by mechanical loading along the x and y directions.

22

L. X. Wang, R. V. N. Melnik

0.2

0.2
0.1

e

e2

2

0.1
0

−0.1

−0.1
−0.2
0

−0.2
0

1
0.5

y

1 0

0.5

x

1 0

y

0.02

0.1
1

1
0.5

0.5

x

0.05

0

u

2

u

0

0
1

−0.02

0.8

1
0.6

0.8
0.6

0.4

1

0.8
0.6

0.2
0 0

x

0.5

0.4

0.4

0.2
y

1

y

0.2

x

0

Fig. 7. Phase combination in a SMA patch induced by mechanical loading along only the x direction

