# Sharp phase-field modeling of isotropic solidification with a super
efficient spatial resolution

Michael Fleck<EMAIL_ADDRESS>Metals and Alloys, University of
Bayreuth, Prof.-Rüdiger-Bormann-Straße 1, 95447 Bayreuth, Bavaria, Germany
Felix Schleifer Metals and Alloys, University of Bayreuth, Prof.-Rüdiger-
Bormann-Straße 1, 95447 Bayreuth, Bavaria, Germany

###### Abstract

The numerical resolution efficiency of phase-field models is limited by grid
friction, grid anisotropy and pinning. The 1D sharp phase-field model
eliminates grid friction and pinning by a global restoration of Translational
Invariance (TI) in the discretized phase-field equation (Phys. Rev. Lett. 121,
025501, 2018). In 3D global TI restricts the beneficial modeling properties to
a finite number of fixed interface orientations. We propose an accurate scheme
to restore TI locally in the local interface normal direction. At one-grid-
point interface resolutions, the new model captures the formation of isotropic
seaweed structures without spurious dendritic selection by grid anisotropy.

Diffuse interface descriptions, such as phase-field models, are widely used
for the microscopic modeling of solidification as well as related
microstructure evolution problems (KurzRappazTrievedi2021PartII, ;
TourretLiuLLorca2021, ; TonksAagesen2019, ). Quantitative simulations require
a proper numerical resolution of the diffuse solid/liquid interface, i.e. the
diffuse interface profile has to be resolved by a certain minimal amount of
grid points. In case of numerical under-resolution, the simulation is
subjected to spurious grid anisotropy as well as grid friction, which in the
worst case leads to the “pinning” of the diffuse interface on the
computational grid. In conventional phase-field models the minimal number of
grid points used to resolve the profile is about $4$
(JokisaariVoorheesGuyerWarrenHeinonen2018, ). In our notation this corresponds
to the dimensionless grid resolution number $\tilde{0pt}=2$. However,
depending on the accuracy demands of the simulation, the double, triple or
even quadruple amount of grid points can be required.

Recently, Finel et al. found a striking new way to eliminate grid friction and
pinning in one dimension, called the sharp phase-field model
(FinelLeBouarDabasAppolairYamada2018, ). This method is conceptually related
to other techniques to improve the performance of phase-field models based on
the phase-field profile function (Glasner2001, ; Weiser2009, ; Eiken2012, ;
DebierreGuerinKassner2016, ; ShenXuYang2019, ; JIMolaviTabriziKarma2022, ).
The 1D sharp phase-field model operates at one-grid-point profile resolutions
($\tilde{0pt}=0.5$) and below without the occurrence of grid pinning!

However, beside the profile resolution, there is one other important aspect
that limits the spatial resolution efficiency of phase-field models in
general: They cannot operate at arbitrarily small interface energy densities
$\Gamma$. Consider an interface between two phases at different bulk free
energy density levels. The latter, also called the driving force $\mu$,
induces an interface motion lowering the total free energy of the system. For
too small interface energies or too large driving forces either the high
energy phase turns unstable (phase stability limit) or the phase-field profile
is spuriously altered. The alternation is accompanied by strong grid friction
effects. We define the dimensionless driving force $\tilde{\mu}=\mu\Delta
x/\Gamma$, which relates to the spatial resolution of the simulation via the
grid spacing $\Delta x$. Imposing constant driving forces, we consider the
simulation of stationary interface motion in 1D at different dimensionless
spatial resolution numbers $\tilde{0pt},\tilde{\mu}$. Reasonable model
operation at the resolution $\tilde{0pt},\tilde{\mu}$ is said to require phase
stability and less than $10\%$ relative deviations from the energetically
exact interface velocity. Further details on this study are given in the
supplementary material. In Fig. 1, we compare the resulting parameter windows
of reasonable model operation for the most frequently used conventional phase-
field model (blue) and the sharp phase-field model (green). The elimination of
spurious grid friction in the sharp phase-field model allows for orders of
magnitude more efficient simulations than possible with the conventional
phase-field model.

$1$$10^{-4}$$10^{-2}$$10^{2}$$10^{4}$$1/3$$1/2$$1$$2$$4$$8$

dimensionless driving force $\tilde{\mu}$

dimensionless profile resolution $\tilde{0pt}$Sharp phase-field model

Phase stability limit

Pinning

Grid friction

conven-tionalmodel$10^{3}\times$more efficientspatial resolution3D
simulation:isotropic densebranching Figure 1: Comparison of the parameter
windows of reasonable model operation (stationary interface motion with
relative errors $<0.1$) for two different phase-field models: The most
frequently used conventional model (blue area) and the new sharp phase-field
model (green area).

During diffusion limited solidification the complex evolution of the
solid/liquid interface undergoes a branching instability (MullinsSekerka1964,
). In a fully isotropic system, this leads to the self-organized formation of
so-called isotropic dense branching or seaweed microstructures
(BrenerHMKTemkinAbel1998, ), as visible in the inset of Fig. 1. The structure
shows a characteristic distance between branches, which nontrivially relates
to the atomistically small capillary length $d_{0}$, that is proportional to
the interface energy density $\Gamma$ (KurzRappazTrievedi2021PartII, ). A
fundamental challenge in solidification modeling is the fact that the
microscopic distance between branches is typically several orders of magnitude
larger than a central aspect of its cause, i.e. the atomistically small
capillary length. If, however, the phase-field model is able to stably operate
at a certain small interface energy or, in other words, a certain large
dimensionless driving force, then the grid spacing $\Delta x$ can exceed
$d_{0}$ in a respective proportion (FleckQuerfurthGlatzel2017, ).

Here, we propose a new sharp phase-field model, which captures the 3D
formation of isotropic dense branching even at one-grid-point profile
resolutions ($\tilde{0pt}=0.5$), see Fig. 1! The absence of any spurious
dendritic selection by the computational grid indicates quite high degrees of
isotropy (Ihle2000, ; BragardKarmaLee2002, ). In this simulation the driving
forces are largely inhomogeneous. In Fig. 1, we visualize the respective
driving force distribution by a boxplot with whiskers to the maximal and
minimal value. In this work, we show that the sharp phase-field model provides
quantitative interface velocities within the full range of different driving
forces! To achieve a comparable accuracy over a similarly wide range of
driving forces, the conventional phase-field model would require profile
resolutions of $\tilde{0pt}=5$, as shown in Fig. 1. In this regard, the new
sharp phase-field model allows for 3D simulations of isotropic solidification
with a $10^{3}\times$more efficient spatial resolution.

### The new sharp phase-field model

The derivation of the new sharp phase-field formulation is started from a
discrete Helmholtz free energy functional
$F\left[\phi_{\mathbf{p}}\right]=\sum_{\mathbf{p}}f_{\mathbf{p}}\Delta x^{3}$,
where $\mathbf{p}$ denotes the locations of the grid points within the simple
cubic 3D numerical lattice with a grid spacing $\Delta x$. The discrete
Helmholtz free energy density $f_{\mathbf{p}}$ associated with the grid point
$\mathbf{p}$ is

$\displaystyle f_{\mathbf{p}}\\!=$
$\displaystyle\frac{\Gamma}{C_{\Gamma}0pt}\sum_{j,k}\\!\gamma_{j}\nu_{j}\Big{(}\frac{0pt^{2}}{2}(\partial_{k}^{+}\phi_{\mathbf{p}})^{2}\\!+\\!g_{k}(\phi_{\mathbf{p}})\Big{)}\\!+\\!\mu_{\mathbf{p}}h(\phi_{\mathbf{p}}).$
(1)

