# Frictionless motion of marginally resolved diffuse interfaces in phase-field
modeling

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  Uwe Glatzel  Metals and
Alloys, University of Bayreuth, Prof.-Rüdiger-Bormann-Straße 1, 95447
Bayreuth, Bavaria, Germany

###### Abstract

We investigate the influence of artificial grid friction in different phase-
field models by considering the stationary motion of an interface between two
phases at different bulk free energy levels. Following the striking idea of a
Sharp Phase-Field Model (SPFM) from Finel et al.
(FinelLeBouarDabasAppolairYamada2018, ), we proof that restoring translational
invariance indeed eliminates artificial grid friction effects during
stationary interface propagation. Over a largely extended range of possible
driving forces the theoretic interface velocities are reproduced by orders of
magnitude more accurately, even if the diffuse interface profile is only
marginally resolved by just one grid point. We propose a new variant of the
SPFM, which restores translational invariance locally in the direction of
interface motion. It is shown that, even for marginally resolved interface-
profiles, the new SPFM variant provides frictionless motion for arbitrarily
oriented planar interfaces. Finally, by considering thermal diffusion limited
solidification, we demonstrate the capability of the SPFM-approach to also
deal with inhomogeneous driving forces using a one-gird point interface
resolution.

Diffuse interface models provide an elegant description of microstructure
evolution involving phase or domain boundary motion. As compared to sharp
interface descriptions, the difficult problem of explicit interface tracking
is avoided, which allows for any topological evolution of the phase or domain
structures, such as interface instabilities, shape bifurcations, coagulation
events, particle nucleation or dissolution. Phase-field methods are
extensively used in the simulation of complex evolution problems, such as
solidification (Boettinger2002, ), solid-state transformations
(AstaBeckermKarma012009, ), crack propagation (SpatschBrenerKarma012010, ;
PonsKarma2010, ; LubomirskyChenKarmaBouchbinder2018, ; ChenBouchbiKarma082017,
), dislocation dynamics (WangLi012010, ; RuffiniLeBouarFinel2017, ;
MianroodiShanthrajKontisCormierGaultSvendsenRaabe2019, ), ferro-electric
domain evolution (YadavNguyenChenSalahuddin2019, ), grain growth
(MoelansBlanpaiWollant072008, ;
DarvishiKamachaliAbbondandoloSiburgSteinbach2015, ), as well as many other
(TonksAagesen2019, ).

Quite often, the width of the diffuse interface appears to be the smallest
length-scale in the system. Obviously, in all these cases one is interested in
choosing the smallest possible width, while still keeping the benefits from
the diffuse interface description. So far, artificial gird friction or grid
pinning in the phase-field equation, which naturally results from the
numerical discretization, has been the major limiting factor in this regard.

Recently Finel et al. found a striking new way to eliminate these effects in
the 1D phase-field equation. The key is to restore the translational
invariance (TI), that has been broken by the numerical discretization. Then a
one-grid point interface resolution is possible without grid pinning
(FinelLeBouarDabasAppolairYamada2018, ). The method is related to previous
suggestions to improve the numerical performance of phase-field solvers via a
nonlinear scaling of the phase-field equation (Glasner2001, ; Weiser2009, ).
Similar formulations involving the section wise defined sinus like phase-field
profile have been independently proposed by J. Eiken (Eiken2012, ).

We consider the stationary motion of a planar interface between two phases at
different bulk free energy levels. The bulk free energy difference drives the
growth of the low energy phase on the expanse of the high energy phase. The
coupling of the phase-field to the bulk energy density difference is
prototypical for many advanced phase-field models, where the driving forces
are calculated from local concentrations or strains (Kassner2001, ;
EchebarFolchKarma122004, ; Steinba2009, ; SteinbaShchygl022011, ;
PilipenkoFleckEmmerich2011, ; GaubertJouiadCormierLeBouarGhighi2015, ;
MushongFleckQuerfur032015, ; JokisaariNaghaviWolvertonVoorheesHeinonen2017, ;
FleckSchleiferHolzingerGlatzel2018, ; AagesenGaoSchwenAhmed2018, ).

Figure 1: The simulation configuration: Stationary motion of a planar
interface.

Keeping the simple cubic computational grid with grid spacing $h$ as fixed, we
investigate the influence of varying interface orientations on the stationary
motion, see Fig. 1. Deviations of the resulting interface velocity from
theoretic expectations reveal the effect of artificial grid friction and
pinning.

### Sharp Phase-Field Modeling

The derivation of the phase-field models is stated from a Helmholtz free
energy functional

$\displaystyle F\left[\phi\right]$
$\displaystyle=\int_{V}\underbrace{\left(\frac{0pt\Gamma\left(\nabla\phi\right)^{2}}{2C_{\Gamma}}+\frac{\Gamma\;g(\phi,\tilde{0pt})}{0ptC_{\Gamma}}+\mu\;p\left(\phi\right)\right)}_{f(\phi,\nabla\phi)}\mathrm{d}V,$
(1)

