Comptes Rendus
Mathématique

Sergio Conti, Martin Lenz, Nora Lüthen, Martin Rumpf and Barbara
Zwicknagl

**Geometry of martensite needles in shape memory alloys**

Volume 358, issue 9-10 (2020), p. 1047-1057

Published online: 5 January 2021

https://doi.org/10.5802/crmath.120

<figure><img src="image_1.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_2.png" /></figure>

Les Comptes Rendus. Mathématique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1778-3569
Numerical Analysis / *Analyse numérique*

# Geometry of martensite needles in shape
memory alloys

*Géométries des aiguilles de martensite dans les alliages à
mémoire de forme*

Sergio Conti<sup>a, ∗</sup>, Martin Lenz<sup>∗, b</sup>, Nora Lüthen<sup>c</sup>, Martin Rumpf<sup>a, b</sup>
and Barbara Zwicknagl<sup>b, d</sup>

**Abstract.** We study the geometry of needle-shaped domains in shape-memory alloys. Needle-shaped domains are ubiquitously found in martensites around macroscopic interfaces between regions which are laminated in different directions, or close to macroscopic austenite/twinned-martensite interfaces. Their geometry results from the interplay of the local nonconvexity of the effective energy density with long-range (linear) interactions mediated by the elastic strain field, and is up to now poorly understood. We present a two-dimensional shape optimization model based on finite elasticity and discuss its numerical solution. Our results indicate that the tapering profile of the needles can be understood within finite elasticity, but not with linearized elasticity. The resulting tapering and bending reproduce the main features of experimental observations on Ni₆₅Al₃₅.

**Résumé.** Nous étudions la géométrie des domaines en forme d'aiguille dans les alliages à mémoire de forme. Les domaines en forme d'aiguille sont omniprésents dans les martensites près des interfaces macroscopiques entre régions laminées dans des directions différentes, ou près d'interfaces macroscopiques entre austénite et martensites jumelées. Leur géométrie résulte de l'influence relative de la non-convexité locale de la densité d'énergie effective et des interactions à longue portée (linéaires) engendrées par le champ de déformation élastique, et est pour le moment assez mal comprise. Nous présentons un modèle d'optimisation de forme bi-dimensionnel basé sur l'élasticité non-linéaire et étudions son approximation numérique. Nos résultats

<footnote><sup>a</sup> Institute for Applied Mathematics, University of Bonn, Endenicher Allee 60, 53115 Bonn, Germany</footnote>

<footnote><sup>b</sup> Institute for Numerical Simulation, University of Bonn, Endenicher Allee 60, 53115 Bonn, Germany</footnote>

<footnote><sup>c</sup> Chair of Risk, Safety and Uncertainty Quantification, ETH Zürich, Stefano-Francini-Platz 5, 8093 Zürich, Switzerland</footnote>

<footnote><sup>d</sup> Humboldt-Universität zu Berlin, Departement of Mathematics, Unter den Linden 6,
10099 Berlin, Germany</footnote>

<footnote>E-mails: sergio.conti@uni-bonn.de, martin.lenz@ins.uni-bonn.de, luethen@ibk.baug.ethz.ch, martin.rumpf@ins.uni-bonn.de, barbara.zwicknagl@math.hu-berlin.de</footnote>

<footnote>* Corresponding author.</footnote>
montrent que le profil effilé des aiguilles peut être expliqué dans le cadre de l'élasticité non-linéaire, mais pas dans le cadre linéarisé. L'incéclessment et la flexion qui en résultent reproduisent les caractéristiques principales observées expérimentalement sur le Ni₆₅Al₃₅.

**Funding.** This work was partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via Project 211504053 – SFB 1060.

Manuscript received 5th November 2019, revised 16th September 2020, accepted 22nd September 2020.

# 1. Introduction

Shape-memory alloys couple complex macroscopic material behavior with specific microstructures [2, 18]. Low-energy interfaces between different variants of martensite, or between martensite and austenite, are possible only under specific conditions, and only with a few special orientations. Indeed, in the martensitic phase one often observes finely twinned laminates with a crystallographically-determined orientation. Most experimental situations require deviations from those orientations and lead to complex structures, which involve curved or complex interfaces.

One typical situation is the appearance of so-called “needle-like” domains [4, 18, 19], see Figure 1 for an illustration. They often appear close to macrointerfaces between austenite and finely twinned martensite or between regions where the martensite is twinned in different directions. Needle-like domains are believed to be crucial for the macroscopic energetics and for the hysteresis of the phase transition, as they determine the cost of the transition state [22–24]. They are thin domains of one martensite variant, which taper approaching a tip. The bending of the needle was related to its tapering [4], but the tapering mechanism is up to now not understood. In particular it is not clear what determines the length scale of the needles.