We restrict the interaction between grid points to the first three neighboring
shells $j\\!=\\!1,2,3$, with
$\left|\mathbf{r}_{k}\right|_{j}\\!=\\!\sqrt{j}\Delta x$ and $\mathbf{r}_{k}$
being a numerical lattice vector that connects two neighboring grid points
along the direction $k$. $\partial_{k}^{+}\phi_{\mathbf{p}}$ denotes the
discrete directional derivative, which is approximated by the forward finite
difference expression
$\partial_{k}^{+}\phi_{\mathbf{p}}\\!\equiv\\!(\phi_{\mathbf{p}+\mathbf{r}_{k}}\\!-\\!\phi_{\mathbf{p}})/\left|\mathbf{r}_{k}\right|$.
For a given neighboring shell with $m_{j}$ neighboring nodes, the coefficients
$\nu_{j}\\!=\\!3/m_{j}$ correct for the multiplicity of the shell. Similar to
(FinelLeBouarDabasAppolairYamada2018, ), the ponderation coefficients
$\gamma_{j}$ are chosen to get best possible energetic equality of differently
oriented ideal interfaces.

The equilibrium potentials $g_{k}(\phi)$ are minimal at $\phi\\!=\\!0$ and
$\phi\\!=\\!1$, which corresponds to the two distinct phases of the system.
$0pt$ denotes the width of the diffuse interface, $\Gamma$ is the interface
energy density, and $C_{\Gamma}$ is the interface energy calibration
parameter. A positive bulk free energy density difference $\mu_{\mathbf{p}}$
favors the growth of phase $\phi\\!=\\!0$ on the expanse of phase
$\phi\\!=\\!1$. Concerning the interpolation function $h(\phi)$, we focus on
the natural interpolation $h_{3}\\!=\\!\phi^{2}(3\\!-\\!2\phi)$
(FinelLeBouarDabasAppolairYamada2018, ) and the most frequently used
polynomial $h_{\mathrm{5}}\\!=\\!\phi^{3}(10\\!-\\!15\phi\\!+\\!6\phi^{2})$
(Plapp092011, ; OhnoTakakiShibuta2017, ; AagesenGaoSchwenAhmed2018, ;
GreenwoodShampurOforiOpokuPinomaaGurevichProvatas2018, ;
GranasyTothWarrenTegzeRatkaiPusztai2019, ; KimShermanAagesenVoorhees2020, ).

The functional phase-field derivative of the discrete Helmholtz free energy is
given by
$\delta_{\phi}F\\!=\\!\partial_{\phi}f_{\mathbf{p}}\\!-\\!\sum_{j,k}\partial_{k}^{-}(\partial_{\left(\partial_{k}\phi\right)}^{+}f_{\mathbf{p}}),$
where the second directional derivative $\partial_{k}^{-}$ is approximated by
$\partial_{k}^{-}\left(\partial
f_{\mathbf{p}}\right)\\!\equiv\\!\left(\partial f_{\mathbf{p}}\\!-\\!\partial
f_{\mathbf{p}-\mathbf{r}_{k}}\right)/\left|\mathbf{r}_{k}\right|$. The phase-
field evolution equation demands that the time derivative
$\partial_{t}\phi_{\mathbf{p}}$ is proportional to $-\delta_{\phi}F$. We
write, $30pt\Gamma\partial_{t}\phi_{\mathbf{p}}\\!=\\!-2M\delta_{\phi}F$,
where $M$ is a kinetic coefficient with the dimension
$\left[M\right]\\!=\\!\mathrm{m}^{2}\mathrm{s}^{-1}$(FleckFedermannPogorelov2018,
). During stationary interface motion, driven by a constant $\mu$, total
energy conservation demands $v_{\mathrm{th}}\\!=\\!-M\mu/\Gamma$. The phase-
field profile function is

$\displaystyle\phi_{\mathrm{\mathbf{p}}}$ $\displaystyle=\left(1-\tanh
2\left(\mathbf{p}\cdot\mathbf{n}-c_{n}\right)/0pt\right)/2,$ (2)

which is an analytic solution of the continuum phase-field equation, if
$g(\phi)\\!=\\!\sum_{j,k}g_{k}\\!\equiv\\!8\phi^{2}\left(1\\!-\\!\phi\right)^{2}$and
$h(\phi)\\!=\\!h_{3}\\!=\\!\phi^{2}(3\\!-\\!2\phi)$. $\mathbf{n}$ is the unit
normal interface vector and $c_{n}\\!=\\!v_{\mathrm{th}}t$ denotes the central
interface position, moving with the velocity $v_{\mathrm{th}}$. The profile
width of $20pt$ is understood as $96.4\%$ of the total transition from
$\phi\\!=\\!0$ to $\phi\\!=\\!1$ ($\tanh 2\\!\simeq\\!0.964$)
(DimokratiLeBouarBenyoucefFinel2020, ).

For vanishing driving forces $\mu\\!=\\!0$ and no phase-field motion
$\partial_{t}\phi\\!=\\!0$ and the phase-field equation reduces to

$\displaystyle{\textstyle\sum_{j,k}}\gamma_{j}\nu_{j}\\{0pt^{2}(\phi_{\mathbf{p}+\mathbf{r}_{k}}\\!\\!-\\!2\phi_{\mathbf{p}}\\!+\\!\phi_{\mathbf{p}-\mathbf{r}_{k}})/\mathbf{r}_{k}^{2}-\partial_{\phi}g_{k}\\}$
$\displaystyle=0,$ (3)

where $\partial_{\phi}\\!=\\!\partial/\partial\phi$ denotes the partial phase-
field derivative. The condition holds, if all individual $k-$components are
simultaneously satisfied. Those can be satisfied at any real time during the
propagation of the interface using the addition property of the hyperbolic
tangent profile (2)
$\phi_{\mathbf{p}\pm\mathbf{r}_{k}}\\!=\\!(1\\!\pm\\!a_{k})\phi_{\mathbf{p}}/(1\\!\pm\\!(2\phi_{\mathbf{p}}\\!-\\!1)a_{k}),$
where the grid coupling parameters $a_{k}\\!\left(\mathbf{n}\right)$ are
defined as
$a_{k}\\!=\\!\tanh\left(2\mathbf{r}_{k}\\!\cdot\\!\mathbf{n}/0pt\right).$
Inserting this property into the phase-field equilibrium condition, we obtain
the $k-$th component of the modified equilibrium potential

$\displaystyle
g_{k}(\phi)\frac{\mathbf{r}_{k}^{2}}{0pt^{2}}=\phi(1\\!-\\!\phi)\,+\,$
$\displaystyle\frac{1\\!-\\!a_{k}^{2}}{4a_{k}^{2}}\ln\bigg{(}\frac{1\\!-\\!a_{k}^{2}}{1\\!-\\!a_{k}^{2}\left(1\\!-\\!2\phi\right)^{2}}\bigg{)},$
(4)