where $V$ is the volume of the system, and $\nabla$ denotes the gradient. The
equilibrium potential $g(\phi,\tilde{0pt})$, depends on the phase-field $\phi$
and and on $\tilde{0pt}=0pt/h$. It has to vanish at the two local minima at
$\phi=0$ and $\phi=1$, which correspond to the distinct phases of the system.
$0pt$ denotes the width of the diffuse interface, $\Gamma$ is the interface
energy density, and $C_{\Gamma}$ is a dimensionless calibration factor. The
calibration factor $C_{\Gamma}$ is calculated via $C_{\Gamma}=\int
f(\phi_{0},\nabla\phi_{0})_{\mu=0}\mathrm{d}\mathbf{n}/\Gamma$, where
$\mathbf{n}$ denotes the interface normal direction, and the phase-field
$\phi_{0}$ is given by Eq. (3). For further information, we refer to the
supplementary material. A positive bulk free energy density difference $\mu$
favors the growth of phase $\phi=0$ on the expanse of phase $\phi=1$. The
interpolation function $p(\phi)$ satisfies $p(0)=0$ and $p(1)=1$ and further
provides vanishing slope at $\phi=0$ and $\phi=1$.

The time derivative of the phase-field, $\partial_{t}\phi$, is proportional to
the variational derivative of the free energy functional (1),
$-\delta_{\phi}F$. We write the equation of motion as
$30pt\Gamma\cdot\partial_{t}\phi=-2M\delta_{\phi}F$
(FleckMushongPilipen102011, ), and obtain

$\displaystyle\partial_{t}\phi=$
$\displaystyle\frac{2M}{3C_{\Gamma}}\left(\nabla^{2}\phi-\frac{1}{0pt^{2}}\partial_{\phi}g(\phi,\tilde{0pt})\right)-\frac{2M}{30pt\Gamma}\mu\;\partial_{\phi}p(\phi),$
(2)

where $\nabla^{2}$ denotes the Laplace operator, $M$ is a kinetic coefficient
comparable to a diffusion coefficient with dimension
$\left[M\right]=\mathrm{m}^{2}\mathrm{s}^{-1}$, and
$\partial_{\phi}=\partial/\partial\phi$ abbreviates the partial derivative
with respect to the phase-field $\phi$.

During stationary interface motion a constant amount of energy per unit time
interval dissipates via the proceeding phase transformation. Thus, the phase
transformation rate and interface velocity is exactly determined by the diving
force $\mu$ via $v_{\mathrm{th}}=-M\mu/\Gamma$. Choosing
$g(\phi,\infty)=8\phi^{2}\left(1-\phi\right)^{2}$ and
$p(\phi)=p_{3}(\phi)=\phi^{2}(3-2\phi)$, the following special analytic
solution for Eq. (2) can be found

$\displaystyle\phi_{\mathrm{0}}(\mathbf{r},t)$
$\displaystyle=\left(1-\tanh\left(2\left(\mathbf{r}\cdot\mathbf{n}-r_{n}\right)/0pt\right)\right)/2.$
(3)

Here $\mathbf{n}$ is the unit normal interface vector, $\mathbf{r}$ is the
position vector and $r_{n}=v_{\mathrm{th}}t$ denotes the central interface
position, moving with the velocity $v_{\mathrm{th}}$.

For vanishing driving forces $\mu=0$ and vanishing phase-field motion
$\partial_{t}\phi=0$, Eq. (2) reduces to the equilibrium condition, which
reads in the discrete form as

$\displaystyle{\textstyle\sum_{k=1}^{3}}0pt^{2}/h^{2}\left(\phi_{\mathbf{i}+\mathbf{u}_{k}}-2\phi_{\mathbf{i}}+\phi_{\mathbf{i}-\mathbf{u}_{k}}\right)$
$\displaystyle-\partial_{\phi}g(\phi,\tilde{0pt})=0,$ (4)

where $\mathbf{i}$ denotes the grid points of the simple cubic 3D numerical
lattice and the lattice vector $\mathbf{u}_{k}$ connects two neighboring grid
points along the direction $\mathbf{e}_{k}$ with a spacing
$h=\left|\mathbf{u}_{k}\right|$. Similar to
(FinelLeBouarDabasAppolairYamada2018, ), the equilibrium potential is
decomposed into $g(\phi,\tilde{0pt})=\sum_{k=1}^{3}g_{k}(\phi,\tilde{0pt})$.
The 3D equilibrium condition (4) is satisfied, if all three 1D cartesian
contributions are independently satisfied at the same time. Interestingly, the
$k-$th component of the condition (4) can be satisfied at any real time during
the propagation of the interface by using the following addition property of
the hyperbolic tangent profile (3) (FinelLeBouarDabasAppolairYamada2018, ):

$\displaystyle\phi_{\mathbf{i}\pm\mathbf{u}_{k}}$
$\displaystyle=\frac{\left(1\pm
a_{k}\right)\phi_{\mathbf{i}}}{1\pm\left(2\phi_{\mathbf{i}}-1\right)a_{k}},$
(5)