A related (and competing) microstructure is the so-called branching pattern, first studied in this context by Kohn and Müller [13], which can be seen as a laminate which refines close to the macrointerface. In some models branching extends down to scale zero, leading to asymptotically self-similar deformation patterns [6]. If interfacial energy penalizes interfaces with all orientations, however, branching is expected to stop at some point, and, close to the interface, to be replaced by a single-scale interpolation. The latter might correspond to needles, see for example [7]. The detailed geometry of branching patterns was studied in [8]. The role of needles can be further strengthened by evolutionary aspects, as the gain in energy made possible by branching patterns is not necessarily dynamically accessible to the material during the phase transformation. We do not investigate branching patterns further here.

The geometry of needles was studied analytically in a geometrically linear setting in [3, 4, 19]. This permitted in particular to find a simple relation between the bending of the thin needles and their tapering, but the tapering profile itself could not be predicted. Numerical studies with geometrically linear models also did not lead to a stable prediction of the tapering profile [16]. There has been a considerable effort to reproduce the formation of needles numerically, with models based on nonlinear elasticity and many different discretizations schemes. Simulations based on a phase field method [9] have shown that a geometrically nonlinear model is needed for the modeling of polytwinned microstructures. Needle formation, bending, tapering and branching have been observed in studies based on finite elements with a nonconvex energy density [15]. Needles close to a free surface have been studied with an atomistic discretization [17], showing that they have a strong tendency to either retract or to transform into a laminate. An analysis based on three-dimensional FEM with parametrized boundaries was presented, for the case of Cu-Al-Ni,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Experimental observation of needle-shaped domains via transmission electron microscopy in a thin sample of Ni₆₅Al₃₅. Reprinted with permission from [4, Figure 5], available under https://doi.org/10.1103/PhysRevB.64.144105. Copyright (2001) by the American Physical Society. (b) Sketch of the geometry in the reference configuration. The red line outlines the fundamental cell, which is repeated periodically in the vertical direction. (c) Sketch of the geometry in matrix space. The macroscopic deformation gradient V is a weighted average of U¹ and U², and is itself compatible (along a different direction) with QU¹ for some Q ∈ SO(2).</figcaption></figure>

in [21]. This paper focused on energy concentration around the needle tip. The resulting needle geometry was qualitatively similar to the experimental one, but there remained significant differences [21, Figure 8 and 9]. A phase-field approach, also with a geometrically nonlinear model, has been presented in [14] and also shows bending and tapering of needle-like structures (see in particular Figure 3 there). In [5] the stored elastic energy was numerically minimized with respect to variation in the period of the twinned microstructure. The formation of needles was also observed using a Fourier-space discretization of a viscoelastic model [20]. Whereas all these studies agree in showing that models based on nonlinear elasticity can predict the formation of needles which qualitatively resemble the observed ones, the detailed shape, and its dependence on the parameters of the problem, was not explored.

In this work, we use a two-dimensional nonlinear elasticity model to show numerically that the tapering profile can be understood as a consequence of geometric nonlinearity. Our results indicate that the effective length of the needle is proportional to 1/δ, where δ is the order parameter in the bulk. Correspondingly, a geometrically linear version of the model predicts an infinite length of the needles. For definiteness, we focus on needles completely in the martensitic phase of Ni₆₅Al₃₅ and on the TEM observations reported in [4], but our model can be easily generalized to other materials.

## 2. Model

We work with nonlinear elasticity, following [1], in two spatial dimensions. We assume that there are two martensitic variants, with eigenstrains

$$ U^1 := \begin{pmatrix} 1 & \delta \\ 0 & 1 \end{pmatrix} \quad \text{and} \quad U^2 := \begin{pmatrix} 1 & -\delta \\ 0 & 1 \end{pmatrix} \qquad (1) $$

