# Towards a sharper phase-field method: a hybrid diffuse–semisharp approach
for microstructure evolution problems†

Jędrzej Dobrzański<EMAIL_ADDRESS>Stanisław Stupkiewicz
<EMAIL_ADDRESS>Institute of Fundamental Technological Research, Polish
Academy of Sciences,
Pawińskiego 5B, 02-106 Warsaw, Poland.

###### Abstract

A new approach is developed for computational modelling of microstructure
evolution problems. The approach combines the phase-field method with the
recently-developed laminated element technique (LET) which is a simple and
efficient method to model weak discontinuities using non-conforming finite-
element meshes. The essence of LET is in treating the elements that are cut by
an interface as simple laminates of the two phases, and this idea is here
extended to propagating interfaces so that the volume fraction of the phases
and the lamination orientation vary accordingly. In the proposed LET-PF
approach, the phase-field variable (order parameter), which is governed by an
evolution equation of the Ginzburg–Landau type, plays the role of a level-set
function that implicitly defines the position of the (sharp) interface. The
mechanical equilibrium subproblem is then solved using the semisharp LET
technique. Performance of LET-PF is illustrated by numerical examples. In
particular, it is shown that, for the problems studied, LET-PF exhibits higher
accuracy than the conventional phase-field method so that, for instance,
qualitatively correct results can be obtained using a significantly coarser
mesh, and thus at a lower computational cost. 22footnotetext: Published in
_Comput. Methods Appl. Mech. Engrg._ , 423, 116841, 2024, doi:
10.1016/j.cma.2024.116841

###### keywords:

microstructure evolution , interfaces , laminate , phase-field method , finite
element method

††journal: Computer Methods in Applied Mechanics and Engineering

##  1 Introduction

The phase-field method is a well-established, versatile and powerful approach
to computational modelling of microstructure evolution problems [1, 2, 3, 4,
5]. Microstructure evolution is inherently associated with propagation of
interfaces that separate individual phases, and this constitutes the main
challenge in developing computational approaches for the corresponding
problems. The essence of the phase-field method is in treating the interfaces
as diffuse, rather than sharp, which is achieved by introducing a phase-field
variable (so-called order parameter) that differentiates the phases and varies
in a continuous manner within the diffuse interfaces. Accordingly, propagation
of interfaces can be simulated on a fixed finite-element mesh (this work
focuses on the finite element method). Tracking of interfaces is thus avoided,
which would otherwise require, for instance, adaptive remeshing or enhancement
techniques such as X-FEM.

The phase-field method is highly versatile because the diffuse-interface
framework can be combined with virtually any physics and, in fact, has been
used in numerous contexts, for instance, solidification [6], solid-state
transformations [1, 7, 8], including martensitic transformation [9, 10, 11,
12] and twinning [13, 14, 15], fracture [16, 17], corrosion [18], and many
more.

Another important feature of the phase-field method is that the order
parameter carries information about the energy of the interfaces. The
interfacial energy is naturally included in the energy balance and thus
contributes to the driving force for interface propagation, which is
particularly important at smaller scales at which the interfacial energy
effects are more pronounced. It is also the interfacial energy term in the
free energy function that sets the thickness of the diffuse interfaces through
the term involving the gradient of the order parameter. Moreover, considering
that the contribution of interfaces to the total energy of a system is size-
dependent, the phase-field method is able to model the size effects [19, 20,
21].

However, the beneficial features of the phase-field method, as discussed
above, come at a price of a high computational cost. This is because a
sufficiently fine computational grid (mesh) is needed to correctly represent
the profile of the order parameter across the diffuse interface and to
correctly represent the corresponding interfacial energy. This sets a severe
constraint on the maximum allowable element size, and the associated high
computational cost limits the maximum physical dimension of the domain that
can be considered in the computations.

With the aim to overcome the limitations discussed above, the so-called ‘sharp
phase-field method’ has been developed by Finel et al. [22], see also [23, 24,
25]. The method allows the computational grid to be larger than the
theoretical interface thickness, hence significantly coarser meshes can be
used compared to the conventional phase-field method, and thus larger physical
domains can be effectively simulated. However, the method relies on the notion
of the discrete gradient, which is limited to the finite difference method and
does not generalize to other discretization techniques, in particular, to the
finite-element method.

In this work, we address the mentioned limitations of the phase-field method
from a different perspective. Focusing on microstructure evolution problems
involving elastic interactions and on the finite element method for spatial
discretization, we propose a hybrid diffuse–semisharp approach. The
microstructure, including diffuse interfaces and their evolution, is described
in a phase-field-like fashion using a continuous order parameter, and the
order parameter plays the role of a level-set function that implicitly defines
the position of the sharp interfaces. Considering that the position of the
interfaces is arbitrary (and changes in a continuous manner), a fixed finite-
element mesh does not conform to those sharp interfaces. To solve the
corresponding mechanical equilibrium subproblem, we employ a recently-
developed semisharp approach, namely the laminated element technique (LET)
[26]. In this approach, the elements that are cut by an interface are treated
as laminates of the two involved phases (hence the name of the method), so
that the interface-affected elements create a thin (one element thick) layer
of elements between the two phases, and the method can thus be classified as
semisharp. The method is inspired by the composite voxel technique that has
been developed for FFT-based solvers in computational homogenization [27, 28,
29, 30, 31]. In the present context, it is an important feature of LET that
the response is a continuous function of the position and orientation of the
interface [26], hence the method is suitable for problems involving moving
interfaces.

In the proposed LET-PF method, the two phases are mixed only within a thin
layer of laminated elements (of the thickness of one element) along the sharp
interface, while in the conventional phase-field method the phases are mixed
within the entire diffuse interface that spans several elements (at least 3–4,
say). Note that the mixing based on laminated microstructures, thus in a sense
similar to LET-PF, is employed in some formulations of the phase-field method
[32, 33, 34, 35], however, consistent with the phase-field framework, it is
applied within the entire diffuse interface.

The proposed approach bears some similarity to the level-set method [36]
combined with X-FEM [37, 38, 39]. That approach appears attractive because
X-FEM is a sharp-interface technique so that the interfaces can be modelled
with a high accuracy using a non-matching mesh. However, compared to the
phase-field method, the related applications for microstructure evolution
problems are much more scarce, see e.g. [40, 41, 42]. The LET approach used in
LET-PF is not as accurate as X-FEM [26], but it is much simpler, in
particular, its implementation is carried out solely at the element level and
does not require intervention into the structure of the finite-element code.
Another pronounced difference between the two approaches is in the way the
interfacial energy is treated. The level-set function in the level-set method
is a signed distance function and its role is only to (implicitly) define the
position of the sharp interface. The interfacial energy (density per unit
area) is then assigned directly to the sharp interface. In LET-PF, a diffuse-
interface description of the interfacial energy is adopted as in the
conventional phase-field method so that the interfacial energy contribution to
the total energy is a bulk-like energy (density per unit volume) distributed
in the vicinity of the interfaces.

In the following, the presentation is limited to the small-strain framework,
however, extension to finite deformations is immediate. The implementation and
the examples are performed in 2D, and, again, extension to 3D is immediate. On
the other hand, the proposed LET-PF approach is limited to the case of two
phases, but generalization to a multi-phase framework is not immediate, see
Remark 2 in Section 3.1.

The remainder of the paper is organized as follows. In Section 2, the
reference sharp-interface microstructure evolution problem is formulated and
its conventional phase-field counterpart is described, including the finite-
element treatment of the latter. The new hybrid diffuse–semisharp approach
(LET-PF) is introduced in Section 3. Numerical examples illustrating the
performance of LET-PF are reported in Section 4.

##  2 Microstructure evolution problem

This section briefly summarizes the classical formulation of the
microstructure evolution problem in the sharp-interface framework (Section
2.1) and in the diffuse-interface framework of the phase-field method (Section
2.2). Also, the finite-element treatment of the simple phase-field model
introduced in Section 2.2 is discussed in Section 2.3.

###  2.1 Reference sharp-interface problem

Let us consider domain $\Omega$ with two subdomains $\Omega_{1}$ and
$\Omega_{2}$ separated by an interface $\Gamma$. Each subdomain $\Omega_{i}$
($i=1,2$) is occupied by an elastic phase characterized by the (isothermal)
Helmholtz free energy density $\psi_{i}$ of the form

$\psi_{i}(\bm{\varepsilon})=\psi^{0}_{i}+\frac{1}{2}(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}_{i}):\mathbb{L}_{i}:(\bm{\varepsilon}-\bm{\varepsilon}^{\rm t}_{i}),$ (1)

so that the local constitutive equation reads

$\bm{\sigma}=\frac{\partial\psi_{i}}{\partial\bm{\varepsilon}}=\mathbb{L}_{i}:(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}_{i}).$ (2)

Here, $\bm{\sigma}$ is the Cauchy stress tensor, $\bm{\varepsilon}=\nabla_{\rm
s}\bm{u}=\frac{1}{2}(\nabla\bm{u}+(\nabla\bm{u})^{\rm T})$ is the strain
tensor, $\bm{\varepsilon}^{\rm t}_{i}$ is the transformation strain (known and
constant within each phase), $\psi^{0}_{i}$ is the free energy in the stress-
free state (so-called chemical energy), and $\mathbb{L}_{i}$ is a positive-
definite fourth-order elastic moduli tensor possessing the usual symmetries.

The displacement field is assumed continuous at $\Gamma$, and $\Gamma$ is
assumed to be smooth with $\bm{n}$ denoting the unit normal pointing from
phase 1 to phase 2. Then, the classical compatibility conditions at $\Gamma$
read

$[\\![\bm{\varepsilon}]\\!]=\frac{1}{2}(\bm{c}\otimes\bm{n}+\bm{n}\otimes\bm{c}),\qquad[\\![\bm{\sigma}]\\!]\cdot\bm{n}=\bm{0},$
(3)

where
$[\\![\mbox{\footnotesize$\square$}]\\!]=(\mbox{\footnotesize$\square$})_{2}-(\mbox{\footnotesize$\square$})_{1}$
denotes the jump at the interface and $\bm{c}$ is a vector. The compatibility
condition (3)1 is a consequence of the continuity of the displacement at
$\Gamma$. On the other hand, mechanical equilibrium implies that Eq. (3)2
holds, in addition to the usual equilibrium equation,

$\operatorname{div}\bm{\sigma}=\bm{0}\quad\mbox{in}\;\,\Omega\setminus\Gamma,$
(4)

supplemented by the boundary conditions,

$\bm{u}=\bm{u}^{\ast}\quad\mbox{on}\;\,\partial\Omega_{\bm{u}},\qquad\bm{\sigma}\cdot\bm{\nu}=\bm{t}^{\ast}\quad\mbox{on}\;\,\partial\Omega_{\bm{t}},$
(5)

where $\bm{u}^{\ast}$ and $\bm{t}^{\ast}$ denote the displacement and the
surface traction prescribed on $\partial\Omega_{\bm{u}}$ and
$\partial\Omega_{\bm{t}}$, respectively, and $\bm{\nu}$ denotes the unit outer
normal to $\Omega$.

Consider now $\Gamma$ to be a phase transformation front, so that $\Gamma$ is
a moving interface, and assume the following (viscous) kinetic law for the
interface [43]

$\hat{v}_{n}=\hat{m}\hat{f},\qquad\hat{f}=[\\![\psi]\\!]-\tilde{\bm{\sigma}}:[\\![\bm{\varepsilon}]\\!]+\gamma\kappa,$
(6)

where $\hat{v}_{n}$ is the interface speed in the direction $\bm{n}$,
$\hat{m}$ is the interface mobility, $\hat{f}$ is the thermodynamic driving
force, $\gamma$ is the interfacial energy (energy density per unit area of the
interface), and $\kappa$ is the total curvature (i.e., twice the mean
curvature) of $\Gamma$. Moreover,
$\tilde{\bm{\sigma}}=\frac{1}{2}(\bm{\sigma}_{1}+\bm{\sigma}_{2})$ denotes the
average stress at the interface; however, in view of the compatibility
conditions (3), we have
$\tilde{\bm{\sigma}}:[\\![\bm{\varepsilon}]\\!]=\bm{\sigma}_{1}:[\\![\bm{\varepsilon}]\\!]=\bm{\sigma}_{2}:[\\![\bm{\varepsilon}]\\!]$,
hence $\tilde{\bm{\sigma}}$ in Eq. (6)2 can be replaced by $\bm{\sigma}_{1}$
or $\bm{\sigma}_{2}$. Quantities referring to the sharp interface are here
denoted by a symbol with a hat.

Direct computational (e.g., finite-element) treatment of the microstructure
evolution problem within the sharp-interface framework is not straightforward
because it would require remeshing to ensure that the mesh matches the
evolving microstructure at each time instant (and possibly at each iteration),
see [44] for a recent development towards this goal. A feasible approach to
model moving interfaces is to adopt a diffuse-interface approximation
independent of the spatial discretization, which is discussed next.