which further satisfies $g_{k}(\phi\\!=\\!0,\\!1)\\!=\\!0$, to allow an easy
calculation of the system’s total interface energy by
$F_{\mathrm{int}}(\phi_{\mathbf{p}})\\!=\\!\sum_{\mathbf{p}}f_{\mu=0}$ using
an arbitrary phase-field (SchleiferHolzingerLinGlatzelFleck2019, ;
SchleiferFleckHolzingerLinGlatzel_superalloys2020, ). In the continuum limit
$\left|\mathbf{r}_{k}\right|\\!\rightarrow\\!0$, Eq. (4) converges to the
conventional Continuum Field (CF) potential
$g_{k}^{\infty}\\!=\\!8\phi^{2}\left(1\\!-\\!\phi\right)^{2}$.

$1$$10^{-12}$$10^{-08}$$10^{-04}$$0$$20$$40$$60$$20$$40$$60$

amplitude $A$

$\vartheta_{[001]}$CFTI⟨110⟩$\vartheta_{[011]}$TI⟨n⟩TI⟨111⟩TI⟨100⟩

Figure 2: Test of Translational Invariance (TI) of the ideal profile (2)
within the equilibrium condition (3). We plot the oscillation amplitude $A$ of
the system integral over Eq. (3) during profile motion for different interface
orientation angles $\vartheta_{[001]}$ and $\vartheta_{[011]}$. Profile
resolution $\tilde{0pt}\\!=\\!0pt/\Delta x\\!=\\!0.5$; system size
$300\\!\times\\!1\\!\times\\!1$.

Translational Invariance (TI) in the phase-field equation is restored based on
properties of the profile function (2). Without TI, the system integral over
Eq. (3) oscillates, when the ideal profile (2) is moved on the grid, as
plotted in Fig. 2. To determine the grid coupling parameters
$a_{k}\\!\left(\mathbf{n}\right)$ Finel et al. proposed to represent the
interface normal vector $\mathbf{n}$ by a constant unit vector $\mathbf{u}$,
perpendicular to a properly chosen lattice plane
(FinelLeBouarDabasAppolairYamada2018, ). This globally restores TI for
interface orientations that agree to one of the equivalent lattice
orientations $\left\langle\mathbf{u}\right\rangle$, as shown in Fig. 2 for
different TI⟨u⟩-models. Further details are given in the supplementary
material. The newly proposed TI⟨n⟩-model (green curve) uses grid coupling
parameters calculated from the local interface normal direction, leading to
very small oscillations regardless of the interface orientation.

For the sufficiently accurate determination of
$a_{k}\\!\left(\mathbf{n}\right)$, we proceed as follows: First, we calculate
preliminary grid coupling parameters by
$\hat{a}_{k}\\!=\\!(\hat{a}_{k}^{+}\\!+\\!\hat{a}_{k}^{-})/2$, where

$\displaystyle\hat{a}_{k}^{\pm}$
$\displaystyle=\frac{\pm\left(\phi_{\mathbf{p}\pm\mathbf{r}_{k}}-\phi_{\mathbf{p}}\right)}{\phi_{\mathbf{p}}-2\phi_{\mathbf{p}\pm\mathbf{r}_{k}}\phi_{\mathbf{p}}+\phi_{\mathbf{p}\pm\mathbf{r}_{k}}}.$
(5)

Using the modified equilibrium potentials the explicit dependence of the
phase-field equation on the profile width $0pt$ cancels out. Then, $0pt$ is
solely controlled by the preliminary grid coupling parameters, which also
contain the a priori unknown interface normal vector $\hat{\mathbf{n}}$. Thus,
without length control of $\hat{\mathbf{n}}$ the profile width $0pt$ wouldn’t
be defined in the model. Thus, we locally calculate all components of the
interface normal vector
$\hat{n}_{k}\\!=\\!0pt\mathrm{arctanh}(\hat{a}_{k})/\left|2\mathbf{r}_{k}\right|$,
restore unit length via
$\mathbf{n}\\!=\\!\hat{\mathbf{n}}/\left|\hat{\mathbf{n}}\right|$ and
calculate corrected grid coupling parameters $a_{k}\left(\mathbf{n}\right)$.

The advancing solidification is accompanied by a release of latent heat at the
solid/liquid interface (KassnerGuerinDucouss082010, ). Thus, the dimensionless
temperature field $U_{\mathbf{p}}\\!=\\!C(T_{\mathbf{p}}\\!-\\!T_{M})/L$ is
introduced, where $T_{M}$, $L$ and $C$ denote the melting temperature, latent
heat and heat capacity, respectively (FleckHuterPilipen012010, ). The driving
force for solidification is given by
$\mu_{\mathbf{p}}\\!=\\!-U_{\mathbf{p}}\Gamma/d_{0}$, where
$d_{0}\\!=\\!\Gamma T_{M}C/L^{2}$ denotes the capillary length. The
temperature obeys a diffusion equation,
$\partial_{t}U_{\mathbf{p}}\\!=\\!D\nabla^{2}U_{\mathbf{p}}\\!+\\!R(\phi_{\mathbf{p}})\partial_{\phi}h\partial_{t}\phi_{\mathbf{p}}$,
with equal diffusion coefficients $D$ in the solid and liquid phase. For small
phase-field widths, $0pt/\Delta x\\!\leq\\!2$, and $R\\!=\\!1$, we observe
spuriously inhomogeneous releases of latent heat, whenever a grid point is
close to the interface center. The spurious heat release provides oscillations
in the solidification velocity as well as some degree of kinetic anisotropy.
Therefore, we propose the regularization $R(\phi_{\mathbf{p}})$ in the
diffusion equation

$\displaystyle R(\phi_{\mathbf{p}})$
$\displaystyle=\frac{3C_{R}}{4}\frac{a_{\left\langle
100\right\rangle}0pt}{\Delta x}\left(1-a_{\left\langle
100\right\rangle}^{2}\left(1-2\phi_{\mathbf{p}}\right)^{2}\right)^{-2},$ (6)

where the grid coupling parameter is $a_{\left\langle
100\right\rangle}\\!=\\!\tanh 2\Delta x/0pt$ and $C_{R}$ denotes a calibration
constant, which is required to maintain total energy conservation during
solidification. The dependence of $C_{R}$ as a function of the dimensionless
profile resolution is plotted in Fig. 3.

### Model calibration

The interface energy calibration $C_{\Gamma}$ is calculated via
$C_{\Gamma}\\!=\\!\sum_{\mathbf{p}_{[100]}}\\!\mathbf{e}_{[100]}\\!\cdot\\!\mathbf{n}f(\phi_{\mathbf{p}}(\mathbf{n}))_{\mu=0}/\Gamma$,
where $\mathbf{e}_{[100]}$ denotes a unit vector pointing in the
$[100]-$direction of the computational grid, $\sum_{\mathbf{p}_{[100]}}$
denotes the sum in the $[100]-$direction, $\mathbf{n}$ is again the direction
normal to the interface, and the phase-field values
$\phi_{\mathbf{p}}(\mathbf{n})$ are given by the ideal profile (2) with
orientation $\mathbf{n}$. For the determination of $C_{\Gamma}$ we chose the
$[100]-$direction as interface orientation. The determination of the energy
calibration factor is independent from the choice of the ponderation
coefficients. Fig. 3 shows the phase-field width dependence of the different
calibration factors. The continuum limit for the calibration factor,
$C_{\Gamma}^{\infty}\\!=\\!2/3$, is indicated by the solid black line in Fig.
3. For sharp diffuse interfaces with a phase-field width below
$\tilde{0pt}\\!<\\!2$, we obtain substantially smaller values for the
calibration line integral as compared to the limiting value.

