Simulation of an Austenite-twinned-Martensite interface.

Developing numerical methods for predicting microstructure in materials is a large and important research area. Two examples of material microstructures are Austenite and Martensite. Austenite is a microscopic phase with simple crystallographic structure while Martensite is one with a more complex structure. One important task in materials science is the development of numerical procedures which accurately predict microstructures in Martensite. In this paper we present a method for simulating material microstructure close to an Austenite-Martensite interface. The method combines a quasi-Newton optimization algorithm and a nonconforming finite element scheme that successfully minimizes an approximation to the total stored energy near the interface of interest. Preliminary results suggest that the minimizers of this energy functional located by the developed numerical algorithm appear to display the desired characteristics.

Key words: Austenite-Martensite interface; finite element method; quasi-Newton method.

**********

1. Introduction

Simulation of microstructures is important for determining the behavior of complex materials [1, 2, 3, 4, 5, 6]. Applying macroscopic loads to some materials can reveal their mechanical properties. These properties are associated with the microscopic structures created due to the material deformation. Microscopic structures of interest are often found near internal surfaces which, in turn are found in systems of liquid-liquid, solid-liquid and solid-solid boundaries [7].

An Austenite-twinned-Martensite interface is an example of an internal surface. This interface conjoins two phases: Austenite and Martensite. Figure 1 shows an example of the twinned-Martensitic microstructures. These microscopic structures are the result of thermal or mechanical loading of the material. When a material is deformed, the structural organization of the atoms is rearranged. At the microscopic level, this rearrangement can create the crystallographic patterns shown in Figure 1, for example. These patterns represent 3-dimensional microscopic structures. There are many different ways to rearrange the atoms that make up a material and each variant corresponds to a particular arrangement. We focus on the two specific variants of Martensite distinguished by the alternating shaded areas in Figure 1.

One important mechanical property associated with an Austenite-twinned-Martensite interface is the Shape Memory Effect (SME) characteristic of Shape Memory Alloys (SMAs). SMAs are materials that, when deformed for an indefinite period of time, return to their original shape. The deformed state is a metastable state [8]. In metastability, the total stored energy is at a local minimum, thus requiring energy to induce a transformation. In SMAs, heat will trigger the SME which in turn, returus the SMA to its original shape. An Austenite-twinned-Martensite interface is observed in materials when they are in a metastable state. Therefore, to simulate Martensite, we minimize the total stored energy near an internal surface. In this paper, we present a numerical technique used to obtain a solution which represents the spatial structure of the twinned Martensite. The technique employs a [Q.sub.1] finite element discretization of a function representing total energy and a limited memory quasi-Newton method to minimize the resulting approximate total energy function. The gradient of the function is calculated by computing the partial derivatives using finite differences. The result is an apparently robust method for approximating what are usually difficult minimizers to locate.

[FIGURE 1 OMITTED]

1.1 Total Stored Energy

A myriad of research methods can be found which focus on the simulation of twinned Martensite, and related microstructures, e.g. Refs. [9, 10, 11, 12, 13, 14]). Previous numerical work includes minimization of the total stored energy using the conjugate gradient methods [15], the method of steepest descent [15], and a descent method [16]. Further references using and studying the use of finite element methods to simulate Martensitic microstructures can be found in [17, 18, 19, 20, 21, 22, 23, 24]. We are interested in simulating the microstructures in the metastable state, thus our work will focus on the static problem.

The functional representing the total stored energy is taken from work by Kohn-Muller in [10, 11]. The Martensite region is represented by the domain [OMEGA] = (0, L) X (0, K). Let x = (x, y) [member of] [OMEGA] [subset] [R.sup.2], u : [R.sup.2] [right arrow] R and u = u(x, y). The double-well energy function is given by