###  2.2 Diffuse-interface framework: the phase-field method

The sharp-interface problem discussed in the preceding section can be
approximated by adopting a diffuse-interface description using the
conventional phase-field method [1, 2, 10, 13, 45]. To this end, a continuous
order parameter $\phi$, $0\leq\phi\leq 1$, is introduced such that $\phi=0$
corresponds to phase 1 and $\phi=1$ to phase 2, while the intermediate values
correspond to the diffuse interface. Accordingly, the Helmholtz free energy
$\psi$, characterizing both the individual phases and the diffuse interface,
is assumed to be a function of the strain $\bm{\varepsilon}$ and of the order
parameter $\phi$ and its spatial gradient $\nabla\phi$, and comprises the bulk
and interfacial energy contributions, $\psi_{\rm bulk}$ and $\psi_{\rm int}$,
respectively,

$\psi(\bm{\varepsilon},\phi,\nabla\phi)=\psi_{\rm
bulk}(\bm{\varepsilon},\phi)+\psi_{\rm int}(\phi,\nabla\phi).$ (7)

The bulk energy contribution $\psi_{\rm bulk}$ is specified as

$\psi_{\rm
bulk}(\bm{\varepsilon},\phi)=\psi^{0}(\phi)+\frac{1}{2}(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}(\phi)):\mathbb{L}(\phi):(\bm{\varepsilon}-\bm{\varepsilon}^{\rm t}(\phi)),$
(8)

where

$\psi^{0}(\phi)=(1-h(\phi))\psi^{0}_{1}+h(\phi)\psi^{0}_{2},\quad\bm{\varepsilon}^{\rm
t}(\phi)=(1-h(\phi))\bm{\varepsilon}^{\rm t}_{1}+h(\phi)\bm{\varepsilon}^{\rm
t}_{2},\quad\mathbb{L}(\phi)=(1-h(\phi))\mathbb{L}_{1}+h(\phi)\mathbb{L}_{2}.$
(9)

Here, the chemical energy $\psi^{0}$, the transformation strain
$\bm{\varepsilon}^{\rm t}$ and the elastic stiffness tensor $\mathbb{L}$ are
interpolated between the respective parameters of the pure phases and depend
on the order parameter $\phi$ through the weighting function $h(\phi)$ such
that $h(0)=0$, $h(1)=1$, and $h^{\prime}(0)=h^{\prime}(1)=0$. We adopt here
$h(\phi)=3\phi^{2}-2\phi^{3}$, i.e. the popular 2-3-4 polynomial (with
$\alpha=3$), see [10, 2].

The interfacial energy contribution $\psi_{\rm int}$ is adopted in the form of
the so-called double-well potential [2],

$\psi_{\rm
int}(\phi,\nabla\phi)=\frac{6\gamma}{\ell}\left(\phi^{2}(1-\phi)^{2}+\frac{\ell^{2}}{4}\nabla\phi\cdot\nabla\phi\right),$
(10)

where $\gamma$ is the interfacial energy, see Eq. (6)2, and $\ell$ is the
interface thickness parameter.

Note that, when $\phi=0$ , the bulk energy $\psi_{\rm bulk}$, Eq. (8), reduces
to $\psi_{1}$, i.e. to the free energy specified by Eq. (1) for the pure phase
1, and likewise for the free energy $\psi_{2}$ of the pure phase 2 when
$\phi=1$.

In the phase-field method, the interfaces are treated as diffuse hence the
discontinuities at the interfaces are smeared out, and the compatibility
conditions (3) need not be considered. Since the strains and stresses are
continuous, the mechanical equilibrium equation holds in the entire domain
$\Omega$, i.e. also within the diffuse interface,

$\operatorname{div}\bm{\sigma}=\bm{0}\qquad\mbox{in}\;\,\Omega,$ (11)

with
$\bm{\sigma}=\partial\psi/\partial\bm{\varepsilon}=\mathbb{L}(\phi):(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}(\phi))$, subject to the boundary conditions (5).

Evolution of the microstructure is governed by the time-dependent
Ginzburg–Landau equation [46, 1],

$\dot{\phi}=-m\frac{\delta\Psi}{\delta\phi},\qquad\frac{\delta\Psi}{\delta\phi}=\frac{\partial\psi}{\partial\phi}-\nabla\cdot\frac{\partial\psi}{\partial\nabla\phi},$
(12)

where $\delta\Psi/\delta\phi$ is the variational derivative of the total free
energy functional $\Psi=\int_{\Omega}\psi{\rm{d}}V$, and $m$ is the mobility
parameter. Homogeneous Neumann boundary conditions are assumed on the entire
boundary,

$\nabla\phi\cdot\bm{\nu}=0\quad\mbox{on}\;\,\partial\Omega.$ (13)

It can be shown that the equilibrium profile of the diffuse interface has the
following form [45],

$\phi(\xi)=\frac{1}{2}\tanh\left(\frac{\xi-\xi_{0}}{\ell}\right)+\frac{1}{2},$
(14)

which can be obtained by minimizing the total interfacial energy $\Psi_{\rm
int}=\int_{\Omega}\psi_{\rm int}{\rm{d}}V$ of a planar interface located at
$\xi=\xi_{0}$ (i.e. $\phi(\xi_{0})=\frac{1}{2}$), where $\xi$ denotes the
coordinate in the direction normal to the interface. This profile is only an
approximation of the actual profile when the interface is non-planar or when
the interface propagates due to a mechanical driving force that is not
constant within the interface. In fact, special conditions must be met so that
a travelling wave solution is obtained with the interface profile independent
of the velocity [2].

Adopting the equilibrium profile (14) as an approximation of the actual
profile, the mobility parameter $m$ in the evolution equation (12) can be
related to the mobility parameter $\hat{m}$ in the sharp-interface
description, cf. Eq. (6)1. By equating the energy dissipated at the
propagating sharp interface to that dissipated within the diffuse interface
propagating with the same velocity, thus
$\hat{f}\hat{v}_{n}=\hat{v}_{n}^{2}/\hat{m}=\int_{-\infty}^{+\infty}(\dot{\phi}^{2}/m){\rm{d}}\xi$,
where $\xi_{0}=\hat{v}_{n}t$, the following relationship is found,

$\hat{m}=3m\ell.$ (15)

It follows that when the interface parameter $\ell$ is varied, then the
mobility parameter $m$ must be scaled according to Eq. (15) so that the
effective mobility $\hat{m}$ of the interface is not affected (to the first
order, in view of the approximate profile of the diffuse interface).

###  2.3 Finite-element treatment of the phase-field model

In order to arrive at the finite-element formulation, the mechanical
equilibrium equation (11) is expressed in a weak form by following the
standard procedure,

$\int_{\Omega}\bm{\sigma}:\delta\bm{\varepsilon}\,{\rm{d}}V-\int_{\partial\Omega_{\bm{t}}}\bm{t}^{\ast}\cdot\delta\bm{u}\,{\rm{d}}S=0\qquad\forall\,\delta\bm{u},$
(16)

where $\delta\bm{\varepsilon}=\nabla_{\rm s}\delta\bm{u}$, and $\delta\bm{u}$
is the virtual displacement (test function) that vanishes,
$\delta\bm{u}=\bm{0}$, on $\partial\Omega_{\bm{u}}$. Likewise, the evolution
equation (12) is expressed in a weak form,

$\int_{\Omega}\Big{[}\Big{(}\frac{\dot{\phi}}{m}+\frac{\partial\psi_{\rm
bulk}}{\partial\phi}+\frac{\partial\psi_{\rm
int}}{\partial\phi}\Big{)}\delta\phi+\frac{\partial\psi_{\rm
int}}{\partial\nabla\phi}\cdot\nabla\delta\phi\Big{]}{\rm{d}}V=0\qquad\forall\,\delta\phi,$
(17)

where $\delta\phi$ is the respective test function and

$\frac{\partial\psi_{\rm
bulk}}{\partial\phi}=\Big{(}\Delta\phi^{0}-\bm{\sigma}:\Delta\bm{\varepsilon}^{\rm
t}+\frac{1}{2}(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}(\phi)):\Delta\mathbb{L}:(\bm{\varepsilon}-\bm{\varepsilon}^{\rm
t}(\phi))\Big{)}h^{\prime}(\phi),$ (18) $\frac{\partial\psi_{\rm
int}}{\partial\phi}=\frac{12\gamma}{\ell}\phi(1-\phi)(1-2\phi),\qquad\frac{\partial\psi_{\rm
int}}{\partial\nabla\phi}=3\gamma\ell\,\nabla\phi,$ (19)

where
$\Delta(\mbox{\footnotesize$\square$})=(\mbox{\footnotesize$\square$})_{2}-(\mbox{\footnotesize$\square$})_{1}$.

The incremental (time-discrete) form of the evolution equation (17) is
obtained by applying the implicit backward-Euler time integration scheme.
Considering a typical time increment $t_{n}\to t_{n+1}=t_{n}+\tau$, the rate
of $\phi$ is thus approximated by the finite difference,

$\dot{\phi}\approx\frac{1}{\tau}(\phi-\phi_{n}),$ (20)

where $\phi_{n}$ is the known value at the previous time step $t_{n}$, and
$\phi=\phi_{n+1}$ is the unknown value at the current time step. Here and
below, to make the notation more compact, the subscript $n+1$ is omitted for
the quantities evaluated at the current time step $t_{n+1}$. Adopting this
notation, the time-discrete problem at $t_{n+1}$ is given directly by the weak
forms (16) and (17) with $\dot{\phi}$ in Eq. 17 replaced by its finite-
difference approximation (20).

Spatial discretization is performed using the finite element method. The
fields of displacement $\bm{u}$ and order parameter $\phi$ are approximated
using the respective basis functions $N_{k}^{(\bm{u})}$ and $N_{k}^{(\phi)}$,

$\bm{u}^{h}=\sum_{k}N_{k}^{(\bm{u})}\bm{u}_{k},\qquad\phi^{h}=\sum_{k}N_{k}^{(\phi)}\phi_{k},$
(21)

where $\bm{u}_{k}$ and $\phi_{k}$ denote the respective nodal values with
$\textbf{{p}}{}_{\bm{u}}=\\{\bm{u}_{k}\\}$ and
$\textbf{{p}}{}_{\phi}=\\{\phi_{k}\\}$ denoting the global vectors. Following
the standard approach, $\bm{u}^{h}$ and $\phi^{h}$ are introduced into the
weak forms (16) and (17), and integration is performed over individual finite
elements $\omega$ of the triangulation ${\cal T}$ of the domain $\Omega$,

$\sum_{\omega\in{\cal
T}}\int_{\omega}\bm{\sigma}:\delta\bm{\varepsilon}^{h}\,{\rm{d}}V-\sum_{\partial\omega\in{\cal
S}_{\bm{t}}}\int_{\partial\omega}\bm{t}^{\ast}\cdot\delta\bm{u}^{h}\,{\rm{d}}S=0\qquad\forall\,\delta\textbf{{p}}{}_{\bm{u}},$
(22) $\sum_{\omega\in{\cal
T}}\int_{\omega}\Big{[}\Big{(}\frac{\phi^{h}-\phi^{h}_{n}}{m\tau}+\frac{\partial\psi_{\rm
bulk}}{\partial\phi}+\frac{\partial\psi_{\rm
int}}{\partial\phi}\Big{)}\delta\phi^{h}+\frac{\partial\psi_{\rm
int}}{\partial\nabla\phi}\cdot\nabla\delta\phi^{h}\Big{]}{\rm{d}}V=0\qquad\forall\,\delta\textbf{{p}}{}_{\phi},$
(23)

where $\delta\bm{\varepsilon}^{h}=\nabla_{\rm s}\delta\bm{u}^{h}$, and ${\cal
S}_{\bm{t}}$ denotes the triangulation of the boundary
$\partial\Omega_{\bm{t}}$ into surface segments $\partial\omega$, consistent
with the triangulation ${\cal T}$ of the bulk.

The discretized weak forms (22) and (23) define the set of coupled nonlinear
equations that can be written in the following residual form,

$\textbf{{R}}_{\bm{u}}(\textbf{{p}}{}_{\bm{u}};\textbf{{p}}{}_{\phi})=\textbf{{0}},\qquad\textbf{{R}}_{\phi}(\textbf{{p}}{}_{\phi};\textbf{{p}}{}_{\bm{u}})=\textbf{{0}}.$
(24)

At each time step of the incremental procedure, these equations are solved in
a monolithic manner, thus $\textbf{{R}}(\textbf{{p}})=\textbf{{0}}$ with
$\textbf{{R}}=\\{\textbf{{R}}_{\bm{u}},\textbf{{R}}_{\phi}\\}$ and
$\textbf{{p}}=\\{\textbf{{p}}{}_{\bm{u}},\textbf{{p}}{}_{\phi}\\}$, using the
Newton method.

##  3 LET-PF: hybrid diffuse–semisharp treatment of propagating interfaces