$0.2$$0.4$$0.6$$0$$1$$2$$3$$4$$5$

calibration parameter

dimensionless profile resolution
$\tilde{0pt{}}$$C_{\Gamma{}}^{\mathrm{TI}}$$\gamma_{2}$$\gamma_{3}$TI⟨100⟩TI⟨n⟩continuum
limit$C_{\Gamma{}}^{\mathrm{CF}}$$C_{R}$ Figure 3:  Plot of the different
calibration parameters $C_{\Gamma}$ (solid green), $C_{\Gamma}^{\mathrm{CF}}$
(dashed green), $C_{R}$ (red) and the ponderation coefficients $\gamma_{2}$
(violet) and $\gamma_{3}$ (blue) as a function of the phase-field width
$\tilde{0pt}$. $\gamma_{1}\\!=\\!1\\!-\\!\gamma_{2}\\!-\\!\gamma_{3}$

For the determination of the ponderation coefficients an optimization
procedure similar to the one proposed by Finel et al.
(FinelLeBouarDabasAppolairYamada2018, ) has been developed. The ponderation
coefficients $\gamma_{j}$ should be chosen such that the interface energy
becomes as isotropic as possible, i.e. the discrete interface energy integral
$\Gamma(\mathbf{n})\\!=\\!\big{<}\sum_{\mathbf{p}_{[100]}}\\!\mathbf{e}_{[100]}\\!\cdot\\!\mathbf{n}f(\phi_{\mathbf{p}}(\mathbf{n}))_{\mu=0}\big{>}_{r_{n}}$
should dependent on the interface orientation $\mathbf{n}$ as little as
possible. Since at least some of these line integrals may not be
Translationally Invariant, we further average over a number of different
values obtained for different positions $r_{n}$ of the interface center, as
denoted by the angle brackets with index $r_{n}$. Given a starting set for the
ponderation coefficients $\gamma_{j}$, we calculate the following three
different interface energy densities: $\Gamma([100])\\!=\\!\Gamma_{j=1}$,
$\Gamma([110])\\!=\\!\Gamma_{j=2}$ and $\Gamma([111])\\!=\\!\Gamma_{j=3}$. As
a measure for interface energy isotropy and as the minimization target, the
square root of the sum of the deviations from the average interface energy
value in square of these three energy densities is chosen. i.e.

$\displaystyle\left\\{\gamma_{j}\right\\}:$
$\displaystyle\min\sqrt{{\textstyle\sum_{j}}(\overline{\Gamma}-\Gamma_{j})^{2}},$
(7)

with $\overline{\Gamma}\\!=\\!\sum_{j}\Gamma_{j}/3$. The optimal choice for
the ponderation coefficients $\left\\{\gamma_{j}\right\\}$, with respect to
this minimization target and under the constraint
$\sum_{j}\gamma_{j}\\!=\\!1$, has been calculated by a simple steepest decent
algorithm. In Fig. 3, the optimal ponderation coefficients are plotted as
function of the phase-field width for the TI⟨n⟩-model (solid curves) as well
as for the TI⟨100⟩-model (dashed curves). The ponderation coefficients
obtained for the CF-model are nearly identical to those of the TI⟨100⟩-model.

$1.0$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$-6$$-4$$-2$$0$$2$$4$$6$$0.0$$0.5$$1.0$

heat release

$\phi(x,t)$

coordinate $\tilde{x}=x/\Delta x$$S(\phi)$:
$\tilde{0pt}=2$$R(\phi)S(\phi)$:$\tilde{0pt}=0.4$ Figure 4: Comparison of the
heat release range of the source term
$S(\phi)\\!=\\!\partial_{\phi}h\partial_{t}\phi\\!\sim\\!\partial_{\phi}h^{2}(\phi)$
with ($\tilde{0pt}\\!=\\!0.4$) and without ($\tilde{0pt}\\!=\\!2$) the
regularization factor $R(\phi)$ Eq. (6). The comparison is based on the ideal
phase-field profile function (2)

The idea behind the source term regularization $R(\phi)$ Eq. (6) is to
distribute the latent heat release over a slightly enlarged range, involving
more than just a single grid point. The different ranges of heat release of
different source term variants are compared in Fig. 4. The regularization
requires a phase-field width dependent calibration procedure. For a given
phase-field width the calibration parameter $C_{R}$ has to take a specific
value in order to ensure the conservation of the total energy in the system.
Using some arbitrary starting value for $C_{R}$, we perform a long term
simulations of solidification until quasi two phase equilibrium in a small,
thermally isolated, one-dimensional system with an initially homogeneous
undercooling temperature of $U_{0}\\!=\\!-0.7$. Then, based of the deviation
of the measured solid phase fraction from the expected outcome of $0.7$, we
successively optimize the $C_{R}$ value.

### Quantitative stationary solidification

$20$$30$$40$$50$$60$TI$+h_{3}$CF$+h_{5}$$-2.0$$-1.5$$-1.0$

solidification

direction

coordinate $\tilde{x}$

$U(x,t)$

solidTI$+h_{3}+R$$0.0$$0.5$$1.0$$25$$30$$0.01$$0.1$meanvalueoscillationamplitudea)Error
description

${v}/v_{\mathrm{th}}$

interface center $\tilde{c}_{n}(t)$

$|1-v_{\mathrm{fit}}/v_{\mathrm{th}}|$

$0.01$$0.1$$1$$1$$100$$10^{4}$b)

$|1-v_{\mathrm{fit}}/v_{\mathrm{th}}|$

dimensionless driving force $\tilde{\mu{}}_{\mathrm{min}}$TI$+h_{3}$:
$\tilde{0pt}=0.4$TI$+h_{3}+R$CF$+h_{5}$: $\tilde{0pt}=2.0$$\tilde{0pt}=0.4$

Figure 5: Stationary solidification using (i) the Continuum Field model
(CF$+h_{5}$) for $\tilde{0pt}\\!=\\!2$ in blue, (ii) the Translationally
Invariant model (TI$+h_{3}$) for $\tilde{0pt}\\!=\\!0.4$ in red and (iii) the
TI-model with regularization (TI$+h_{3}\\!+\\!R$) in green. a) Exemplary
simulation results and a plot of the velocity as function of the interface
center ($\tilde{\mu}_{\mathrm{int}}\\!=\\!100$). The temperature $U$ is given
by colored lines and the phase-field values by black full symbols. b) Plot of
the interface velocity error as function of the dimensionless driving force
$\tilde{\mu}_{\mathrm{int}}\\!=\\!\mu_{\mathrm{int}}\Delta x/\Gamma$.

In Fig. 5a), the configuration of stationary solidification is shown. An
animation of this figure is provided in the supplementary material. Far in
front of the solid/liquid interface the temperature is $U(L)\\!=\\!-2.0$. When
the system reaches a stationary state, the solid phase is found at the minimal
undercooling temperature of $U_{\mathrm{int}}\\!=\\!-1.0$. Then, the
theoretically expected solidification velocity is given by
$v_{\mathrm{th}}\\!=\\!MU_{\mathrm{int}}/d_{0},$ where $M$ denotes the kinetic
coefficient, and $d_{0}$ is the capillary length (Caginalp1989, ). We restrict
to the comparison with the sharp interface equation and omit more
sophisticated thin interface corrections (KarmaRappel041998, ). The ratio
between the total system length and the theoretic stationary diffusion length
$l_{D}\\!=\\!2D/v_{\mathrm{th}}$ is chosen to be $L/l_{D}\\!=\\!5$. The system
is resolved by $200$ grid points, i.e. $L/\Delta x\\!=\\!200$, with a solid
phase fraction of $12\%$. The fraction is kept constant by incremental
shifting of the whole system (FleckBrenerSpatsch2010, ).