J(u) = [[integral].sub.[OMEGA]]([[partial derivative].sub.x]u)[.sup.2] + (([[partial derivative].sub.y]u)[.sup.2] dx + [epsilon][[integral].sub.[OMEGA]]|[[partial derivative].sub.yy]u| dx (1)

where u equals 0 at x = 0. The x = 0 boundary corresponds to the internal surface. The function u = 0 at x = 0, due to the constraint of elastic compatibility [10, 11]. The first integral is the elastic energy and the last integral is the surface energy. We seek a function u [member of] [W.sup.1.4] ([OMEGA]) which minimizes Eq. (1) where

[W.sup.1.4]([OMEGA]) = {v|v [member of] [L.sup.4] ([OMEGA]), [[partial derivative].sub.x] v [member of] [L.sup.4] ([OMEGA]), [[partial derivative].sub.y] v [member of] [L.sup.4] ([OMEGA])

and v = 0 at x = 0.}

1.2 Elastic Energy

The 2-D scalar model of elastic energy, in Eq. (1) is

[[integral].sub.[OMEGA]][([[partial derivative].sub.x]u)[.sup.2] + (([[partial derivative].sub.y]u)[.sup.2] - 1)[.sup.2]]dx. (2)

This term is minimized by a function u such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (3)

Each of the gradients in Eq. (3) corresponds to a stress-free state in one of the two distinct variants of Martensite [10,11]. Examples of functions satisfying Eq.(3) are plotted in Fig.2.

In Fig. 2 functions are piecewise linear in the y-direction and constant along the x-direction. Requiring that u equal zero at x = 0 creates more oscillations. This is similar to the behavior of Martensite near an internal surface. The function with Amplitude 1 is closest to zero, on average, than those with Amplitudes 2 and 3.

The number of oscillations is directly related to the number of discontinuities in [[partial derivative].sub.y]u. We note that the function with Amplitude 1 also has the largest number of discontinuities, occurring at the peaks of u, where [[partial derivative].sub.y]u changes from +1 to -1 (or vice-versa). Therefore, a desirable minimizer u has the characteristics of the function with Amplitude 1. The consequence of u satisfying Eq. (2) and u = 0 at x = 0 is the occurrence of more discontinuities in [[partial derivative].sub.y]u. At the minimizer u, however, the surface energy penalizes these discontinuities.

1.3 Surface Energy

We employ the surface energy presented by Kohn and Muller in Ref. [10,11],

[FIGURE 2 OMITTED]

[epsilon][[integral].sub.[OMEGA]]|[[partial derivative].sub.yy]u| dx, (4)

where [epsilon] is a constant.

The definition proposed by Kohn and Muller replaces the integral of Eq. (4) in the y-direction with the total variation of [[partial derivative].sub.y]u over (0, K) [10,11]. This change designates the role of the surface energy at the minimizer u as a "counter" because [1/2][[integral].sub.0.sup.K] |[[partial derivative].sub.yy]u|dy counts the number of discontinuities in [[partial derivative].sub.y]u, where [[partial derivative].sub.y] u [member of] {[+ or -]1}. Thus, at the minimizer u of Eq. (1), the surface energy Eq. (4) exhibits opposite behavior to Eq.(2) since we are minimizing Eq. (1). To illustrate this counting role of Eq. (4), we briefly review functions of bounded variation.

1.4 Functions of Bounded Variation

The definition of functions of bounded variation is taken from Refs. [25,26, 27]. Let

P = {[y.sub.0], [y.sub.1],..., [y.sub.t]}

be any partition of (0, K) and let

|P| = [max.[j]]|[y.sub.j+1] - [y.sub.j]|.

The partitioning P is a collection of points [y.sub.j], j = 0, 1,..., l, such that [y.sub.0] = 0 and [y.sub.l] = K. For each partitioning P, let

[S.sub.P]([bar.x]) = [l-1.summation over (i=0)]|[[partial derivative].sub.y]u([bar.x], [y.sub.j+1]) - [[partial derivative].sub.y]u([bar.x], [y.sub.j])|, (5)

for a fixed [bar.x] [member of] (0, L). The total variation of [[partial derivative].sub.y]u over (0, K) is

[sup.[P]][S.sub.P] ([bar.x]), (6)

where the supremum is taken over all countable partitions P of (0, K). If 0 [less than or equal to] [S.sub.P]([bar.x]) [less than or equal to] +[infinity], then 0 [less than or equal to] [epsilon][[integral].sub.[OMEGA]]|[[partial derivative].sub.yy]u | dx [less than or equal to] +[infinity]. If [epsilon][[integral].sub.[OMEGA]]|[[partial derivative].sub.yy]u | dx < +[infinity], then [[partial derivative].sub.y] u is a function of bounded variation.

Next, we describe the substitution of [[integral].sub.0.sup.K]|[[partial derivative].sub.yy]u | dy by Eq. (6). We assume that [[partial derivative].sub.y]u [member of] [C.sup.1]([OMEGA]). Using Eq. (5) and the Mean Value Theorem, we have for some [[zeta].sub.j] [member of] ([y.sub.j], [y.sub.j+1]), j = 1, 2,..., l that

[S.sub.p]([bar.x]) = [l-1.summation over (j=0)]|[[partial derivative].sub.y]u([bar.x], [y.sub.j+1]) - [[partial derivative].sub.y]u([bar.x], [y.sub.j])| (7)

= [l-1.summation over (j=0)]|[[partial derivative].sub.yy]u([bar.x],[[zeta].sub.j])|([y.sub.j+1] - [y.sub.j]), (8)

and we obtain by Theorem 2.9 in Ref. [25]:

[sub.[P]][S.sub.P]([bar.x]) = [lim.[|P|[right arrow]0]][S.sub.P]([bar.x]) = [lim.[|P|[right arrow]0]][l-1.summation over (j=0)]|[[partial derivative].sub.yy]u([bar.x],[[zeta].sub.j])|([y.sub.j+1] - [y.sub.j]) = [[integral].sub.0.sup.K]|[[partial derivative].sub.yy]u|d y (9)

The assumption that u [member of] [C.sup.1]([OMEGA]) is stronger than the assumtion that u [member of] [W.sup.1,4]([OMEGA]); however, the equivalence in Eq. (9) is possible since the changes in [[partial derivative].sub.y]u do not occur in subintervals of measure zero, [9, 10, 11, 28]. Figure 3 shows that [1/2][[integral].sub.0.sup.K]|[[partial derivative].sub.yy]u| dy counts the number of sign changes in [[partial derivative].sub.y]u. Let (0, K) = (0, 1).

In Fig. 3a, the piecewise linear function u([bar.x], y) is plotted over (0, 1). In Fig. 3b, the function [[partial derivative].sub.y]u([bar.x], y) is plotted over (0, 1). Clearly, the number of discontinuities in [[partial derivative].sub.y]u is 5 from Fig. 3b. For a discontinuity at the partition point [y.sub.j], we adopt the following convention:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

Now consider ([y.sub.j-1], [y.sub.j+1]) = ([y.sub.j-1], [y.sub.j])[union]([y.sub.j], [y.sub.j+1]) in Fig. 3b. Evaluating at the endpoint of the two subintervals gives

|[[partial derivative].sub.y]u([bar.x], [y.sub.j]) - [[partial derivative].sub.y]u([bar.x], [y.sub.j-1])| = |1 - (-1)| = 2 and (11)

|[[partial derivative].sub.y]u([bar.x], [y.sub.j+1]) - [[partial derivative].sub.y]u([bar.x], [y.sub.j])| = |1 - 1| = 0. (12)

[FIGURE 3 OMITTED]

A jump occurred in ([y.sub.j-1], [y.sub.j]) because [[partial derivative].sub.y]u is discontinuous at a point in this subinterval. Based on our conventional definition of [[partial derivative].sub.y]u([bar.x], [y.sub.j]) at a point [y.sub.j] in Eq. (10), we say a "jump" occurs at [y.sub.j] if [[partial derivative].sub.y]u is discontinuous at this point. This "jump" refers to the change in values of [[partial derivative].sub.y]u at [y.sub.j] from +1 to -1 or vice-versa.

This jump gives the nonzero value in Eq. (11). In the subinterval ([y.sub.j], [y.sub.j+1]), no jump occurs in [[partial derivative].sub.y]u, therefore the difference in Eq. (12) is zero. We compute Eq. (9):

[[integral].sub.0.sup.1]|[[partial derivative].sub.yy]u([bar.x], y)|dy = 2 X 5 = 10,

where the magnitude of the jump is 2 and the number of discontinuities is 5. Table 1 shows the number of discontinuities for the functions presented in Fig. 2.

Therefore, we have:

[[integral].sub.0.sup.1]|[[partial derivative].sub.yy]u([bar.x], y)| dy = 2 X (number of discontinuities of [[partial derivative].sub.y]u),

and the surface energy succeeds in serving as a "counter." The next section describes the numerical approximation to Eq. (1).

2. Description of Method

The total stored energy Eq. (1) consists of elastic plus surface energy. First, we describe the implementation of the elastic energy followed by a discussion of the surface energy implementation. We computed Eq. (2) using affine finite elements, [29]. We chose the finite element space:

[Q.sub.1] = span{1, x, y, xy}, (13)

because functions in this polynomial space lie in [W.sup.1,4] ([OMEGA]), which satisfy the criteria in Eq. (3) for a minimizer u of Eq. (1). The elements are rectangles with the function values at the vertices as degrees of freedom. The parent element is [^.Q] = (0, 1) X (0, 1) and for some k [member of] N, the reference element is

[Q.sub.k] = (i[h.sub.1], (i + 1)[h.sub.1]) X (j[h.sub.2], (j + 1) [h.sub.2]) [subset] [OMEGA] = [[union].p][Q.sub.p],

where [h.sub.1] is the mesh size along the x-axis and [h.sub.2] is the mesh size along the y-axis. The degrees of freedom of the elements are the function values at the vertices, and for [^.Q] we label them u([a.sub.1]), u([a.sub.2], u([a.sub.3]), and u([a.sub.4]). Figure 4 shows the parent element with the vertices also denoted by concentric circles.

2.1 Affine Finite Elements

Let F: [^.Q][right arrow][Q.sub.k] be given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (14)

We simplify, by letting F = F([^.x], [^.y]). We use the affine map Eq. (14) to evaluate the function u, to compute [nable]u and to integrate. The function u is spanned by the four basis functions defined in each element [Q.sub.k]. Using Eq. (14), we can determine all basis functions by the four basis functions in [^.Q]. Thus, the approximation to u is:

u(x, y) = [u.sub.1][[phi].sub.1]([^.x], [^.y]) + [u.sub.2][[phi].sub.2]([^.x], [^.y]) + [u.sub.3][[phi].sub.3]([^.x],[^.y]) + [u.sub.4][[phi].sub.4]([^.x], [^.y]), (15)

[FIGURE 4 OMITTED]

where x = [^.x](x, y), y = [^.y](x, y). The nth basis function in Q is

[[phi].sub.n]([^.x], [^.y]) = [[alpha].sub.0.sup.(n)] + [[alpha].sub.1.sup.(n)][^.x] + [[alpha].sub.2.sup.(n)][^.y] + [[alpha].sub.3.sup.(n)] [^.x][^.y],

and it satisfies the property that [[phi].sub.n]([a.sub.m]) = [[delta].sub.mn], for n, m = 1, 2, 3, 4, where [[delta].sub.mn] is the Kronecker delta function.

To compute the gradient of u and to integrate u, we need: [nable]F, det[nable]F and [nable][F.sup.-1]. The gradient of F is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

and, we have

[nable][[phi].sub.n](x, y) = [nable][F.sup.-1][nable][[phi].sub.n]([^.x], [^.y]). (17)

Using Eqs. (15), (16), and (17) we have

[nable]u(x, y)=[nable][F.sup.-1][4.summation over (n=1)][u.sub.n][nable][[phi].sub.n]([^.x], [^.y]).

To integrate the elastic energy over [Q.sub.k] we use Eqs. (15)-(17) and obtain

[[integral].sub.Q.sub.k]([[partial derivative].sub.x]u)[.sup.2]dxdy = [[integral].sub.[^.Q]]([1/[h.sub.1]][4.summation over (n=1)][u.sub.n][[partial derivative].sub.[^.x]][[phi].sub.n]([^.x], [^.y]))[.sup.2] det [nable] Fd[^.x]d[^.y],

=[[integral].sub.[^.Q]]([1/[h.sub.1]][4.summation over (n=1)][u.sub.n][[partial derivative].sub.[^.x]][[phi].sub.n]([^.x], [^.y]))[.sup.2] [h.sub.1][h.sub.2]d[^.x]d[^.y]. (18)

The procedure is the same for

[[integral].sub.Q.sub.k](([[partial derivative].sub.y]u)[.sup.2] - 1)[.sup.2] dxdy.

2.2 Computation of Surface Energy

We evaluate the surface energy term in Eq. (4) using Eq. (6), approximating [[partial derivative].sub.y]u using backward finite differences,

[[partial derivative].sub.y]u([bar.x], [y.sub.j]) [approximately equal to] [u([bar.x], [y.sub.j]) - u([bar.x], [y.sub.j-1])]/[h.sub.2]

and

[[partial derivative].sub.y]u([bar.x], [y.sub.j+1]) [approximately equal to] [u[bar.x], [y.sub.j+1]) - u([bar.x], [y.sub.j])]/[h.sub.2]

Here,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

where N([h.sub.2]) is the number of nodes along the y-direction. We can integrate along the x-direction to complete the approximation of Eq. (4).

3. Numerical Implementation

The computation of the integral in Eq. (18) is done exactly in the parent element [^.Q]. The symbolic package Maple (1) [30] was employed to integrate this term exactly. The algebraic expression consisting of nodes [u.sub.j] and mesh sizes [h.sub.1] and [h.sub.2], is evaluated over each element.

The gradient of Eq. (1) is approximated using cellcentered finite differences. Let u [member of] [W.sup.1,4] ([OMEGA]). The Gateaux derivative is:

G(u) = [lim.[t[right arrow]0]][1/t][J(u + t[phi]) - J(u)] (20)

where [phi] [member of] [W.sup.1,4] ([OMEGA]) [31]. The linear operation G is the directional derivative of J in the direction [phi]. From the Euler-Lagrange equations, we see that the gradient of J is difficult to compute. Considering the Euler-Lagrange equation we have,

G(u) = [d/dt]J(u + t[phi])|[.sub.t=0] = [[integral].sub.[OMEGA]][[[partial derivative]J(u)]/[[partial derivative][nabla]u]][nabla][phi]dx,

= -[[integral].sub.[OMEGA]]div[[[partial derivative]J]/[[partial derivative][nabla]u]][phi]dx,

where the gradient is div ([partial derivative]J/[partial derivative][nabla]u). Let J(u) be the discretized representation of the total energy Eq. (1) for u [member of] [R.sup.N(h)], where N(h) is the total number of nodes in [OMEGA].

With this numerical approximation to [nabla]J(u), we use the limited memory BFGS method [32] to numerically minimize J(u). This inexpensive quasi-Newton method seeks to build an approximation [H.sup.(k)] to the inverse of the Hessian matrix of second derivatives of J at the point [u.sup.(k)]. These inexpensive approximations are constructed from a small number of vectors updated at every iteration.

The k + 1 step of the limited memory BFGS method is given by:

[u.sup.(k+1)] = [u.sup.(k)] - [[alpha].sup.(k)][H.sup.(k)][nabla]J([u.sup.(k)]).

where the parameter [[alpha].sup.(k)] is a line search parameter and is computed at each iteration using a line-search procedure.

The Hessian matrix [B.sup.(k)] satisfies the secant equation:

[B.sup.(k)][s.sup.(k)] = [y.sup.(k)],

where

[s.sup.(k)] = ([u.sup.(k)] - [u.sup.(k-1)]), and [y.sup.(k)] = [nabla]J([u.sup.(k)]) - [nabla]J([u.sup.(k-1)]).

We employ a limited memory BFGS method because the Hessian matrix [B.sup.(k)] is expensive to compute since we are solving a large scale optimization problem [32,33]. This is done by storing a certain number of vector pairs {[s.sup.(k)], [y.sup.(k)]}. At each new iterate, the oldest vector pair is deleted and replaced by the new vector pair, thus preserving curvature information only from the most recent iterations. The formula for updating the inverse of the Hessian matrix [H.sup.(k)] = ([B.sup.(k)])[.sup.-1] is:

[H.sup.(k + 1)] = ([V.sup.(k - m)] ... [V.sup.(k-1)])[.sup.T] ([H.sup.(k)])[.sup.0] ([V.sup.(k - m)] ... [V.sup.(k - 1)]) + [[rho].sup.(k - m)]([V.sup.(k - m + 1)] ... [V.sup.(k - 1)])[.sup.T][s.sup.(k - m)] ([s.sup.(k - m)])[.sup.T]([V.sup.(k - m + 1)] ... [V.sup.(k - 1)]) + [[rho].sup.(k - m + 1)]([V.sup.(k - m + 2)] ... [V.sup.(k - 1)])[.sup.T][s.sup.(k - m + 1)]([s.sup.(k - m + 1)])[.sup.T] ([V.sup.(k - m + 2)] ... [V.sup.(k - 1)]) +...+ [[rho].sup.(k - 1)][s.sup.(k - 1)]([s.sup.(k - 1)])[.sup.T].

where ([H.sup.(k)])[.sup.0] is some initial Hessian approximation,

[V.sup.(k)] = I - [[rho].sup.(k)][y.sup.(k)]([s.sup.(k)])[.sup.T],

and

[rho] = 1/[([y.sup.(k)])[.sup.T][s.sup.(k)]].

4. Results

In this section, results computed on the various grids, 30 X 30, 60 X 60, and 120 X 120, are presented. In all minimization, the initial iterate [u.sup.0] is a small ([approximately equal to] [10.sup.-3]) perturbation of pure Austenite, which we take to be u = 0. Density plots of the minimizer u are shown in Figs. 5-10 and demonstrate the existence of the desired microstructures in the twinned-Martensite phase. These figures also demonstrate the role of the surface energy as a penalization term of Eq. (1). We seem to be able to control the number of microstructures by varying the values of [epsilon], with a larger value of [epsilon] resulting in a smaller number of microstructures in Martensite and viceversa. Figures 11, 12, and 13 show plots of the profiles of minimizer u at x = 4. Figures 11, 12, and 13 exhibit a larger number of discontinuities in [[partial derivative].sub.y]u for [epsilon] = 2 X [10.sup.-6] than for [epsilon] = 2 X [10.sup.-2]. In Figs. 5-9 one sees the effect of grid-size comparing results on 60 X 60 grid with 120 X 120 grid for [epsilon] = 2 X [10.sup.-6] through [epsilon] = 2 X [10.sup.-2]. The black and white stripes correspond directly to the saw-toothed behavior of u at x = 4. For small [epsilon] values, the number of discontinuities is close to the number of partitions along the y-direction. Tables 2-4 show the energy values for the minimizers shown in Figs. 5-10. The column with E(u) corresponds to the elastic energy values of the minimizer while the column with S(u) corresponds to the surface energy values of the minimizer. The energy values presented in Tables 3 and 4 appear consistent with these conclusions. The larger the [epsilon], the larger the value of the surface energy at the minimizer. In Table 3 the value of the surface energy increases by a larger order than that of the elastic energy.

One final and curious observation is the appearance of a "diagonal band" structure in Figs. 5-10. It would seem that this band, which is apparent for solutions on various grid sizes, is related to the competing roles of the surface and elastic energies is and may be associated to the "equipartitioning of energy" principle proposed in Kohn and Muller in Ref. [10,11]. Future research plans include a deeper investigation of this diagonal-band structure.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

Table 1. Number of "jumps" in [[partial derivative].sub.y]u

Function     [[integral].sub.0.sup.1]|          Discontinuities
             [[partial derivative].sub.yy]u|dy

Amplitude 1          24                                12
Amplitude 2          12                                 6
Amplitude 3           6                                 3

Table 2. Energy values for various [epsilon] on 120 X 120 grid

[epsilon]          E(u)   [epsilon]S(u)  Total energy

2 X [10.sup.-6]   6.3014  0.5691e-2       6.3071

2 X [10.sup.-5]   6.1122  0.5698e-1       6.1692
2 X [10.sup.-4]   7.0859  0.5474          7.6334
2 X [10.sup.-3]   8.2222  4.5200         12.7423
2 X [10.sup.-2]  23.3218  1.2559         24.5770

Table 3. Energy values for various [epsilon] on 60 X 60 grid

[epsilon]        E(u)    [epsilon]S(u)  Total energy

2 X [10.sup.-6]  4.3684   0.3120e-2       4.3716
2 X [10.sup.-5]  6.3005   0.2808e-1       6.3286
2 X [10.sup.-4]  6.3010   0.2798          6.5809
2 X [10.sup.-3]  5.7738   2.7356          8.5095
2 X [10.sup.-2]  3.7945  14.6825         18.477

Table 4. Energy values for various [epsilon] on 30 X 30 grid

[epsilon]         E(u)   [epsilon]S(u)  Total energy

2 X [10.sup.-6]  4.6093  0.0015          4.6108
2 X [10.sup.-5]  4.9832  0.0145          4.9977
2 X [10.sup.-4]  6.1073  0.1326          6.1499
2 X [10.sup.-3]  5.5335  1.3246          6.8581
2 X [10.sup.-2]  8.0650  9.6770         17.7420

Accepted: January 6, 2004

Available online: http://www.nist.gov/jres

(1) Certain commercial equipment, instruments, or materials are identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose.

5. References

[1] M. Chipot and D. Kinderlehrer, Equilibrium configurations of crystals, Arch. Rational Mech. Anal. 103, 237-227 (1988).

[2] J. L. Erickssen, Some phase transitions in crystals, Arch. Rational Mech. Anal. 73, 99-124 (1980).

[3] J. L. Erickssen, Twinning of Crystals I, Metastability and Incompletely Posed Problem, IMA Vol. Math. Appl. 3, Springer. New York (1987) pp. 77-96.

[4] R. D. James and D. Kinderlehrer, Frustration in ferromagnetic materials, Contin. Mech. Thermodyn. 2, 215-239 (1990).

[5] R. D. James and D. Kinderlehrer, Mathematical approaches to the study of smart materials, in Proc. 1993 N. Amer. Conference on Smart Structures and Materials, H. I. Banks, ed., SPIE 1919, 2-18 (1993).

[6] R. V. Kohn, Mathematics in Materials Science, Opportunities for the Mathematical Sciences, Division of Mathematical Sciences, National Science Foundation Workshop, June 26-27, 2000.

[7] R. G. Linford, The Derivation of Thermodynamic Equations for Solid Surfaces, Chemical Reviews, Vol. 78, No. 2, April 1978.

[8] C. R. Barret, W. D. Nix, and A. S. Tetelman, The Principles of Engineering Materials, Prentice Hall (1973).

[9] R. Choksi and R. V. Kohn, Bounds on the Micromagnetic Energy of a Uniaxial Ferromagnet, Comm. on Pure and Appl. Math. LI, 259-289 (1998).

[10] R. V. Kohn and S. Muller, Branching of twins near an austenite-twinned-martensite interface, Phil. Mag. A 66 (5), 697-715 (1992).

[11] R. V. Kohn and S. Muller, Surface Energy and Microstructure in Coherent Phase Transitions, Comm. Pure Appl. Math. 47 (4), 405-435 (1994).

[12] L. Landau and E. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies, Phy. Z. Sowjet. 8, 153 (1935).

[13] I. A. Privorotskii, Theory of Domain Structure of Uniaxial Ferromagnets, Soviet Physics JETP 32 (5), 964-970 (1971).

[14] I. A. Privorotskii, Thermodynamic Theory of Domain Structures, Israeli University Press (1976).

[15] C. Collins and M. Luskin, The Computation of the Austenitic-Martensitic Phase Transition, Partial Differential Equations and Continuum Models of Phase Transition, Lecture Notes in Physics, Vol. 344, Springer-Verlag, New York (1989) pp. 34-50.

[16] P. Kloucek and L. A. Melara, Computational Modelling of Branching Fine Structures in Constrained Crystals, J. Computational Phys. 183, 623-651 (2002).

[17] B. Li and M. Luskin, Finite element analysis of microstructure for the cubic to tetragonal transformation, SIAM J. Numer. Anal. 35 (1), 376-392 (1998).

[18] B. Li and M. Luskin, Nonconforming finite element approximation of crystalline microstructure, Math. Comp. 67 (223), 917-946 (1998).

[19] B. Li and M. Luskin, Approximation of a martensitic laminate with varying volume fractions, [M.sup.2] AN 33 (1), 67-87 (1999).

[20] B. Li and M. Luskin, Theory and computation for the microstructure near the interface between twinned layers and a pure variant of martensite, Materials Sci. Eng. A, 273-275:0, 237-240 (1999).

[21] B. Li and Z. Zhang, Analysis of a class of superconvergence patch recovery techniques for linear and bilinear finite elements, Numer. Meth. PDEs 15, 151-167 (1999).

[22] K. Bhattacharya, B. Li, and M. Luskin, The simply laminated microstructure in martensitic crystals that undergo a cubic to orthorhombic phase transformation, Arch. Rational Mech. Anal. 149 (2), 123-154 (1999).

[23] B. Li, Approximation of martensitic microstructure with general homogeneous boundary data, J. Math. Anal. Appl. 266, 451-467 (2002).

[24] R. A. Nicolaides, N. Walkington, and H. Wang, Numerical Methods for a Nonconvex Optimization Problem Modeling Martensitic Microstructure, SIAM J. Sci. Comput. 8 (4), 1122-1141 (1997).

[25] R. L. Wheeden and A. Zygmund, Measure and Integral, An Introduction to Real Analysis. Marcel Dekker, Inc. (1997).

[26] H. L. Royden, Real Analysis. Third Edition, Prentice Hall (1988).

[27] F. Jones, Lebesgue Integration on Euclidean Space, Jones and Bartlett Publishers (1993).

[28] W. P. Ziemer, Weakly differentiable functions, Graduate Texts in Math., Springer-Verlag, New York (1989).

[29] P. G. Ciarlet, Basic Error Estimates for Elliptic Problems, Handbook of Numerical Analysis, Vol. II, Elsevier Science Publishers (1991).

[30] M. B. Monagan, K. O. Geddes, K. M. Heal, G. Labahn, S. M. Vorkoetter, and J. McCarron, Maple 6: The Standard for Analytical Computation, Waterloo Maple Inc. (2000).

[31] L. V. Kantorovich, Functional Analysis in Normed Linear Spaces, MacMillan (1964).

[32] D. C. Liu and J. Nocedal, On the Limited Memory BFGS Method for Large Scale Optimization, Math. Prog. 45, 503-528 (1989).

[33] J. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Math. of Computation 35, 151, 773-782 (1980).

A.J. Kearsley and L. A. Melara Jr.

National Institute of Standards and Technology, Gaithersburg, MD 20899-0001

[e-mail]

[e-mail]

About the authors: Anthony J. Kearsley is a member of the Mathematical and Computational Sciences Division at NIST. His research involves developing numerical algorithms. In 1998, he was awarded a Presidential Early Career Award for his work on Numerical Methods for Optimization. Luis A. Melara is a National Research Council Postdoctoral Associate in the Mathematical and Computational Sciences Division at NIST. The National Institute of Standards and Technology is an agency of the Technology Administration, U.S. Department of Commerce.