where the grid coupling parameters $a_{k}$ are defined as

$\displaystyle a_{k}$
$\displaystyle=\tanh\left(2\mathbf{u}_{k}\cdot\mathbf{n}/0pt\right)=\tanh\left(2hn_{k}/0pt\right).$
(6)

$n_{k}$ denotes the projection of $\mathbf{n}$ onto the direction
$\mathbf{e}_{k}$. Using Eq. (5) and (4) we obtain the $k-$th component of the
derivative of the equilibrium potential

$\displaystyle\partial_{\phi}g_{k}(\phi,\tilde{0pt})$
$\displaystyle=a_{k}^{2}\tilde{0pt}^{2}\frac{4\phi\left(1-\phi\right)\left(1-2\phi\right)}{1-a_{k}^{2}\left(1-2\phi\right)^{2}}.$
(7)

The actual potential is given by the integral,
$g(\phi,\tilde{0pt})=\int_{0}^{\phi}\sum_{k}\partial_{\phi}g_{k}(\phi^{\prime},\tilde{0pt})\mathrm{d}\phi^{\prime}$
(also see (FinelLeBouarDabasAppolairYamada2018, )). In the continuum limit
$h\rightarrow 0$, this potential converges the original quartic double well
potential, since $\lim_{h\rightarrow 0}\tilde{0pt}^{2}a_{k}^{2}=4$.

Figure 2: Logarithmic plot of the oscillation amplitude $A$ of the system
integral over the 3D equilibrium condition (4), using the ideal phase-field
profile function (3), for different interface orientation angles
$\vartheta_{[001]}$ and $\vartheta_{[011]}$ (see Fig. 5). The width is
$\tilde{0pt}=1$ and the system size is $300\times 1\times 1$ plus a one
stencil boundary.

For an accurate determination of the local grid coupling parameters $a_{k}$,
we propose the following scheme based on the addition property Eq. (5),
$a_{k}=\left(a_{k}^{+}+a_{k}^{-}\right)/2$, where

$\displaystyle a_{k}^{\pm}$
$\displaystyle=\frac{\pm\left(\phi_{\mathbf{i}\pm\mathbf{u}_{k}}-\phi_{\mathbf{i}}\right)}{\phi_{\mathbf{i}}-2\phi_{\mathbf{i}\pm\mathbf{u}_{k}}\phi_{\mathbf{i}}+\phi_{\mathbf{i}\pm\mathbf{u}_{k}}}.$
(8)

By inserting the derivative (7) into the equilibrium condition (4) the
explicit dependence on the interface width $0pt$ cancels out. Then, the phase-
field width is solely controlled by the grid coupling parameters (8). In
contrast to the original formulation by Finel et al.
(FinelLeBouarDabasAppolairYamada2018, ), here, these parameters depend on the
interface normal vector, which is calculated locally by (8). Consequently,
without additional length control of $\mathbf{n}$ the width of the diffuse
interface wouldn’t be defined in the model. Therefore, we locally calculate
the components $n_{k}=\tilde{0pt}\;\mathrm{arctanh}\;a_{k}/2$, restore unit
length via $\mathbf{n}^{\mathrm{new}}=\mathbf{n}/\left|\mathbf{n}\right|$ and
re-calculate corrected grid coupling parameters
$a_{k}^{\mathrm{new}}\left(\mathbf{n}^{\mathrm{new}}\right)$ using (6).

Testing the translational invariance, we calculate the system integral over
the equilibrium condition (4) using the discrete 3D phase-field as given by
the ideal profile function (3). When the ideal profile is moved across several
grid points, the discrete force integral oscillates. In Fig. 2, we plot the
force oscillation amplitude $A$ for different interface orientations. Large
oscillation amplitudes indicate broken translational invariance (TI), as for
the classical no TI formulation shown by the black curve. Fixed grid coupling
parameters, as proposed by Finel et al. (FinelLeBouarDabasAppolairYamada2018,
), restore translational invariance for fixed interface orientations (fixed
TI). The new TI model (green curve) restores translational invariance locally
in the interface normal direction, leading to evenly small force oscillation
amplitudes $A$, regardless of the interface orientation. On the other hand,
interface orientation dependent grid coupling parameters lead to a
substantially higher algorithmic complexity, than the original model proposed
by A. Finel.

### Possible range of driving forces

Grid friction and pinning during stationary motion has been studied earlier
(KarmaRappel041998, ; BoeschHMKShochet1995, ). It leads to a reduction of the
average interface width and velocity. Both, the velocity and the width
oscillate as the center of the interface passes one grid point after the
other. With decreasing phase-field width, we obtain increasingly larger drops
of the average values as well as increasing oscillation amplitudes. This
culminates in a vanishing velocity, where the phase-field is pinned the to
computational grid.

Figure 3: Logarithmic plot of the relative error in the stationary interface
velocity as a function of the dimensionless driving force $\tilde{\mu}=\mu
h/\Gamma$. The transparently colored areas indicate the oscillation range of
the velocity error.