In Fig. 5b) the relative error in the solidification velocity is plotted as
function of the dimensionless driving force
$\tilde{\mu}_{\mathrm{int}}\\!=\\!\mu_{\mathrm{int}}\Delta x/\Gamma$. The
CF$+h_{5}$-model (blue color) is subjected to strong spurious grid friction
for both small as well as large dimensionless driving forces. In case of
$\tilde{\mu}_{\mathrm{int}}\\!=\\!100$, the observed solidification velocity
is $90\>\%$ smaller than the expectation. The TI-models are limited by phase
stability only. This limit is indicated by the vertical dashed line in Fig.
5b). The TI$+h_{3}$-model (red curve) provides large oscillations in the
interface velocity. These result from spuriously inhomogeneous heat release at
the solid/liquid interface, as visible in Fig. 5a). It can be avoided by
employing the newly proposed source term regularization $R$ Eq. (6), see the
green curves in Fig. 5.

### Diffusion limited solidification

initial state
|  time $tM/\Delta x^{2}$: | $5$ | $50$ | $\geq 100$
---|---|---|---|---
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}\\!+\\!R$ |  |  |
$\tilde{0pt}=0.5$ |  |  |  |
|  TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}$ |  |  |
$\tilde{0pt}=0.5$ |  |  |  |
|  TI${}_{\left\langle 100\right\rangle}\\!+\\!h_{3}$ |  |  |
$\tilde{0pt}=0.5$ |  |  |  |
|  CF$\\!+\\!h_{5}$ |  |  |
$\tilde{0pt}=2.0$ |  |  |  |
Figure 6:  Time series of phase-field simulations of diffusion limited
solidification using four different models: The
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}$ model (i) with and
(ii) without regularization $R$, (iii) the TI${}_{\left\langle
100\right\rangle}\\!+\\!h_{3}$ model each with $\tilde{0pt}\\!=\\!0.5$, and
(iv) the CF$+h_{5}$-model with $\tilde{0pt}\\!=\\!2.0$. The temperature $U$ is
visualized by the coloring and the phase-field is represented by the
$\phi\\!=\\!1/2-$contour. Further parameters: $d_{0}/\Delta x\\!=\\!2\cdot
10^{-3}$, $D/M\\!=\\!5\\!\cdot\\!10^{-3}$, domain size
$120\\!\times\\!60\\!\times\\!60$.

For dimensionless undercooling temperatures $U$ smaller than unity, we obtain
diffusion limited solidification. Four comparable simulations are performed
using four different phase-field models, as shown in Fig. 6. An animation
showing the full courses of all four simulations is provided in the
supplementary material. The simulations are started from the same initial
state at $U=-0.3$. All boundaries are thermally insulating, except for the
boundary at the $\left[100\right]-$end of the simulation domain on the right
hand side, which is held at $U_{\mathrm{max}}=-0.3$. The initial quasi planar
solid/liquid interface has small bumps at regular intervals of $10\Delta x$.
In the beginning, the interface develops the Mullins-Sekerka instability
(MullinsSekerka1964, ), since the dimensionless capillary length is chosen to
be sufficiently small $\tilde{d_{0}}\\!=\\!0.002$
($\tilde{\mu}_{\mathrm{max}}\\!=\\!U_{\mathrm{max}}/\tilde{d_{0}}\\!=\\!150$).
As soon as the most advanced point of the solid/liquid interface exceeds the
fraction of $0.7$ of the simulation domain along the
$\left[100\right]-$direction, the whole system is shifted back by one grid
point (FleckBrenerSpatsch2010, ).

In later stages, the disordered seaweed or dense-branching morphology develops
(IhleMuel1994, ; BrenerHMKTemkinAbel1998, ; UtterBodenschatz2005, ), if the
residual grid anisotropy is sufficiently small. For super efficient one-grid-
point profile resolutions of $\tilde{0pt}\\!=\\!0.5$, this requires the local
restoration of TI in the local interface normal direction as well as the
inclusion of the source term regularization Eq. (6), as shown in first row in
Fig. 6. Without regularization the simulation shows a spurious dendritic
selection in the $\left\langle 110\right\rangle-$directions of the
computational grid, which originates from the inhomogeneous temperature
release via the source term in the diffusion equation. The simulations using
the TI⟨100⟩\- and CF- model show spurious dendritic selection in the
$\left\langle 100\right\rangle-$directions. For the TI⟨100⟩ model, the
selection originates from anisotropic interface kinetics (Ihle2000, ;
BragardKarmaLee2002, ), which result from residual grid friction for interface
orientations that differ from the $\left\langle 100\right\rangle-$directions
(FleckSchleifer2021, ). In case of the CF-model, it results from strong grid
friction.

### Conclusion

A new sharp phase-field model is proposed: Instead of using global grid
dependent equilibrium potentials (4), that restore the Translational
Invariance (TI) for a finite amount of fixed interface orientations, the newly
proposed model restores TI locally for the local interface normal direction
$\mathbf{n}$. Furthermore, we propose a source term regularization Eq. (6) to
effectively suppress spurious inhomogeneous temperature releasees by diffuse
interfaces as sharp as $\tilde{0pt}\\!=\\!0.4$, see Fig. 5. Compared to the
conventional phase-field model with the resolution limits $\tilde{0pt}>2.0$
and $\tilde{\mu}<1.0$, the sharp phase-field model allows for super efficient
quantitative simulations of stationary solidification with phase-field profile
resolutions of $\tilde{0pt}\\!=\\!0.4$ and dimensionless driving forces up to
$\tilde{\mu}\\!=\\!7200$! The new sharp phase-field model with source term
regularization
(TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}\\!+\\!R$) provides
extremely high degrees of isotropy. It provides the expected isotropic seaweed
or dense-branching morphology using extraordinary efficient spatial
resolutions: $\tilde{0pt}\\!=\\!0.5$ and
$\tilde{\mu}_{\mathrm{max}}\\!=\\!U_{\mathrm{max}}/\tilde{d_{0}}\\!=\\!150$!

###### Acknowledgements.

We thank B. Böttger and J. Eiken from ACCESS, Aachen, Germany as well as A.
Finel from ONERA, Châtillon, France for fruit-full discussions on this issue.
The work is funded by the Deutsche Forschungsgemeinschaft (DFG) – 431968427.

© 2022. This version of the manuscript is made available under the CC-BY 4.0
license http://creativecommons.org/licenses/by/4.0/