for some $\delta > 0$; in the simulations we use $\delta = 0.1$ which is a typical value for Ni₆₅Al₃₅ (see [4, Section IV/A]). We use, as customary, the austenite as reference configuration, scaled to have the same density as the martensite (so that $U^1$ and $U^2$ have unit determinant). The minimizers of the free-energy density are then the deformation gradients of the form $RU^1$ or $RU^2$, for any $R \in SO(2)$. If there is a continuous interface between a region with deformation gradient $RU^1$ and one with deformation gradient $R'U^2$, then necessarily $RU^1\tau = R'U^2\tau$ for any tangent vector $\tau$.
This equation has two pairs of solutions for $\tau \in S^1$, the first one being $\tau = \pm e_1$ (with $R' = R$); and the second one being $\tau = \pm e_2$ (with $R, R'$ related by $R(\delta, 1) = R'(-\delta, 1)$).

We assume that the left part of the sample is in variant $U^1$, and the right part is in a mixture of variants $U^1$ and $U^2$ as displayed in Figure 1(b).

We denote by $\Omega^1$ (or $\Omega^2$) the part of the domain in which variant $U^1$ (or $U^2$) is used, and minimize the energy

$$E[\phi, \Omega^1, \Omega^2] := \int_{\Omega^1} W^1(D\phi) \, dx + \int_{\Omega^2} W^2(D\phi) \, dx, \quad (2)$$

where $\phi: \Omega^1 \cup \Omega^2 \to \mathbb{R}^2$ is the elastic deformation, and $W^i$ is the free-energy density of variant $i$. A detailed exposition of the physical configuration will be given in Section 4. Here we confine to the discussion of the relevant parameters of the domains $\Omega^1$ and $\Omega^2$ to describe the important geometric characteristics of the needle configuration.

We do not include a surface-energy term, proportional to the length of the interface between $\Omega^1$ and $\Omega^2$. This contribution is typically much smaller than the elastic energy. The surface energy is crucial in determining the length-scale of the $U^1/U^2$ mixture on the right of the interface, which determines the number of interfaces and hence their total length, an effect we do not study here. The bending and tapering of the needle instead have only a minor effect on the total length of the interface, whereas they have a significant effect on the long-range elastic incompatibility. Therefore surface energy is not crucial in this situation, and for simplicity we do not include it in the model.

We impose boundary conditions on the two domains $\Omega^1, \Omega^2$ corresponding to the geometry illustrated in Figure 1(b), and determine both the elastic deformation $\phi$ and the two sets $\Omega^1, \Omega^2$ by minimizing the energy $E$, resulting in a shape-optimization problem. As we do not optimize over the topology, we opt for a reparametrization scheme (discussed below) for the optimization of the two shapes. The domains $\Omega^1, \Omega^2$ are parametrized via low-order polynomials.

We assume that on the far right a periodic mixture of the $U^1$ and $U^2$ variants is present, with periodicity $H$, and volume fractions $\theta$ and $1-\theta$, $\theta \in (0,1)$ (cf. Figure 2). This is geometrically possible, since $U^1 e_1 = U^2 e_1$. The average deformation gradient in this region is then $V := \theta U^1 + (1-\theta) U^2 = \text{Id} + (2\theta - 1) \delta e_1 \otimes e_2$ and is compatible to a different rotation of the $U^1$ variant, $QU^1 e_2^\delta = Ve_2^\delta$, where $e_2^\delta := \frac{(-\delta\theta, 1)}{\sqrt{1+\delta^2\theta^2}}$. The rotation $Q$ can be easily computed in terms of $\delta$ and $\theta$. We are interested in resolving the structure around this macrointerface, which is oriented along $e_2^\delta$, assuming the periodicity is not changed.

We use a generic polyconvex nonlinear elastic energy density with the cubic symmetry of the austenite [11],

$$W(F) := a_1 \operatorname{tr}(F^T F)^2 + a_2 \det(F^T F) - a_3 \log(\det(F^T F)) + a_4 \left( (F^T F)_{11}^2 + (F^T F)_{22}^2 \right), \quad (3)$$

and extend it to the two martensite variants by setting $W^1(F) := W(F(U^1)^{-1})$, $W^2(F) := W(F(U^2)^{-1})$. If $4a_1 + a_2 - a_3 + 2a_4 = 0$ and suitable inequalities hold, then $W$ is minimized on SO(2). The remaining three conditions for $a_1, a_2, a_3, a_4$ depend on the three elasticity constants of a cubic material. We choose $c_{11} = 200 \text{ GPa}, c_{12} = 130 \text{ GPa}, c_{44} = 110 \text{ GPa}$, which are appropriate for the B2 austenite phase of Ni₆₅Al₃₅, see for example [10, Table I]. Then $d_{11} := c_{11} - c_{12}^2/c_{11}$, $d_{12} := c_{12} - c_{12}^2/c_{11}$ and $d_{44} := c_{44}$ are the effective elastic coefficients for a two-dimensional plane stress reduction. The three conditions $\frac{1}{2}d_{11} = 16a_1 + 2a_2 + 2a_3 + 12a_4$, $d_{12} = 16a_1 + 8a_2$, and $2d_{44} = 16a_1 - 4a_2 + 4a_3 + 8a_4$ ensure the Hessian of $W$ coincides with the two-dimensional elasticity tensor given by these elastic coefficients. We solve this system for $a_1, a_2, a_3, a_4$ and obtain $a_1 = 11.56 \text{ GPa}, a_2 = -17.44 \text{ GPa}, a_3 = 10.04 \text{ GPa}, a_4 = -9.38 \text{ GPa}$. We used the plane-stress reduction to obtain the two-dimensional elastic constants $d_{ij}$ from the three-dimensional constants $c_{ij}$ since the experiments in [4] are performed in thin films. Using a plane-strain reduction
would lead to different values for the elastic constants, but not to changes in the model or the
functional form of W in (3). The model studied in this paper is purely two dimensional.

3. Results

We have numerically optimized the total elastic energy $E[\phi, \Omega^1, \Omega^2]$ both in the deformation $\phi$
and in the shape of the domains $\Omega^1$ and $\Omega^2$. We use a finite-element discretization for the de-
formation and minimize it out for any fixed domain geometry. The resulting effective functional,
which depends only on the domains, is then minimized in the shape of the needles, which is
parametrized in terms of the total needle length $L^n$ and two quadratic curves $\gamma^t$ and $\gamma^b$ describ-
ing the top and bottom needle profiles in the reference configuration, respectively. Since we as-
sume periodicity along the interface, only one needle needs to be resolved numerically. For a de-
tailed explanation of the concrete parametrization of the needle geometry and the shape opti-
mization method we refer to Section 4.

The resulting shape of the needles is illustrated in Figure 2(a). The needle tapering naturally
occurs on a length scale which is of the order of 5 to 10 needle spacings, in agreement with
experiment (see Figure 1(a)). This tapering effect was not resolved in the geometrically linear
model in [4] and is discussed in more detail below. The tapering in turn generates a bending of the
needle, as apparent both in the numerical results of Figure 2(a) and in the experimental image in
Figure 1(a). In the reference configuration the bending of the needle is almost absent, and only a
small asymmetry between the two boundaries is apparent (see Figure 1(b) and discussion below).
This confirms that the experimentally observed bending results from the fact that compatibility
across a tapering needle requires a nontrivial rotation, as predicted by the geometrically linear
analysis [4]. In the following we discuss the tapering length and the asymmetry in more detail.

One important result of our computations is that the tapering length can be explained by a
two-dimensional model based on finite elasticity, but not by the geometrically linear version.
This can be qualitatively understood as a consequence of the well-known differences between
linear and nonlinear elasticity in the theory of elastic plates. Indeed, the key mechanism driving
the deformation in the needle region is the need for the region $\Omega^2$ to accommodate two different
orientations: the deformation gradient is close to $U^2$ in the part around the right boundary, and
close to $R_0U^2$ around the vertical interface, where $R_0 \in SO(2)$ is such that $R_0U^2 - QU^1$ is a rank-
one matrix, as illustrated in Figure 1. A detailed analysis shows that the leading-order energy
contribution in finite kinematics is different from the one in the linearized theory. This will be
discussed in a forthcoming publication.

Needle length

Figure 2 shows the optimal needle shape obtained for $\delta = 0.1$ and volume fraction $\theta = 0.25$. The right panel shows an enlargement of the shape of a single needle, which is almost, but not exactly, symmetric in the reference configuration (see below).

Figure 3 shows that the general structure with an almost symmetric tapering in the reference
configuration and a substantial bending in the deformed configuration is not specific to $\delta = 0.1$.
At the same time the effective tapering length $L^n$ (defined precisely via the formulas for $\gamma^t$ and
$\gamma^b$ given in Section 4 below) depends strongly on $\delta$ (defined in (1)). By scaling one easily sees
that necessarily $L^n = Hf(\delta)$ for the present scale-invariant model. The right panel shows that the
dependence is very well approximated by $L^n = cH/\delta$, predicting an infinite tapering length in the
limit $\delta \to 0$. This is precisely the limit in which linearized elasticity applies.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Periodic needle pattern (in the deformed configuration) as obtained by numerical minimization of the functional (2) with $\delta = 0.1$ and $\theta = 0.25$, to be compared with the experimental observation in Figure 1(a). We also show the needle shape in the reference (b) and deformed (c) configuration, scaled by a factor of 4 vertically in order to better illustrate the shape of the needle. The shape in the reference configuration is almost symmetric, but in the deformed configuration a substantial bending appears.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Optimal needle shape for different values of $\delta$ in the (a) reference and (b) deformed configuration, scaled for clarity by a factor of 4 vertically as in Figure 2(b) and (c). The right panel shows the needle length $L^n$ as a function of $\delta$ (the filled circles mark the configurations on the left). The black line shows the slope of $1/\delta$.</figcaption></figure>

## Linearized elasticity

We also considered a linearization of the present model, in which the energy density defined in (3) is replaced by $W_{\text{Lin}}(\epsilon) := \frac{1}{2}d_{11}(\epsilon_{11}^2 + \epsilon_{22}^2) + d_{12}\epsilon_{11}\epsilon_{22} + 2d_{44}\epsilon_{12}^2$ and $W_{\text{Lin}}^{1,2}(\epsilon) := W_{\text{Lin}}(\epsilon \pm \frac{1}{2}\delta(e_1 \otimes e_2 + e_2 \otimes e_1))$, with $\epsilon(\phi) := \frac{1}{2}(D\phi + D\phi^T) - \text{Id}$. Correspondingly, the macrointerface between the average displacement gradient $(2\theta - 1)\delta e_1 \otimes e_2$ and the variant $\delta e_1 \otimes e_2$ is oriented along $e_2$, which is the leading-order linearization of the vector $e_2^{\delta\theta}$ introduced above. Therefore the appropriate boundary condition in the linear context is $\phi(x_1, H) = \phi(x_1, 0) + He_2 + (2\theta - 1)\delta He_1$. The numerics (cf. Figure 4) show that the stored elastic energy, after optimizing with respect to the remaining degrees of freedom, decays for increasing $L^n$. In fact, the elastic energy appears to be proportional to $1/L^n$, so that the optimal value of $L^n$ seems to be infinity.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Linearized elastic energy of the optimal needle shape for different values of $L^n$, demonstrating that there is no fixed length scale in the linearized case. The black line shows the slope of $1/L^n$.</figcaption></figure>

### Nonsymmetric, curved needle geometry

The optimal needle geometry appears almost, but not exactly, symmetric. In fact, the needle tip is shifted upwards by ≈ 0.003 with respect to the centerline of the corresponding laminate layer in the reference configuration. Closer inspecting the optimized, parametric curves $\gamma^t$ and $\gamma^b$, one observes a tangential alignment of the needle geometry and the corresponding martensitic laminate layer of constant thickness at the right boundary. For our finest mesh computation the slope of $\gamma^t$ is $1.9 \cdot 10^{-4}$ and the slope of $\gamma^b$ is $5.4 \cdot 10^{-5}$ at the point where the needle merges into the laminate (in the reference configuration).

## 4. Numerical scheme for the optimization of the needle geometry

In order to numerically optimize over the geometry we use a change-of-variables technique related to configurational mechanics, see, e.g., [12]. Precisely, we let the physical domain $\Omega$ be the image of a (fixed) computational domain $\hat{\Omega}$ via a map $\psi[\alpha]$ depending on a finite set of design parameters $\alpha$. The physical reference configuration represents the undeformed state of the fundamental cell of the martensitic twin microstructure and is deformed by the actual elastic deformation $\phi$. Computational domain and physical reference configuration are sketched in Figure 5. In the physical reference configuration the top and bottom boundaries of the needle are

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Sketch of the transformation of the fundamental cell from the computational domain to the physical reference configuration, cf. Figure 1.</figcaption></figure>
described by curves $\gamma^t$ and $\gamma^b$, which in turn are parametrized via $\psi[\alpha]$ over two (fixed) piecewise affine curves $\hat{\gamma}^t$ and $\hat{\gamma}^b$ at corresponding positions in the computational domain. The horizontal length $L$ of the fundamental cell is fixed and splits up into the length of left $U^1$ variant $L^l$, which is fixed as well, the length of the needle tip $L^n$ and the length of the laminate component on the right $L^r$ (ideally this is an infinite domain, we make it finite only for computational reasons). The needle length $L^n$ in the physical reference configuration is one of the design parameters and consequently $L^r = L - L^l - L^n$. The overall height for the fundamental cell is denoted by $H$, which splits up into fractions $\theta$ and $(1-\theta)$ corresponding to the two martensitic variants. In our computations we chose $H=1$. The values of $L$ and $L^l$ were chosen depending on the expected range of $L^n$, i.e. $L=14.5$, $L^l=2.5$ for Figure 2 and Figure 6, $L=22.5$, $L^l=2.5$ for Figure 3, and $L=16$, $L^l=2$ for Figure 4.

Another design parameter is the width $\Delta$ of the needle cross section at the end of the needle tip above the $x_1$ axis along the sheared direction $e_2^{\delta\theta}$. In particular the composite component with $L^n < x < L^n + L^r$ is shifted in this direction. Correspondingly, the width below the $x_1$ axis in the physical configuration is $\theta H - \Delta$. The parametrization $\psi[\alpha]$ is composed of three transformations: (1) a linear stretching of the needle tip interval of length $\hat{L}^n$ along the $x_1$ axis by a factor $\frac{L^n}{L^n}$ and a corresponding linear compression of the laminate component to the right, (2) a piecewise affine deformation on slices in the direction of $e_2$ within the needle interval $[0, L^n]$ which maps $\hat{\gamma}^b$ to $\gamma^b$, $\hat{\gamma}^t$ to $\gamma^t$, and $\hat{\gamma}^b + He_2$ to $\gamma^b + He_2$, and (3) a shearing deformation keeping $e_1$ fixed and mapping $e_2$ onto $e_2^{\delta\theta}$. To describe the curved needle geometries we consider simple quadratic curves $\gamma^t = \{(x, a^t x^2 + (\frac{\Delta}{L^n} - a^t L^n)x) : x \in [0, L^n]\}$, $\gamma^b = \{(x, a^b x^2 + (\frac{\Delta - \theta H}{L^n} - a^b L^n)x) : x \in [0, L^n]\}$ where the coefficients $a^t$ and $a^b$ are additional degrees of freedom. Altogether we consider as the vector of design variables $\alpha := (\Delta, L^n, a^b, a^t)$.

Given this geometric configuration we aim at minimizing the shape functional $\alpha \mapsto \mathcal{J}[\alpha]$, where $\mathcal{J}[\alpha] := \mathcal{E}[\alpha, \phi[\alpha]]$ with $\phi[\alpha]$ being the elastic deformation minimizing the stored elastic energy (2), which we assume to be unique. Taking into account the parametrizations $\Omega^{1,2} := \psi[\alpha](\hat{\Omega}^{1,2})$ of the two physical domains $\Omega^1$ and $\Omega^2$, we obtain for the stored elastic energy for a deformation $\hat{\phi} := \phi \circ \psi[\alpha]$

$$ \begin{aligned} \hat{\mathcal{E}}[\alpha, \hat{\phi}] :={}& E[\phi, \Omega^1, \Omega^2] = \int_{\psi[\alpha](\hat{\Omega}^1)} W^1(D\phi) \,dx + \int_{\psi[\alpha](\hat{\Omega}^2)} W^2(D\phi) \,dx \\ ={}& \sum_{m=1,2} \int_{\hat{\Omega}^m} W^m (\hat{D}\hat{\phi}(\hat{x})(\hat{D}\psi[\alpha](\hat{x}))^{-1}) \det \hat{D}\psi[\alpha](\hat{x}) \,d\hat{x}. \end{aligned} $$

Here, we assume that the $\psi[\alpha]$ is bijective and used the chain rule $D\phi(x) = \hat{D}\hat{\phi}[\alpha](\hat{x})(\hat{D}\psi[\alpha](\hat{x}))^{-1}$ where $\hat{D}$ is the differential in reference coordinates $\hat{x}$ with $x = \psi[\alpha](\hat{x})$. Finally, we consider periodic boundary conditions on the top and bottom boundary of our fundamental domain and natural boundary conditions on the left and on the right.

A necessary condition for a reference deformation $\hat{\phi}$ to minimize the elastic energy $\hat{\mathcal{E}}[\alpha, \hat{\phi}]$ for fixed $\alpha$ is given by the Euler-Lagrange equation $0 = \hat{\mathcal{E}}_{,\hat{\phi}}[\alpha, \hat{\phi}]$, where $\hat{\mathcal{E}}_{,\hat{\phi}}$ denotes the Fréchet derivative with respect to $\hat{\phi}$. For the spatial discretization we consider an admissible simplicial finite element mesh for the reference domains $\hat{\Omega}^{1,2}$ with coinciding nodes on the common interface. Denoting by $\hat{V}_h$ the associated space of piecewise affine and continuous functions from $\hat{\Omega}$ to $\mathbb{R}^2$ we consider the conforming finite element discretization of our variational problem with $\hat{\phi}$ in $\hat{V}_h$. The solution $\hat{\phi}[\alpha]$ is computed using Newton's scheme for the state equation $\mathcal{E}_{,\hat{\phi}}[\alpha, \hat{\phi}] = 0$ where a small number of iterations (up to 11) are sufficient to obtain a residual in the $L^2$ norm of $10^{-6}$. This scheme can easily be extended via a suitable projection to incorporate translational invariance $\int_{\hat{\Omega}} (\phi - \text{Id}) \,dx = 0$.

To minimize the cost functional $\mathcal{J}[\alpha] = \hat{\mathcal{E}}[\alpha, \hat{\phi}[\alpha]]$ we consider a cascadic nonlinear conjugate gradient scheme. For the linesearch with step size control in the descent we proceed as follows.
Given a descent direction $q$ in the design variable $\alpha$ we define $f(t) := \mathcal{J}(\alpha + tq)$ and let $p$ be a quadratic polynomial solving the interpolation problem $p(0) = f(0)$, $p'(0) = f'(0) = \mathcal{J}_{,\alpha}[\alpha](q)$, and $p(\lambda) = f(\lambda)$ for a suitable $\lambda > 0$. The minimum of $p$ is attained for the time step $\lambda^* = -\frac{f'(0)\lambda^2}{2(f(\lambda)-f'(0)\lambda-f(0))}$. This time step is then used as input for Armijo's time step control applied to $t \mapsto f(t)$. We descend until $|\mathcal{J}_{,\alpha}[\alpha]| < 10^{-5}$, with at most 101 descent steps on each level of the cascadic scheme in our numerical experiments.