In Fig. 3, we plot the relative deviation of the measured interface velocity
$v$ from its theoretic expectation $v_{\mathrm{th}}$ as function of the
driving force $\tilde{\mu}=\mu h/\Gamma$. The mean error is depicted by the
thick solid/dashed lines. The continuous lines indicate a constant time
resolution of $M\mu\Delta t/(\Gamma h)=8\cdot 10^{-8}$, where $\Delta t$ is
the time step. For all these driving forces the interface center has passes at
least four grid points. The dashed lines indicate shorter simulations, where
the interface center has not been advanced by the distance of a full grid
point. The difference between the oscillation amplitude and the mean value is
plotted as transparently colored area. When the colored area is found above
the mean value, we have the “healthy” situation that the measured interface
velocity oscillates around the theoretic expectation. In contrast, colored
areas below the mean value denote the “unhealthy” case, when the theoretic
expectation is located somewhere outside the oscillation interval.

Without restoring of translational invariance (no TI), grid friction leads to
substantial errors in the interface velocity. Even for interfaces as wide as
$\tilde{0pt}=5$, grid pinning occurs for driving forces below $\mu
h/4\Gamma<10^{-7}$. In contrast, the 1D SPFM is not subjected to this
limitation, even if the phase-field width is chosen as small as
$\tilde{0pt}=0.6$. For the SPFM the error in the interface velocity only
depends on the error in the time discretization!

Conceptually allowed values for driving forces have to be small enough to
guarantee the meta-stability of the high energy phase. In other words, the two
local minima of the potential energy at $\phi=0$ and $\phi=1$ have to be
separated by a maximum. Considering the potentials, $g(\phi,\infty)$ and
$p_{3}(\phi)$, this requires the absolute value of the driving force to be
below $\tilde{\left|\mu\right|}<4/0pt$, with $\tilde{\mu}=\mu h/\Gamma$
denoting the dimensionless driving force. For a fixed interface energy the
range of possible driving forces is smaller the wider the interface is! This
circumstance can be quite restricting as the appearance of interface
instabilities, such as the diffusional Mullins-Sekerka (MullinsSekerka1964, )
or the elastic Asaro-Tiller-Grinfeld instability (AsaroTiller1972, ;
Grinfeld1993, ), require the interface energy to be comparably small
(FleckQuerfurthGlatzel2017, ; SpatschekFleck2007, ). Using the SPFM potentials
$p_{3}(\phi)$ and $g(\phi,\tilde{0pt}=0.6)$ together with the energy
calibration $C_{\Gamma}=0.238$, we calculate the range of possible driving
forces for an interface as sharp as $\tilde{0pt}=0.6$ to be
$\tilde{\left|\mu\right|}\lesssim 230$! The upper bounds for the driving force
for the different cases are indicated by respectively colored dash-dotted
vertical lines in Fig. 3.

Figure 4: Comparison of the fitted interface width $0pt_{\mathrm{fit}}$ as a
function the dimensionless driving force $\tilde{\mu}=\mu h/\Gamma$. The 1D
SPFM with $\tilde{0pt}=2$ is used.

Using interpolation functions of higher polynomial order can also extend the
theoretic range of possible driving forces. However, the effect of using
interpolation functions other than the natural function,
$\partial_{\phi}p_{\mathrm{3}}=6\phi(1-\phi)$, for which a steady state
solution of Eq. (2) exists, is that the nonequilibrium phase-field profile is
altered, see Fig. 4. The steeper interpolation function,
$\partial_{\phi}p_{\mathrm{5}}=30\phi^{2}\left(1-\phi\right)^{2}$
(KarmaRappel041998, ; EchebarFolchKarma122004, ), provide decreasing phase-
field widths with increasing driving forces. Lower order interpolation
functions, such as
$\partial_{\phi}p_{\mathrm{\mathrm{sin}}}=8\sqrt{\phi(1-\phi)}/\pi$
(Steinba2009, ), lead to an increase in phase-field width as the driving force
increases. The 1D SPFM in conjunction with $p_{5}$ provides significant
oscillations in the measured interface width for $\tilde{\mu}\geq 1$, as
indicated by the dark-orange colored area in Fig. 4. A further extension of
the range of possible driving forces within the 1D SPFM by using $p_{5}(\phi)$
requires to account for the altered nonequilibrium profile.

### Frictionless interfaces motion in 3D

In 3D the direction of interface motion appears as an important additional
degree of freedom. We consider interface orientations $\mathbf{n}$ resulting
from rotations around the two different axes $\left[001\right]$ and
$\left[011\right]$, as sketched in Fig. 5. The simulation of interface
propagation in directions other than the $\left\langle
100\right\rangle-$directions of the computational grid requires special
wetting-angle boundary conditions for the phase-field (GuoSteinba042010, ;
BenSaidSelzerNestlerBraunGreinerGarcke2014, ;
DiewaldLautenschlaegerStephanLangenbachKuhnSecklerBungartzHasseMueller2020, ).
These have to account for the contact angle $\alpha$ between the boundary
direction $\mathbf{e}{}_{k}$ and the interface normal $\mathbf{n}$. The
boundary value at $\mathbf{b}$ is calculated as