## References

  * (1) W. Kurz, M. Rappaz, and R. Trivedi.  Progress in modelling solidification microstructures in metals and alloys. part ii: dendrites from 2001 to 2018.  _Int. Mater. Rev._ 66, 30–76 (2021), [link].
  * (2) D. Tourret, H. Liu, and J. LLorca.  Phase-field modeling of microstructure evolution: Recent applications, perspectives and challenges.  _Prog. Mater. Sci._ 100810 (2021), [link].
  * (3) M.R. Tonks and L.K. Aagesen.  The phase field method: Mesoscale simulation aiding material discovery.  _Annu. Rev. Mater. Res._ 49, 79–102 (2019), [link].
  * (4) A.M. Jokisaari, P. W. Voorhees, J. E. Guyer, J. A. Warren, and O. Heinonen.  Phase field benchmark problems for dendritic growth and linear elasticity.  _Comp. Mater. Sci._ 149, 336–347 (2018), [link].
  * (5) A. Finel, Y. Le Bouar, B. Dabas, B. Appolaire, Y. Yamada, and T. Mohri.  Sharp phase field method.  _Phys. Rev. Lett._ 121, 025501 (2018), [link].
  * (6) K. Glasner.  Nonlinear preconditioning for diffuse interfaces.  _J. Comp. Phys._ 174, 695–711 (2001), [link].
  * (7) M. Weiser.  Pointwise nonlinear scaling for reaction–diffusion equations.  _Appl. Num. Math._ 59, 1858–1869 (2009), [link].
  * (8) J. Eiken.  Numerical solution of the phase-field equation with minimized discretization error.  _IOP Conf. Ser.: Mater. Sci. Eng._ 33, 012105 (2012), [link].
  * (9) Jean-Marc Debierre, Rahma Guérin, and Klaus Kassner.  Phase-field study of crystal growth in three-dimensional capillaries: Effects of crystalline anisotropy.  _Phys. Rev. E_ 94, 013001 (2016), [link].
  * (10) J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  _SIAM Review_ 61, 474–506 (2019), [link].
  * (11) K. Ji, A. Molavi Tabrizi, and A. Karma.  Isotropic finite-difference approximations for phase-field simulations of polycrystalline alloy solidification.  _J Comput. Phys._ 111069 (2022), [link].
  * (12) W.W. Mullins and R.F. Sekerka.  Stability of a planar interface during solidification of a dilute binary alloy.  _J. Appl. Phys._ 35, 444–451 (1964), [link].
  * (13) E.A. Brener, H. Müller-Krumbhaar, D.E. Temkin, and T. Abel.  Morphology diagram of possible structures in diffusional growth.  _Physica A_ 249, 73–81 (1998), [link].
  * (14) M. Fleck, F. Querfurth, and U. Glatzel.  Phase field modeling of solidification in multi-component alloys with a case study on the Inconel 718 alloy.  _J. Mater. Res._ 32, 4605–4615 (2017), [link].
  * (15) T. Ihle.  Competition between kinetic and surface tension anisotropy in dendritic growth.  _Euro. Phys. J B_ 16, 337–344 (2000), [link].
  * (16) J. Bragard, A. Karma, Y. H. Lee, and M. Plapp.  Linking Phase-Field and Atomistic Simulations to Model Dendritic Solidification in Highly Undercooled Melts.  _Interf. Sci._ 10, 121 (2002), [link].
  * (17) M. Plapp.  Unified derivation of phase-field models for alloy solidification from a grand-potential functional.  _Phys. Rev. E_ 84, 031601 (2011), [link].
  * (18) M. Ohno, T. Takaki, and Y. Shibuta.  Variational formulation of a quantitative phase-field model for nonisothermal solidification in a multicomponent alloy.  _Phys. Rev. E_ 96, 033311 (2017), [link].
  * (19) L.K. Aagesen, Y. Gao, D. Schwen, and K. Ahmed.  Grand-potential-based phase-field model for multiple phases, grains, and chemical components.  _Phys. Rev. E_ 98, 023309 (2018), [link].
  * (20) M. Greenwood, K. N. Shampur, N. Ofori-Opoku, T. Pinomaa, L. Wang, S. Gurevich, and N. Provatas.  Quantitative 3d phase field modelling of solidification using next-generation adaptive mesh refinement.  _Comp. Mater. Sci._ 142, 153 (2018), [link].
  * (21) L. Gránásy, G. I. Tóth, J. A. Warren, F. Podmaniczky, G. Tegze, L. Rátkai, and T. Pusztai.  Phase-field modeling of crystal nucleation in undercooled liquids – a review.  _Progress in Materials Science_ 106, 100569 (2019), [link].
  * (22) K. Kim, Q.C. Sherman, L.K. Aagesen, and P.W. Voorhees.  Phase-field model of oxidation: Kinetics.  _Phys. Rev. E_ 101, 022802 (2020), [link].
  * (23) M. Fleck, H. Federmann, and E. Pogorelov.  Phase-field modeling of li-insertion kinetics in single LiFePO4-nano-particles for rechargeable li-ion battery application.  _Comp. Mater. Sci._ 153, 288–296 (2018), [link].
  * (24) A. Dimokrati, Y. Le Bouar, M. Benyoucef, and A. Finel.  S-pfm model for ideal grain growth.  _Acta Mater._ 201, 147–157 (2020), [link].
  * (25) F. Schleifer, M. Holzinger, Y.-Y. Lin, U. Glatzel, and M. Fleck.  Phase-field modeling of a $\gamma$/$\gamma^{\prime\prime}$ microstructure in nickel-base superalloys with high $\gamma^{\prime\prime}$ volume fraction.  _Intermetallics_ 120, 106745 (2020), [link].
  * (26) F. Schleifer, M. Fleck, M. Holzinger, Y.-Y. Lin, and U. Glatzel.  _Phase-field modeling of $\gamma^{\prime}$ and $\gamma^{\prime\prime}$ precipitate size evolution during heat treatment of Ni-base superalloys_, chap. 49, 500–508.  Superalloys 2020 (Springer International Publishing, Cham, 2020), [link].
  * (27) K. Kassner, R. Guérin, T. Ducousso, and J.-M. Debierre.  Phase-field study of solidification in three-dimensional channels.  _Phys. Rev. E_ 82, 021606 (2010), [link].
  * (28) M. Fleck, C. Hüter, D. Pilipenko, R. Spatschek, and E. A. Brener.  Pattern formation during diffusion limited transformations in solids.  _Phil. Mag._ 90, 265 (2010), [link].
  * (29) G. Caginalp.  Stefan and hele-shaw type models as asymptotic limits of the phase-field equations.  _Phys. Rev. A_ 39, 5887–5896 (1989), [link].
  * (30) A. Karma and W.-J. Rappel.  Quantitative phase-field modeling of dendritic growth in two and three dimensions.  _Phys. Rev. E_ 57, 4323–4349 (1998), [link].
  * (31) M. Fleck, E. A. Brener, R. Spatschek, and B. Eidel.  Elastic and plastic effects on solid-state transformations: A phase field study.  _Int. J. Mater. Res._ 4, 462 (2010), [link].
  * (32) T. Ihle and H. Müller-Krumbhaar.  Fractal and compact growth morphologies in phase transitions with diffusion transport.  _Phys. Rev. E_ 49, 2972–2991 (1994), [link].
  * (33) B. Utter and E. Bodenschatz.  Double dendrite growth in solidification.  _Phys. Rev. E_ 72, 011601 (2005), [link].
  * (34) M. Fleck and F. Schleifer.  Frictionless motion of marginally diffuse interfaces by sharp phase-field modelling (2022), [link].  Revision in preparation.

Supplementary material for: Sharp phase-field modeling of isotropic
solidification with a super efficient spatial resolution”

In the supplementary material, we provide additional technical information
about the modeling and the simulations presented in the manuscript.

##  .1 Description of the supplementary animations

  1. 1.