We have experimentally studied the convergence of the elastic energy and the geometry degrees of freedom $\alpha$. For the shape shown in Figure 2 convergence plots for a spatial grid size $h = 2^{-n}$ for $n = 2, \dots, 7$ are presented in Figure 6. As a robustness check we have also considered different types of finite element meshes, i.e. meshes with diagonal faces in the $(1, -1)$ (⊠) and $(1, 1)$ (⊠) directions, respectively. All computations in Section 3 have been performed on a mesh with $h = 2^{-6}$ with a ⊠ grid, oriented as in Figure 6. The results have also been verified to be robust with respect to the choice of initial conditions for the optimization parameters.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: Coarsest finite element mesh, $h = 2^{-2}$, ⊠ type. Bottom: Convergence plots for the total energy and the four optimization parameters (the two curvatures $-a^t$ and $a^b$ in the second plot). In each plot the results for ⊠ (blue x markers, corresponding to the mesh depicted on top) and ⊠ (green circle markers) grids are shown for successive refinement by quadrisection. A symmetric needle would have $\Delta = 0.125$.</figcaption></figure>
## Légendes des figures

Figure 1. (a) Observation expérimentale de domaines en forme d'aiguille par microscopie électronique en transmission dans un échantillon mince de Ni₆₅Al₃₅, réimprimé avec la permission de [4, Figure 5], disponible sous https://doi.org/10.1103/PhysRevB.64.144105, Copyright (2001) American Physical Society. (b) Esquisse de la géométrie dans la configuration de référence. La ligne rouge délimite la cellule fondamentale, qui est répétée périodiquement dans le sens vertical. (c) Esquisse de la géométrie dans l'espace matriciel. Le gradient de déformation macroscopique V est une moyenne pondérée de $U¹$ et $U²$, et est lui-même compatible (dans une direction différente) avec $QU¹$ pour certains $Q \in SO(2)$.