In Section 3.1, we briefly present the laminated element technique (LET) that
has been recently developed by Dobrzański et al. [26] for an efficient
treatment of material (i.e. non-propagating) interfaces. The presentation is
limited to the simplest case of small-strain framework and linear elasticity,
but the approach is more general [26]. Next, in Section 3.2, LET is combined
with the phase-field method thus leading to a hybrid diffuse–semisharp
treatment of propagating interfaces.

###  3.1 Laminated element technique (LET) for material interfaces

In LET, a non-matching (non-conforming) finite-element discretization is
considered, i.e., the mesh is not aligned with the interface. Each element
that is cut by the interface is treated as a laminate composed of the two
phases involved, and the corresponding laminated microstructure is fully
characterized by the volume fraction of the phases and by the lamination
orientation which are determined as discussed below. At the same time, no
treatment is applied to the remaining elements that belong to one of the
phases. Considering that the interface affects (through lamination) only one
layer of elements along the interface, see the sketch in Fig. 1, LET can be
classified as a _semisharp_ approach, as opposed to the sharp-interface
approach, characteristic for matching-mesh discretizations and for XFEM-type
approaches, and as opposed to the diffuse-interface approach, characteristic
for the phase-field method.

(a) (b)

Figure 1: Laminated element technique (LET): (a) the elements cut by an
interface are treated as simple laminates; (b) the laminated elements are
grouped in subset ${\cal T}_{\rm int}$, the remaining elements belong to
subset ${\cal T}_{1}$ or ${\cal T}_{2}$.

The interface is implicitly defined by a level set function $\phi$ defined
such that $\phi<\frac{1}{2}$ corresponds to phase $1$, $\phi>\frac{1}{2}$
corresponds to phase $2$, and the level set $\phi=\frac{1}{2}$ specifies the
interface. In the discretized setting, the level set function
$\phi\approx\phi^{h}$ is defined in terms of the basis functions
$N_{i}^{(\phi)}$ and respective nodal values $\phi_{i}$, as in Eq. (21)2.
Note, however, that the level set function is here prescribed and does not
evolve since it describes a fixed geometry. This assumption will be relaxed in
Section 3.2.

For each element $\omega\in{\cal T}$, the volume fraction
$\eta^{(\omega)}=\eta^{(\omega)}_{2}=1-\eta^{(\omega)}_{1}$ of phase $2$
within the element is determined according to

$\eta^{(\omega)}=\frac{\sum_{k=1}^{N_{\rm
n}}\langle\phi_{k}^{(\omega)}-\frac{1}{2}\rangle}{\sum_{k=1}^{N_{\rm
n}}|\phi_{k}^{(\omega)}-\frac{1}{2}|},$ (25)

where
$\langle\mbox{\footnotesize$\square$}\rangle=\frac{1}{2}(\mbox{\footnotesize$\square$}+|\mbox{\footnotesize$\square$}|)$
denotes the Macauley brackets, and $N_{n}$ denotes the number of nodes per
element. The above definition is applicable to four-node quadrilateral
elements in 2D and to eight-node hexahedral elements in 3D. For an element cut
by the interface, we have $0<\eta^{(\omega)}<1$. Otherwise, for an element
that is not cut by the interface and thus fully belongs to one of the phases,
we have either $\eta^{(\omega)}=0$ or $\eta^{(\omega)}=1$.

The lamination orientation is specified by the normal to the interface and is
defined as the normalized gradient of the level-set function evaluated at the
element centre $\bm{x}_{0}^{(\omega)}$,

$\bm{n}^{(\omega)}=\frac{\nabla\phi^{h}(\bm{x}_{0}^{(\omega)})}{\|\nabla\phi^{h}(\bm{x}_{0}^{(\omega)})\|}.$
(26)

Now, the set of all elements can be split into three disjoint subsets, ${\cal
T}={\cal T}_{1}\cup{\cal T}_{2}\cup{\cal T}_{\rm int}$, according to the value
of the volume fraction $\eta^{(\omega)}$ in each element,

${\cal T}_{1}=\\{\omega:\;\eta^{(\omega)}=0\\},\qquad{\cal
T}_{2}=\\{\omega:\;\eta^{(\omega)}=1\\},\qquad{\cal T}_{\rm
int}=\\{\omega:\;0<\eta^{(\omega)}<1\\},$ (27)

so that the elements that are cut by the interface are grouped in ${\cal
T}_{\rm int}$, while sets ${\cal T}_{1}$ and ${\cal T}_{2}$ correspond to pure
phases 1 and 2, respectively, see Fig. 1(b).

For a fixed microstructure specified by a prescribed level-set function
$\phi^{h}$, the problem at hand is a linear elasticity problem governed by the
following weak form,

$\sum_{i=1}^{2}\left(\sum_{\omega\in{\cal
T}_{i}}\int_{\omega}\bm{\sigma}_{i}:\delta\bm{\varepsilon}^{h}\,{\rm{d}}V\right)+\sum_{\omega\in{\cal
T}_{\rm
int}}\int_{\omega}\bar{\bm{\sigma}}:\delta\bm{\varepsilon}^{h}\,{\rm{d}}V-\sum_{\partial\omega\in{\cal
S}_{\bm{t}}}\int_{\partial\omega}\bm{t}^{\ast}\cdot\delta\bm{u}^{h}\,{\rm{d}}S=0\qquad\forall\,\delta\textbf{{p}}{}_{\bm{u}},$
(28)

where $\bm{\sigma}_{i}$ denotes the stress in phase $i$ ($i=1,2$),

$\bm{\sigma}_{i}=\bm{\sigma}_{i}(\bm{\varepsilon}^{h})=\left.\frac{\partial\psi_{i}(\bm{\varepsilon})}{\partial\bm{\varepsilon}}\right|_{\bm{\varepsilon}=\bm{\varepsilon}^{h}}=\mathbb{L}_{i}:(\bm{\varepsilon}^{h}-\bm{\varepsilon}_{i}^{\rm
t}),$ (29)

and $\bar{\bm{\sigma}}$ is the overall (macroscopic) stress in the laminated
microstructure in the elements containing the interface,

$\bar{\bm{\sigma}}=\bar{\bm{\sigma}}(\bm{\varepsilon}^{h},\eta^{(\omega)},\bm{n}^{(\omega)})=\left.\frac{\partial\bar{\psi}_{\rm
bulk}(\bm{\varepsilon},\eta^{(\omega)},\bm{n}^{(\omega)})}{\partial\bm{\varepsilon}}\right|_{\bm{\varepsilon}=\bm{\varepsilon}^{h}}=\bar{\mathbb{L}}:(\bm{\varepsilon}^{h}-\bar{\bm{\varepsilon}}^{\rm
t}).$ (30)

Here, $\bar{\psi}_{\rm bulk}=\bar{\psi}_{\rm
bulk}(\bm{\varepsilon},\eta^{(\omega)},\bm{n}^{(\omega)})$ is the overall free
energy of the laminate characterized by $\eta^{(\omega)}$ and
$\bm{n}^{(\omega)}$,
$\bar{\mathbb{L}}=\bar{\mathbb{L}}(\eta^{(\omega)},\bm{n}^{(\omega)})$ is the
overall stiffness tensor, and $\bar{\bm{\varepsilon}}^{\rm
t}=\bar{\bm{\varepsilon}}^{\rm t}(\eta^{(\omega)},\bm{n}^{(\omega)})$ is the
overall inelastic (transformation) strain. Recall that $\eta^{(\omega)}$ and
$\bm{n}^{(\omega)}$ depend on the level-set function $\phi^{h}$, see Eqs. (25)
and (26). The micro-macro transition for the laminated microstructure is
discussed in A.

##### Remark 1

The presentation above is limited to the constitutive framework of linear
elasticity with eigenstrain. However, the approach is more general,
essentially with no restrictions on the constitutive models of the phases. For
instance, an application to finite-strain elastoplasticity is presented in
[26].

##### Remark 2

The LET approach, as presented above, is limited to the case of only two
phases and at most one interface within a single finite element. The latter
assumption can be relaxed, for instance, by introducing a separate level-set
function for each interface, as discussed in Section 4.4 in [26], although
that approach may not directly generalize to the case of propagating
interfaces, discussed in the next section. Consideration of several phases
would, in particular, require an adequate treatment of triple points. Possible
scenarios to treat such cases, for instance, based on higher-rank laminates,
are left for future studies.

###  3.2 Combining the phase-field method with the laminated element technique

The idea of the proposed hybrid diffuse–semisharp approach is to combine the
phase-field method with LET, thus the method will be referred to as LET-PF.
Specifically, diffuse interfaces, their propagation and the corresponding
microstructure evolution are modelled in a manner similar to the phase-field
method. At the same time, the order parameter $\phi$ plays the role of the
level-set function in the semisharp LET-based treatment of the mechanical
subproblem, see Fig. 2.

Figure 2: LET-PF compared to the conventional phase-field method. The profile
of the order parameter is sketched in the upper figures. Shading of the mesh
in the bottom figures indicates the volume fraction of the phases that governs
the bulk energy. In the conventional phase-field method, the transition layer
is diffuse and spans several elements (left), while in the semisharp LET-PF
method it is localized to only one layer of elements (right).

As discussed in Section 2.3, the weak formulation of the discretized
conventional phase-field method is specified by Eqs. (22) and (23). In LET-PF,
the virtual work principle (22) is replaced by that of the LET method, Eq.
(28), which here refers to the current time instant $t=t_{n+1}$ with the
microstructure specified by the order parameter field $\phi^{h}$ at the
current time $t_{n+1}$. The time-discrete evolution of the order parameter
$\phi^{h}$ is governed by the following weak form,

$\sum_{\omega\in{\cal
T}}\int_{\omega}\Big{[}\Big{(}\frac{\phi^{h}-\phi^{h}_{n}}{m\tau}+\frac{\partial\psi_{\rm
int}}{\partial\phi}\Big{)}\delta\phi^{h}+\frac{\partial\psi_{\rm
int}}{\partial\nabla\phi}\cdot\nabla\delta\phi^{h}\Big{]}{\rm{d}}V+\sum_{\omega\in{\cal
T}_{\rm int}}\int_{\omega}\sum_{i=1}^{N_{\rm n}}\frac{\partial\bar{\psi}_{\rm
bulk}}{\partial\phi_{i}^{(\omega)}}\,\delta\phi_{i}^{(\omega)}{\rm{d}}V=0\qquad\forall\,\delta\textbf{{p}}{}_{\phi},$
(31)

which replaces the weak form (23) of the conventional phase-field model. The
first term in Eq. (31) includes the interfacial energy contribution and the
viscous evolution term and is identical to the respective part of the weak
form (23). The second term in Eq. (31) describes the driving force resulting
from the bulk energy $\bar{\psi}_{\rm bulk}$ of the laminated elements along
the interface. This term is the counterpart of the respective term in the LET
form of the virtual work principle, see Eqs. (28) and (30). Note that, in LET-
PF, the bulk energy depends on the order parameter field $\phi^{h}$ only in
the laminated elements $\omega\in{\cal T}_{\rm int}$ (through
$\eta^{(\omega)}$ and $\bm{n}^{(\omega)}$). Actually, also the set ${\cal
T}_{\rm int}$ itself depends on $\phi^{h}$ (likewise, ${\cal T}_{1}$ and
${\cal T}_{2}$) and may change during a time step and even between the
individual iterations of the global Newton scheme. Accordingly, to improve the
convergence behaviour, a regularization of the formula for $\eta^{(\omega)}$
in Eq. (25) is employed, as discussed in Section 3.3.

Summarizing, the difference with respect to the conventional phase-field
method is that, in LET-PF, the two sub-models are coupled only within the
laminated interface elements in ${\cal T}_{\rm int}$, see Figs. 1 and 2. This
is, in particular, apparent in Eq. (28), where $\bar{\bm{\sigma}}$ is the only
term that depends on $\phi^{h}$, in addition to the set ${\cal T}_{\rm int}$
itself. Likewise, in Eq. (31), the coupling occurs only through the overall
bulk energy $\bar{\psi}_{\rm bulk}$ of the laminates that is defined only
within the laminated interface elements in ${\cal T}_{\rm int}$. At the same
time, in the conventional phase-field method, the coupling occurs within the
whole volume of the diffuse interfaces, i.e., whenever $0<\phi<1$, see Fig. 2.
It follows that the mechanical equilibrium subproblem, which is governed by
the bulk contribution to the free energy, is treated in LET-PF in a
significantly sharper manner than in the conventional phase-field method, and
hence coarser meshes can be used in LET-PF to achieve a similar accuracy, as
illustrated by the numerical examples below.

###  3.3 Finite-element treatment of the LET-PF model