Supplementary_material_1_stationary_solidification.mpg: This movie is an
animation showing the different simulations of stationary solidification
presented in Figure 5. It shows the time evolution of three different
simulations of stationary solidification, for a dimensionless undercooling
temperature of $U_{\mathrm{max.}}\\!=\\!-2$ and a dimensionless driving force
of $\tilde{\mu}_{\mathrm{min}}\\!=\\!100$.

  2. 2.

Supplementary_material_2_diffusion_limited_solidification.mpg: This animation
shows the four different simulations of diffusion limited solidification
presented in Figure 6 using the four different models: (i)
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}\\!+\\!R$ with
$\tilde{0pt}\\!=\\!0.5$, (ii)
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}$ with
$\tilde{0pt}\\!=\\!0.5$, (iii) TI${}_{\left\langle
100\right\rangle}\\!+\\!h_{3}$ with $\tilde{0pt}\\!=\\!0.5$ and (iv)
CF$+\\!h_{5}$ with $\tilde{0pt}\\!=\\!2.0$. The temperature $U$ is indicated
by the coloring and the phase-field is visualized via a
$\phi\\!=\\!1/2-$contour plot. The capillary length is $d_{0}/\Delta
x\\!=\\!0.002$, undercooling: $U(\mathbf{x},0)\\!=\\!-0.3$, kinetics:
$D/M\\!=\\!0.01$.

##  .2 Stationary interface motion

velocity $|1-v_{\mathrm{fit}}/v_{\mathrm{th}}|$

$0.01$$1$$10^{-6}$$10^{-4}$

$1$$10^{-3}$$10^{3}$TI$+h_{3}$$h_{5}$$\tilde{0pt{}}=4.0$

$1$$10^{-3}$$10^{3}$$\tilde{0pt{}}=3.0$

$1$$10^{-3}$$10^{3}$$\tilde{0pt{}}=2.5$

$1$$10^{-3}$$10^{3}$$\tilde{0pt{}}=0.5$

$0.01$$1$$10^{-6}$$10^{-4}$

width $|1-0pt_{\mathrm{fit}}/0pt|$

$0.01$$1$$10^{-6}$$10^{-4}$

$1$$10^{-3}$$10^{3}$TI$+h_{3}$$h_{5}$

$1$$10^{-3}$$10^{3}$

$1$$10^{-3}$$10^{3}$

$1$$10^{-3}$$10^{3}$

$0.01$$1$$10^{-6}$$10^{-4}$

$0.1$$1$$10$$10^{-4}$$10^{-3}$$10^{-2}$$10^{2}$$10^{3}$$10^{4}$parameter
windows ofresonable model
operation$\tilde{0pt{}}=0.4$$\tilde{0pt{}}=0.5$$\tilde{0pt{}}=2.5$$\tilde{0pt{}}=3.0$$\tilde{0pt{}}=4.0$

dimensionless driving force $\tilde{\mu}$

Figure 7: Errors plots of the stationary interface velocity (top row) and the
fitted interface width (middle row) as a function of the dimensionless driving
force $\tilde{\mu}=\mu\Delta x/\Gamma$, for different phase-field widths:
$\tilde{0pt}=0pt/\Delta x=4.0,\;3.0,\;2.5,\;0.5$. Two models are compared: (i)
Continuum Field (CF) model with $h_{5}$ (blue) and (ii) the sharp phase-field
model with Translational Invariance (TI$+h_{3}$) (green). Solid lines denote
the mean relative errors and the oscillations are indicated as transparently
colored areas. The time resolution is $M\mu\Delta t/(\Gamma\Delta x)=1.6\cdot
10^{-7}$.

We consider the constantly driven stationary motion of a planar interface in
one dimension. In Fig. 7, we compare mean errors in the interface velocities
and widths (solid lines) as well as their relative oscillation amplitudes
(colored areas) for different models. As illustrated in the lower panel of
Fig. 5a), the colored areas start from the oscillation amplitude value and end
at the mean value. When the colored area is found above the mean value, we
have the “healthy” situation that the measured value oscillates around the
theoretic expectation. In contrast, colored areas below the mean value denote
the “unhealthy” case, when the theoretic expectation is located outside the
oscillation interval. While the conventional Continuum Field (CF) model is
subjected to pinning, the sharp phase-field model allows for arbitrarily small
driving forces.

The condition of phase stability demands the driving force to be small enough
to guarantee meta-stability of the high energy phase: The two local minima of
the potential energy density at $\phi=0,1$ have to be separated by a maximum.
The TI$+h_{3}$-model provides a phase-field width dependent stability limit,
which can be surprisingly high. For instance, imposing the phase-field width
$\tilde{0pt}=0.4$, then the limiting driving force is
$\tilde{\left|\mu\right|}_{\tilde{0pt}=0.4}\lesssim 7200$! The theoretic
stability limits for the different profile resolutions $\tilde{0pt}=0pt/\Delta
x=4.0,\;3.0,\;2.5,\;0.5$ have been indicated by the vertical dashed green
lines in Fig.7. These theoretical limits nicely reflect the behavior of the
sharp phase-field model.

Switching the interpolation function changes the phase stability limits. The
most common choice for the interpolation function is
$h_{\mathrm{5}}=\phi^{3}(10-15\phi+6\phi^{2})$. The CF$+h_{5}$-model provides
phase stability for infinitely large driving forces! However, using
interpolation functions other than the natural one leads to altered
nonequilibrium phase-field profiles. The resulting deviation of fitted phase-
field width $0pt_{\mathrm{fit}}$ from the theoretic expectation $0pt$ is
plotted in the middle row of Fig. 7. The profile alternation increases with
increasing driving force. Increasingly stronger alternations lead to
increasingly stronger grid friction effects. Consider the phase-field width
$\tilde{0pt}=3.0$ and $\tilde{\mu}=100$, then the diffuse interface is
compressed down to $22\%$ of its original width. Grid friction drops the
interface velocity down to about $5\%$ of the theoretic expectation. Thus, for
large dimensionless driving forces the CF-model $h_{5}$ is effectively limited
by spurious grid friction. In the lower part of Fig. 7, we plot the parameter
window of reasonable model operation. We define the range of reasonable
operation to end when the relative velocity error exceeds $0.1$.

##  .3 Translational Invariance of the ideal profile

Testing the Translational Invariance (TI), we calculate the system integral
over the equilibrium condition (3). We consider a discrete 3D system with a
phase-field as represented by an array of 64bit floating point numbers, each
associated with a grid point within the simple cubic numerical lattice of size
$300\times 1\times 1$ (excluding the one stencil boundary halo). The phase-
field values are initialized according to the ideal profile function (2), such
that the interface is sitting in the middle of the system. Then the total grid
friction forces are defined as the system integral over (3). This integral
value may oscillate, when the ideal profile is moved in such a way that the
interface center $c_{n}$ passes several grid points. In Fig. 2, we plot the
oscillation amplitude $A$ of these forces for different interface
orientations. Large oscillation amplitudes indicate broken Translational
Invariance (TI). The conventional model without restoration of TI is shown by
the black curve. Restoring TI using constant grid coupling parameters
$a_{k}\\!\left(\mathbf{u}\right)$ based on a globally constant unit vector
$\mathbf{u}$, as proposed by Finel et al., provides vanishing force
oscillations for those interface orientations, that match with one of the
equivalent numerical lattice directions $\left\langle\mathbf{u}\right\rangle$:
When the grid coupling parameters are, for instance, chosen based on the unit
lattice vector $\mathbf{u}$ parallel to the $\left[110\right]-$direction
(TI⟨110⟩, see dark blue curve in Fig. 2), then vanishing force amplitudes are
found for interface orientations with normal vectors pointing in all the
$\left\langle 110\right\rangle-$directions. The vanishing force amplitudes are
restricted to very sharp interface orientation windows, as visible in Fig. 2.
The new TI⟨n⟩-model (green curve) uses grid coupling parameters, that are
determined by means of the local interface normal direction. This leads to
very small oscillation amplitudes, regardless of the interface orientation.