Figure 2. (a) Modèle d'aiguille périodique (dans la configuration déformée) obtenu par minimisation numérique de la fonctionnelle (2) avec $\delta = 0.1$ et $\theta = 0.25$, à comparer avec l'observation expérimentale de Figure 1(a). Nous montrons également la forme de l'aiguille dans la configuration de référence (b) et déformée (c), mis à l'échelle verticalement d'un facteur 4 afin de mieux illustrer la forme de l'aiguille. La forme dans la configuration de référence est presque symétrique, mais dans la configuration déformée, une flexion importante apparaît.

Figure 3. Forme optimale de l'aiguille $L^n$ pour différentes valeurs de $\delta$, dans la configuration (a) de référence et (b) déformée, mis à l'échelle pour la clarté d'un facteur 4 verticalement, comme dans la figure 2(b) et (c). Le panneau de droite montre la longueur de l'aiguille $L^n$ en fonction de $\delta$ (les cercles pleins indiquent les configurations à gauche). La ligne noire indique la pente de $1/\delta$.

Figure 4. Énergie élastique linéarisée de la forme optimale de l'aiguille pour différentes valeurs de $L^n$, démontrant qu'il n'y a pas d'échelle de longueur fixe dans le cas linéarisé. La ligne noire indique la pente de $1/L^n$.