$\displaystyle\phi_{\mathbf{b}}$
$\displaystyle=\frac{\left(1-a_{n}\right)\phi_{\mathbf{b}+\mathbf{u}_{k}}}{1-a_{n}\left(2\phi_{\mathbf{b}+\mathbf{u}_{k}}-1\right)},$
(9)

with $a_{n}=\tanh(2h\sin\alpha/0pt)$. The idea is to calculate the boundary
value at $\mathbf{b}$ from the neighboring phase-field value at
$\mathbf{b}+\mathbf{u}_{k}$ using the addition theorem for the hyperbolic
tangent profile (5), and imposing a profile-shift by the length
$s_{n}=h_{k}\sin\alpha$ along the interface normal direction.

Figure 5: Plot of the relative error in the $x-$component of the interface
velocity for different orientations of the interface. The driving force
decreases with increasing angles, i.e. $\mu=\mu_{0}\cos\left(\vartheta\right)$
with $h\mu{}_{0}/\Gamma=1.0$. The time discretization is $M\mu\Delta t/(\Gamma
h)=5\cdot 10^{-5}$.

In Fig. 5, we compare the orientation dependent error in the interface
velocity for two different models. The classical model with no restored
translational invariance (no TI) and $\tilde{0pt}=5$ is shown in blue color.
The green curve in Fig. 5 shows the behavior of the new model, where
translational invariance has been restored in the direction of interface
motion. It provides very accurate interface velocities for all orientations,
even if the phase-field width is as small as $\tilde{0pt}=0.6$. As for the 1D
SPFM, the resulting velocity error only depends on the error in the time
discretization. Thus, the new model indeed provides frictionless motion of
planar interfaces with arbitrary orientations! This study is partly comparable
to the sphericity study of a growing sphere as presented in Fig. 2 in
(FinelLeBouarDabasAppolairYamada2018, ) showing the respective behavior of the
original 3D SPFM.

### Diffusion limited solidification

The feasibility to deal with time dependent and inhomogeneous driving forces
$\mu$ is demonstrated by considering the thermal diffusion limited
solidification of a pure substances within a thermally isolated system. We
introduce the dimensionless temperature field
$U(\mathbf{x},t)=C\left(T(\mathbf{x},t)-T_{M}\right)/L$, where $T_{M}$, $L$
and $C$ denote melting temperature, latent heat and heat capacity,
respectively. $\mu$ is proportional to the temperature field by
$\mu=U\Gamma/d_{0}$, where $d_{0}=\Gamma T_{M}C/L^{2}$ denotes the capillary
length, which scales with the interface energy density $\Gamma$. The
temperature obeys a diffusion equation,
$\partial_{t}U=D\nabla^{2}U+\partial_{\phi}p\partial_{t}\phi$, with the
thermal diffusion coefficient $D$ (Nestler2005221, ; FleckHuterPilipen012010,
; JokisaariVoorheesGuyerWarrenHeinonen2018, ). The new model is capable to
simulate diffusion limited solidification with interfaces as sharp as
$\tilde{0pt}=0.6$! For an initial undercooling temperature of
$U(\mathbf{x},0)=-0.3$ the smallest possible capillary length is
$d_{0}/h=0.003$. This value is about a factor 30 smaller than the smallest
value used in (KarmaRappel041998, ). It is small enough to resolve a full
dendrite with side-arms within a $60\times 30\times 30-$system (see Fig. 6).
Such a small capillary length to grid spacing ratio relates to a maximal
dimensionless driving force of
$\tilde{\mu}_{\mathrm{max}}=U(\mathbf{x},0)/d_{0}=100$. Therefore, the SPFM
related extension of the range of possible driving forces is accessible in
practical 3D simulations of thermal diffusion limited solidification!

Figure 6:  3D Simulation of thermal diffusion limited solidification using the
new TI model with $\tilde{0pt}=0.6$, within a thermally insulating $60\times
30\times 30-$system. The temperature $U$ is indicated by the coloring and the
phase-field is visualized via $\phi=1/2-$contour. Further parameters are
$d_{0}/h=0.003$, $U(\mathbf{x},0)=0.3$ and $D/M=0.01$.

The interface energy anisotropy, necessary for the dendritic selection
(BrenerBoussinHuter2009, ) visible in Fig. 6, results from the underlying
cubic grid. Equilibrium particle shape calculations with conserved phase
volumes (KarmaRappel041998, ; NestlerWendlerSelzer072008, ;
BhadakSankarasubramanianChoudhury2018, ; HolzingerSchleiferGlatzelFleck2019, )
show that for $\tilde{0pt}=0.6$ the new TI model provides a lattice anisotropy
corresponding to
$\epsilon_{4}=\left(R_{\mathrm{max}}-R_{\mathrm{min}}\right)/2R_{\mathrm{max}}=0.064$,
or $R_{\mathrm{min}}/R_{\mathrm{max}}=0.881$, where, $R_{\mathrm{min}}$ and
$R_{\mathrm{max}}$ denote the smallest and largest distance between the
nonlinear $\phi=1/2-$contour and the particle’s center, respectively.