Translational Invariance (TI) of the ideal planar front solution can also be
tested with regard to oscillations in the total interface energy. Then the
system integral over the interface energy density (1) has to be evaluated
instead. However, in contrast to the forces, theoretically, only the total
interface energy, i.e. the density integral over the full, infinite profile,
provides a TI value. The interface energy density alone does not need to show
this property. For a single direction $k$, the new sharp phase-field
formulation provides TI total interface energies for arbitrarily oriented
ideal phase-field profiles, as long as the full profile function is evaluated
in that direction.

##  .4 Construction of the models

Here, we explain how the different models are constructed from the given
finite difference equations. A overview over all the different models is given
in Tab. 1. The models differ by different choices for the equilibrium
potentials $g_{k}(\phi)$ and for the interpolation function $h(\phi)$.
Further, the source term regularization factor $R(\phi)$ can be either imposed
or otherwise set to unity. All models are separately calibrated. Thus, the
imposed calibration parameters, $C_{\Gamma}$, $\gamma_{j}$, can be different
for the different models. The Continuum Field (CF) model is obtained in the
limit $\lim_{\left|\mathbf{u}_{k}\right|\rightarrow 0}$. In this limit the
equilibrium potentials (4) converge to the classical quartic double-well
potential. For the CF-model, we impose the equilibrium potentials
$g_{k}^{\infty}\\!=\\!\bar{\nu}8\phi^{2}(1\\!-\\!\phi)^{2}$, where the
multiplicity correction $\bar{\nu}\\!=\\!1/3$ equilibrates for the
overweighting by the sum in the equilibrium potentials within each neighboring
shell $j$.

Table 1: The construction of the different models. model |  | equilibrium potential |  | interpolation function |  | regularization |  | calibration
---|---|---|---|---|---|---|---|---
CF$+\\!h_{5}$ |  | $g_{k}^{\infty}\\!=\\!\bar{\nu}8\phi^{2}(1\\!-\\!\phi)^{2}$ |  | $h_{\mathrm{5}}\\!=\\!\phi^{3}(10\\!-\\!15\phi\\!+\\!6\phi^{2})$ |  | – |  | $C_{\Gamma}^{\mathrm{CF}}$, $\gamma_{j}^{\mathrm{CF}}$
TI${}_{\left\langle 100\right\rangle}\\!+\\!h_{3}$ |  |  $g_{k}$: Eq. (4), $a_{k}(\mathbf{u}_{\left\langle 100\right\rangle})$ |  | $h_{3}\\!=\\!\phi^{2}(3\\!-\\!2\phi)$ |  | – |  |  $C_{\Gamma}$, $\gamma_{j}^{\mathrm{TI}_{\left\langle 100\right\rangle}}$
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}$ |  |  $g_{k}$: Eq. (4), $a_{k}\left(\mathbf{n}\right)$ |  | $h_{3}\\!=\\!\phi^{2}(3\\!-\\!2\phi)$ |  | – |  |  $C_{\Gamma}$, $\gamma_{j}^{\mathrm{TI}_{\left\langle\mathbf{n}\right\rangle}}$
TI${}_{\left\langle\mathbf{n}\right\rangle}\\!+\\!h_{3}\\!+\\!R$ |  |  $g_{k}$: Eq. (4), $a_{k}\left(\mathbf{n}\right)$ |  | $h_{3}\\!=\\!\phi^{2}(3\\!-\\!2\phi)$ |  |  $R:$ Eq. (6) |  |  $C_{\Gamma}$, $\gamma_{j}^{\mathrm{TI}_{\left\langle\mathbf{n}\right\rangle}}$, $C_{R}$

Translational Invariance (TI) is obtained when the new equilibrium potentials
Eqs. (4) are imposed in conjunction with the natural interpolation function
$h_{3}$. When all the grid coupling parameters $a_{k}$ in the equilibrium
potentials are set as fixed, based on the globally fixed lattice vector
$\mathbf{u}\\!=\\!\left[100\right]$, then TI is restored for all equivalent
$\left\langle 100\right\rangle-$directions of the computational grid. This
model is denoted as TI${}_{\left\langle 100\right\rangle}+h_{3}$. A
combination of the new equilibrium potentials with the other interpolation
function is not useful, because the nonequilibrium phase-field profile
alternation destroys the carefully restored TI again. In case of the
TI⟨n⟩-models, the locally calculated and length corrected grid coupling
parameters $a_{k}\left(\mathbf{n}\right)$ (see Eq. (5) ff.) are used in the
equilibrium potentials $g_{k}(\phi)$ Eq. (4).

##  .5 Additional information on the modeling

The phase-field equation of motion is given as

$\displaystyle\partial_{t}\phi_{\mathbf{p}}=$
$\displaystyle\frac{2M}{3C_{\Gamma}}\sum_{j,k}\gamma_{j}\nu_{j}\left(\partial_{k}^{-}\left(\partial_{k}^{+}\phi_{\mathbf{p}}\right)-\frac{1}{0pt^{2}}\partial_{\phi}g_{k}(\phi_{\mathbf{p}})\right)$
$\displaystyle-\frac{2M}{30pt\Gamma}\mu\;\partial_{\phi}h(\phi_{\mathbf{p}}),$
(8)

where $M$ is a kinetic coefficient comparable to a diffusion coefficient with
dimension $\left[M\right]=\mathrm{m}^{2}\mathrm{s}^{-1}$. We know that the
phase-field equation promotes solution of the form of Eq. (2). The hyperbolic
tangent function provides the following addition property,

$\displaystyle\tanh\left(p\pm q\right)$
$\displaystyle=\frac{\tanh\left(p\right)\pm\tanh\left(q\right)}{1\pm\tanh\left(p\right)\tanh\left(q\right)}.$
(9)

This property can be reformulated in terms of the phase-field profile
function, and we obtain the relation

$\displaystyle\phi_{\mathbf{p}\pm\mathbf{r}_{k}}$
$\displaystyle=\frac{\left(1\pm
a_{k}\right)\phi_{\mathbf{p}}}{1\pm\left(2\phi_{\mathbf{p}}-1\right)a_{k}},$
(10)

where the grid coupling parameter $a_{k}$ has been introduced as
$a_{k}=\tanh\left(2\mathbf{r}_{k}\cdot\mathbf{n}/0pt\right).$The equilibrium
condition (3) holds, if all 1D $k-$components are simultaneously satisfied.
The individual $k-$component can be satisfied at any real time during the
propagation of the interface by using the addition property of the hyperbolic
tangent profile. Inserting (10) into the $k-$th component of the equilibrium
condition (3) yields

$\displaystyle\partial_{\phi}g_{k}(\phi)$
$\displaystyle=a_{k}^{2}\frac{0pt^{2}}{\mathbf{r}_{k}^{2}}\frac{4\phi\left(1-\phi\right)\left(1-2\phi\right)}{1-a_{k}^{2}\left(1-2\phi\right)^{2}}.$
(11)