Figure 5. Esquisse de la transformation de la cellule fondamentale du domaine de calcul à la configuration de référence physique, cf. figure 1.

Figure 6. En haut: le maillage par éléments finis le plus grossier, $h = 2^{-2}$, type $\square$. En bas: graphique de convergence pour l'énergie totale et les quatre paramètres d'optimisation (les deux courbures $-a^t$ et $a^b$ dans le deuxième tracé). Dans chaque graphique, les résultats pour les maillages $\square$ (marqueurs x bleus, correspondant au maillage représenté en haut) et $\square$ (marqueurs circulaires verts) s'affichent pour un raffinement successif par quadrisection. Une aiguille symétrique aurait $\Delta = 0.125$.

## References

[1] J. M. Ball, R. D. James, "Fine phase mixtures as minimizers of the energy", *Arch. Ration. Mech. Analysis* **100** (1987), p. 13-52.

[2] K. Bhattacharya, *Microstructure of martensite: Why it forms and how it gives rise to the shape-memory effect*, Oxford Series in Materials Modelling, vol. 2, Oxford University Press, 2003.

[3] P. Boullay, D. Schryvers, J. M. Ball, "Nano-structures at martensite macrotwin interfaces in Ni₆₅Al₃₅", *Acta Materialia* **51** (2003), no. 5, p. 1421-1436.