### Conclusion

The effect of artificial grid friction on the stationary motion of planar
interfaces is investigated in a classical phase-field model as well as
variants of the Sharp Phase-Field Model (SPFM)
(FinelLeBouarDabasAppolairYamada2018, )

  * •

With a one-grid point interface resolution the SPFM provides a largely
extended range of possible driving forces.

  * •

The extended parameter range is accessible in 3D simulations of thermal
diffusion limited solidification.

  * •

The newly proposed model with translational invariance in the direction of
interface motion provides frictionless motion of marginally resolved diffuse
interfaces for arbitrary interface orientations.

###### Acknowledgements.

We thank the Federal Ministry for Economics and Energy, Germany (BMWi), for
financially supporting the project COORETEC: ISar 03ET7047D. Further, we thank
the Federal Ministry of Education and Research (BMBF) for financial support in
the project ParaPhase 01IH15005B. The work is also funded by the Deutsche
Forschungsgemeinschaft (DFG) 387117768 (assigned to the priority program 1713)
and 431968427. We are very thankful to 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.

## References

  * (1) A. Finel, Y. Le Bouar, B. Dabas, B. Appolaire, Y. Yamada, and T. Mohri.  _Phys. Rev. Lett._ 121, 025501 (2018), [link].
  * (2) W. Boettinger, J. Warren, C. Beckermann, and A. Karma.  _Annu. Rev. Mater. Res._ 32, 163 (2002), [link].
  * (3) M. Asta, C. Beckermann, A. Karma, W. Kurz, R. Napolitano, M. Plapp, G. Purdy, M. Rappaz, and R. Trivedi.  _Acta Mater._ 57, 941 (2009), [link].
  * (4) R. Spatschek, E. A. Brener, and A. Karma.  _Phil. Mag._ 91, 75–95 (2011), [link].
  * (5) A. J. Pons and A. Karma.  _Nature_ 464, 85 (2010), [link].
  * (6) Y. Lubomirsky, C.-H. Chen, A. Karma, and E. Bouchbinder.  _Phys. Rev. Lett._ 121, 134301 (2018), [link].
  * (7) C.-H. Chen, E. Bouchbinder, and A. Karma.  _Nature Physics_ 13, 1186–1190 (2017), [link].
  * (8) Y. Wang and J. Li.  _Acta Mater._ 58, 1212 (2010), [link].
  * (9) A. Ruffini, Y. Le Bouar, and A. Finel.  _J. Mech. Phys. Sol._ 105, 95–115 (2017), [link].
  * (10) J.R. Mianroodi, P. Shanthraj, P. Kontis, J. Cormier, G. Gault, B. Svendsen, and D. Raabe.  _Acta Mater._ (2019), [link].
  * (11) A.K. Yadav, K.X. Nguyen, Z. Hong, P. García-Fernández, P. Aguado-Puente, C.T. Nelson, S. Das, B. Prasad, D. Kwon, S. Cheema, A.I. Khan, C. Hu, J. Íñiguez, J. Junquera, L.-Q. Chen, D. A. Muller, R. Ramesh, and S. Salahuddin.  _Nature_ 565, 468–471 (2019), [link].
  * (12) N. Moelans, B. Blanpain, and P. Wollants.  _Phys. Rev. Lett._ 101, 025502 (2008), [link].
  * (13) R. Darvishi Kamachali, A. Abbondandolo, K.F. Siburg, and I. Steinbach.  _Acta Mater._ 90, 252–258 (2015), [link].
  * (14) M.R. Tonks and L.K. Aagesen.  _Annu. Rev. Mater. Res._ 49, 79–102 (2019), [link].
  * (15) K. Glasner.  _J. Comp. Phys._ 174, 695–711 (2001), [link].
  * (16) M. Weiser.  _Appl. Num. Math._ 59, 1858–1869 (2009), [link].
  * (17) J. Eiken.  _IOP Conf. Ser.: Mater. Sci. Eng._ 33, 012105 (2012), [link].
  * (18) K. Kassner, C. Misbah, J. Müller, J. Kappey, and P. Kohlert.  _Phys. Rev. E_ 63, 036117 (2001), [link].
  * (19) B. Echebarria, R. Folch, A. Karma, and M. Plapp.  _Phys. Rev. E_ 70, 061604 (2004), [link].
  * (20) I. Steinbach.  _Modell. Sim. Mater. Sci. Eng._ 17, 073001 (2009), [link].
  * (21) I. Steinbach and O. Shchyglo.  _Curr. Opin. Solid State Mater. Sci._ 15, 87–92 (2011), [link].
  * (22) D. Pilipenko, M. Fleck, and H. Emmerich.  _Eur. Phys. J. Plus_ 126, 100 (2011), [link].
  * (23) A. Gaubert, M. Jouiad, J. Cormier, Y. Le Bouar, and J. Ghighi.  _Acta Mater._ 84, 237–255 (2015), [link].
  * (24) L. T. Mushongera, M. Fleck, J. Kundin, F. Querfurth, and H. Emmerich.  _Adv. Eng. Mater._ 17, 1149–1157 (2015), [link].
  * (25) A.M. Jokisaari, S.S. Naghavi, C. Wolverton, P. W. Voorhees, and O. Heinonen.  _Acta Mater._ 141, 273–284 (2017), [link].
  * (26) M. Fleck, F. Schleifer, M. Holzinger, and U. Glatzel.  _Metall. Mater. Trans. A_ 49, 4146–4157 (2018), [link].
  * (27) L.K. Aagesen, Y. Gao, D. Schwen, and K. Ahmed.  _Phys. Rev. E_ 98, 023309 (2018), [link].
  * (28) M. Fleck, L. Mushongera, D. Pilipenko, K. Ankit, and H. Emmerich.  _Eur. Phys. J. Plus_ 126, 95 (2011), [link].
  * (29) A. Karma and W.-J. Rappel.  _Phys. Rev. E_ 57, 4323–4349 (1998), [link].
  * (30) A. Bösch, H. Müller-Krumbhaar, and O. Shochet.  _Z. Phys. B_ 97, 367–377 (1995), [link].
  * (31) W.W. Mullins and R.F. Sekerka.  _J. Appl. Phys._ 35, 444–451 (1964), [link].
  * (32) R.J. Asaro and W.A. Tiller.  _Metall. Mater. Trans. A_ 3, 1789–1796 (1972), [link].
  * (33) M.A. Grinfeld.  _J. Nonlin. Sci._ 3, 35–83 (1993), [link].
  * (34) M. Fleck, F. Querfurth, and U. Glatzel.  _J. Mater. Res._ 32, 4605–4615 (2017), [link].
  * (35) R. Spatschek and M. Fleck.  _Phil. Mag. Lett._ 87, 909 (2007), [link].
  * (36) W. Guo and I. Steinbach.  _Int. J. Mater. Res._ 4, 480 (2010), [link].
  * (37) M. Ben Said, M. Selzer, B. Nestler, D. Braun, C. Greiner, and H. Garcke.  _Langmuir_ 30, 4033–4039 (2014), [link].
  * (38) F. Diewald, M.P. Lautenschlaeger, S. Stephan, K. Langenbach, C. Kuhn, S. Seckler, H.-J. Bungartz, H. Hasse, and R. Müller.  _Comp. Meth. Appl. Mech. Eng._ 361, 112773 (2020), [link].
  * (39) B. Nestler, D. Danilov, and P. Galenko.  _J. Comp. Phys._ 207, 221 (2005), [link].
  * (40) M. Fleck, C. Hüter, D. Pilipenko, R. Spatschek, and E. A. Brener.  _Phil. Mag._ 90, 265 (2010), [link].
  * (41) A.M. Jokisaari, P. W. Voorhees, J. E. Guyer, J. A. Warren, and O. Heinonen.  _Comp. Mater. Sci._ 149, 336–347 (2018), [link].
  * (42) E.A. Brener, G. Boussinot, C. Hüter, M. Fleck, D. Pilipenko, R. Spatschek, and D.E. Temkin.  _J. Phys. Cond. Mat._ 21, 464106 (2009), [link].
  * (43) B. Nestler, F. Wendler, M. Selzer, B. Stinner, and H. Garcke.  _Phys. Rev. E_ 78, 011604 (2008), [link].
  * (44) B. Bhadak, R. Sankarasubramanian, and A. Choudhury.  _Metall. Mater. Trans. A_ 49, 5705–5726 (2018), [link].
  * (45) M. Holzinger, F. Schleifer, U. Glatzel, and M. Fleck.  _Euro. Phys. J B_ 92, 208 (2019), [link].