The governing equations of LET-PF, see Eqs. (28) and (31), define the set of
non-linear equations to be solved at each time step. The unknowns of the
problem are the nodal values of the displacement ($\bm{u}_{k}$) and order
parameter ($\phi_{k}$), both entering the formulation through the finite-
element approximation (21) with $\textbf{{p}}{}_{\bm{u}}$ and
$\textbf{{p}}{}_{\phi}$ denoting the respective global vectors of unknowns.
The general structure of the problem is the same as in the case of the
conventional phase-field model, and the residual form of the governing
equations of the two subproblems can be written as in Eq. (24), of course,
with due differences in the form of the residual vectors
$\textbf{{R}}_{\bm{u}}$ and $\textbf{{R}}_{\phi}$. A monolithic Newton-based
solution scheme is employed to solve these equations, just like in the case of
the conventional phase-field model, see the end of Section 2.3.

In the numerical examples, which are limited here to 2D, four-node
quadrilateral elements are employed with bilinear shape functions used for
both the displacement and the order parameter and with the standard $2\times
2$ Gauss quadrature. Computer implementation has been performed using _AceGen_
, a symbolic code generation system employing the automatic differentiation
(AD) and expression optimization techniques [47, 48], thus resulting in a
computationally efficient computer code. The actual finite-element
computations have been carried out using _AceFEM_ , a finite-element system
integrated with _AceGen_.

Preliminary computations have revealed that the convergence of the Newton
method at the global level is severely affected by the non-differentiable
functions (the absolute value and Macauley brackets) that are involved in the
formula for the volume fraction $\eta^{(\omega)}$ in element $\omega$, see Eq.
(25). In order to improve the robustness of the method, the following
regularization is applied to the absolute value function and Macauley brackets
used in Eq. (25),