[4] P. Boullay, D. Schryvers, R. V. Kohn, "Bending martensite needles in Ni₆₅Al₃₅ investigated by two-dimensional elasticity and high-resolution transmission electron microscopy", *Phys. Rev. B* **64** (2001), article ID 144105.

[5] E. Bronstein, E. Faran, D. Shilo, "Analysis of austenite-martensite phase boundary and twinned microstructure in shape memory alloys: The role of twinning disconnections", *Acta Materialia* **164** (2019), p. 520-529.

[6] S. Conti, "Branched microstructures: scaling and asymptotic self-similarity", *Commun. Pure Appl. Math.* **53** (2000), no. 11, p. 1448-1474.

[7] S. Conti, B. Zwicknagl, "Low volume-fraction microstructures in martensites and crystal plasticity", *Math. Models Methods App. Sci.* **26** (2016), p. 1319-1355.

[8] P. Dondl, B. Heeren, M. Rumpf, "Optimization of the branching pattern in coherent phase transitions", *C. R. Math. Acad. Sci. Paris* **354** (2016), no. 6, p. 639-644.
[9] A. Finel, Y. Le Bouar, A. Gaubert, U. Salman, "Phase field methods: Microstructures, mechanical properties and complexity", *C. R. Physique* **11** (2010), p. 245-256.