Supplementary material for: “Frictionless motion of marginally diffuse
interfaces in phase-field modeling”

Whenever, the energy is evaluated across the diffuse interface region, the
model should provide the exact amount of energy excess, which corresponds to
the imposed value interface energy density $\Gamma$. Therefore, the
calibration factor $C_{\Gamma}$ for the interface energy can be calculated via
the following line integral $C_{\Gamma}=\int
f(\phi_{0},\partial_{i}\phi_{0})_{\mu=0}\mathrm{d}\mathbf{n}/\Gamma$, where
$\mathbf{n}$ denotes the direction normal to the interface, and $\phi_{0}$
denotes a phase-field profile function, which contains one respectively
oriented interface, where the phase-field variable undergoes the full
transition from $\phi=0$ to $\phi=1$. In Fig. 7, we plot the interface energy
calibration factor $C_{\Gamma}$ as function of the dimensionless interface
grid resolution number $0pt/h$. For the marginally resolved diffuse
interfaces, we obtain values for the calibration line integral mentioned above
that substantially differ from the continuum value.

Figure 7: Plot of the interface energy calibration factor $C_{\Gamma}$ as a
function of the interface grid resolution number $0pt/h$.

##  .1 Grid friction

Figure 8: The influence of artificial grid friction on the motion of a planar
interface. Comparison of a conventional phase-field formulation for different
phase-field widths $0pt/h$ with the one dimensional Sharp Phase-Field Model
(1D SPFM), as proposed by (FinelLeBouarDabasAppolairYamada2018, ), for a
marginally diffuse interface (green curves).