$|\phi|_{\rm
reg}=\left\\{\begin{array}[]{ll}|\phi|&\mbox{if}\;|\phi|\geq\phi_{\rm reg},\\\
\frac{1}{2\phi_{\rm reg}}(\phi^{2}+\phi^{2}_{\rm
reg})&\mbox{if}\;|\phi|<\phi_{\rm
reg},\end{array}\right.\qquad\langle\phi\rangle_{\rm
reg}=\frac{1}{2}(\phi+|\phi|_{\rm reg}),$ (32)

where $\phi_{\rm reg}>0$ is a regularization parameter with the typical value
of 0.1 used in the computations reported in Section 4. As illustrated in C,
the above regularization significantly improves the convergence behaviour so
that the analysis can proceed with a larger time step, while it has a small
effect on the results.

##  4 Illustrative examples

###  4.1 Evolving circular inclusion

The aim of this study is to evaluate the overall performance of LET-PF and to
compare it to the conventional phase-field method (abbreviated to PFM in the
sequel) for a 2D problem that has an analytical solution, to be used as a
reference. The sharp-interface problem is specified first, followed by the
description of the computational model and by the results of a comprehensive
study.

####  4.1.1 Sharp-interface benchmark problem

We consider an elastic circular domain $\Omega_{2}$ (matrix) of the radius $R$
and an evolving elastic circular inclusion $\Omega_{1}$ of the radius
$\rho=\rho(t)$, see Fig. 3(a). The phase transformation and thus the evolution
of the inclusion are induced by the internal stresses resulting from the
volumetric eigenstrain $\bm{\varepsilon}_{1}^{\text{t}}=\epsilon\bm{I}$ in the
inclusion and by the energy of the interface $\Gamma$ of the density $\gamma$.
Propagation of the interface is assumed to be governed by the viscous kinetic
law, cf. Eq. (6)1,

$\hat{v}_{n}=-\dot{\rho}=\hat{m}\hat{f},\qquad\hat{f}=\hat{f}_{\rm
bulk}+\hat{f}_{\rm int},$ (33)

where $\hat{f}_{\rm bulk}$ and $\hat{f}_{\rm int}$ denote the local
thermodynamic driving forces originating from the bulk and interfacial energy,
respectively. The outer boundary is assumed to be free.

(a) (b)

Figure 3: Evolving circular inclusion: (a) scheme of the problem; (b)
computational domain with a regular (non-conforming) mesh of quadrilateral
elements. A coarse mesh of $11\times 11$ elements within the internal square
part (element size $h=0.1$) is shown, the actual computations are carried out
using significantly finer meshes.

The sharp-interface problem specified above is a 1D problem and admits an
analytical solution that is derived in B. In the special case of identical
elastic properties of the two phases, the evolution of the inclusion radius is
governed by the following differential equation, cf. Eq. (52),

$\dot{\rho}=-\frac{\hat{m}\gamma}{\rho}\left(1+A\frac{\rho}{\rho_{0}}\left(\frac{1}{2}-\left(\frac{\rho}{R}\right)^{2}\right)\right),\qquad
A=\frac{\rho_{0}E\epsilon^{2}}{(1-\nu^{2})\gamma},$ (34)

where $\rho_{0}$ denotes the initial radius of the inclusion, $E$ is the
Young’s modulus and $\nu$ is the Poisson’s ratio. Note that a dimensionless
parameter $A$ has been introduced which combines all the problem parameters
related to both the elastic strain energy and the interfacial energy (the
chemical energy contribution is not considered, hence the elastic strain
energy is the only contribution to the bulk energy). For small values of $A$,
the evolution is mainly governed by the interfacial energy, and for large
values of $A$ by the elastic strain energy. Exact integration of Eq. (34) is
not possible; integration is thus performed numerically.

The explicit formulae for the driving force contributions $\hat{f}_{\rm bulk}$
and $\hat{f}_{\rm int}$ are given by Eq. (53). Specifically, $\hat{f}_{\rm
int}=\gamma/\rho$ is the driving force originating from the interfacial energy
of the density $\gamma$, and the interface curvature is $1/\rho$. The driving
force $\hat{f}_{\rm bulk}$ originates from the elastic strain energy induced
by the inclusion eigenstrain. It is proportional to the elastic modulus $E$
and squared eigenstrain magnitude $\epsilon^{2}$ and depends on the
geometrical factor $\rho/R$. For $\rho_{0}=R/2$, as assumed in the
computations below, the driving force $\hat{f}=\hat{f}_{\rm bulk}+\hat{f}_{\rm
int}$ is positive and hence $\dot{\rho}$ is negative, which means that the
inclusion decreases in size and ultimately vanishes.

####  4.1.2 Computational model

The two methods (LET-PF and PFM) involve identical material and numerical
parameters, including those related to the finite-element model, and the
specifications provided below apply to both methods. The finite-element
computations are carried in 2D on a computational domain encompassing one
quarter of the circular domain with adequate symmetry conditions imposed along
the horizontal and vertical edges, see Fig. 3(b). In the computations, it is
assumed that the material properties of the inclusion and matrix are identical
with the Young’s modulus $E=1$ and Poisson’s ratio $\nu=0.25$. The radius of
the domain is $R=2$, and the initial radius of the inclusion is $\rho_{0}=1$.
The volumetric eigenstrain of the inclusion is assumed as $\epsilon=0.1$. The
effective mobility parameter is set to $\hat{m}=1$, and the actual mobility
parameter $m$ used in the computations is then specified as a function of
$\ell$ according to Eq. (15).

One of the goals here is to study the performance of LET-PF over a wide range
of the values of the interfacial energy density $\gamma$. To this end, the
values of $\gamma$ have been chosen from within the range
$\gamma\in[0.0001,0.003]$ (12 values have been chosen), so that the
dimensionless parameter $A$ varies within the range ${A}\in[3.6,107.]$. This
range has been chosen such that the extreme values correspond to two regimes
of interest. For $\gamma=0.0001$ ($A=107.$), the bulk contribution
$\hat{f}_{\rm bulk}$ to the driving force is much greater than the interfacial
contribution $\hat{f}_{\rm int}$ over nearly the whole range of inclusion
radii $\rho$ (note that $\hat{f}_{\rm int}\to\infty$ for $\rho\to 0$), see
Fig. 4. Accordingly, in this (elasticity-driven) regime, evolution is mostly
governed by the elastic strain (bulk) energy. On the other hand, for
$\gamma=0.003$ ($A=3.6$), both contributions are initially of the same order,
see Fig. 4. The third regime, in which $\hat{f}_{\rm int}\gg\hat{f}_{\rm
bulk}$, is not interesting here because the evolution is then fully governed
by the interfacial energy and the effect of LET-PF is negligible. Moreover, it
has been checked that the response and the performance of both methods are
then quite similar as for $\hat{f}_{\rm int}\approx\hat{f}_{\rm bulk}$.

Figure 4: Thermodynamic driving forces $\hat{f}_{\rm bulk}$ and $\hat{f}_{\rm
int}$ as a function of the inclusion radius $\rho$, see Eq. (53).

The impact of the interface thickness parameter $\ell$ on the accuracy of the
results is also examined, with the values of $\ell$ specified relative to the
characteristic element size $h$, namely, $\ell/h\in\\{0.75,1,1.5,2\\}$. As
mentioned earlier, for the selected parameter values, the inclusion disappears
completely during the evolution process by gradually reducing its radius to
zero. Therefore, a regular grid of square finite elements is employed in the
central region containing the inclusion, and quadrilateral elements with
somewhat less regular shapes are employed in the remaining part of the domain,
as depicted in Fig. 3(b). The main computations are performed for two mesh
sizes, $h=0.02$ and $h=0.01$, which correspond to, respectively, $55\times 55$
and $110\times 110$ elements in the central region. At the initial time $t=0$,
the initial profile of the order parameter $\phi$ is prescribed according to
formula (14) with $\xi$ replaced by the radius $r$ and $\xi_{0}$ replaced by
the initial inclusion radius $\rho_{0}$.

In all computations, an adaptive time incrementation scheme has been used,
such that the time increment is increased or decreased depending on the
current convergence behaviour. In each case, the maximum time increment
$\Delta t_{\rm max}$ has been prescribed relative to $T_{\rm exact}$, the time
of the complete evolution as resulting from the analytical solution.

The results will be presented in terms of the mean inclusion radius
$\bar{\rho}$ which is defined as the average radius of the inclusion defined
by the $\phi=\frac{1}{2}$ level set. The individual radii are measured along a
number of directions taken every 1 degree.

####  4.1.3 Results

The time evolution of the mean inclusion radius $\bar{\rho}$ is shown in Fig.
5 for the two extreme values of the interfacial energy, $\gamma=0.0001$ and
$\gamma=0.03$. As a reference, the analytical solution is also included in
Fig. 5. The computations have been performed with small time increments,
$\Delta t_{\rm max}=T_{\rm exact}/500$, so that the error introduced by the
time integration is insignificant. It can be seen that, in the elasticity-
driven regime ($\gamma=0.0001$), PFM is highly inaccurate (Fig. 5(b)), while
LET-PF performs much better (Fig. 5(a)). In the interfacial-energy-driven
regime ($\gamma=0.003$), both methods perform similarly well. Fig. 5
illustrates also a significant effect of the interface thickness parameter
$\ell$.

(a) (b) (c) (d)

Figure 5: Time evolution of the mean inclusion radius $\bar{\rho}$ predicted
by LET-PF (a,c) and PFM (b,d) for $\gamma=0.0001$ (a,b) and for $\gamma=0.003$
(c,d). The computations are performed using a coarse mesh, $h=0.02$, and small
time increment, $\Delta t_{\rm max}=T_{\rm exact}/500$. As a reference, the
analytical solution is depicted by a solid black line.

The effect of mesh density and time increment is illustrated in Fig. 6 for an
intermediate value of the interfacial energy, $\gamma=0.0008$. In addition to
the small time increment ($\Delta t_{\rm max}=T_{\rm exact}/500$) and thus
highly accurate time integration, see the top figures in Fig. 6, the analysis
has also been performed with no restriction on the time increment ($\Delta
t_{\rm max}=T_{\rm exact}$) so that the actual time increments result from the
current convergence behaviour and are significantly larger, see the bottom
figures in Fig. 6. It can be seen that the effect of time step is significant.
As discussed later, the number of time steps needed to complete the simulation
can be (and will be) used as a measure of the robustness of the model.
Concerning the effect of mesh density, it is more pronounced in the case of
PFM than LET-PF. Again, this effect will be studied in more detail later.

(a) (b) (c) (d)

Figure 6: Time evolution of the mean inclusion radius $\bar{\rho}$ for
$\gamma=0.0008$: (a) $\ell=0.75h$, (b) $\ell=h$, (c) $\ell=1.5h$, (d)
$\ell=2h$. The results shown in the top row are obtained for a small time
increment (adaptive time stepping with $\Delta t_{\rm max}=T_{\rm
exact}/500$), and those in the bottom row for significantly larger time
increments ($\Delta t_{\rm max}=T_{\rm exact}$).

Fig. 7 shows the maps of the order parameter $\phi$, radial stress
$\sigma_{rr}$ and hoop stress $\sigma_{\theta\theta}$ at the instant when the
inclusion radius reaches $\bar{\rho}\approx 0.85$. While the order parameter
$\phi$ and the hoop stress $\sigma_{\theta\theta}$ exhibit no substantial
differences between LET-PF and PFM, the radial stress $\sigma_{rr}$, which
according to the analytical solution should be continuous at the interface, is
better represented by LET-PF. This is particularly evident along the
horizontal and vertical edges of the computational domain, as the interface is
then approximately parallel to the element edges and the LET technique is then
highly accurate [26]. A more profound insight into these fields behaviours is
presented in Fig. 8, showing the profiles of the radial and hoop stresses
along the horizontal direction. Regarding the radial stress, as shown in
Figure 8(a), a significant distinction between the two methods is readily
apparent. PFM exhibits fluctuations within the interface region, while LET-PF
produces results much closer to the analytical solution. As for the hoop
stress, as illustrated in Figure 8(b), the situation appears somewhat
different—both methods seem to yield a similar response. However, it is
crucial to note that, when utilizing LET-PF, the information is available also
about the local strains and stresses resulting from the homogenization process
in the laminated elements. The local stresses are indicated in the insets in
Fig. 8 by markers, thereby revealing that LET-PF provides also the local
response that aligns more closely with the analytical solution.

PFM LET-PF (a) (b) (c)

Figure 7: Maps of the order parameter $\phi$ (a), radial stress $\sigma_{rr}$
(b), and hoop stress $\sigma_{\theta\theta}$ (c) obtained using LET-PF (top)
and PFM (bottom) at the instant when the inclusion radius reaches
$\bar{\rho}\approx 0.85$ ($\gamma=0.0008$, $h=0.02$, $\ell=1.5h$).

(a) (b)

Figure 8: Profiles of the radial stress $\sigma_{rr}$ (a) and hoop stress
$\sigma_{\theta\theta}$ (b) along the bottom edge of the computational domain
($\bar{\rho}\approx 0.85$, $\gamma=0.0008$, $h=0.02$, $\ell=1.5h$). The lines
depict the overall stresses, while the markers in the insets represent the
local stresses within each phase, as predicted by LET-PF.

The results presented in Figs. 5 and 6 suggest that LET-PF outperforms PFM in
terms of accuracy. To examine the accuracy and efficiency of LET-PF in more
detail, a comprehensive study has been performed by varying the governing
material and numerical parameters. The material parameters are fully
characterized by the dimensionless parameter $A$, Eq. (34)2, and this
parameter is controlled by varying the interfacial energy $\gamma$, as
described above, while the remaining material parameters are fixed. Concerning
the numerical parameters, two mesh densities are considered,
$h\in\\{0.01,0.02\\}$, four values of the interface thickness parameter,
$\ell/h\in\\{0.75,1,1.5,2\\}$, and two values of the maximum time increment,
$\Delta t_{\rm max}\in\\{T_{\rm exact},T_{\rm exact}/500\\}$.

As shown in Fig. 4, when the inclusion radius approaches zero, the
thermodynamic driving force, and specifically its interfacial contribution
$\hat{f}_{\rm int}$, increases to infinity. The evolution is then fully
governed by $\hat{f}_{\rm int}$ and the difference between LET-PF and PFM
diminishes. Moreover, the radius of curvature of the interface becomes then
small compared to the element size and interface thickness so that additional
effects influence the solution. Accordingly, in the following, the final part
of the evolution is excluded from the analysis, and only the range of the
inclusion radii between $\rho_{0}$ and $0.15\rho_{0}$ is considered. The
corresponding duration of the evolution process, resulting from the analytical
solution, depends on $A$ and is denoted by $T_{\rm exact}$. The maximum time
increment $\Delta t_{\rm max}$ in the adaptive time incrementation scheme is
then expressed in terms of $T_{\rm exact}$, as specified above.

Fig. 9 presents the relative error as a function of parameter $A$. The
relative error is defined here as

${\rm relative\,error}=\dfrac{\int_{0.15\rho_{0}}^{\rho_{0}}\lvert\tau^{\rm
exact}(r)-\tau^{\rm
num}(r)\rvert\,{\rm{d}}r}{\int_{0.15\rho_{0}}^{\rho_{0}}\tau^{\rm
exact}(r)\,{\rm{d}}r},$ (35)

where $\tau={\rho}^{-1}$ is the inverse function of $\rho(t)$ that describes
the time as a function of the inclusion radius, and $\tau^{\rm exact}$ and
$\tau^{\rm num}$ correspond to the analytical and numerical solutions,
respectively. The use of the inverse function ${\rho}^{-1}$ is here necessary
because the duration of the evolution process in each simulation is, in
general, different than the duration in the analytical solution so that
integration of $\lvert\rho^{\text{exact}}(t)-\rho^{\rm num}(t)\rvert$ would
not be feasible. The integration in Eq. (35) is performed numerically using
the trapezoidal rule. The results in Fig. 9 illustrate that, in terms of
accuracy, LET-PF outperforms PFM in all cases studied. Although both methods
exhibit similar trends for lower values of $A$ (where the interfacial energy
has a greater contribution), LET-PF demonstrates superior performance as the
contribution of the bulk (elastic) energy becomes more prominent. This is
particularly evident in Fig. 9(a,b), where the errors for PFM reach extreme
values of nearly 0.8, while LET-PF performs much better. In general, it can be
seen that the accuracy of LET-PF only weakly depends on $A$.

(a) (b) (c) (d)

Figure 9: Relative error, Eq. (35), evaluated for all considered cases for
LET-PF (solid lines) and PFM (dashed lines). The shaded regions correspond to
the error that is approximately twice higher than the average error of LET-PF,
i.e., exceeding 0.3 for $\Delta t_{\max}=T_{\rm exact}$ (panels (a) and (b))
and exceeding 0.15 for $\Delta t_{\max}=T_{\rm exact}/500$ (panels (c) and
(d)). The data points falling within the shaded regions are marked with
crosses, see also Fig. 10.

Fig. 10 shows the number of time steps needed to complete the simulation for
all cases included in Fig. 9. When there is no restriction on the maximum time
increment ($\Delta t_{\rm max}=T_{\rm exact}$), so that the simulation
proceeds with the time increment dictated by the current convergence
behaviour, the number of time steps may serve as an indicator of the
robustness of the computational scheme (the smaller the number of time steps,
the more robust the scheme). Moreover, the number of time steps quantifies the
computational cost since the total computation time is proportional to the
number of time steps. However, in terms of the cost, the comparison of LET-PF
and PFM is not immediate, because the cost of evaluation of the element
quantities (residual vector and tangent matrix) differs for the two methods.
The corresponding results are presented in Fig. 10(a,b). It can be seen that,
for small $A$, the number of time steps is similar for both methods, PFM being
slightly more efficient. For high values of $A$, the number of time steps
significantly increases for LET-PF, while it slightly decreases for PFM.
However, caution is needed when interpreting these results. As has been shown
in Fig. 9, PFM is highly inaccurate for high $A$, and the corresponding data
points are marked with crosses in Fig. 10, consistent with Fig. 9. It follows
that PFM is indeed cheaper and more robust for high $A$, but it is then
unacceptable in terms of accuracy. On the other hand, for high $A$, LET-PF can
only proceed with small time increments, but it delivers results of similar
accuracy as for small $A$.

(a) (b) (c) (d)

Figure 10: The number of time steps needed to complete the simulation for
$\Delta t_{\rm max}=T_{\rm exact}$ (a,b) and for $\Delta t_{\rm max}=T_{\rm
exact}/500$ (c,d). Two mesh sizes are used, namely $h=0.02$ (a,c) and $h=0.01$
(b,d). The results correspond to the cases shown in Fig. 9, and the cross
markers indicate the data points of exceedingly high error, see Fig. 9.

When the time increment is small ($\Delta t_{\rm max}=T_{\rm exact}/500$), for
small and intermediate $A$, the simulation proceeds with the prescribed
maximum time increment, and the number of time steps is close to 500 for both
LET-PF and PFM, see Fig. 10(c,d). For high $A$, in the case of LET-PF, the
number of time steps increases as in the case of $\Delta t_{\rm max}=T_{\rm
exact}$. On the other hand, in the case of PFM, the number of time steps
decreases for high $A$, and this is related to the high error of the method
and to the related reduction of the duration of the evolution process, see,
for instance, Fig. 5(b). Note that the corresponding data points are marked by
crosses, as in Fig. 9.

The results reported in Fig. 9 show that the relative error significantly
depends on parameter $A$, on the maximum time increment $\Delta t_{\rm max}$,
and on the interface thickness parameter $\ell$. The effect of mesh density,
although present, is hard to perceive as the other dependencies prevail.
Accordingly, additional computations have been performed in which four mesh
densities have been used, $h\in\\{0.05,0.02,0.01,0.005\\}$, and two scenarios
concerning the interface thickness parameter $\ell$ have been considered. In
the first scenario, $\ell$ is proportional to the mesh size $h$ and is set to
$\ell=h$. In the second scenario, $\ell$ is set to a constant value
$\ell=0.02$ so that $\ell/h\in\\{0.4,1,2,4$). The simulations are performed
for $\gamma=0.0008$ and $\Delta t_{\rm max}=T_{\rm exact}/200$. Fig. 11 shows
the total elastic strain energy $\Psi_{\rm el}$ as a function of the mean
inclusion radius $\bar{\rho}$.

The energy is normalized by the area of the inclusion, $\Psi_{\rm
el}/(\pi\bar{\rho}^{2})$, so that the differences between the individual
curves are better visible. It can be seen that LET-PF outperforms PFM in both
scenarios. In the fist scenario ($\ell=h$), LET-PF converges to the analytical
solution visibly faster than PFM, see Fig. 11(a). This is because LET-PF
treats the interface in a semisharp manner, and it is the diffuseness of the
interface in PFM that is the source of the observed error. This is illustrated
in Fig. 11(b), where the interface thickness is kept constant ($\ell=0.02$)
and PFM converges to a solution that is significantly different from the
correct one, particularly at small $\bar{\rho}$, when the interface thickness
is large relative to the inclusion radius. At the same time, the large
thickness of the interface does not affect the performance of LET-PF because
the $\phi=\frac{1}{2}$ level set, which specifies the actual interface in the
elasticity problem, only weakly depends on the thickness of the diffuse
interface.

(a) (b)

Figure 11: Normalized total elastic strain energy $\Psi_{\rm
el}/(\pi\bar{\rho}^{2})$ as a function of the mean radius $\bar{\rho}$ for
increasing mesh density (element size $h$ decreases from left to right) and
for two scenarios: (a) $\ell=h$, (b) $\ell=0.02$ ($\gamma=0.0008$, $\Delta
t_{\rm max}=T_{\rm exact}/200$).

Since the reference sharp-interface problem is axisymmetric and the interface
$\Gamma$ is then a circle, any deviation from the circular shape can be
treated as a measure of the numerical solution error. As a quantitative
measure of this deviation, the coefficient of variation of the radius can be
calculated at each time instant,

$\mathrm{CV}_{\rho}=\frac{\sigma_{\rho}}{\bar{\rho}},$ (36)

where $\sigma_{\rho}$ is the standard deviation of the orientation-dependent
inclusion radius $\rho(\theta)$. Some representative results presenting
$\mathrm{CV}_{\rho}$ as a function of the mean inclusion radius $\bar{\rho}$
are shown in Fig. 12. It follows that the deviation from the circular shape
depends on the model parameters and varies during the evolution process,
including small oscillations resulting from the interface traversing the
finite-element mesh. Anyway, the deviation is relatively small in all cases
(typically $\mathrm{CV}_{\rho}<0.02$).

(a) (b) (c)

Figure 12: Deviation of the inclusion shape from a circular one, quantified by
the coefficient of variation of the radius, $\mathrm{CV}_{\rho}$, Eq. (36), as
a function of the mean inclusion radius $\bar{\rho}$ for three representative
cases: (a) $\gamma=0.003$, $\ell=0.75h$; (b) $\gamma=0.0001$, $\ell=0.75h$;
(c) $\gamma=0.0001$, $\ell=1.5h$ ($h=0.01$, $\Delta t_{\rm max}=T_{\rm
exact}/500$).

In order to comprehensively investigate this characteristic, the mean value of
$\mathrm{CV}_{\rho}$ has been calculated for each case,

$\overline{\mathrm{CV}_{\rho}}=\frac{1}{0.85\rho_{0}}\int\limits_{0.15\rho_{0}}^{\rho_{0}}\mathrm{CV}_{\rho}\,\mathrm{d}\rho,$
(37)

and the aggregated results, separately for the two mesh sizes, are shown in
Fig. 13. This analysis reveals that both LET-PF and PFM result in a low mean
value of the coefficient of variation of the radius with the maximum value
less than 2%. However, it can be seen that, for the interface thickness
$\ell=0.75h$, for both methods the deviation from the circular shape is higher
than for larger interface thicknesses, which can be attributed to the
interface pinning which is expected to diminish with the interface thickness
increasing relative to the element size.

(a) (b)

Figure 13: Mean coefficient of variation of the inclusion radius
$\overline{\mathrm{CV}_{\rho}}$, Eq. (37), for (a) $h=0.02$ and (b) $h=0.01$
($\Delta t_{\rm max}=T_{\rm exact}/500$).

###  4.2 Single inclusion in a constrained domain

In this example, a 2D square domain with a single evolving inclusion is
considered. To ensure a non-trivial steady-state solution, the displacements
are fully constrained on the boundary and a volumetric eigenstrain of the same
value but with opposite signs is introduced in both phases, namely
$\bm{\varepsilon}^{\rm t}_{2}=-\bm{\varepsilon}^{\rm t}_{1}=\epsilon\bm{I}$
with $\epsilon=0.1$. Accordingly, considering a viscous evolution towards the
minimum of the elastic strain energy, a steady-state solution is expected with
the overall volume fraction of both phases close to 0.5, because this leads to
a null overall eigenstrain in the domain, thus effectively minimizing the
total elastic strain energy (the chemical energy is not considered).

The dimensions of the domain are $1\times 1$, the inclusion is initially
located in the domain centre and its initial radius is $\rho_{0}=0.1$.
Homogeneous elastic properties are assumed, specifically, the Young’s modulus
$E=1$ and the Poisson’s ratio $\nu=0.25$. Further, the following interface
parameters are assumed: the interfacial energy density $\gamma=0.0003$, the
effective mobility $\hat{m}=1$, and the interface thickness parameter
$\ell=1.5h$.

Preliminary computations have been performed on the mesh of $100\times 100$
elements. Fig. 14 illustrates the evolution of the inclusion shape for both
methods, starting from $t=0$ and ending in a steady state at
$t\rightarrow\infty$ (a sufficiently large final time has been adopted).
Discernible differences between LET-PF and PFM become evident almost
instantaneously. The final shape of the inclusion is different, with the
corners being more acute for PFM. Also, the evolution towards a square-like
shape proceeds faster for PFM, compare the snapshots corresponding to $t=15$.
Moreover, for PFM, artifacts manifest in the corners of the domain at a
relatively early stage, a phenomenon not observed for LET-PF. Finally, in the
case of PFM, the interface thickness increases with respect to the theoretical
one, which is prescribed at $t=0$, while no visible difference in the
interface thickness is observed for LET-PF.

Figure 14: Evolution of the inclusion within a constrained domain for LET-PF
(top) and PFM (bottom) obtained for the mesh of $100\times 100$ elements. The
colour maps depict the order parameter $\phi$. The black contours indicate the
position of the $\phi=\frac{1}{2}$ level set.

To get a better insight, the computations have been performed for a range of
mesh densities with $N\times N$ elements where
$N\in\\{26,50,100,200,400,800,1600\\}$. Fig. 15 shows the steady-state shapes
for selected mesh densities (the results for two finest meshes are omitted, as
no differences are visible with further mesh refinement). Based on the results
obtained, a clear advantage of LET-PF over PFM can be observed. For coarse
meshes, PFM fails, yielding a trivial solution with $\phi=0.5$ in the entire
domain. This stands in contrast to LET-PF which produces qualitatively correct
results already for the coarsest mesh. Secondly, while the final shape seems
to converge to the same solution for both methods, the convergence is faster
for LET-PF. This is also apparent in Fig. 16 which shows the total elastic
strain energy $\Psi_{\rm el}$ and the total interfacial energy $\Psi_{\rm
int}$ in the steady state as a function of mesh density.

Figure 15: Steady-state shapes for LET-PF (top) and PFM (bottom) as a function
of mesh density (increasing from left to right).

(a) (b)

Figure 16: Convergence of the total elastic strain energy $\Psi_{\rm el}$ (a)
and of the total interfacial energy $\Psi_{\rm int}$ (b) with mesh refinement.
The points corresponding to the trivial solution obtained for PFM for $N=26$
and $N=50$ are not included, the dashed lines show that these points lie well
outside the plot range.

To examine the evolution of shape quantitatively, the distance from the centre
to the interface (level set $\phi=0.5$), measured in the horizontal direction
and along the diagonal, is shown in Fig. 17 for three mesh densities. For the
finest mesh ($N=1600$), the results obtained for LET-PF and PFM are very close
one to the other, which again suggests that both methods converge to the same
solution. However, consistent with Fig. 16, the convergence rate is visibly
different. For PFM, the difference between the results obtained for $N=200$
and $N=1600$ is significantly larger than in the case of LEF-PF. For the
coarsest mesh ($N=26$), a trivial solution is obtained for PFM, while LET-PF
delivers a reasonably accurate solution.

(a) (b)

Figure 17: Evolution of the inclusion size expressed by the distance from the
centre to the propagating interface measured in the horizontal (a) and
diagonal (b) directions, shown for selected mesh densities. Cross markers
indicate the instant when the trivial solution is achieved for PFM with
$\phi=0.5$ in the entire domain.

###  4.3 Three inclusions in a constrained domain

As the last example, a problem similar to that of Section 4.2 is studied with
the only difference that three inclusions are considered in the initial
conditions. The initial inclusion radii are 0.1, 0.15 and 0.2, and the
corresponding positions of the inclusion centres are (0.25,0.25), (0.75,0.30)
and (0.35,0.75), respectively, for the origin of the coordinate system located
in the lower-left corner of the domain. The remaining parameters of the
problem are the same as in Section 4.2.

Fig. 18 shows selected snapshots of the evolution process for four mesh
densities. Consider first the case of the finest mesh of $400\times 400$
elements. Here, the evolution is very similar for LET-PF and PFM. At the
beginning, all inclusions start to grow. Then two inclusions coalesce and grow
further, while the third inclusion (that in the lower-left corner) shrinks and
ultimately vanishes. In the case of LET-PF, this scenario is obtained also for
two coarser meshes, although the two inclusions coalesce at an earlier stage
(at $t=29$ for $N=200$ and at $t=22$ for $N=100$), and only for the coarsest
mesh ($N=50$) a different scenario is obtained in which all inclusions
coalesce. However, in the case of PFM, the correct scenario is obtained for
$N=200$, even if the evolution proceeds then significantly faster (see the
snapshot at $t=29$). At the same time, an incorrect scenario is obtained for
$N=100$ (three inclusions coalesce), while for $N=50$ the microstructure
evolves towards the trivial solution, as in Section 4.2.

Figure 18: Evolution of three inclusions within a constrained domain obtained
for four mesh densities and for LET-PF (upper snapshots) and PFM (lower
snapshots). The colour maps depict the order parameter $\phi$. The black
contours indicate the position of the $\phi=0.5$ level set.

The final shapes correspond, in general, to those obtained for the single-
inclusion problem in Section 4.2. Note that, for LET-PF and $N=100$, the
shrinking inclusion does not vanish completely, rather a small remnant
inclusion persists in the steady state, see the lower-left corner of the
domain. Also, the final shape obtained for LET=PF for $N=50$ is not perfectly
symmetric with the upper-right corner of the inclusion being somewhat more
rounded than the other ones. Both effects suggest that LET-PF may be prone to
mesh pinning, and this effect may require further attention in the future.

Nevertheless, the present numerical example confirms that LET-PF performs
significantly better than PFM for the class of problems considered here, in
particular, it is capable of capturing the correct evolution scenario for a
coarser mesh.

##  5 Conclusion

The laminated element technique (LET), which has been recently developed for
an efficient treatment of weak discontinuities on non-matching meshes [26],
has been combined with the phase-field method and thus extended to problems
involving moving interfaces and, specifically, to microstructure evolution
problems. The essence of LET is in treating an element that is cut by an
interface as a simple laminate with the lamination orientation and phase
volume fraction depending on the orientation and position of the interface
within the element. The position of a propagating interface is then determined
using a continuous order parameter which plays the role of a level set
function and, at the same time, carries the information about the interfacial
energy. Evolution of the phase field order parameter is governed by a
Ginzburg–Landau type evolution equation.

The resulting LET-PF method has been implemented in a finite-element code, and
its performance has been illustrated by a comprehensive set of numerical
examples, including a comparison to the conventional phase-field method (PFM).
The implementation and the computations are limited to 2D problems in the
framework of small-strain elasticity, but generalization to 3D and finite
deformations is immediate. The results show that, when the interfacial energy
has a significant contribution to the total thermodynamic driving force
(interfacial-energy-driven evolution), the two methods behave similarly in
terms of both accuracy and robustness. However, when the elastic strain energy
contribution dominates (elasticity-driven evolution), LET-PF is found
significantly more accurate than PFM. As a result, a coarser mesh (along with
a larger interface thickness) can be used in LET-PF to get results of similar
accuracy, thus reducing the computational cost.

Overall, LET-PF seems to be a promising approach for problems involving
propagating interfaces. Further work is needed to improve its robustness and
to examine its performance for a wider class of problems. While generalization
to a multi-phase framework is an open problem, it can be readily applied to a
wide range of two-phase evolution problems.

### Acknowledgement

This work has been partially supported by the National Science Centre (NCN) in
Poland through Grant No. 2018/29/B/ST8/00729.

##  Appendix A Micro-macro transition for a simple laminate

In this section, a laminate composed of two phases is considered in the
framework of linear elasticity with eigenstrain, and the micro-macro
transition is discussed for this microstructure. A simple laminate is a
microstructure in which layers of two phases are separated by parallel planar
interfaces, and the microstructure is defined by the interface normal $\bm{n}$
and volume fraction $\eta$. Under the assumption of separation of scales,
strains and stresses are homogeneous within each layer and are identical in
the layers of the same phase. Accordingly, the compatibility conditions
formulated in Eq. (3) for the interface jumps hold for the strains
$\bm{\varepsilon}_{i}$ and stresses $\bm{\sigma}_{i}$ within individual
phases,

$\Delta\bm{\varepsilon}=\frac{1}{2}(\bm{c}\otimes\bm{n}+\bm{n}\otimes\bm{c}),\qquad\Delta\bm{\sigma}\cdot\bm{n}=\bm{0},$
(38)

with
$\Delta(\mbox{\footnotesize$\square$})=(\mbox{\footnotesize$\square$})_{2}-(\mbox{\footnotesize$\square$})_{1}$,
as defined in Section 2.3.

The overall quantities, denoted by a superimposed bar, are obtained by simple
averaging,

$\bar{\psi}=\\{\psi\\},\qquad\bar{\bm{\varepsilon}}=\\{\bm{\varepsilon}\\},\qquad\bar{\bm{\sigma}}=\\{\bm{\sigma}\\},$
(39)

where
$\\{\mbox{\footnotesize$\square$}\\}=(1-\eta)(\mbox{\footnotesize$\square$})_{1}+\eta(\mbox{\footnotesize$\square$})_{2}$.
Note that the overall free energy $\bar{\psi}$ of the laminate constitutes the
bulk free energy in the LET setting, Eqs. (30) and (31), hence
$\bar{\psi}_{\rm bulk}=\bar{\psi}$, and the subscript ‘bulk’ is used in
Section 3 to differentiate $\bar{\psi}_{\rm bulk}$ from the interfacial energy
contribution, as in the conventional phase-field model of Section 2.2.

Using the compatibility condition (38)1 and the averaging rule (39)2, the
local strains $\bm{\varepsilon}_{i}$ can be expressed as a function of the
overall strain $\bar{\bm{\varepsilon}}$ and vector $\bm{c}$,

$\bm{\varepsilon}_{1}=\bar{\bm{\varepsilon}}-\eta\Delta\bm{\varepsilon},\qquad\bm{\varepsilon}_{2}=\bar{\bm{\varepsilon}}+(1-\eta)\Delta\bm{\varepsilon}.$
(40)

The unknown vector $\bm{c}$ is then determined by solving the compatibility
condition (38)2 with the local stresses governed by the constitutive equation
$\bm{\sigma}_{i}=\partial\psi_{i}/\partial\bm{\varepsilon}_{i}=\mathbb{L}_{i}:(\bm{\varepsilon}_{i}-\bm{\varepsilon}^{\rm
t}_{i})$. In the present case of linear elasticity, this yields a linear
equation for vector $\bm{c}$ (note that the transformation strains
$\bm{\varepsilon}^{\rm t}_{i}$ are here constant and known).

In view of linearity, the complete micro-macro transition scheme can be solved
in closed form. In particular, the overall free energy $\bar{\psi}$ can be
expressed in the following form,

$\bar{\psi}=\bar{\psi}(\bar{\bm{\varepsilon}},\eta,\bm{n})=\bar{\psi}{}^{0}+\frac{1}{2}(\bar{\bm{\varepsilon}}-\bar{\bm{\varepsilon}}^{\rm
t}):\bar{\mathbb{L}}:(\bar{\bm{\varepsilon}}-\bar{\bm{\varepsilon}}^{\rm t}),$
(41)

where

$\bar{\psi}{}^{0}=(1-\eta)\psi^{0}_{1}+\eta\psi^{0}_{2}+\frac{1}{2}\eta(1-\eta)\Delta\bm{\varepsilon}^{\rm
t}:\mathbb{S}:\Delta\bm{\varepsilon}^{\rm t},\qquad\bar{\bm{\varepsilon}}^{\rm
t}=\\{\mathbb{B}^{\rm T}:\bm{\varepsilon}^{\rm t}\\}.$ (42)

Here, $\bar{\mathbb{L}}$ is the overall elastic stiffness tensor,
$\mathbb{B}_{1}$ and $\mathbb{B}_{2}$ are the fourth-order stress
concentration tensors, and $\mathbb{S}$ is a fourth-order tensor. Closed-form
formulae for $\bar{\mathbb{L}}$, $\mathbb{B}_{1}$, $\mathbb{B}_{2}$ and
$\mathbb{S}$ in the general case of distinct anisotropy of the phases can be
found, for instance, in [49, 50]. In the case of homogeneous elastic
properties, $\mathbb{L}=\mathbb{L}_{1}=\mathbb{L}_{2}$, those formulae
simplify significantly, namely $\bar{\mathbb{L}}=\mathbb{L}$ and
$\mathbb{B}_{1}=\mathbb{B}_{2}=\mathbb{I}$, where $\mathbb{I}$ denotes the
fourth-order identity tensor, so that $\bar{\bm{\varepsilon}}^{\rm
t}=\\{\bm{\varepsilon}^{\rm t}\\}$, while $\mathbb{S}$ is then given by,

$\mathbb{S}=\mathbb{L}-\mathbb{L}:(\bm{n}\otimes(\bm{n}\cdot\mathbb{L}\cdot\bm{n})^{-1}\otimes\bm{n}):\mathbb{L}.$
(43)

The explicit expression (41) for the the overall free energy $\bar{\psi}$ is
provided above in order to illustrate the general structure of the micro-macro
transition for a simple laminate, and the details are omitted here. In the
actual computer implementation of LET-PF, these closed-form formulae are not
used. Rather, the micro-macro transition equations (38)–(40) along with the
local constitutive equations are solved directly on demand. The computational
scheme relies on the use of the automatic differentiation (AD) technique
implemented in _AceGen_ [47, 48], which facilitates consistent linearization
and computation of the derivatives with respect to $\eta$ and $\bm{n}$ that
are needed in LET-PF, cf. Eq. (31). The computational scheme is a special case
of a more general one that is discussed in Appendix A in [26] in the context
of elastoplastic laminates at finite deformation. The reader is referred to
[26] for details.

##  Appendix B Analytical solution for the evolving circular inclusion

For a prescribed inclusion radius $\rho$, the elasticity problem corresponding
to the setup defined in Section 4.1, see Fig. 3(a), can be solved in a closed
form by using the classical Lamé solution of a thick-walled cylinder under
internal pressure. Enforcing the continuity of the displacement and traction
at the interface between the elastic cylinder and the elastic inclusion with
eigenstrain $\bm{\varepsilon}_{1}^{\text{t}}=\epsilon\bm{I}$, the radial
displacement $u_{r}$ can be determined as

$u_{r}(r)=\begin{cases}(1-\epsilon^{\ast})\epsilon\,r&0\leqslant
r\leqslant\rho,\\\\[8.61108pt]
\displaystyle\frac{\left(\lambda_{1}+\mu_{1}\right)\left(\mu_{2}r^{2}+\left(\lambda_{2}+\mu_{2}\right)R^{2}\right)\rho^{2}\epsilon^{\ast}\epsilon}{\mu_{2}\left(\lambda_{2}+\mu_{2}\right)\left(R^{2}-\rho^{2}\right)r}\quad&\rho<r\leqslant
R,\end{cases}$ (44)

while the circumferential displacement is zero, $u_{\theta}=0$. Here,
$\epsilon^{\ast}$ denotes the fraction of the eigenstrain that is recovered
elastically in the inclusion,

$\epsilon^{\ast}=\frac{\mu_{2}\left(\lambda_{2}+\mu_{2}\right)\left(R^{2}-\rho^{2}\right)}{\mu_{2}\left(\lambda_{1}-\lambda_{2}+\mu_{1}-\mu_{2}\right)\rho^{2}+\left(\lambda_{2}+\mu_{2}\right)\left(\lambda_{1}+\mu_{1}+\mu_{2}\right)R^{2}},$
(45)

so that $(1-\epsilon^{\ast})\epsilon$ is the total strain in the inclusion,
and $\lambda_{i}$ and $\mu_{i}$ are the Lamé constants of the phases, $i=1,2$.
The total elastic strain energy of the system is then given by

$\hat{\Psi}_{\rm
el}=2\pi\left(\lambda_{1}+\mu_{1}\right)\rho^{2}\epsilon^{\ast}\epsilon^{2}.$
(46)

The total free energy is a sum of the elastic strain energy $\hat{\Psi}_{\rm
el}$ and interfacial energy $\hat{\Psi}_{\rm int}$ (the chemical energy
contribution is not considered in this example),

$\hat{\Psi}=\hat{\Psi}_{\rm el}+\hat{\Psi}_{\rm int},\qquad\hat{\Psi}_{\rm
int}=\int_{\Gamma}\gamma\,{\rm{d}}S=2\pi\rho\gamma,$ (47)

where $\gamma$ is the interfacial energy density, assumed isotropic.

Now consider that the interface $\Gamma$ is evolving, i.e., the inclusion
radius is a function of time, $\rho=\rho(t)$. The global rate-potential
$\hat{\Pi}$ is thus formulated as a sum of the rate of the total free energy,
$\dot{\hat{\Psi}}={\rm{d}}\hat{\Psi}/{\rm{d}}t$, and the global dissipation
potential $\hat{\mathcal{D}}$,

$\hat{\Pi}=\dot{\hat{\Psi}}+\hat{\mathcal{D}}.$ (48)

Note that no external loading is applied so that the potential energy of the
loading is equal to zero and does not contribute to $\hat{\Pi}$. The
dissipation is here assumed to be of purely viscous nature, i.e., the local
dissipation potential $\hat{D}$ is quadratic in terms of the interface speed
$\hat{v}_{n}=-\dot{\rho}$,

$\hat{\mathcal{D}}=\int_{\Gamma}\hat{D}\left(\dot{\rho}\right)\,{\rm{d}}S=2\pi\rho\hat{D}\left(\dot{\rho}\right)=\frac{\pi\rho\dot{\rho}^{2}}{\hat{m}},\qquad\hat{D}(\hat{v}_{n})=\frac{\hat{v}_{n}^{2}}{2\hat{m}},$
(49)

where $\hat{m}$ is the local interface mobility parameter.

The evolution problem can now be formulated as a minimization problem for the
global rate-potential $\hat{\Pi}$,

$\dot{\rho}=\arg\min_{\dot{\rho}}\hat{\Pi}(\dot{\rho}).$ (50)

Since $\hat{\Pi}$ is a quadratic function of $\dot{\rho}$, the condition of
stationarity of $\hat{\Pi}$, ${\rm{d}}\hat{\Pi}/{\rm{d}}\dot{\rho}=0$,
delivers a linear equation for $\dot{\rho}$ that can be easily solved. In the
general case of distinct elastic properties, the resulting expression for
$\dot{\rho}$ is rather lengthy and is not provided here for brevity.

Assuming identical elastic properties of the two phases,
$\lambda_{1}=\lambda_{2}=\lambda$ and $\mu_{1}=\mu_{2}=\mu$, the formula for
$\hat{\Psi}_{\rm el}$ simplifies considerably,

$\hat{\Psi}_{\rm
el}=\frac{2\pi\mu(\lambda+\mu)(R^{2}-\rho^{2})\rho^{2}\epsilon^{2}}{(\lambda+2\mu)R^{2}},$
(51)

and so does the resulting evolution equation,

$\dot{\rho}=-\hat{m}\left(\left(\frac{1}{2}-\left(\frac{\rho}{R}\right)^{2}\right)\frac{E\epsilon^{2}}{1-\nu^{2}}+\frac{\gamma}{\rho}\right),$
(52)

where the Lamé constants have been converted into the Young’s modulus $E$ and
Poisson’s ratio $\nu$. The above equation can be rewritten in a more concise
format,

$\dot{\rho}=-\hat{m}\left(\hat{f}_{\rm bulk}+\hat{f}_{\rm
int}\right),\qquad\hat{f}_{\rm
bulk}=\frac{E\epsilon^{2}}{1-\nu^{2}}\left(\frac{1}{2}-\left(\frac{\rho}{R}\right)^{2}\right),\qquad\hat{f}_{\rm
int}=\frac{\gamma}{\rho},$ (53)

where $\hat{f}_{\rm bulk}$ and $\hat{f}_{\rm int}$ denote the local
thermodynamic driving forces related to the elastic and interfacial energy,
respectively. It can be checked that the thermodynamic driving force
$\hat{f}=\hat{f}_{\rm bulk}+\hat{f}_{\rm int}$ derived above from the total
energy balance can be equivalently obtained from the local definition, Eq.
(6)2 (to have a consistent sign of the driving force, the interface normal
must point into the inclusion so that the inclusion is treated as phase 2 in
the notation of Section 2.1).

##  Appendix C Regularization of the phase volume fraction in laminated
elements

As discussed in Section 3.3, the formula (25) for the phase volume fraction
$\eta^{(\omega)}$ in the laminated elements is regularized according to Eq.
(32). In this appendix, we illustrate the effect of the regularization
parameter $\phi_{\rm reg}$ on the efficiency and accuracy of LET-PF.

Additional simulations of the problem studied in Section 4.1 have thus been
performed for three values of $\phi_{\rm reg}\in\\{0.001,0.01,0.1\\}$ and for
a wide range of model parameters. Fig. 1 presents representative results
concerning the efficiency and robustness (quantified by the total number of
time steps when $\Delta t_{\rm max}=T_{\rm exact}$) and concerning the
accuracy (expressed by the relative error, Eq. (35)). In the latter case,
small time increments have been used ($\Delta t_{\rm max}=T_{\rm exact}/500$)
to isolate the error related to the regularization from that resulting from
time integration. It has been observed that the dependence of both indicators
on the interfacial energy $\gamma$ is (approximately) monotonic, hence only
the extreme values $\gamma=0.0001$ and $\gamma=0.003$ are included in Fig. 1.
The computations have been performed for the coarse mesh ($h=0.02$) and for
two representative values of the interface thickness parameter
$\ell\in\\{h,1.5h\\}$.

(a) (b)

Figure 1: The effect of regularization parameter $\phi_{\rm reg}$: (a) on the
total number of time steps (for $\Delta t_{\rm max}=T_{\rm exact}$) and (b) on
the relative errer (for $\Delta t_{\rm max}=T_{\rm exact}/500$).

Fig. 1(a) shows that the number of time steps needed to complete the
simulation decreases with increasing $\phi_{\rm reg}$ for $\gamma=0.0001$, and
it only weakly depends on $\phi_{\rm reg}$ for $\gamma=0.003$. This confirms
that the robustness and efficiency of the method increases with increasing
$\phi_{\rm reg}$. Note that, in many cases, the simulations could not be
completed without regularization, i.e., for $\phi_{\rm reg}=0$. Concerning the
relative error, Fig. 1(b), the effect of $\phi_{\rm reg}$ is weak, except one
case ($\ell=1.5h$, $\gamma=0.0001$) when the error visibly increases with
increasing $\phi_{\rm reg}$.

Concluding, the results presented above, as well as other results not reported
here, do not give a definite answer on the choice of the regularization
parameter. Higher values of $\phi_{\rm reg}$ are, in general, beneficial for
robustness and computational efficiency. At the same time, the effect on
accuracy is in general weak.

## References

  * Chen [2002] L.-Q. Chen.  Phase-field models for microstructure evolution.  _Annual Review of Materials Research_ , 32(1):113–140, 2002.  doi: 10.1146/annurev.matsci.32.112001.132041.
  * Steinbach [2009] I. Steinbach.  Phase-field models in materials science.  _Modelling and Simulation in Materials Science and Engineering_ , 17(7):073001, 2009.  doi: 10.1088/0965-0393/17/7/073001.
  * Provatas and Elder [2010] N. Provatas and K. Elder.  _Phase-Field Methods in Materials Science and Engineering_.  Wiley, 2010.  doi: 10.1002/9783527631520.
  * Wang and Li [2010] Y. Wang and J. Li.  Phase field modeling of defects and deformation.  _Acta Materialia_ , 58(4):1212–1235, 2010.  doi: 10.1016/j.actamat.2009.10.041.
  * Tourret et al. [2022] D. Tourret, H. Liu, and J. LLorca.  Phase-field modeling of microstructure evolution: Recent applications, perspectives and challenges.  _Progress in Materials Science_ , 123:100810, 2022.  doi: 10.1016/j.pmatsci.2021.100810.
  * Ode et al. [2001] M. Ode, S. G. Kim, and T. Suzuki.  Mathematical modeling of iron and steel making processes. recent advances in the phase-field model for solidification.  _ISIJ International_ , 41(10):1076–1082, 2001.  doi: 10.2355/isijinternational.41.1076.
  * Ubachs et al. [2004] R. Ubachs, P. Schreurs, and M. Geers.  A nonlocal diffuse interface model for microstructure evolution of tin–lead solder.  _Journal of the Mechanics and Physics of Solids_ , 52(8):1763–1792, 2004.  doi: 10.1016/j.jmps.2004.02.002.
  * Guin and Kochmann [2023] L. Guin and D. M. Kochmann.  A phase-field model for ferroelectrics with general kinetics, part i: Model formulation.  _Journal of the Mechanics and Physics of Solids_ , 176:105301, 2023.  doi: 10.1016/j.jmps.2023.105301.
  * Wang and Khachaturyan [1997] Y. Wang and A. Khachaturyan.  Three-dimensional field model and computer modeling of martensitic transformations.  _Acta Materialia_ , 45(2):759–773, 1997.  doi: 10.1016/s1359-6454(96)00180-2.
  * Levitas and Preston [2002] V. I. Levitas and D. L. Preston.  Three-dimensional Landau theory for multivariant stress-induced martensitic phase transformations. I. Austenite$\leftrightarrow$martensite.  _Phys. Rev. B_ , 66:134206, 2002.  doi: 10.1103/PhysRevB.66.134206.
  * Xu et al. [2020] B. Xu, G. Kang, Q. Kan, C. Yu, and X. Xie.  Phase field simulation on the cyclic degeneration of one-way shape memory effect of NiTi shape memory alloy single crystal.  _International Journal of Mechanical Sciences_ , 168:105303, 2020.  doi: 10.1016/j.ijmecsci.2019.105303.
  * Tůma et al. [2021] K. Tůma, M. Rezaee-Hajidehi, J. Hron, P. Farrell, and S. Stupkiewicz.  Phase-field modeling of multivariant martensitic transformation at finite-strain: Computational aspects and large-scale finite-element simulations.  _Computer Methods in Applied Mechanics and Engineering_ , 377:113705, 2021.  doi: 10.1016/j.cma.2021.113705.
  * Clayton and Knap [2011] J. Clayton and J. Knap.  A phase field model of deformation twinning: Nonlinear theory and numerical simulations.  _Physica D: Nonlinear Phenomena_ , 240(9-10):841–858, 2011.  doi: 10.1016/j.physd.2010.12.012.
  * Liu et al. [2018] C. Liu, P. Shanthraj, M. Diehl, F. Roters, S. Dong, J. Dong, W. Ding, and D. Raabe.  An integrated crystal plasticity–phase field model for spatially resolved twin nucleation, propagation, and growth in hexagonal materials.  _International Journal of Plasticity_ , 106:203–227, 2018.  doi: 10.1016/j.ijplas.2018.03.009.
  * Rezaee-Hajidehi et al. [2022] M. Rezaee-Hajidehi, P. Sadowski, and S. Stupkiewicz.  Deformation twinning as a displacive transformation: Finite-strain phase-field model of coupled twinning and crystal plasticity.  _Journal of the Mechanics and Physics of Solids_ , 163:104855, 2022.  doi: 10.1016/j.jmps.2022.104855.
  * Bourdin et al. [2000] B. Bourdin, G. Francfort, and J.-J. Marigo.  Numerical experiments in revisited brittle fracture.  _Journal of the Mechanics and Physics of Solids_ , 48(4):797–826, 2000.  doi: 10.1016/s0022-5096(99)00028-9.
  * Ambati et al. [2014] M. Ambati, T. Gerasimov, and L. D. Lorenzis.  A review on phase-field models of brittle fracture and a new fast hybrid formulation.  _Computational Mechanics_ , 55(2):383–405, 2014.  doi: 10.1007/s00466-014-1109-y.
  * Cui et al. [2021] C. Cui, R. Ma, and E. Martínez-Pañeda.  A phase field formulation for dissolution-driven stress corrosion cracking.  _Journal of the Mechanics and Physics of Solids_ , 147:104254, 2021.  doi: 10.1016/j.jmps.2020.104254.
  * Tůma et al. [2016] K. Tůma, S. Stupkiewicz, and H. Petryk.  Size effects in martensitic microstructures: Finite-strain phase field model versus sharp-interface approach.  _Journal of the Mechanics and Physics of Solids_ , 95:284–307, 2016.  doi: 10.1016/j.jmps.2016.04.013.
  * Yeddu [2018] H. K. Yeddu.  Phase-field modeling of austenite grain size effect on martensitic transformation in stainless steels.  _Computational Materials Science_ , 154:75–83, 2018.  doi: 10.1016/j.commatsci.2018.07.040.
  * Rezaee-Hajidehi and Stupkiewicz [2020] M. Rezaee-Hajidehi and S. Stupkiewicz.  Phase-field modeling of multivariant martensitic microstructures and size effects in nano-indentation.  _Mechanics of Materials_ , 141:103267, 2020.  doi: 10.1016/j.mechmat.2019.103267.
  * Finel et al. [2018] A. Finel, Y. L. Bouar, B. Dabas, B. Appolaire, Y. Yamada, and T. Mohri.  Sharp phase field method.  _Physical Review Letters_ , 121(2), 2018.  doi: 10.1103/physrevlett.121.025501.
  * Dimokrati et al. [2020] A. Dimokrati, Y. L. Bouar, M. Benyoucef, and A. Finel.  S-PFM model for ideal grain growth.  _Acta Materialia_ , 201:147–157, 2020.  doi: 10.1016/j.actamat.2020.09.073.
  * Fleck et al. [2022] M. Fleck, F. Schleifer, and P. Zimbrod.  Frictionless motion of diffuse interfaces by sharp phase-field modeling.  _Crystals_ , 12(10):1496, 2022.  doi: 10.3390/cryst12101496.
  * Fleck and Schleifer [2022] M. Fleck and F. Schleifer.  Sharp phase-field modeling of isotropic solidification with a super efficient spatial resolution.  _Engineering with Computers_ , 39(3):1699–1709, 2022.  doi: 10.1007/s00366-022-01729-z.
  * Dobrzański et al. [2024] J. Dobrzański, K. Wojtacki, and S. Stupkiewicz.  Lamination-based efficient treatment of weak discontinuities for non-conforming finite element meshes.  _Computers & Structures_, 291:107209, 2024.  doi: 10.1016/j.compstruc.2023.107209.
  * Gélébart and Ouaki [2015] L. Gélébart and F. Ouaki.  Filtering material properties to improve FFT-based methods for numerical homogenization.  _Journal of Computational Physics_ , 294:90–95, 2015.  doi: 10.1016/j.jcp.2015.03.048.
  * Kabel et al. [2015] M. Kabel, D. Merkert, and M. Schneider.  Use of composite voxels in FFT-based homogenization.  _Computer Methods in Applied Mechanics and Engineering_ , 294:168–188, 2015.  doi: 10.1016/j.cma.2015.06.003.
  * Kabel et al. [2017] M. Kabel, A. Fink, and M. Schneider.  The composite voxel technique for inelastic problems.  _Computer Methods in Applied Mechanics and Engineering_ , 322:396–418, 2017.  doi: 10.1016/j.cma.2017.04.025.
  * Mareau and Robert [2017] C. Mareau and C. Robert.  Different composite voxel methods for the numerical homogenization of heterogeneous inelastic materials with FFT-based techniques.  _Mechanics of Materials_ , 105:157–165, 2017.  doi: 10.1016/j.mechmat.2016.12.002.
  * Keshav et al. [2022] S. Keshav, F. Fritzen, and M. Kabel.  FFT-based homogenization at finite strains using composite boxels (ComBo).  _Computational Mechanics_ , 2022.  doi: 10.1007/s00466-022-02232-4.
  * Durga et al. [2013] A. Durga, P. Wollants, and N. Moelans.  Evaluation of interfacial excess contributions in different phase-field models for elastically inhomogeneous systems.  _Modelling and Simulation in Materials Science and Engineering_ , 21(5):055018, 2013.  doi: 10.1088/0965-0393/21/5/055018.
  * Mosler et al. [2014] J. Mosler, O. Shchyglo, and H. M. Hojjat.  A novel homogenization method for phase field approaches based on partial rank-one relaxation.  _Journal of the Mechanics and Physics of Solids_ , 68:251–266, 2014.  doi: 10.1016/j.jmps.2014.04.002.
  * Schneider et al. [2015] D. Schneider, O. Tschukin, A. Choudhury, M. Selzer, T. Böhlke, and B. Nestler.  Phase-field elasticity model based on mechanical jump conditions.  _Computational Mechanics_ , 55(5):887–901, 2015.  doi: 10.1007/s00466-015-1141-6.
  * Bartels and Mosler [2017] A. Bartels and J. Mosler.  Efficient variational constitutive updates for Allen–Cahn-type phase field theory coupled to continuum mechanics.  _Computer Methods in Applied Mechanics and Engineering_ , 317:55–83, 2017.  doi: 10.1016/j.cma.2016.11.024.
  * Sethian [1999] J. A. Sethian.  _Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science_.  Cambridge Monographs on Applied and Computational Mathematics (No. 3). Cambridge University Press, Cambridge, England, 2nd edition, 1999.
  * Moës et al. [1999] N. Moës, J. Dolbow, and T. Belytschko.  A finite element method for crack growth without remeshing.  _International Journal for Numerical Methods in Engineering_ , 46(1):131–150, 1999.  doi: 10.1002/(sici)1097-0207(19990910)46:1<131::aid-nme726>3.0.co;2-j.
  * Sukumar et al. [2001] N. Sukumar, D. Chopp, N. Moës, and T. Belytschko.  Modeling holes and inclusions by level sets in the extended finite-element method.  _Computer Methods in Applied Mechanics and Engineering_ , 190(46-47):6183–6200, 2001.  doi: 10.1016/s0045-7825(01)00215-8.
  * Moës et al. [2003] N. Moës, M. Cloirec, P. Cartraud, and J.-F. Remacle.  A computational approach to handle complex microstructure geometries.  _Computer Methods in Applied Mechanics and Engineering_ , 192(28-30):3163–3177, 2003.  doi: 10.1016/s0045-7825(03)00346-3.
  * Ji et al. [2002] H. Ji, D. Chopp, and J. E. Dolbow.  A hybrid extended finite element/level set method for modeling phase transformations.  _International Journal for Numerical Methods in Engineering_ , 54(8):1209–1233, 2002.  doi: 10.1002/nme.468.
  * Duddu et al. [2011] R. Duddu, D. L. Chopp, P. Voorhees, and B. Moran.  Diffusional evolution of precipitates in elastic media using the extended finite element and the level set methods.  _Journal of Computational Physics_ , 230(4):1249–1264, 2011.  doi: 10.1016/j.jcp.2010.11.002.
  * Munk et al. [2022] L. Munk, S. Reschka, H. J. Maier, P. Wriggers, and S. Löhnert.  A sharp-interface model of the diffusive phase transformation in a nickel-based superalloy.  _Metals_ , 12(8):1261, 2022.  doi: 10.3390/met12081261.
  * Gurtin [2000] M. E. Gurtin.  _Configurational Forces as Basic Concepts of Continuum Physics_.  Springer New York, 2000.  doi: 10.1007/b97847.
  * Moës et al. [2023] N. Moës, J.-F. Remacle, J. Lambrechts, B. Lé, and N. Chevaugeon.  The eXtreme Mesh deformation approach (X-MESH) for the Stefan phase change model.  _Journal of Computational Physics_ , 477:111878, 2023.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2022.111878.
  * Hildebrand and Miehe [2012] F. Hildebrand and C. Miehe.  A phase field model for the formation and evolution of martensitic laminate microstructure at finite strains.  _Philosophical Magazine_ , 92(34):4250–4290, 2012.  doi: 10.1080/14786435.2012.705039.
  * Penrose and Fife [1990] O. Penrose and P. C. Fife.  Thermodynamically consistent models of phase-field type for the kinetic of phase transitions.  _Physica D: Nonlinear Phenomena_ , 43(1):44–62, 1990.  doi: 10.1016/0167-2789(90)90015-h.
  * Korelc [2009] J. Korelc.  Automation of primal and sensitivity analysis of transient coupled problems.  _Computational Mechanics_ , 44(5):631–649, 2009.  doi: 10.1007/s00466-009-0395-2.
  * Korelc and Wriggers [2016] J. Korelc and P. Wriggers.  _Automation of Finite Element Methods_.  Springer International Publishing, 2016.  doi: 10.1007/978-3-319-39005-5.
  * Stupkiewicz and Petryk [2002] S. Stupkiewicz and H. Petryk.  Modelling of laminated microstructures in stress-induced martensitic transformations.  _Journal of the Mechanics and Physics of Solids_ , 50(11):2303–2331, 2002.  doi: 10.1016/s0022-5096(02)00029-7.
  * Stupkiewicz [2007] S. Stupkiewicz.  _Micromechanics of Contact and Interphase Layers_.  Springer Berlin Heidelberg, 2007.  doi: 10.1007/978-3-540-49717-2.