[10] X. Huang, I. I. Naumov, K. M. Rabe, "Phonon anomalies and elastic constants of cubic NiAl from first principles", *Phys. Rev. B* **70** (2004), p. 064301.

[11] N. Kambouchev, J. Fernandez, R. Radovitzky, "A polyconvex model for materials with cubic symmetry", *Modelling and Simulation in Materials Science and Engineering* **15** (2007), no. 5, p. 451-467.

[12] R. Kienzler, G. A. Maugin (eds.), *Configurational mechanics of materials*, CISM Courses and Lectures, vol. 427, Springer, 2001.

[13] R. V. Kohn, S. Müller, "Surface energy and microstructure in coherent phase transitions", *Commun. Pure Appl. Math.* **47** (1994), p. 405-435.

[14] V. I. Levitas, A. M. Roy, D. L. Preston, "Multiple twinning and variant-variant transformations in martensite: Phase-field approach", *Phys. Rev. B* **88** (2013), p. 054113.

[15] B. Li, M. Luskin, "Theory and computation for the microstructure near the interface between twinned layers and a pure variant of martensite", *Materials Science and Engineering. Ser. A* **273-275** (1999), p. 237-240.

[16] B. K. Muite, O. U. Salman, "Computations of geometrically linear and nonlinear Ginzburg-Landau models for martensitic pattern formation", *ESOMAT* (2009), article ID 03008.

[17] J. Novak, U. Bismayer, E. K. H. Salje, "Simulated equilibrium shapes of ferroelastic needle domains", *Journal of Physics: Condensed Matter* **14** (2002), no. 3, p. 657-664.

[18] E. K. H. Salje, *Phase transitions in ferroelastic and co-elastic crystals*, Cambridge University Press, 1990.

[19] E. K. H. Salje, Y. Ishibashi, "Mesoscopic structures in ferroelastic crystals: needle twins and right-angled domains", *Journal of Physics: Condensed Matter* **8** (1996), no. 44, p. 8477-8495.

[20] O. U. Salman, B. Muite, A. Finel, "Origin of stabilization of macrotwin boundaries in martensites", *The European Physical Journal B* **92** (2019), no. 1, p. 20.

[21] H. Seiner, O. Glatz, M. Landa, "A finite element analysis of the morphology of the twinned-to-detwinned interface observed in microstructure of the Cu-Al-Ni shape memory alloy", *International Journal of Solids and Structures* **48** (2011), no. 13, p. 2005-2014.

[22] Z. Zhang, R. D. James, "A way to search for multiferroic materials with unlikely combinations of physical properties", in *Magnetism and structure in functional materials*, Springer Series in Materials Science, vol. 79, Springer, 2005.

[23] Z. Zhang, R. D. James, S. Müller, "Energy barriers and hysteresis in martensitic phase transformations", *Acta Materialia* **57** (2009), p. 4332-4352.

[24] B. Zwicknagl, "Microstructures in low-hysteresis shape memory alloys: scaling regimes and optimal needle shapes", *Arch. Ration. Mech. Anal.* **213** (2014), no. 2, p. 355-421.