Fig. 8 illustrates the influence from artificial grid friction on the
stationary motion of planar interfaces in one dimension. An animated version
of Fig. 8 is also provided in the supplementary material. We compare the
resulting interface motion for different choices of the ratio between the
phase-field width and the numerical grid spacing $0pt/h$. On top a snapshot of
the phase-field values at different grid points (full symbols) around the
interface are plotted. This is accomplished by a least square fit of the
profile function (Eq. (3) in the manuscript) to the phase-field values, as
indicated by the solid lines. Below, first the energy density and then the
fitting value for the phase-field width are plotted both as function of the
advancing center of the interface $x_{0}=v_{\mathrm{th}}t$.

Due to artificial grid friction, the average values of the measured interface
width and energy are both smaller than the analytic expectations. The squeezed
interface propagates with a clearly smaller average velocity than expected.
Further, we obtain oscillations in the interface energy and the fitted width
as the center of the interface passes one grid point after the other. With
decreasing phase-field width, we obtain increasingly larger drops of the
average values as well as increasing oscillation amplitudes. For the
conventional model this culminates in fully destroyed interface kinetics,
which is commonly referred to as grid pinning. In contrast, the
translationally invariant one dimensional Sharp Phase-Field Model (1D SPFM) by
A. Finel et al. (FinelLeBouarDabasAppolairYamada2018, ), as shown by the
green-curves in Fig. 8, is not subjected to this effects even-though the
phase-field width has been chosen as small as $0pt/h=0.5$.

##  .2 Measure of the interface position

It should be mentioned that for marginally diffuse interfaces the conventional
measure of the interface position by a linear interpolation of the $\phi=1/2$
contour-position is not the best choice to approximate the real center of the
phase-field profile. Especially, when the phase-field width is small compared
to the grid spacing, a linear interpolation can provide a quite different
result as compared to the interface position resulting from a least square fit
using the analytic profile function. On the other hand performing a least
square-fit only to get a measure for the local interface position is most
often not appropriate.

Therefore, we propose a nonlinear interpolation of the interface position
based on the analytic phase-field profile. The interpolated position, where
the phase-field profile takes the arbitrary contour-level $l$ ($0<l<1$) can be
calculated as

$\displaystyle x_{k}^{\mathrm{int}}\left(\mathbf{i}\right)$
$\displaystyle=\mathbf{i}h+\frac{\lambda}{2n_{k}}\left|\mathsf{arctanh}\frac{l-\phi_{\mathbf{i}}}{2l\phi_{\mathbf{i}}-\phi_{\mathbf{i}}-l}\right|,$
(10)

where $n_{k}$ denotes the projection of the unit normal interface orientation
vector onto the $k-$th cartesian direction, which also corresponds to one of
the principle directions of the simple cubic numerical lattice in this case.
In order to further regularize the resulting interface position as function of
time, we impose an averaging between the two slightly different interface
positions $x_{k}^{\mathrm{int}}\left(\mathbf{i}\right)$ and
$x_{k}^{\mathrm{int}}\left(\mathbf{i}+\mathbf{u}_{k}\right)$ of the two
neighboring grid points, connected by the principle lattice vector
$\mathbf{u}_{k}$, next to the position of the contour-level $l$ , i.e. we have
the condition
$\left(\phi_{\mathbf{i}}-l\right)\cdot\left(\phi_{\mathbf{i}+\mathbf{u}_{k}}-l\right)\leq
0$.

We attach three movies to this supplemental material:

  1. 1.

Supplementary_material_2_Pinning_animation.mpg: This movie is an animated
version of Figure 1 in the manuscript. Thus it illustrates the influence of
artificial grid-friction on the motion of a planar interface in one dimension.
We compare the behavior of the conventional phase-field formulation for
different phase-field widths $0pt/h$ with the behavior of the Sharp Phase-
Field Model (1D SPFM), as proposed by (FinelLeBouarDabasAppolairYamada2018, ),
for a marginally diffuse interface (green curves)

  2. 2.

Supplementary_material_3_Steady_state_interface_motion.mpg: This movie is an
animated visualization of the evolution of the phase-field during the
simulation of stationary motion of a planar interface with propagates under an
angle of 45° with respect to the computational grid.

  3. 3.

Supplementary_material_4_ContactAngleBD.mpg: This movie illustrates the
function of the newly proposed boundary conditions. They enforce a finite
contact angle between the interface normal and the boundary plane. In this
movie, we show a simulation of the shape-evolution of an initially cubic
particle toward its spherical equilibrium shape under conserved phase volume.
The particle is in contact with the bottom boundary, where a contact angle of
80° with respect to the boundary plane is enforced.