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.

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 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.

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-Müller in [10,11]. The Martensite region is represented by the domain Ω = (0, L) × (0, K). Let x = (x, y) ∈ Ω ⊂ and u = u(x, y). The double-well energy function is given by (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 which minimizes Eq. (1) where

Elastic Energy
The 2-D scalar model of elastic energy, in Eq. (1) is (2) This term is minimized by a function u such that Each of the gradients in Eq. (3) corresponds to a stressfree 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 ∂ y u. We note that the function with Amplitude 1 also has the largest number of discontinuities, occurring at the peaks of u, where ∂ 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 ∂ y u. At the minimizer u, however, the surface energy penalizes these discontinuities.

Surface Energy
We employ the surface energy presented by Kohn and Müller in Ref. [10,11], (4) where ε is a constant.
The definition proposed by Kohn and Müller replaces the integral of Eq. (4) in the y-direction with the total variation of ∂ y u over (0, K) [10,11]. This change designates the role of the surface energy at the minimizer u as a "counter" because counts the number of discontinuities in ∂ y u, where ∂ y u ∈ {±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.

Functions of Bounded Variation
The definition of functions of bounded variation is taken from Refs. [25,26,27]. Let be any partition of (0, K) and let   The partitioning P is a collection of points y j , j = 0, 1, …, l, such that y 0 = 0 and y l = K. For each partitioning P, let (5) for a fixed ∈ (0, L). The total variation of ∂ y u over (0, K) is (6) where the supremum is taken over all countable partitions P of (0, K).
In Fig. 3a, the piecewise linear function u(x, y) is plotted over (0, 1). In Fig. 3b, the function ∂ y u(x, y) is plotted over (0, 1). Clearly, the number of discontinuities in ∂ y u is 5 from Fig. 3b. For a discontinuity at the partition point y j , we adopt the following convention: (10) Now consider (y j-1 , y j+1 ) = (y j-1 , y j ) (y j , y j+1 ) in Fig. 3b. Evaluating at the endpoint of the two subintervals gives A jump occurred in (y j-1 , y j ) because ∂ y u is discontinuous at a point in this subinterval. Based on our conventional definition of ∂ y u(x, y j ) at a point y j in Eq. (10), we say a "jump" occurs at y j if ∂ y u is discontinuous at this point. This "jump" refers to the change in values of ∂ y u at y j from +1 to -1 or vice-versa.
This jump gives the nonzero value in Eq. (11). In the subinterval (y j , y j+1 ), no jump occurs in ∂ y u, therefore the difference in Eq. (12) is zero. We compute Eq. (9): 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.
Volume 108, Number 6, November-December 2003 Journal of Research of the National Institute of Standards and Technology Surface energy term as a counter. Pictorial example illustrating the role of surface energy as a counter term. The notation u(·, y) means that the graph shows the y-dependence only.  and the surface energy succeeds in serving as a "counter." The next section describes the numerical approximation to Eq. (1).

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: (13) because functions in this polynomial space lie in , 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) × (0, 1) and for some , the reference element is where h 1 is the mesh size along the x-axis and h 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 1 ), u(a 2 ), u(a 3 ), and u(a 4 ). Figure  4 shows the parent element with the vertices also denoted by concentric circles.

Affine Finite Elements
Let F : Q→Q k be given by (14) We simplify, by letting F = F(x, y). We use the affine map Eq. (14) to evaluate the function u, to compute ∇u and to integrate. The function u is spanned by the four basis functions defined in each element Q k . Using Eq.
, we can determine all basis functions by the four basis functions in Q. Thus, the approximation to u is: where x = x(x, y), y = y(x, y). The nth basis function in Q is and it satisfies the property that φ n (a m ) = δ mn , for n, m = 1, 2, 3, 4, where δ mn is the Kronecker delta function.
To compute the gradient of u and to integrate u, we need: ∇F, det∇F and ∇F -1 . The gradient of F is (16) and, we have .
The procedure is the same for

Computation of Surface Energy
We evaluate the surface energy term in Eq. (4) using Eq. (6), approximating ∂ y u using backward finite differences, and Here, (19) where N(h 2 ) is the number of nodes along the y-direction. We can integrate along the x-direction to complete the approximation of Eq. (4).

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 j and mesh sizes h 1 and h 2 , is evaluated over each element.
The gradient of Eq. (1) is approximated using cellcentered finite differences. Let u ∈ (Ω). The Gateaux derivative is: (20) where ϕ ∈ (Ω) [31]. The linear operation G is the directional derivative of J in the direction ϕ. From the Euler-Lagrange equations, we see that the gradient of J is difficult to compute. Considering the Euler-Lagrange equation we have, where the gradient is div (∂J/∂∇u). Let J(u) be the discretized representation of the total energy Eq. (1) for u ∈ ∈ , where N(h) is the total number of nodes in Ω. With this numerical approximation to ∇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 (k) to the inverse of the Hessian matrix of second derivatives of J at the point u (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: where the parameter α (k) is a line search parameter and is computed at each iteration using a line-search procedure.
The Hessian matrix B (k) satisfies the secant equation: where We employ a limited memory BFGS method because the Hessian matrix B (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 (k) , y (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 (k) = (B (k) ) -1 is: 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.

Results
In this section, results computed on the various grids, 30 × 30, 60 × 60, and 120 × 120, are presented. In all minimization, the initial iterate u 0 is a small (≈10 -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 ε, with a larger value of ε 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 ∂ y u for ε = 2 × 10 -6 than for ε = 2 × 10 -2 . In Figs. 5-9 one sees the effect of grid-size comparing results on 60 × 60 grid with 120 × 120 grid for ε = 2 × 10 -6 through ε = 2 × 10 -2 . The black and white stripes correspond directly to the saw-toothed behavior of u at x = 4. For small ε values, the number of discontinuities is close to the number of partitions along the y-direction. Tables 2-4 Tables 3 and  4 appear consistent with these conclusions. The larger the ε, 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 Müller in Ref. [10,11]. Future research plans include a deeper investigation of this diagonal-band structure.