Geometry of Needle-Like Microstructures in Shape-Memory Alloys

Needle-like microstructures are often observed in shape memory alloys near macro-interfaces that separate regions with different laminate orientation. We study their shape with a two-dimensional model based on nonlinear elasticity, that contains an explicit parametrization of the needle profiles. Energy minimization leads to specific predictions for the geometry of needle-like domains. Our simulations are based on shape optimization of the needle interfaces, using a polyconvex energy density with cubic symmetry for the elastic problem, and a numerical implementation via finite elements on a dynamically changing grid.


Introduction
Martensitic materials are well known to exhibit rich microstructures that are associated with the solid-solid phase transformation [1,2].It is widely believed that the formation of those patterns leads to an energy barrier which plays a crucial role in the hysteresis and reversibility of the phase transformation itself [3,4].Mathematically, these microstructures are typically modelled within the framework of the phenomenological theory of martensite [5], i.e., as minimizers of continuum elasticity functionals of the form Z X Wðru; TÞ dx; ð1:1Þ where X denotes the reference domain and W is a free elastic energy density which reflects the crystallography of the phase transformation and depends on the deformation gradient and the temperature.Here we work at fixed temperature below the transformation temperature, and we drop it from the notation.This modelling ansatz has proven successful to predict certain properties of the microstructures such as orientations of interfaces.However, it is often not sufficient to predict length scales of multiscale patterns.Therefore, one often introduces a regularizing surface energy term that penalizes interfaces between regions of different martensitic variants.Starting with the work by Kohn and Mu ¨ller in the early 90s [6,7], there has been a huge body of literature on the understanding of minimizers of such functionals in simplified settings, in particular scaling laws for the minimal energy (see for example [8][9][10][11][12] and references therein).These results generally indicate that for certain materials one expects rather uniform structures close to interfaces, while for others, complicated multiscale branching patterns are predicted.This is in accordance with experimental results (see e.g.[13]).However, only in very few and restricted situations, there are results on more quantitative properties of minimizers, such as periodicity or self-similarity (see [14][15][16]).
In this note, we report results on the modelling and quantitative study of special microstructures, so called needle-type patterns, building in particular on the recent works [17][18][19].Needle-like patterns consist of laminated structures where one variant thins out close to the interface.These patterns are often found in experiments, and appear to be rather stable.It is remarkable that they occur in very different materials, including shape-memory alloys (for instance, NiAl) and perovskites (for instance, YBaCuO), at very different length-scales ranging from a few nanometers to tenths of micrometers, and at both, martensite/martensite and martensite/austenite interfaces (see e.g.[13,[20][21][22]).
It has been found that models based on finite elasticity are well-suited to describe such needle patterns, while linearized elasticity is not able to predict the relevant length scales.This has been observed in various numerical experiments (see [17] and the references discussed there) and has been made rigorous in terms of an energy scaling law for the minimal energy [18].Therefore, we deal only with numerical experiments in the context of finite elasticity.We focus on three materials in which needle-like patterns have often been observed, namely NiAl, CuAlNi and YBaCuO.We study the dependence of the needle shape on various parameters, in particular the elastic constants (the anisotropy ratio and the Poisson ratio), the eigenstrain, and the volume fraction of the laminate.To obtain quantitative results, a careful and physically reasonable choice of the free energy density W is essential.Therefore, after setting the notation in section ''Geometry of Needles'', we will discuss our choice of the free energy density in section ''Energy Density''.The numerical implementation is described in section ''Shape Optimization, Numerical Implementation'', and the results are presented in section ''Results''.

Geometry of Needles
We focus on the structure of needle-type interfaces, and make an ansatz which in particular excludes branching of the needles.It is widely believed that branched patterns of two martensitic variants are favored if the transition layer is sufficiently wide compared to the height of the laminate, see e.g.[13].This can be seen already in the following back-of-the-envelope computation for scalar-valued toy models [7]: Assume that a planar microstructure is mainly described by the out-of plane component of the displacement u : R !R in a rectangular transition layer R :¼ ð0; LÞ Â ð0; HÞ, where we assume the interface to be at the left edge of the rectangle (at fx ¼ 0g).Within this ansatz, an energy based on the phenomenological theory from [5] is essentially estimated by the Kohn-Mu ¨ller-type energy (see e.g.[23] for a computation) where the preferred gradients ð0; hÞ and ð0; À1 þ hÞ represent the two variants of martensite, and the second term is to be understood distributionally.

ðh; 1Þ:
& There is a large body of literature on this model, and the scaling law of the minimal energy in various settings is well-understood.It has been shown that many properties of more general vectorial models (see e.g.[8,10,11,24] and the references therein) are already captured by such toy models.We refer to [25] for a recent discussion in relation with experimental images.
The test functions that have been used to prove the upper bound in the energy scaling laws quantify the expectation that for small values of the surface energy e, branched patterns are expected, while for large values of e, optimal structures are more uniform.We focus here on the case of small e where complex pattern formation is expected.More precisely, within the model described above, for branching-type patterns with periodic boundary conditions (see the sketch in Fig. 1) the energy scales like eL þ h 2 H 3 L (see e.g.[10, Lemma 5.1] with g ¼ hH for a building block and the branching-type constructions described thereafter, as well as the references given there).On the other hand, simple needle-type structures have an energy scaling like eL þ h 2 H 3 ' þ h 2 'H, where ' is the tapering length of the needle, or equivalently of the transition region (see also [4,23] for an analogous computation in the context of linearized elasticity).Therefore, in terms of energy scaling, needles are competitive compared to branching if H $ ' $ L, and can thus be even energetically preferred, depending on the precise constants.
In this note, we continue and complement the recent studies in [17,19] and consider a martensitic macrotwin where a laminate of two variants of martensite of volume fractions h and 1 À h meets a rotated martensitic variant at a straight interface (see Fig. 1).The model is vectorial and geometrically nonlinear.The two variants are given by eigenstrains U 1 and U 2 , which we can assume to be of the form (see [26,Section 5]) where d measures the strain; relevant values are given in Table 1.These variants can form a laminate with (vertical) normal e 2 , and a macrotwin as in Fig. 1 with a normal vector parallel to e 1 þ dhe 2 , i.e., a vector rotated by an angle of order dh with respect to the (horizontal) vector e 1 (for details see e.g.[18, Lemma 2.1]).
To determine the optimal needle geometry within this ansatz, we perform a two-level optimization: Given the boundary curves of the needle, we minimize the elastic energy in the transition layer [see (4.1)], the result only depends on the shape of the needle (see Fig. 2).The optimal needle structure is then the one for which this minimal elastic energy is itself minimized.The specific boundary conditions and the class of ansatz functions are presented in section ''Shape Optimization, Numerical Implementation''.An essential ingredient here is the free energy density, that is discussed in section ''Energy Density''.
As mentioned above, we restrict ourselves to the setting of nonlinear elasticity.Although the linear theory has proven useful to explain higher order effects (see e.g.[27] and the references therein), it has been found that approximations in terms of linearized elasticity are not appropriate to understand the qualitative behaviour in terms of the tapering length scale (see [18] for a heuristic as well as a rigorous argument).Roughly speaking, the main factor determining the length scale of the pattern seems to be the change in rotation: Close to the macrotwin interface the relevant rotation is the one which is determined by the other rank-one connection between the martensitic variants, which differs from the one used in the laminates by an angle of order dh.

Energy Density
We shall write the energy density W i : R 2Â2 !½0; 1 of the i-th martensitic variant, with eigenstrain U i as in (2.1), in terms of the energy density W A of the austenite, which we take as reference configuration, as Therefore, it suffices to identify a suitable function W A .We start listing our requirements on this function, then discuss some proposals present in the literature and motivate our choice.
Requirements on W A The energy density W A needs to satisfy the invariance property Fig. 1 Sketch of a laminate ending in a needle (left) and of a branching pattern (right).As both patterns are periodic, in the following we focus only on one period.The macrointerface is drawn vertical for simplicity, the true slope has an angle of order dh with the vertical for some convex and lower-semicontinuous function g : R 2Â2 Â R ! ½0; 1.We recall that polyconvexity implies lower semicontinuity of variational problems of the form R X ðW A ðDuÞ À f Á uÞdx, and is therefore crucial in proving existence of minimizers.Further, polyconvexity implies the Legendre-Hadamard condition, which in turn is equivalent to rank-one convexity, i.e., to the fact that the function We refer to [28] for a mathematical background on these concepts.Finally, we require W A to reproduce the experimentally known elastic constants of austenite around the identity, in the sense that where C 2 R 2Â2Â2Â2 is the elasticity tensor.In the cubic case this reduces to where the elastic constants C 12 , C 44 and C 0 :¼ ðC 11 À C 12 Þ=2 are all positive in the relevant situation (see Table 1).The cubic anisotropy of the material is normally quantified by the parameter A defined by one can easily check that for isotropic materials A ¼ 1.In summary, our aim is an expression for W A which has the invariances stated in (3.2), is minimized on rotations as in (3.4), is polyconvex in the sense of (3.5), linearizes to (3.8) around the identity, and has sufficient growth and regularity to permit easy and stable computations.
Some candidates from the literature In the literature, one often uses a linear elastic energy of the form 1 2

P
ijkl C ijkl ij kl evaluated on the nonlinear elastic strain 1 2 ðF T F À IdÞ, resulting in expressions of the type  However, the polynomial pðtÞ :¼ ðt 2 þ 2tÞ 2 ¼ ðt þ 2Þ 2 t 2 is not convex.Therefore W C is not rank-one convex, and not polyconvex, for any sensible choice of the elastic coefficients C. One expects that this lack of convexity leads to numerical instabilities if this part of the space of strains (which is relatively far from the identity, as one can see from det F À1=2 ¼ 1=2) is explored in a numerical simulation; apparently this was not the case in the cited papers.
A polyconvex expression with cubic symmetry was proposed by Kambouchev et al. in [31, Eq. ( 9)].In 3D, they define and show that if the four coefficients are nonnegative then Comparing with the cubic linearization (3.8) they obtain a 3 ¼ ðC 0 À C 44 Þ=4 and a 4 ¼ ð2C 44 À C 0 Þ=2, which are nonnegative only if the anisotropy ratio fulfills 1 2 A 1 ð3:14Þ (see [31, Eq. ( 14)]).This condition is fulfilled by many metals but not by the shape-memory alloys of interest here, which have larger values of A (see Table 1).In two dimensions, one would equivalently write with similar results.We show that the relevant condition A 1 is indeed necessary for polyconvexity of (3.15).It suffices to consider matrices of the form which give For large t the leading-order term is a 3 t 4 , which is convex only if a 3 !0. Therefore the condition a 3 !0, which implies A 1, is necessary for rank-one convexity and, hence, for polyconvexity.
A variant of this formula was used in [17], precisely, The change was that the first term is fourth-order and not quadratic, so that it matches the order of the second and the fourth one.The expression in (3.16) is obviously polyconvex if the four coefficients are nonnegative.Even more, one can show that it is polyconvex provided a 1 !0, a 2 !0, a 3 !0, and a 4 !À 2 3 a 1 .To see this, it suffices to first check (for example, computing the Hessian) that the function g : R 2 !R, gðx; yÞ :¼ x 4 þ y 4 þ 6x 2 y 2 , is convex, and then to observe that by the chain rule this implies convexity of F7 !gðjFe 1 j; jFe 2 jÞ ¼ 3jFj 4 À 2jFe 1 j 4 À 2jFe 2 j 4 .
The difference from (3.15) to (3.16) however improves the range of admissible values of A. In fact, condition (3.4) implies 4a 1 þ a 2 À a 3 þ 2a 4 ¼ 0. Linearizing (3.16), we obtain that (3.8) is equivalent to ð3:17Þ which readily shows that the largest value of A that can be attained in the polyconvex range is A ¼ 2. This improvement over (3.14) is not sufficient for the present purposes.Indeed, in [17], using values for the elastic constants appropriate for NiAl listed in Table 1, the coefficients a 1 ¼ 11:56 GPa, a 2 ¼ À17:44 GPa, a 3 ¼ 10:04 GPa, a 4 ¼ À9:38 GPa were obtained (up to an irrelevant global factor of 2).Clearly they are not in the range in which the expression above is guaranteed to be polyconvex.In the numerical computations presented in [17], as in the case of the papers based on (3.10), this did not cause any difficulty.
In closing this review we show that the expression in (3.16) is not rank-one convex whenever a 4 \ À 2 3 a 1 , which is the same as A [ 2. Consider for some b [ 0 the rankone line Clearly det F t ¼ 1 for all t and b, so that the two terms depending on the determinant are irrelevant, and W 2020 ðF t Þ is a fourth-order polynomial in t.One computes 3 a 1 then the first coefficient is negative, and choosing b sufficiently small leads to a negative second derivative in a rank-one direction.
Construction of W A In [19] we developed a new expression that is able to obtain polyconvexity for arbitrarily large values of the cubic anisotropy parameter A. It depends on four independent small parameters q 1 , q 2 , q 3 , q 4 2 ð0; 1Þ via the functions h q : ½0; 1Þ !R, h q ðxÞ :¼ 0; if x\1 À q; ðx À ð1 À qÞÞ 2 ; if x ! 1 À q; & ð3:20Þ g q : ½0; 1Þ 2 !R, g q ðx; yÞ :¼ ð1 À ð1 À qÞ 2 Þx 2 ; if y\ð1 À qÞx; ð3:21Þ and f q : R ! ½0; 1, One can easily check that all three functions are convex, continuously differentiable (whenever they are finite), lower semicontinuous, and that the first two are nondecreasing.Finally, for any q 2 ð0; 1Þ 4 and b 2 ð0; 1Þ 4 we set ð3:23Þ This function is polyconvex, has the desired invariance and minimizers.Computing the Taylor series to second order and comparing to (3.8) one can see that where

ð3:24Þ
The matrix Id þ M q is, for q small, invertible.Further, for any positive value of the three elastic constants one obtains that if the q's are chosen sufficiently small then the resulting b's are nonegative, and W A is polyconvex.The smallness of q depends, however, on C. Examples are given in Table 2 below.

Parametrization of Subdomains
In what follows we will decompose the energy in the two phase components with fixed topology using the domain decomposition approach developed in [17,19].To this end we consider a physical reference configuration X ¼ X 1 [ X 2 where the actual domains of the martensitic variants X 1 , X 2 are parametrized over subdomains X1 , X2 of a computational reference configuration X via a parameter-dependent bijective transformation w½a : X !X with X i ¼ w½að Xi Þ, where a is a small set of parameters (Fig. 3).
Given the parametrization w½a of the physical domains the energy can be rewritten as We minimize Ê with respect to a and with respect to /. Minimizing over / for fixed a and thus for fixed geometry reflects the solution of the elastic state equation in reference coordinates x.Let us denote the minimizer by /½a.Minimizing J ½a :¼ Ê½a; /½a ð4:3Þ over a turns into solving the actual free boundary problem.Thus, in the presented form the free boundary problems can be considered as a joint minimization with respect to a and /.The solution to this elastic shape optimization problem will be described below in section ''Optimization''.The physical domain parametrization w½a primarily parametrizes the geometry of the phase boundaries.The parametrization of the interior of the phases can be chosen freely as long as the mapping w½a is bijective and for implementational purposes sufficiently regular.In what follows, we explicitly describe w½a on o X1 which in particular also determines the transformation of o X2 .
In the computational reference configuration o X1 is assumed to be polygonal.For each edge ê from â ¼ ð â1 ; â2 Þ to b ¼ ð b1 ; b2 Þ we consider w½aj ê as a graph over x1 (which reflects that the needles are horizontally stretched).
Using quadratic polynomials for the parametrization of the upper and lower edge of the needle tip we obtain The above parameters a and b are given functions of the needle length ' and the shift parameter D.Besides ' and D the parameter vector a includes the parameters j for the upper and the lower edge.
To define w½a not only on o X1 and o X but also in the interior of the phases (and thus on the whole domain X) we use a harmonic extension approach.We define w½a : X !X by If w½a À Id is small and sufficiently regular on o X1 , the extended w½a is ensured to be bijective.

Finite Element Discretization
We consider piecewise affine finite elements on the computational domain X, where o X1 is composed of edges of the triangulation T h with vertices xh , the index h indicating the grid size.The corresponding (vector valued) finite element space is denoted by V h .
The extension (4.4) of the phase boundary parametrization on o X1 then reads as follows: ð4:6Þ into a standard conforming finite element discretization of the Euler-Lagrange-equation of (4.2) and requires the solution of a nonlinear system of equations in V h , performed here via Newton's method.For fixed a, the minimizer of Êh ½a; Á is then denoted by /h ½a and an evaluation of Êh then immediately yields the value of the discrete shape functional J h ½a :¼ Êh ½a; /h ½a.

Optimization
So far, the free boundary problem is described as an optimization problem of the parametrization w½a of the phase boundary as a function of the parameter vector a subject to the constraint that / minimizes the elastic energy given in (4.2).As described above, the finite element discretization of this optimization problem turns into a minimization problem with nonlinear constraints.For the solution of the discrete state equation we use the finite element library FEniCS.For details we refer to [36][37][38][39][40][41][42][43][44].
To evaluate the discrete shape gradient o a J h ½a via the adjoint calculus we use the additional package Dolfin Adjoint [45][46][47].As input this only requires the solution of the forward map a7 !/ h ½a via solving the state equation.In fact, Dolfin Adjoint tracks all steps of this computation and automatically configures the adjoint derivative.Based on this shape gradient, a limited memory BFGS method with constraints [48] is used for the actual optimization (with SciPy minimization parameters ftol ¼ 10 À16 , gtol ¼ 10 À12 and maxiter ¼ 200).

Results
In this section, we present numerical results for the optimal shapes in different materials.
Figure 4 shows the needle geometries for the three different materials NiAl, CuAlNi and YBaCuO in the physical reference configuration.One clearly observes a significant influence of the material constants on the needle length.If one rescales the optimal needle geometries so that they Fig. 3 The computational domain X (left) is transformed by w½a into the physical reference domain X (middle), which is then deformed elastically by / to the deformed configuration /ðXÞ (right) have the same length, they become very similar.In the rest of this section we discuss the effect of various parameters on the length of the needle.
Influence of the order parameter d In Fig. 5, we show the dependence of the needle length on the order parameter d for the elastic parameters of NiAl.It is apparent that the needle length has a very strong dependence on the order parameter d and it diverges as d !0, as was already discussed in [17,18].Indeed, this dependence of the needle length on the order parameter explains most of the variability observed in Fig. 4 The appropriate b i are then calculated via Eq.(3.24).The results are shown in Fig. 6.
Influence of the volume fraction h Figure 7 shows that needle tips are shorter for larger values of the volume fraction h, again using the elastic parameters of NiAl.Even though for varying values of h the underlying finite element mesh is updated to avoid over stretching of elements the general profile of the dependence of ' on h is clearly visible.The needles become longer as h decreases, but-at variance with what has been discussed above in the case of d-we do not expect the length of the needle to diverge as h !0.
Acknowledgements This work was partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through project 211504053/SFB1060 and project 441211072/SPP2256.J. Verhu ¨lsdonk was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2047/1 -390685813 and -EXC2151 -390873048.We would like to thank Nora Lu ¨then for helpful discussions.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the  source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

2tÞ 2
Cb b b b: ð3:12Þ Since C defines a quadratic form which is strictly positive on nonzero symmetric matrices, Cb b b b [ 0.
xÞ Dw½að xÞ À Á À1 det Dw½að xÞd x; ð4:2Þ where x are computational reference coordinates representing the change of variables x ¼ w½að xÞ and D is the Jacobian in these coordinates, and the associated elastic deformation / ¼ / w½a in reference coordinates.For the differentiation we have applied the chain rule D /ð xÞ ¼ D/ðw½að xÞÞ Dw½að xÞ.
. The asymmetry reflected by the shift D of the needle tip gets larger in absolute value for increasing d.Influence of material constants The energy density depends on the three material constants C 11 , C 12 and C 44 .Neglecting the overall scaling there are two remaining degrees of freedom.In what follows, we investigate the influence of the anisotropy ratio A ¼ 2C 44 =ðC 11 À C 12 Þ and the Poisson ratio m ¼ C 12 =C 11 via modulation of the parameters A and m starting from the material constants C 11 ¼ 115:5, C 12 ¼ 45:5 and C 44 ¼ 110 of NiAl.Here, we choose fixed q 1

Fig. 4 Fig. 5
Fig. 4 Top: optimal needle geometries for NiAl, CuAlNi and YBaCuO (from top to bottom, scaled anisotropically with a factor of 0.25 in x-direction) are plotted in physical reference configuration.The length of the needle tip is indicated by triangle pointers and a dotted line.The right boundary of the computational domain is marked by a small white gap with a horizontal continuation on the right.Bottom: Fixing the tip of all needles at the same position and scaling them all to the same length, the different needle contours (NiAl solid red, CuAlNi-dashed green, YBaCuO dotted blue) are compared

Fig. 6 Fig. 7
Fig. 6 Optimal needle length ' as a function of the anisotropy parameter A (left) and of the Poisson ratio m (right) Roughly speaking, for functions u with o 2 u 2 fÀ1 þ h; hg almost everywhere, this term can be computed by integrating Nðx 1 Þ over (0, L), where Nðx 1 Þ denotes the number of interfaces on the slice fx 1 g Â ð0; HÞ (the number of times that o 2 u flips from h to

Table 1
Material parameters in three dimensions and after plane-stress reduction (rounded) Sketch of the shape-optimization problem in a periodic cell.The thick white curves are the boundaries of the needle which are optimized, together with the length ' and the shift D of the needleW A ðQR T FRÞ ¼ W A ðFÞ whenever Q 2 SOð2Þ; F 2 R 2Â2 ; R 2 P AOð2Þis the point group of the austenite.In our case P A is the cubic group.We require the growth condition We refer to[19,Table1] for details.The elastic constants are in GPa.The experimental data are from [32, Table 4], [33, Table I], [34, Section 3.1], [35, Table 7.1] Fig. 2 A ðFÞ ¼ gðF; det FÞ for all F 2 R 2Â2 ð3:5Þ [30,kl ; ð3:10Þ see for example[29, Eq. (40)] or[30, Section 2].However, these expressions are not rank-one convex, and therefore not polyconvex.To see this, we fix any unit vector b and consider the rank-one line

Table 2
Choices of admissible parameters entering the energy (3.23) for the materials of interest here [19,coefficients b are in GPa.Table from[19, Table 2] In detail, fixing C 11 and C 12 we set C 44 ¼ 0:5AðC 11 À C 12 Þ for varying A. For varying m and fixed A, we hold C 11 fixed and set C 12 ¼ mC 11 and