# A Cartesian grid-based boundary integral method for moving interface
problems

Han Zhou<EMAIL_ADDRESS>Wenjun Ying<EMAIL_ADDRESS>School of
Mathematical Sciences, MOE-LSC and Institute of Natural Sciences, Shanghai
Jiao Tong University, Minhang, Shanghai 200240, PR China

###### Abstract

This paper proposes a Cartesian grid-based boundary integral method for
efficiently and stably solving two representative moving interface problems,
the Hele-Shaw flow and the Stefan problem. Elliptic and parabolic partial
differential equations (PDEs) are reformulated into boundary integral
equations and are then solved with the matrix-free generalized minimal
residual (GMRES) method. The evaluation of boundary integrals is performed by
solving equivalent and simple interface problems with finite difference
methods, allowing the use of fast PDE solvers, such as fast Fourier transform
(FFT) and geometric multigrid methods. The interface curve is evolved
utilizing the $\theta-L$ variables instead of the more commonly used $x-y$
variables. This choice simplifies the preservation of mesh quality during the
interface evolution. In addition, the $\theta-L$ approach enables the design
of efficient and stable time-stepping schemes to remove the stiffness that
arises from the curvature term. Ample numerical examples, including
simulations of complex viscous fingering and dendritic solidification
problems, are presented to showcase the capability of the proposed method to
handle challenging moving interface problems.

###### keywords:

Hele-Shaw flow; The Stefan problem; Cartesian grid; Boundary integral
equations; Kernel-free boundary integral method; Small scale decomposition

##  1 Introduction

Moving interface problems are encountered in various fields of natural
sciences and industrial applications, ranging from mathematics [1, 2, 3] and
fluid mechanics [4, 5, 6] to material sciences [7, 8, 9, 10] and imaging
sciences [11, 12]. In these problems, interfaces are present, dividing the
surrounding region into different sub-regions where the underlying physics is
governed by PDEs.

When the motion of the interface is not known in advance and needs to be
determined as part of solving the entire problem, it is referred to as a free
boundary problem. Free boundary problems are inherently nonlinear and pose
computational challenges due to the coupling between the dynamics of the
interface and the underlying PDEs.

In this paper, we focus on two representative free boundary problems: the
Hele-Shaw flow [13, 14, 15] and the Stefan problem [9, 16, 17]. The Hele-Shaw
problem describes the flow of a viscous fluid in a thin gap between two
parallel plates and is an elliptic type problem. The Stefan problem, on the
other hand, models the solidification or melting of a material with a moving
interface and is a parabolic type problem. These problems have been
extensively studied in terms of their numerical methods and analysis, due to
their significance in applications and the computational difficulties they
pose.

Developing accurate and efficient numerical methods for moving interface
problems poses several challenges. The first challenge lies in accurately
approximating a complex and evolving interface. Although various numerical
methods have been developed, such as the front-tracking method [18, 19, 16],
the level set method [20, 21, 22, 23, 24], the volume of fluid method [25,
26], and phase-field methods [27, 28, 29], obtaining both simplicity and
accuracy in approximations remains challenging. The second challenge is
related to solving PDEs on complex and time-varying domains. Due to the
evolving interface, methods like the classical finite element method, which
relies on body-fitted meshes, require frequent re-meshing procedures to
maintain mesh quality. This not only leads to computational inefficiencies but
also requires additional implementation efforts. Additionally, conventional
boundary integral equation methods [17, 14, 30, 31, 32, 15] are efficient for
homogeneous PDEs but have limitations such as the need for dense matrix-vector
multiplications and evaluations of singular integrals. In recent decades,
Cartesian grid-based methods have gained popularity for moving interface
problems. These methods, including the immersed boundary method (IBM) [33, 34,
35], the immersed interface method (IIM) [36, 37, 38, 39], and the ghost fluid
method (GFM) [40, 41, 42, 43], involve immersing the moving interface into a
fixed background mesh, typically a Cartesian grid. This approach simplifies
the algorithm and improves computational efficiency.

When considering the effect of surface tension on the interface, the problem
formulation includes the Laplace-Young equation or the Gibbs-Thomson relation
to account for local curvature. However, the presence of high-order
derivatives in the curvature introduces stiffness into the evolution problem
and imposes strict stability constraints on the time step when using explicit
time-stepping schemes. Conversely, using a straightforward implicit scheme
becomes complicated and computationally expensive due to the nonlinear and
nonlocal nature of the interface velocity as a function of the interface
position. To address these challenges, Hou et al. developed the small scale
decomposition (SSD) method [14]. This method successfully removes stiffness
induced by surface tension by employing a special $\theta-L$ formulation of
the interface and an implicit discretization for the stiff but linear part of
the evolution equation. As a result, the SSD method allows for the application
of large time steps, improving computational efficiency. The SSD method has
been adapted for solving multiple moving interface problems, including
microstructural evolution in inhomogeneous elastic media [44], elastic
membranes in viscous flows [45, 46], solid tumor growth [47], and crystal
growth [48, 49, 50], among others.

The main objective of this paper is to develop an efficient and stable
Cartesian grid-based boundary integral method for solving the Hele-Shaw flow
and the Stefan problem. By reformulating the PDEs as boundary integral
equations, we use a finite difference analogue of the boundary integral method
known as the kernel-free boundary integral (KFBI) method. The KFBI method is
based on potential theory and is specifically designed to solve boundary value
problems of elliptic equations on irregular domains [51]. It takes advantage
of fast PDE solvers on Cartesian grids and eliminates the need to evaluate
singular and nearly singular integrals. The method has been extended to
higher-order versions and has been successfully applied to solve various
problems in the past [52, 53, 54]. To accurately track the moving interface,
the $\theta-L$ approach is utilized to represent the Jordan curve and achieve
equidistant interface discretization. In addition, the SSD method is employed
to enhance the efficiency and stability of numerical schemes for interface
evolutions in both the Hele-Shaw flow and the Stefan problem.

The remainder of this paper is organized as follows: In Section 2, we outline
the governing equations of the Hele-Shaw problem and the Stefan problem.
Section 3 introduces time discretization methods for time-dependent PDEs and
formulations of boundary integral equations. The kernel-free boundary integral
method is elaborated upon in Section 4. Numerical approaches for interface
evolution are described in Section 5. Numerical results are presented in
Section 6. Finally, in Section 7, we briefly discuss the advantages of the
proposed method and potential areas for improvement.

##  2 Moving interface problems

Let $\Gamma:[0,2\pi)\times[0,T]\rightarrow\mathbb{R}^{2}$ be a time-dependent
and closed curve separating a domain $\mathcal{U}\subset\mathbb{R}^{2}$ into
an interior domain $\Omega^{+}$ and an exterior domain $\Omega^{-}$.
Generally, the domain $\mathcal{U}$ can be bounded or unbounded. We shall
assume that the interior domain $\Omega^{+}$ is bounded and can be covered by
a bounding box $\mathcal{B}$, see Fig 1. In moving interface problems, the
physical quantities defined in $\Omega^{+}$ and $\Omega^{-}$ satisfy certain
PDEs. Interface or boundary conditions are also prescribed on the moving
interface $\Gamma$ for the PDEs.

Figure 1: A Schematic of a moving interface problem.

###  2.1 The Hele-Shaw flow

We consider an exterior moving interface problem of the Hele-Shaw flow that
describes the motion of an air bubble in a radial Hele-Shaw cell. The moving
interface $\Gamma$ separates the full space $\mathbb{R}^{2}$ into an air
domain $\Omega^{+}$ and an oil domain $\Omega^{-}$. In the air domain, the
pressure is assumed to be constant, which we may take to be zero. In the oil
domain, the fluid velocity $\mathbf{u}$ and pressure $p$ satisfy Darcy’s law,
together with the incompressibility constraint:

$\mathbf{u}=-M\nabla p,\quad\nabla\cdot\mathbf{u}=0,\quad\text{in
}\Omega^{-},$ (1)

where $M=\frac{b^{2}}{12\mu}$ is the mobility of the fluid, $b$ is the gap
width of the Hele-Shaw cell, and $\mu$ is the viscosity. Consider constant
injection of air at the origin,

$\nabla\cdot\mathbf{u}=2\pi J\delta(\mathbf{x}),\quad\text{in }\Omega^{-},$
(2)

where $J\geq 0$ is a constant injection rate and $\delta$ is the Dirac delta
function modeling the air injection. Combining equations (1) and (2), and
assuming $M=1$, we obtain the Poisson equation for the pressure

$\Delta p=-2\pi J\delta(\mathbf{x}),\quad\text{in }\Omega^{-}.$ (3)

Although the Dirac delta function in the right-hand side of (3) vanishes in
$\Omega^{-}$, it prescribes the behavior of the solution at infinity
$p=-J\ln|\mathbf{x}|+C+o(1)$ as $|\mathbf{x}|\rightarrow\infty$, where $C$ is
a constant ambient pressure. On the moving interface $\Gamma$, the pressure is
given by the Laplace-Young condition

$p=-\sigma\kappa,\quad\text{on }\Gamma,$ (4)

where $\sigma>0$ is the surface tension coefficient and $\kappa$ is the local
curvature of $\Gamma$. In addition, the motion of the moving interface
$\Gamma$ follows the kinematic condition

$\dfrac{d\mathbf{x}}{dt}=\mathbf{u},\quad\text{for }\mathbf{x}\in\Gamma.$ (5)

As a matter of fact, since we are only interested in the shape of $\Gamma$, it
suffices only to prescribe the normal velocity $V_{n}$,

$V_{n}=\dfrac{d\mathbf{x}}{dt}\cdot\mathbf{n}=\mathbf{u}\cdot\mathbf{n}=-\partial_{\mathbf{n}}p,\quad\text{for
}\mathbf{x}\in\Gamma.$ (6)

where $\mathbf{n}$ is the unit outward normal of $\Gamma$.

###  2.2 The Stefan problem

We also consider the Stefan problem that models diffusion-driven phase changes
between solid and liquid phases. Here, the moving solid-liquid interface
$\Gamma$ separates the box domain $\mathcal{B}$ into the solid region
$\Omega^{+}$ and the liquid region $\Omega^{-}$. In the classical Stefan
problem, the temperature field $T$ satisfies the heat equation in both solid
and liquid regions,

$\displaystyle\partial_{t}T=\Delta T,\quad\text{in }\Omega^{+}\cup\Omega^{-}.$
(7)

On the solid-liquid interface, the temperature is continuous and is coupled
with the surface tension and the molecular kinematic effect through the Gibbs-
Thomson relation

$\displaystyle
T+\varepsilon_{C}(\mathbf{n})\kappa+\varepsilon_{V}(\mathbf{n})V_{n}=0,\quad\text{on
}\Gamma,$ (8)

where $\varepsilon_{C}(\mathbf{n})$ and $\varepsilon_{V}(\mathbf{n})$ are
surface tension and molecular kinetic coefficients, which are non-negative and
may depend on the orientation of the interface. The normal velocity of
$\Gamma$ is determined by the Stefan equation.

$\displaystyle V_{n}=[\partial_{\mathbf{n}}T].$ (9)

Here, the notation $[\cdot]$ means the jump value of a quantity across the
interface. For example, given a piece-wise continuous function $q$, we have

$[q](\mathbf{x})=q^{+}(\mathbf{x})-q^{-}(\mathbf{x})=\lim_{\mathbf{y}\in\Omega^{+},\mathbf{y}\rightarrow\mathbf{x}}q(\mathbf{y})-\lim_{\mathbf{y}\in\Omega^{-},\mathbf{y}\rightarrow\mathbf{x}}q(\mathbf{y}),\quad\text{for
}\mathbf{x}\in\Gamma.$ (10)

A suitable boundary condition should also be prescribed on the outer boundary
$\partial\mathcal{B}$, for which we choose a no-flux boundary condition
$\partial_{\mathbf{n}}T=0$. For modeling solidification problems, initially, a
solid seed is placed in an undercooled surrounding liquid. The temperature in
the solid seed is assumed to be equal to the melt temperature $T_{m}$, and the
temperature in the undercooled liquid is given by $T_{\infty}$. The degree of
undercooling is described by the Stefan number $St=T_{\infty}-T_{m}$.

In a more general setting, the temperature differences in the liquid phase
lead to changes in the specific volume of fluid parcels and, hence, the fluid
density. The density changes further lead to buoyancy force-driven convection
of the fluid. With the Boussinesq approximation [55], we assume the fluid is
incompressible, and the effect of density changes only appears in the buoyancy
force. Incorporating the natural convection effect, the Stefan problem is
governed by the following system

$\displaystyle\partial_{t}T=\Delta T,$ $\displaystyle\quad\text{in
}\Omega^{+},$ (11) $\displaystyle\partial_{t}T+\mathbf{u}\cdot\nabla T=\Delta
T,$ $\displaystyle\quad\text{in }\Omega^{-},$ (12)
$\displaystyle\partial_{t}\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}=\Delta\mathbf{u}-\nabla
p+\mathbf{G},$ $\displaystyle\quad\text{in }\Omega^{-},$ (13)
$\displaystyle\nabla\cdot\mathbf{u}=0,$ $\displaystyle\quad\text{in
}\Omega^{-},$ (14)

where $\mathbf{u}$ is fluid velocity, $p$ is pressure. Here, we have assumed
that the density and viscosity of the fluid are $1$. The buoyancy force
$\mathbf{G}$ under the Boussinesq approximation is linearly proportional to
the temperature difference,

$\mathbf{G}=-g\beta(T-T_{\infty})\mathbf{j},$ (15)

where $g$ represents the gravity acceleration, $\beta$ denotes the thermal
expansion coefficient, $\mathbf{j}$ denotes the unit vector in the vertical
direction. On the solid-liquid interface $\Gamma$, the Gibbs-Thomson relation
(8) and the Stefan equation (9) are given for the temperature field, and the
no-slip boundary condition is specified for the fluid. On the outer boundary
$\partial\mathcal{B}$, we set $\partial_{\mathbf{n}}T=0$ and
$\mathbf{u}=\mathbf{u}_{b}$ where $\mathbf{u}_{b}$ is the boundary data
describing in/out-flow or no-slip boundary conditions. It is worth mentioning
that the model is a two-way coupling of the temperature and the velocity field
through the fluid convection and buoyancy force.

##  3 Time discretization and boundary integral equations

The PDEs in the last section are all solved in a boundary integral equation
framework. Time-dependent parabolic PDEs are discretized in time to obtain a
sequence of elliptic PDEs, for which we can formulate equivalent boundary
integral equations. Let $t_{n}=n\tau$, $n=0,1,\cdots,N_{T}$ be the uniform
temporal mesh where $\tau=T/N_{T}$ is the time step. For a function $f$,
denote by $f^{n}$ the numerical approximation of $f(t_{n})$.

###  3.1 The Hele-Shaw flow

In the Hele-Shaw flow, the solution to the Poisson equation (3) with the
boundary condition (4) can be represented as the sum of two functions $v$ and
$w$, where

$\displaystyle v(\mathbf{x})=-J\ln|\mathbf{x}|.$ (16)

comes from the point-source term and $w(\mathbf{x})$ satisfies the exterior
problem of the Laplace equation

$\displaystyle\Delta w=0,$ $\displaystyle\quad\text{in }\Omega^{-},$ (17)
$\displaystyle w=-\sigma\kappa-v,$ $\displaystyle\quad\text{on }\Gamma.$ (18)

Let $G_{0}(\mathbf{y},\mathbf{x})=(1/2\pi)\ln|\mathbf{y}-\mathbf{x}|$ be the
free-space Green function associated with the Laplacian $\Delta$, the solution
$w(\mathbf{x})$ can be represented as a modified double-layer potential

$w(\mathbf{x})=(D\varphi)(\mathbf{x})+\int_{\Gamma}\varphi(\mathbf{y})\,d\mathbf{s}_{\mathbf{y}}=\int_{\Gamma}\varphi(\mathbf{y})\left(\dfrac{\partial
G_{0}(\mathbf{y},\mathbf{x})}{\partial\mathbf{n}_{\mathbf{y}}}+1\right)\,d\mathbf{s}_{\mathbf{y}},$
(19)

where $\varphi$ is an unknown dipole density function defined on $\Gamma$. The
boundary integral formulation of $w$ naturally matches the boundary condition
at infinity. Restricting (19) on $\Gamma$ and using the boundary condition
(18), we obtain a boundary integral equation for the density function
$\varphi$,

$-\dfrac{1}{2}\varphi(\mathbf{x})+\int_{\Gamma}\varphi(\mathbf{y})\left(\dfrac{\partial
G_{0}(\mathbf{y},\mathbf{x})}{\partial\mathbf{n}_{\mathbf{y}}}+1\right)\,d\mathbf{s}_{\mathbf{y}}=-\sigma\kappa-v(\mathbf{x}),\quad\text{for
}\mathbf{x}\in\Gamma.$ (20)

The boundary integral equation is a Fredholm integral equation of the second
kind and is well-conditioned.

###  3.2 The classical Stefan problem

We first discretize the time-dependent heat equation in time to reduce the
problem into solving an elliptic equation in each time step. For better
accuracy and stability, the second-order backward differentiation formula
(BDF2) is employed for the time discretization

$\displaystyle\dfrac{3T^{n+1}-4T^{n}+T^{n-1}}{2\tau}=\Delta
T^{n+1},\quad\text{in }\Omega^{+}\cup\Omega^{-}.$ (21)

It leads to a modified Helmholtz equation for $T^{n+1}$

$\displaystyle(\Delta-\dfrac{3}{2\tau})T^{n+1}=\dfrac{T^{n-1}-4T^{n}}{2\tau},\quad\text{in
}\Omega^{+}\cup\Omega^{-},$ (22)

The modified Helmholtz equation is also subject to jump conditions
$[T^{n+1}]=0$ and $[\partial_{\mathbf{n}}T^{n+1}]=V_{n}$ on $\Gamma$ and an
outer boundary condition $\partial_{\mathbf{n}}T^{n+1}=0$ on
$\partial\mathcal{B}$. Let $c=\sqrt{3/(2\tau)}$. We split $T^{n+1}$ into two
parts $T^{n+1}=T_{1}+T_{2}$, where $T_{1}$ is the solution to the modified
Helmholtz equation with an inhomogeneous right-hand side,

$\displaystyle\Delta T_{1}-c^{2}T_{1}=\dfrac{T^{n-1}-4T^{n}}{2\tau},\quad$
$\displaystyle\text{in }\mathcal{B},$ (23)
$\displaystyle\partial_{\mathbf{n}}T_{1}=0,\quad$ $\displaystyle\text{on
}\partial\mathcal{B},$ (24)

and $T_{2}$ is the solution to the interface problem with a homogeneous right-
hand side,

$\displaystyle\Delta T_{2}-c^{2}T_{2}=0,\quad$ $\displaystyle\text{in
}\Omega^{+}\cup\Omega^{-},$ (25) $\displaystyle[T_{2}]=0,\quad$
$\displaystyle\text{on }\Gamma,$ (26)
$\displaystyle[\partial_{\mathbf{n}}T_{2}]=V_{n},\quad$ $\displaystyle\text{on
}\Gamma,$ (27) $\displaystyle\partial_{\mathbf{n}}T_{2}=0,\quad$
$\displaystyle\text{on }\partial\mathcal{B}.$ (28)

Since the temperature field $T$ is continuous across the interface $\Gamma$,
the right-hand side of (23) is also continuous. Therefore, the part $T_{1}$ is
much more regular than $T_{2}$, and the problem (23) and (24) can be solved
with a standard finite difference method. For each fixed
$\mathbf{x}\in\mathcal{B}$, let $G_{c}(\mathbf{y},\mathbf{x})$ be the Green
function that satisfies

$\displaystyle(\Delta-c^{2})G_{c}(\mathbf{y},\mathbf{x})=\delta(\mathbf{y}-\mathbf{x}),\quad$
$\displaystyle\text{for }\mathbf{y}\in\mathcal{B},$ (29)
$\displaystyle\partial_{\mathbf{n}}G_{c}(\mathbf{y},\mathbf{x})=0,\quad$
$\displaystyle\text{for }\mathbf{y}\in\partial\mathcal{B}.$ (30)

The part $T_{2}$ can be expressed as a single layer potential

$\displaystyle
T_{2}(\mathbf{x})=-(S\psi)(\mathbf{x})=-\int_{\Gamma}G_{c}(\mathbf{y},\mathbf{x})\psi(\mathbf{y})\,d\mathbf{s}_{y},\quad\text{for
}\mathbf{x}\in\mathcal{B}.$ (31)

where $\psi$ is an unknown density function defined on $\Gamma$. Note that we
also have $\psi=[\partial_{\mathbf{n}}T_{2}]=V_{n}$. Further, with the Gibbs-
Thomson relation (8), a boundary integral equation can be obtained for the
density function $\psi$,

$\displaystyle\varepsilon_{V}\psi(\mathbf{x})-\int_{\Gamma}G_{c}(\mathbf{y},\mathbf{x})\psi(\mathbf{y})\,d\mathbf{s}_{y}=-\varepsilon_{C}\kappa(\mathbf{x})-T_{1}(\mathbf{x}),\quad\text{for
}\mathbf{x}\in\Gamma.$ (32)

For a nonzero $\varepsilon_{V}$, the boundary integral equation (32) is a
Fredholm integral equation of the second kind. It degenerates to a first-kind
Fredholm integral equation if $\varepsilon_{V}$ vanishes.

###  3.3 The Stefan problem with natural convection

Similar to the heat equation, we first discretize time derivatives in the
time-dependent advection-diffusion equation (12) and the Navier-Stokes
equation (13), (14), and reduce the problems into elliptic equations in each
time step. To avoid treating the nonlinear advection term, we advance
$T,\mathbf{u},p$ with a semi-Lagrangian method in time. In the semi-Lagrangian
discretization, the advection terms are incorporated into material
derivatives, namely,

$\displaystyle\dfrac{dT}{dt}=\Delta T,$ $\displaystyle\quad\text{in
}\Omega^{+}\cup\Omega^{-},$ (33)
$\displaystyle\dfrac{d\mathbf{u}}{dt}=\Delta\mathbf{u}-\nabla p+\mathbf{G},$
$\displaystyle\quad\text{in }\Omega^{-},$ (34)
$\displaystyle\nabla\cdot\mathbf{u}=0,$ $\displaystyle\quad\text{in
}\Omega^{-},$ (35)

where $\frac{d}{dt}=\partial_{t}+\mathbf{u}\cdot\nabla$ is the material
derivative. Here, we assume $\mathbf{u}$ equals zero in the solid domain, so
the temperature in both domains satisfies the same equation (33). A second-
order semi-implicit is used to discretize equations 33, 34 and 35,

$\displaystyle\dfrac{3T^{n+1}-4\tilde{T}^{n}+\tilde{T}^{n-1}}{2\tau}=\Delta
T^{n+1},$ $\displaystyle\quad\text{in }\Omega^{+}\cup\Omega^{-},$ (36)
$\displaystyle\dfrac{3\mathbf{u}^{n+1}-4\tilde{\mathbf{u}}^{n}+\tilde{\mathbf{u}}^{n-1}}{2\tau}=\Delta\mathbf{u}^{n+1}-\nabla
p^{n+1}+2\mathbf{G}^{n}-\mathbf{G}^{n-1},$ $\displaystyle\quad\text{in
}\Omega^{-},$ (37) $\displaystyle\nabla\cdot\mathbf{u}^{n+1}=0,$
$\displaystyle\quad\text{in }\Omega^{-}.$ (38)

where $\tilde{\mathbf{u}}^{n}$, $\tilde{T}^{n}$, $\tilde{\mathbf{u}}^{n-1}$,
and $\tilde{T}^{n-1}$ are the fluid velocities and temperatures at the
departure points $\mathbf{x}^{n}$ and $\mathbf{x}^{n-1}$, respectively. The
scheme treats the buoyancy term explicitly such that the resulting discrete
system is not coupled and can be solved separately. The departure points
$\mathbf{x}^{n}$ and $\mathbf{x}^{n-1}$ can be found by solving the initial
value problem backward in time,

$\displaystyle\dfrac{d\mathbf{x}(t)}{dt}=\mathbf{u}(\mathbf{x}(t),t),\quad\mathbf{x}(t_{n+1})=\mathbf{x}_{0}.$
(39)

A second-order mid-point method is used for computing the positions of the
departure points,

$\displaystyle\mathbf{x}^{*}=\mathbf{x}_{0}-\frac{\tau}{2}\mathbf{u}\left(\mathbf{x}_{0}-\frac{\tau}{2}\mathbf{u}^{n+\frac{1}{2}},t_{n+\frac{1}{2}}\right),\quad\mathbf{x}^{n}=\mathbf{x}_{0}-\tau\mathbf{u}\left(\mathbf{x}^{*},t_{n+\frac{1}{2}}\right),$
(40)
$\displaystyle\mathbf{x}^{*}=\mathbf{x}_{0}-\tau\mathbf{u}\left(\mathbf{x}_{0}-\tau\mathbf{u}^{n},t_{n}\right),\quad\mathbf{x}^{n-1}=\mathbf{x}_{0}-2\tau\mathbf{u}\left(\mathbf{x}^{*},t_{n}\right).$
(41)

Off-grid values of $\mathbf{u}$ are computed with the cubic Lagrangian
interpolation and the velocity at $t_{n+\frac{1}{2}}$ is computed with a
second-order extrapolation scheme
$\mathbf{u}^{n+\frac{1}{2}}=\frac{3}{2}\mathbf{u}^{n}-\frac{1}{2}\mathbf{u}^{n-1}$.
After making some rearrangements, we arrive at the modified Helmholtz equation

$(\Delta-\dfrac{3}{2\tau})T^{n+1}=\dfrac{\tilde{T}^{n-1}-4\tilde{T}^{n}}{2\tau},\quad\text{in
}\Omega^{+}\cup\Omega^{-},,$ (42)

subject to interface and boundary conditions on $\Gamma$ and
$\partial\mathcal{B}$ as mentioned before, and the modified Stokes equation

$\displaystyle(\Delta-\dfrac{3}{2\tau})\mathbf{u}^{n+1}-\nabla
p^{n+1}=\dfrac{\tilde{\mathbf{u}}^{n-1}-4\tilde{\mathbf{u}}^{n}}{2\tau}+2\mathbf{G}^{n}-\mathbf{G}^{n-1},$
$\displaystyle\quad\text{in }\Omega^{-},$ (43a)
$\displaystyle\nabla\cdot\mathbf{u}^{n+1}=0,$ $\displaystyle\quad\text{in
}\Omega^{-},$ (43b)

subject to boundary conditions

$\displaystyle\mathbf{u}^{n+1}=\mathbf{0},\quad\text{on
}\Gamma,\quad\mathbf{u}^{n+1}=\mathbf{u}_{b},\quad\text{on
}\partial\mathcal{B},$ (44)

The boundary integral formulation for the modified Helmholtz (42) is similar
to that for (22), and we omit it. We shall give the boundary integral equation
for the modified Stokes (43). We split the solution pair into two parts
$(\mathbf{u}^{n+1},p^{n+1})=(\mathbf{u}_{1},p_{1})+(\mathbf{u}_{2},p_{2})$.
The particular solution $(\mathbf{u}_{1},p_{1})$ satisfies the modified Stokes
equation with an inhomogeneous right-hand side

$\displaystyle(\Delta-c^{2})\mathbf{u}_{1}-\nabla p_{1}=\mathbf{f},$
$\displaystyle\quad\text{in }\mathcal{B},$ (45a)
$\displaystyle\nabla\cdot\mathbf{u}_{1}=0,$ $\displaystyle\quad\text{in
}\mathcal{B},$ (45b) $\displaystyle\mathbf{u}_{1}=\mathbf{u}_{b},$
$\displaystyle\quad\text{on }\partial\mathcal{B},$ (45c)

where $c=\sqrt{3/(2\tau)}$ and $\mathbf{f}$ is given by

$\mathbf{f}=\left\\{\begin{aligned}
&\dfrac{\tilde{\mathbf{u}}^{n-1}-4\tilde{\mathbf{u}}^{n}}{2\tau}+2\mathbf{G}^{n}-\mathbf{G}^{n-1},\quad\text{in
}\Omega^{-},\\\ &2\mathbf{G}^{n}-\mathbf{G}^{n-1},\quad\text{in
}\Omega^{+}.\end{aligned}\right.$ (46)

Note that $\mathbf{f}$ is a continuous extension of the right-hand side of
(43a) due to the no-slip boundary condition on $\Gamma$. Then the part
$(\mathbf{u}_{1},p_{1})$ has high regularity, and standard finite difference
methods for the Stokes equation can be applied. The second part
$(\mathbf{u}_{2},p_{2})$ satisfies an exterior Dirichlet boundary value
problem

$\displaystyle(\Delta-c^{2})\mathbf{u}_{2}-\nabla p_{2}=\mathbf{0},$
$\displaystyle\quad\text{in }\Omega^{-},$ (47)
$\displaystyle\nabla\cdot\mathbf{u}_{2}=0,$ $\displaystyle\quad\text{in
}\Omega^{-},$ (48) $\displaystyle\mathbf{u}_{2}=-\mathbf{u}_{1},$
$\displaystyle\quad\text{on }\Gamma,$ (49)
$\displaystyle\mathbf{u}_{2}=\mathbf{0},$ $\displaystyle\quad\text{on
}\partial\mathcal{B}.$ (50)

For each fixed $\mathbf{x}\in\mathcal{B}$, let
$(\mathbf{G}_{\mathbf{u}}(\mathbf{y},\mathbf{x}),\mathbf{G}_{p}(\mathbf{y},\mathbf{x}))$
be the Green function pair that satisfies

$\displaystyle(\Delta-c^{2})\mathbf{G}_{\mathbf{u}}(\mathbf{y},\mathbf{x})-\nabla\mathbf{G}_{p}(\mathbf{y},\mathbf{x})=\delta(\mathbf{y}-\mathbf{x})\mathbf{I},\quad$
$\displaystyle\text{for }\mathbf{y}\in\mathcal{B},$ (51)
$\displaystyle\nabla\cdot\mathbf{G}_{p}(\mathbf{y},\mathbf{x})=0,\quad$
$\displaystyle\text{for }\mathbf{y}\in\mathcal{B},$ (52)
$\displaystyle\mathbf{G}_{\mathbf{u}}(\mathbf{y},\mathbf{x})=\mathbf{0},\quad$
$\displaystyle\text{for }\mathbf{y}\in\partial\mathcal{B}.$ (53)

Then the solution pair $(\mathbf{u}_{2},p_{2})$ can be represented as a
double-layer potential
$D\bm{\varphi}=(D_{\mathbf{u}}\bm{\varphi},D_{p}\bm{\varphi})^{T}$ [56],

$\displaystyle\mathbf{u}_{2}(\mathbf{x})$
$\displaystyle=(D_{\mathbf{u}}\bm{\varphi})(\mathbf{x})=\int_{\Gamma}T(\mathbf{G}_{\mathbf{u}},\mathbf{G}_{p})(\mathbf{y},\mathbf{x})\bm{\varphi}(\mathbf{y})\,d\mathbf{s}_{\mathbf{y}},$
(54) $\displaystyle p_{2}(\mathbf{x})$
$\displaystyle=(D_{p}\bm{\varphi})(\mathbf{x})=2\int_{\Gamma}\dfrac{\partial\mathbf{G}_{p}(\mathbf{y},\mathbf{x})}{\partial\mathbf{n}_{\mathbf{y}}}\bm{\varphi}(\mathbf{y})\,d\mathbf{s}_{\mathbf{y}}.$
(55)

where $\bm{\varphi}=(\varphi_{1},\varphi_{2})^{T}$ is a vector-valued unknown
density function defined on $\Gamma$ and
$T(\mathbf{u},p)=-p\mathbf{n}+(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})\mathbf{n}$
is the traction operator. Restricting (54) to $\Gamma$ leads to the boundary
integral equation

$\displaystyle-\dfrac{1}{2}\bm{\varphi}(\mathbf{x})+\int_{\Gamma}T(\mathbf{G}_{\mathbf{u}},\mathbf{G}_{p})(\mathbf{y},\mathbf{x})\bm{\varphi}(\mathbf{y})\,d\mathbf{s}_{\mathbf{y}}=-\mathbf{u}_{1}(\mathbf{x}),\quad\text{for
}\mathbf{x}\in\Gamma.$ (56)

The boundary integral equation is also a Fredholm integral equation of the
second kind and is well-conditioned.

##  4 Kernel-free boundary integral method

The boundary integral equations (20), (32) and (56) are solved with the
kernel-free boundary integral method. The discrete systems of boundary
integral equations are solved with a Krylov subspace iterative method, the
GMRES method [57]. In each iteration, one only needs to compute the matrix-
vector multiplication operation, which mainly consists of evaluations of
boundary integrals. The procedure can be implemented in a matrix-free manner
to avoid forming the full matrix. The main idea of the KFBI method is to make
use of a Cartesian grid-based PDE solver instead of numerical quadratures for
evaluating the boundary integrals.

###  4.1 Equivalent interface problems

Let $\mathcal{A}$ be an elliptic differential operator that can be the
Laplacian, modified Helmholtz, or modified Stokes operator. The evaluation of
boundary integral operators associated with $\mathcal{A}$ can be described in
the same framework. According to the classical potential theory, the single-
layer potential $S\bm{\psi}(\mathbf{x})$ and the double-layer potential
$D\bm{\varphi}(\mathbf{x})$ satisfy equivalent interface problems [56, 51],
which can be unified as

$\left\\{\begin{aligned} &\mathcal{A}\mathbf{v}=\mathbf{0},&\quad\text{in
}\Omega^{+}\cup\Omega^{-},\\\
&[\bm{\pi}_{D}(\mathbf{v})]=\mathbf{\Phi},&\quad\text{on }\Gamma,\\\
&[\bm{\pi}_{N}(\mathbf{v})]=\mathbf{\Psi},&\quad\text{on }\Gamma,\\\
&\text{some suitable BCs},&\quad\partial\mathcal{B}.\end{aligned}\right.$ (57)

where $(\bm{\pi}_{D},\bm{\pi}_{N})$ is the Cauchy data pair and is specified
in (1).

Table 1: Cauchy data pair for different elliptic differential operators. | $\mathbf{v}$ | $\bm{\pi}_{D}(\mathbf{v})$ | $\bm{\pi}_{N}(\mathbf{v})$
---|---|---|---
$\mathcal{A}=\Delta-c^{2}$ | $u$ | $u$ | $\partial_{\mathbf{n}}u$
$\mathcal{A}=((\Delta-c^{2})\mathbf{I}-\nabla,\nabla\cdot)^{T}$ | $(\mathbf{u},p)^{T}$ | $\mathbf{u}$ | $T(\mathbf{u},p)$

The functions $\bm{\Phi}$ and $\bm{\Psi}$ are given in (2).

Table 2: Jump relations of the single- and double-layer potentials. | $\bm{\Phi}$ | $\bm{\Psi}$
---|---|---
$\mathbf{v}=S\bm{\psi}$ | $\mathbf{0}$ | $-\bm{\psi}$
$\mathbf{v}=D\bm{\varphi}$ | $\bm{\varphi}$ | $\mathbf{0}$

The boundary condition on $\partial\mathcal{B}$ depends on which Green’s
function is used. For the Hele-Shaw flow, since the problem is defined on an
unbounded domain and the box boundary $\partial\mathcal{B}$ is only an
artificial boundary. A natural choice is to directly use the potential value
on $\partial\mathcal{B}$ as a Dirichlet-type boundary condition. The integral
form boundary condition is discretized with the composite trapezoidal rule.
Due to the periodicity of the integrand function on $\Gamma$, the composite
trapezoidal rule is highly accurate. No singular or nearly singular integral
needs to be handled here since the quadrature points are on $\Gamma$ and
target points are on $\partial\mathcal{B}$. In the Stefan problem, boundary
conditions for the temperature or the fluid field are of Dirichlet or Neumann-
type, depending on the problem at hand.

The equivalent interface problem (57) is much simpler than the original
problems since the interface condition can be easily decoupled into derivative
jumps in each direction, which is shown in A. After solving the interface
problem (57) on a Cartesian grid with an efficient PDE solver, boundary
integral values, as well as their normal derivatives, can be extracted from
grid values with an interpolation procedure.

We stress that the analytical expression of Green’s function is only used in
the Hele-Shaw flow to convert the infinite boundary condition to a bounded
one. In other cases, the method is kernel-free in the sense that there is no
need to use any of Green’s functions explicitly. In all cases, evaluations of
singular or nearly singular integrals are avoided. Accelerated by fast PDE
solvers, the computation of boundary integrals with the KFBI method is
efficient and is comparable with the fast multipole method [52].

###  4.2 Corrected finite difference scheme

The interface problem (57) is solved with a finite difference method with
additional corrections near the interface. For simplicity, suppose the
computational domain $\mathcal{B}$ is a unit square $[0,1]^{2}$ and is
uniformly partitioned into $N$ intervals in each spatial direction. The grid
nodes are denoted as $(x_{i},y_{j})$, $x_{i}=ih$, $y_{j}=jh$, $i,j=0,\cdots,N$
where $h=1/N$ is the mesh size.

####  4.2.1 The modified Helmholtz equation

Denote by $u_{i,j}$ the numerical approximation of $u(x_{i},y_{j})$. A
standard second-order five-point central difference scheme for the modified
Helmholtz equation is given by

$\displaystyle\Delta_{h}u_{i,j}-c^{2}u_{i,j}=\dfrac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}}{h^{2}}-c^{2}u_{i,j}=0.$
(58)

A grid node $(x_{i},y_{j})$ is called an irregular node if at least one of the
stencil nodes is on the other side of the interface. Clearly, irregular nodes
are always in the vicinity of $\Gamma$. Since the solution has certain jumps
on the interface $\Gamma$, the local truncation error at an irregular node is
on the order of $\mathcal{O}(h^{-2})$ and is too large to achieve an accurate
approximation. To show this, we assume at an irregular node
$(x_{i},y_{j})\in\Omega^{+}$, a stencil node $(x_{i+1},y_{j})\in\Omega^{-}$ is
on the different side of $\Gamma$ and other stencil nodes are all in
$\Omega^{+}$. Suppose the interface $\Gamma$ intersects the grid lines segment
between $(x_{i},y_{j})$ and $(x_{i+1},y_{j})$ at $(\xi,y_{j})$. By Taylor
expansion, the local truncation error at $(x_{i},y_{j})$ is given by

$\displaystyle
E_{h}(x_{i},y_{j})=(\Delta_{h}-c^{2})u(x_{i},y_{j})=-\dfrac{1}{h^{2}}\\{[u]+(x_{i+1}-\xi)[u_{x}]+\dfrac{1}{2}(x_{i+1}-\xi)^{2}[u_{xx}])\\}+\mathcal{O}(h).$
(59)

To avoid the large local truncation error, one can use the leading term as a
correction term to the right-hand side of (58). It gives the corrected scheme

$\displaystyle\Delta_{h}u_{i,j}-c^{2}u_{i,j}=C_{i,j}=-\dfrac{1}{h^{2}}\\{[u]+(x_{i+1}-\xi)[u_{x}]+\dfrac{1}{2}(x_{i+1}-\xi)^{2}[u_{xx}])\\}.$
(60)

The resulting local truncation error becomes $\mathcal{O}(h)$. Here, the
correction term $C_{i,j}$ is from the contribution of the intersection point
$(\xi,y_{j})$ and is a linear combination of jump values $[u]$, $[u_{x}]$,
and$[u_{yy}]$, which are known before solving the interface problem, see A.
Similarly, one can derive correction terms for different intersection patterns
by Taylor expansions. Evidently, $C_{i,j}$ is non-zero only at irregular nodes
and is always a linear combination of the jump values
$[u],[u_{x}],[u_{y}],\cdots$.

Since irregular nodes are near the interface $\Gamma$, which has co-dimension
one, the overall accuracy can still be second-order [58]. The correction terms
only appear on the right-hand side and the coefficient matrix is not altered,
FFT-based fast Poisson solvers can be applied to solve the resulting linear
system efficiently.

####  4.2.2 The modified Stokes equation

The discretization of the modified Stokes equations is based on the marker and
cell (MAC) scheme on a staggered grid. The pressure is at the cell center, the
x-component of the velocity is at the center of the east and west edges of a
cell, and the y-component of the velocity is at the center of the north and
south edges of a cell. The discrete mesh functions $u_{i,j},v_{i,j},p_{i,j}$
are defined as

$\displaystyle u_{i,j}=u(x_{i},y_{j}-h/2),\quad i=0,1\cdots,N_{x},\quad
j=1,2\cdots,N_{y},$ (61) $\displaystyle v_{i,j}=v(x_{i}-h/2,y_{j}),\quad
i=1,2\cdots,N_{x},\quad j=0,1\cdots,N_{y},$ (62) $\displaystyle
p_{i,j}=p(x_{i}-h/2,y_{j}-h/2),\quad i=1,2\cdots,N_{x},\quad
j=1,2\cdots,N_{y}.$ (63)

The MAC scheme for the modified Stokes equation is given by

$\displaystyle(\Delta_{h}-c^{2})\mathbf{u}_{i,j}-\nabla_{h}p_{i,j}$
$\displaystyle=\mathbf{0},$ (64a) $\displaystyle\nabla_{h}\cdot\mathbf{u}$
$\displaystyle=0,$ (64b)

where $\Delta_{h},\nabla_{h}$ are the discrete approximations for
$\Delta,\nabla$, respectively,

$\displaystyle\Delta_{h}\mathbf{u}_{ij}=\dfrac{\mathbf{u}_{i+1,j}+\mathbf{u}_{i-1,j}+\mathbf{u}_{i,j+1}+\mathbf{u}_{i,j-1}-4\mathbf{u}_{i,j}}{h^{2}},$
(65)
$\displaystyle\nabla_{h}p_{i,j}=(\dfrac{p_{i+1,j}-p_{i,j}}{h},\dfrac{p_{i,j+1}-p_{i,j}}{h})^{T},$
(66)
$\displaystyle\nabla_{h}\cdot\mathbf{u}_{i,j}=\dfrac{u_{i,j}-u_{i-1,j}}{h}+\dfrac{v_{i,j}-v_{i,j-1}}{h}.$
(67)

The Dirichlet boundary condition is discretized with a symmetric approach for
the $x$-component velocity on the top and bottom boundary and the
$y$-component velocity on the right and left boundary. In the presence of
interfaces, one can also derive local truncation errors at irregular nodes for
each finite difference equation in (64). Similarly, by using the leading term
in the local truncation error as a correction term of the right-hand side in
(64), we arrived at the corrected MAC scheme

$\displaystyle(\Delta_{h}-c^{2})\mathbf{u}_{i,j}-\nabla_{h}p_{i,j}$
$\displaystyle=\mathbf{C}_{i,j},$ (68a)
$\displaystyle\nabla_{h}\cdot\mathbf{u}$ $\displaystyle=D_{i,j},$ (68b)

where $\mathbf{C}_{i,j}$ and $D_{i,j}$ are correction terms, which are linear
combinations of jump values $[u]$, $[v]$, $[p]$, $[u_{x}]$, $[u_{y}]$,
$\cdots$. The resulting linear system (68) is solved with an efficient V-cycle
geometric multigrid method with the Distributive Gauss-Seidel (DGS) smoother
[59, 60].

######  Remark 4.1.

Since particular solution parts have higher regularity, classical numerical
approaches, such as finite difference methods, finite element methods, and
spectral methods, can be applied. In this work, the particular parts are also
solved with a five-point central difference scheme and the MAC scheme without
any correction terms near the interface.

###  4.3 Extracting boundary values from the grid data

The single- and double-layer potential functions $S\bm{\psi}$ and
$D\bm{\varphi}$ defined by boundary integrals are not smooth on $\Gamma$.
Define the single- and double-layer boundary integral operators
$\mathcal{L}\bm{\psi}$ and $\mathcal{K}\bm{\varphi}$ as values of the
potential functions $S\bm{\psi}$ and $D\bm{\varphi}$ on $\Gamma$,
respectively. By potential theory, for smooth $\Gamma$, the integral operators
satisfy

$\displaystyle\mathcal{L}\bm{\psi}(\mathbf{x})$
$\displaystyle=\dfrac{1}{2}((S\bm{\psi})^{+}(\mathbf{x})+(S\bm{\psi})^{-}(\mathbf{x})),\quad\mathbf{x}\in\Gamma,$
(69) $\displaystyle\mathcal{K}\bm{\varphi}(\mathbf{x})$
$\displaystyle=\dfrac{1}{2}((D\bm{\varphi})^{+}(\mathbf{x})+(D\bm{\varphi})^{-}(\mathbf{x})),\quad\mathbf{x}\in\Gamma.$
(70)

The normal derivatives of the potential functions
$\partial_{\mathbf{n}}(S\bm{\psi})$ and $\partial_{\mathbf{n}}(D\bm{\varphi})$
have similar expressions

$\displaystyle\partial_{\mathbf{n}}(S\bm{\psi})(\mathbf{x})$
$\displaystyle=\mathbf{n}(\mathbf{x})\cdot\dfrac{1}{2}((\nabla
S\bm{\psi})^{+}(\mathbf{x})+(\nabla
S\bm{\psi})^{-}(\mathbf{x})),\quad\mathbf{x}\in\Gamma,$ (71)
$\displaystyle\partial_{\mathbf{n}}(D\bm{\varphi})(\mathbf{x})$
$\displaystyle=\mathbf{n}(\mathbf{x})\cdot\dfrac{1}{2}((\nabla
D\bm{\varphi})^{+}(\mathbf{x})+(\nabla
D\bm{\varphi})^{-}(\mathbf{x})),\quad\mathbf{x}\in\Gamma.$ (72)

The normal derivatives of potential functions are closely related to adjoint
double-layer integral and hyper-singular integral operators [56]. As long as
the grid-valued single- and double-layer potential functions
$(S\bm{\psi})_{h}$ and $(D\bm{\varphi})_{h}$ are obtained, boundary integral
operators can be computed by extracting boundary values and normal derivatives
of the potential functions from grid data via polynomial interpolation

In doing the interpolation, one should take into account the jump values of
the potential functions. Suppose one needs to evaluate the one-sided limit of
a discontinuous function $v(\mathbf{x})$ at a point
$\mathbf{p}=(p_{1},p_{2})^{T}\in\Gamma$ from the side of $\Omega^{+}$. Denote
by $v^{+}(\mathbf{x})$ and $v^{-}(\mathbf{x})$ two functions that coincide
with $v(\mathbf{x})$ in $\Omega^{+}$ and $\Omega^{-}$, respectively. Given a
set of interpolation points $\\{\mathbf{q}_{m}\\}_{m=1}^{M}$, by Taylor
expansion at $\mathbf{p}$, we have

$\displaystyle v(\mathbf{q}_{m})$
$\displaystyle=v^{+}(\mathbf{p})+(\Delta\mathbf{x}_{m})^{T}\nabla
v^{+}(\mathbf{p})+\dfrac{1}{2}(\Delta\mathbf{x}_{m})^{T}\nabla^{2}v^{+}(\mathbf{p})\Delta\mathbf{x}_{m}+\mathcal{O}(|\Delta\mathbf{x}_{m}|^{3}),\quad\text{
if }\mathbf{q}_{m}\in\Omega^{+},$ (73) $\displaystyle
v(\mathbf{q}_{m})+C(\mathbf{q}_{m})$
$\displaystyle=v^{+}(\mathbf{p})+(\Delta\mathbf{x}_{m})^{T}\nabla
v^{+}(\mathbf{p})+\dfrac{1}{2}(\Delta\mathbf{x}_{m})^{T}\nabla^{2}v^{+}(\mathbf{p})\Delta\mathbf{x}_{m}+\mathcal{O}(|\Delta\mathbf{x}_{m}|^{3}),\quad\text{
if }\mathbf{q}_{m}\in\Omega^{-},$ (74)

where $\Delta\mathbf{x}_{m}=\mathbf{q}_{m}-\mathbf{p},m=1,\cdots,M$ and
$C(\mathbf{p})$ is given by

$\displaystyle
C(\mathbf{q}_{m})=[v](\mathbf{p})+(\Delta\mathbf{x}_{m})^{T}[\nabla
v](\mathbf{p})+\dfrac{1}{2}(\Delta\mathbf{x}_{m})^{T}[\nabla^{2}v](\mathbf{p})\Delta\mathbf{x}_{m}.$
(75)

Dropping the $\mathcal{O}(|\Delta\mathbf{x}_{m}|^{3})$ term, one obtains an
$M\times M$ linear system to be solved for the boundary values $v^{+}$,
$\nabla v^{+}$, $\nabla^{2}v^{+}$. Then, the values $v^{-}$, $\nabla v^{-}$,
$\nabla^{2}v^{-}$ can also be obtained after simple algebraic manipulations.
To select interpolation points, for example, we first find the closest grid
node $(x_{i_{0}},y_{j_{0}})$ to $\mathbf{p}$. Let $d_{1},d_{2}$ be two
integers,

$d_{1}=\left\\{\begin{aligned} 1,\text{ if }p_{1}>x_{i_{0}},\\\ -1,\text{ if
}p_{1}\leq x_{i_{0}},\\\ \end{aligned}\right.\quad
d_{2}=\left\\{\begin{aligned} 1,\text{ if }p_{2}>y_{j_{0}},\\\ -1,\text{ if
}p_{2}\leq y_{j_{0}},\\\ \end{aligned}\right.$ (76)

Then six interpolation points are chosen

$\\{(x_{i_{0}+r},y_{j_{0}+s})\\}_{(r,s)\in\mathcal{I}},\quad\mathcal{I}=\\{(0,0),(1,0),(0,1),(0,-1),(-1,0),(d_{1},d_{2})\\}.$
(77)

This choice of interpolation points results in a $6\times 6$ system, which can
be solved with a direct method, such as the LU decomposition.

##  5 Interface evolution method

###  5.1 The $\theta-L$ approach

Suppose the interface $\Gamma$ is represented as a closed curve given by
$\mathbf{X}(\alpha,t)=(x(\alpha,t),y(\alpha,t))$, where $\alpha\in[0,2\pi)$
parameterizes the curve. The evolution of the curve is given by

$\displaystyle\mathbf{X}_{t}=U\mathbf{n}+V\mathbf{s},$ (78)

where $\mathbf{n}$ is the unit outward normal, $\mathbf{s}$ is the unit
tangent, $U$ and $V$ are the normal and tangential components of the curve
velocity, respectively. In most cases, only the shape of the moving interface
is of practical interest, which can be solely determined by the normal motion
since the tangential motion only changes the parameterization. However, pure
normal motion for tracking a moving interface is not a good choice for
numerical computation due to the clustering or spread of marker points on the
interface. Another issue is that the curvature term introduces numerical
stiffness into the interface evolution and causes strict stability constraints
on the time step. To mitigate these issues, we use the $\theta-L$ approach
instead of (78) for the curve evolution,

$\displaystyle L_{t}$
$\displaystyle=\int_{0}^{2\pi}\theta_{\alpha^{\prime}}U\,d\alpha^{\prime},$
(79) $\displaystyle\theta_{t}$
$\displaystyle=(\dfrac{2\pi}{L})(\theta_{\alpha}V-U_{\alpha}).$ (80)

where $\theta$ is the tangent angle to the curve $\Gamma$ and $L$ is the curve
length. By setting $s_{\alpha}=\sqrt{x_{\alpha}^{2}+y_{\alpha}^{2}}\equiv
L/2\pi$ for equal-arclength parameterization, one obtains the artificial
tangential velocity

$V(\alpha,t)=\dfrac{\alpha}{2\pi}\int_{0}^{2\pi}\theta_{\alpha^{\prime}}U\,d\alpha^{\prime}-\int_{0}^{\alpha}\theta_{\alpha^{\prime}}U\,d\alpha^{\prime},$
(81)

The mapping from $(\theta,L)$ to $\mathbf{X}=(x,y)$ still needs two more
integration constants to determine the position of the curve, for which we
track $\mathbf{\bar{X}}=\frac{1}{2\pi}\int_{0}^{2\pi}\mathbf{X}\,d\alpha$
using the evolution equation

$\displaystyle\mathbf{\bar{X}}_{t}=\dfrac{1}{2\pi}\int_{0}^{2\pi}U\mathbf{n}+V\mathbf{s}\,d\alpha.$
(82)

While it is possible to track a single point, as demonstrated in [14], we find
that using the averaged value provides better rotational symmetry.

###  5.2 Small scale decomposition

In the Hele-Shaw flow, the mapping from curvature $\kappa$ to velocity $U$ is
a one-phase Dirichlet-to-Neumann (DtN) mapping [61], denoted as
$\mathcal{T}_{1}:\sigma\kappa\rightarrow U$. Upon small scale decomposition,
it is found that the highest-order term in $U$, which contributes to the
stiffness, is linear. Specifically, we have

$\displaystyle
U\sim-\dfrac{\sigma}{2}\left(\dfrac{2\pi}{L}\right)^{2}\mathscr{H}[\theta_{\alpha\alpha}],$
(83)

where $\mathscr{H}$ represents the Hilbert transform [14]. This approach
offers a straightforward means to develop efficient semi-implicit schemes for
alleviating the stiffness associated with the curvature term.

For the Stefan problem, the mapping
$\mathcal{T}_{2}:\varepsilon_{C}\kappa\rightarrow U$ is a two-phase DtN
mapping. We shall also derive a small-scale decomposition for the normal
velocity $U$. For simplicity, we assume $\varepsilon_{V}$ and
$\varepsilon_{C}$ are constants. Due to the relation $U=\psi$, the DtN mapping
is defined through the boundary integral equation (32). Since the integral
equation (32) becomes a Fredholm integral equation of the first kind if
$\varepsilon_{V}$ = 0, we need to consider two cases: (1)
$\varepsilon_{V}=\mathcal{O}(1)$, and (2) $\varepsilon_{V}\ll 1$.

####  5.2.1 The first case: $\varepsilon_{V}=\mathcal{O}(1)$

Denote by
$K\psi=(1/\varepsilon_{V})\int_{\Gamma}G_{c}(\mathbf{y},\mathbf{x})\psi(\mathbf{y})\,d\mathbf{s}_{\mathbf{y}}$
and $g=-(\varepsilon_{C}/\varepsilon_{V})\kappa-T_{1}/\varepsilon_{V}$. We
write the boundary integral equation (32) as an operator equation

$(I-K)\psi=g,$ (84)

where the operator $K$ is compact and $I-K$ has a bounded inverse. Then, we
have

$\displaystyle\psi=g+K\psi=g+K[(I-K)^{-1}g].$ (85)

The integral operator $K$ is a pseudo-differential operator of order $-1$ and
defines the mapping $K:C^{\beta}(\Gamma)\rightarrow C^{1+\beta}(\Gamma)$ for a
fix constant $\beta$. Therefore, the term $S[(I-S)^{-1}g]$ is more regular
than $g$. Further, since $T_{1}$ is equivalent to a volume integral, which is
more regular than $\kappa=\theta_{\alpha}/s_{\alpha}=2\pi\theta_{\alpha}/L$,
we obtain

$U\sim-\dfrac{\varepsilon_{C}}{\varepsilon_{V}}\left(\dfrac{2\pi}{L}\right)\theta_{\alpha}.$
(86)

####  5.2.2 The second case: $\varepsilon_{V}\ll 1$

We rewrite $\psi$ as a formal asymptotic expansion
$\psi=\psi_{0}+\varepsilon_{V}\psi_{1}+\varepsilon_{V}^{2}\psi_{2}+\cdots$ and
substitute it into the boundary integral equation (32). The leading order
equation implies that $\psi_{0}$ should satisfy

$\displaystyle\int_{\Gamma}G_{c}(\mathbf{y},\mathbf{x})\psi_{0}(\mathbf{y})\,d\mathbf{s}_{y}=\varepsilon_{C}\kappa(\mathbf{x})+T_{1}(\mathbf{x}),\quad\text{for
}\mathbf{x}\in\Gamma.$ (87)

Since we only need the highest-order term in the small scale decomposition, we
can simply proceed by assuming $\varepsilon_{V}=0$ in (32) and use a similar
technique as that in [30] to obtain the small scale decomposition. Note that
the Green function $G_{c}(\mathbf{y},\mathbf{x})$ associated with the modified
Helmholtz operator $\Delta-c^{2}$ has the same singular behavior as the Green
function $G_{0}(\mathbf{y},\mathbf{x})$ associated with the Laplacian
$\Delta$. A simple singularity subtraction leads to

$\displaystyle\dfrac{1}{2\pi}\int_{\Gamma}\ln|\mathbf{y}-\mathbf{x}|\psi(\mathbf{y})\,d\mathbf{s}_{y}-\int_{\Gamma}(G_{0}(\mathbf{y},\mathbf{x})-G_{c}(\mathbf{y},\mathbf{x}))\psi(\mathbf{y})\,d\mathbf{s}_{y}=\varepsilon_{C}\kappa(\mathbf{x})+T_{1}(\mathbf{x}),\quad\text{for
}\mathbf{x}\in\Gamma,$ (88)

where the kernel $G_{0}(\mathbf{y},\mathbf{x})-G_{c}(\mathbf{y},\mathbf{x})$
no longer has a $log$-type singularity. Hence, in the left-hand side of (88),
the second term is more regular than the first one. For every $2\pi$-periodic
function $v$, define the operators $\mathscr{M}$ and $\mathscr{R}$ by

$\displaystyle\mathscr{M}[v](\alpha)$
$\displaystyle=\dfrac{1}{2\pi}\int_{0}^{2\pi}\ln\left|2\sin\dfrac{\alpha-\alpha^{\prime}}{2}\right|v(\alpha^{\prime})\,d\alpha^{\prime},$
(89) $\displaystyle\mathscr{R}[v](\alpha)$
$\displaystyle=\dfrac{1}{2\pi}\int_{0}^{2\pi}\ln\left|\dfrac{\mathbf{x}(\alpha)-\mathbf{x}(\alpha^{\prime})}{2\sin\dfrac{\alpha-\alpha^{\prime}}{2}}\right|v(\alpha^{\prime})\,d\alpha^{\prime}.$
(90)

Using the relation $s_{\alpha}\equiv L/2\pi$, the equation (88) can be written
as

$\dfrac{L}{2\pi}(\mathscr{M}[\psi]+\mathscr{R}[\psi])-\int_{\Gamma}(G_{0}(\mathbf{y},\mathbf{x})-G_{c}(\mathbf{y},\mathbf{x}))\psi(\mathbf{y})\,d\mathbf{s}_{y}=\varepsilon_{C}\dfrac{2\pi}{L}\theta_{\alpha}+T_{1}.$
(91)

Note that for every $2\pi$-periodic function $v$, the relation
$\mathscr{M}^{-1}[v]=c-2\mathscr{H}[v]$ holds, where $c$ is an arbitrary
constant. Taking $\mathscr{M}^{-1}$ on both sides of (91), it yields

$\displaystyle\psi-\bar{\psi}=-2\varepsilon_{C}\left(\dfrac{2\pi}{L}\right)^{2}\mathscr{H}[\theta_{\alpha\alpha}]-\dfrac{4\pi}{L}\mathscr{H}[\int_{\Gamma}(G_{0}(\mathbf{y},\mathbf{x})-G_{c}(\mathbf{y},\mathbf{x}))\psi(\mathbf{y})\,d\mathbf{s}_{y}+T_{1}-\dfrac{L}{2\pi}\mathscr{R}[\psi]]_{\alpha},$
(92)

where $\bar{\psi}$ is the mean of $\psi$. Finally, since the second term in
the right-hand side of (92) is more regular than the first term, we have

$\displaystyle
U\sim-2\varepsilon_{C}\left(\dfrac{2\pi}{L}\right)^{2}\mathscr{H}[\theta_{\alpha\alpha}],$
(93)

By inserting $U$ into the $\theta$ evolution equation (80) and using the small
scale expressions equations 86 and 93 , we can extract the linear and highest-
order terms

$\theta_{t}=\left\\{\begin{aligned}
&\dfrac{\varepsilon_{C}}{\varepsilon_{V}}\left(\dfrac{2\pi}{L}\right)^{2}\theta_{\alpha\alpha}+\mathscr{N}(\alpha),&\quad\text{if
}\varepsilon_{V}=\mathcal{O}(1),\\\
&2\varepsilon_{C}\left(\dfrac{2\pi}{L}\right)^{3}\mathscr{H}[\theta_{\alpha\alpha\alpha}]+\mathscr{N}(\alpha),&\quad\text{if
}\varepsilon_{V}\ll 1,\end{aligned}\right.$ (94)

where $\mathscr{N}$ consists of the remaining lower-order terms. The case
$\varepsilon_{V}=\mathcal{O}(1)$ is second-order diffusive and is similar to a
heat equation or a mean curvature flow. The stiffness of the second-order
derivative term can be removed by employing an implicit time-stepping scheme
for the stiff term. The second case $\varepsilon_{V}\ll 1$ is third-order
diffusive [14]. Implicit schemes for the case are more difficult due to the
nonlocal Hilbert transform $\mathscr{H}$. Using the fact that $\mathscr{H}$ is
diagonalizable under the Fourier transform, an accurate and efficient semi-
implicit scheme can be devised.

###  5.3 Semi-implicit scheme

The evolution equations (79) and (82) are not stiff, they can be discretized
with an explicit scheme. We use the second-order Adams-Bashforth scheme,

$\displaystyle L^{n+1}=L^{n}+\tau(3M^{n}-M^{n-1}),\quad
M=-\int_{0}^{2\pi}\theta_{\alpha^{\prime}}U\,d\alpha^{\prime},$ (95)
$\displaystyle\mathbf{\bar{X}}^{n+1}=\mathbf{\bar{X}}^{n}+\tau(3Q^{n}-Q^{n-1}),\quad
Q=\dfrac{1}{2\pi}\int_{0}^{2\pi}V_{n}\mathbf{n}+V_{s}\mathbf{s}\,d\alpha.$
(96)

The evolution equation (80) involves a stiff term, but it is hidden in the
normal velocity $U$ through a DtN mapping. In order to both remove the
stiffness and avoid solving nonlinear algebraic systems, only the highest-
order term in (80) is discretized implicitly. We first rewrite equation (80)
as

$\displaystyle\theta_{t}$
$\displaystyle=(\dfrac{2\pi}{L})(U_{\alpha}+\theta_{\alpha}T)+(\mathscr{L}(\alpha)-\mathscr{L}(\alpha))=\mathscr{L}(\alpha)+\mathscr{N}(\alpha),$
(97)

where $\mathscr{L}$ is the highest-order and linear term and $\mathscr{N}$
consists of the remaining lower-order terms. We have
$\mathscr{L}=\mathscr{L}_{1}=\lambda_{1}(2\pi/L)^{2}\theta_{\alpha\alpha}$ for
(86) and
$\mathscr{L}=\mathscr{L}_{2}=\lambda_{2}(2\pi/L)^{3}\mathscr{H}[\theta_{\alpha\alpha\alpha}]$
for (83) and (93), where $\lambda_{1}$ and $\lambda_{2}$ are two constant
parameters. For the Hele-Shaw flow and the Stefan problem with constant
$\varepsilon_{C}$ and $\varepsilon_{V}$, the parameters $\lambda_{1}$ and
$\lambda_{2}$ can be chosen to match the corresponding constant factors in the
highest-order terms. For the Stefan problem with anisotropic $\varepsilon_{C}$
and $\varepsilon_{V}$, $\lambda_{1}$ and $\lambda_{2}$ are regarded as
stabilization parameters, which should be chosen as sufficiently large to
ensure stability. With a frozen coefficient analysis, we choose the parameters
as

$\displaystyle\lambda_{1}=\max_{\alpha\in[0,2\pi)}\left|\dfrac{\varepsilon_{C}(\mathbf{n}(\alpha))}{\varepsilon_{V}(\mathbf{n}(\alpha))}\right|,\quad\lambda_{2}=\max_{\alpha\in[0,2\pi)}|2\varepsilon_{C}(\mathbf{n}(\alpha))|.$
(98)

In the Fourier space, the Hilbert transform $\mathscr{H}$ becomes diagonal,
and (97) simplifies to

$\displaystyle\hat{\theta}_{t}(k)$
$\displaystyle=-\lambda_{1}\left(\dfrac{2\pi}{L}\right)^{2}k^{2}\hat{\theta}(k)+\mathscr{\hat{N}}(k)\quad\text{for
}\mathscr{L}=\mathscr{L}_{1},$ (99) $\displaystyle\hat{\theta}_{t}(k)$
$\displaystyle=-\lambda_{2}\left(\dfrac{2\pi}{L}\right)^{3}|k|^{3}\hat{\theta}(k)+\mathscr{\hat{N}}(k)\quad\text{for
}\mathscr{L}=\mathscr{L}_{2}.$ (100)

A linear propagator and a second-order Adams-Bashforth method are used to
discretize the stiff part and the non-stiff part in (99) and (100),
respectively,

$\displaystyle\hat{\theta}^{n+1}(k)=e_{k}(t_{n},t_{n+1})\hat{\theta}^{n}(k)+\dfrac{\tau}{2}(3e_{k}(t_{n},t_{n+1})\mathscr{\hat{N}}^{n}(k)-e_{k}(t_{n-1},t_{n+1})\mathscr{\hat{N}}^{n-1}(k)),$
(101)

where the factors $e_{k}(t_{n},t_{n+1})$ and $e_{k}(t_{n-1},t_{n+1})$ are
specified as

$\displaystyle e_{k}(t_{n},t_{n+1})=\exp\left(-\frac{\lambda_{1}\tau}{2}(2\pi
k)^{2}\left[\frac{1}{\left(L^{n}\right)^{2}}+\frac{1}{\left(L^{n+1}\right)^{2}}\right]\right),$
(102) $\displaystyle e_{k}(t_{n-1},t_{n+1})=\exp\left(-\lambda_{1}\tau(2\pi
k)^{2}\left[\frac{1}{2\left(L^{n-1}\right)^{2}}+\frac{1}{\left(L^{n}\right)^{2}}+\frac{1}{2\left(L^{n+1}\right)^{2}}\right]\right),$
(103)

for the case $\mathscr{L}=\mathscr{L}_{1}$ and

$\displaystyle
e_{k}(t_{n},t_{n+1})=\exp\left(-\frac{\lambda_{2}\tau}{2}(2\pi|k|)^{3}\left[\frac{1}{\left(L^{n}\right)^{3}}+\frac{1}{\left(L^{n+1}\right)^{3}}\right]\right),$
(104) $\displaystyle
e_{k}(t_{n-1},t_{n+1})=\exp\left(-\lambda_{2}\tau(2\pi|k|)^{3}\left[\frac{1}{2\left(L^{n-1}\right)^{3}}+\frac{1}{\left(L^{n}\right)^{3}}+\frac{1}{2\left(L^{n+1}\right)^{3}}\right]\right).$
(105)

for the case $\mathscr{L}=\mathscr{L}_{2}$.

######  Remark 5.2.

The numerical operations for the evolution of $\theta$ and $L$, such as
differentiation, integration, and solving ODEs, are all performed in the
Fourier space using FFTs. This allows the interface evolution method to have
spectral accuracy, resulting in smaller errors compared to second-order time
integration methods. It is important to note that $\theta$ is periodic but has
an increment of $2\pi$ from $\alpha=0$ to $\alpha=2\pi$. To obtain accurate
results, one can use the auxiliary variable $\eta=\theta-\alpha$ when
performing Fourier transforms.

##  6 Numerical results

In this section, we demonstrate the application of the proposed method in
solving moving interface problems through numerical examples. In the first
part, we focus on the Hele-Shaw flow. We examine the convergence of the method
and utilize it to simulate the various dynamics of the Hele-Shaw bubble. By
using the proposed method, we are able to accurately capture the complex
behavior of the bubble, such as its growth and deformation. Furthermore, we
showcase the method’s ability to handle long-time computations by modeling a
bubble that exhibits intricate finger-like structures. Moving on to the second
part, we tackle the Stefan problem. We evaluate the convergence and stability
of the method using multiple examples. The proposed method proves to be robust
and accurate in solving the Stefan problem, allowing us to accurately track
the solidification interface and capture its evolution over time.
Additionally, we utilize the method to simulate dendritic solidification
problems, both with and without the presence of flow effects in the liquid
region. Through these examples, we demonstrate the versatility of the method
in handling different scenarios and accurately capturing the complex dynamics
of moving boundaries in solidification processes.

###  6.1 The Hele-Shaw flow

####  6.1.1 Convergence test

In this example, we investigate the convergence of the method for solving the
Hele-Shaw flow. The initial shape we consider is a four-fold flower, defined
as

$\displaystyle(x(\alpha,0),y(\alpha,0))=r(\alpha)(\cos\alpha,\sin\alpha),\quad
r(\alpha)=0.8+0.2\cos 4\alpha,\quad\alpha\in[0,2\pi).$ (106)

The surface tension coefficient is chosen to be $0.01$, and the numerical
computation is performed within the bounding box of $[-1.5,1.5]^{2}$. The air
injection rate is set to zero, leading to the area-preserving property of the
Hele-Shaw flow (since the fluid is incompressible). The numerical error is
measured by comparing the enclosed area of the curve at $T=1$ with the initial
area. To study the spatial accuracy, we use a fixed time step size of
$\tau=1\times 10^{-3}$ and vary the mesh size as $\Delta x=3\times 2^{-l}$,
with $l=5,6,\cdots,9$. Similarly, to analyze the temporal accuracy, we use a
fixed mesh size of $\Delta x=3\times 2^{-10}$ and vary the time step size as
$\tau=2^{-l}\times 10^{-3}$, with $l=1,2,\cdots,5$. We compute the enclosed
area of the curve by numerically evaluating the integral
$A[\Gamma]=\frac{1}{2}\int_{\Gamma}\mathbf{x}\cdot\mathbf{n}\,d\mathbf{s}$.
The numerical results are shown in Figure 2, and it can be observed that the
method exhibits second-order accuracy both spatially and temporally.

Figure 2: (a) shows the spatial accuracy of the method with different mesh
sizes $\Delta x=3\times 2^{-l},l=5,6,\cdots,9$. (b) shows the temporal
accuracy of the method with different time steps $\tau=2^{-l}\times
10^{-3},l=1,2,\cdots,5$.

####  6.1.2 Bubble relaxation

The Hele-Shaw flow without air injection presents interesting characteristics
due to the combined effects of surface tension and incompressibility. These
effects result in a curve-shortening and area-preserving behavior [62]. In
order to investigate this phenomenon, we consider an initial curve in the form
of a six-fold flower, described by the parametric equation:

$\displaystyle(x(\alpha,0),y(\alpha,0))=r(\alpha)(\cos\alpha,\sin\alpha),\quad
r(\alpha)=0.8+0.2\cos 6\alpha,\quad\alpha\in[0,2\pi).$ (107)

A surface tension coefficient of $0.01$ is used, and the evolution is computed
within a computational domain of $[-1.5,1.5]^{2}$. The numerical simulation is
conducted on a $512\times 512$ grid, with a time step size of $\tau=0.001$.
The evolution of the curve is illustrated in Figure 3. Additionally, the area
and length profiles of the curve, as well as the iteration numbers of the
GMRES method, are presented in Figure 4. As a result of the stabilizing
influence of surface tension, the initially irregular interface gradually
relaxes and approaches a circular shape over time. Notably, the enclosed area
remains constant throughout the evolution, while the length of the curve
decreases over time and eventually converges. These observations align with
the theoretical understanding of Hele-Shaw flow without air injection. The
GMRES iteration number remains relatively stable and decreases as the curve
approaches a circular shape.

Figure 3: Morphologies of the interface from $t=0$ to $t=1.2$ with a time
increment of $0.2$.

Figure 4: (a) shows the time evolutions of the enclosed area and the length of
the curve. (b) shows the iteration number of the GMRES method.

####  6.1.3 Unstable viscous fingering

In this example, the initial curve is given by a three-fold flower,

$\displaystyle(x(\alpha,0),y(\alpha,0))=r(\alpha)(\cos\alpha,\sin\alpha),\quad
r(\alpha)=0.8+0.2\cos 3\alpha,\quad\alpha\in[0,2\pi).$ (108)

The air injection rate is chosen as $J=1$ such that the air bubble
demonstrates an unstable growth. The surface tension coefficient varies from
$1\times 10^{-2}$ to $5\times 10^{-4}$. The computational domain is chosen as
$[-4,4]^{2}$. We solve the example on a $512\times 512$ grid with a time step
$\tau=0.005$. The morphologies of the evolutionary curve are shown in Figure 5
in time increments of $0.2$ to the final time $T=3$. The competition between
the stabilizing effect of the surface tension and the destabilizing effect of
the driven force due to the injection leads to the viscous fingering feature
of the Hele-Shaw flow. With small surface tension, which has a less
stabilizing effect, the growing bubble develops more branches. It also can be
observed that for small surface tensions, the symmetry of the interface is
difficult to capture since grid-induced anisotropy dominates the interface
evolution.

Figure 5: Detailed morphologies of the curve with different surface tension
coefficients: (a) $1\times 10^{-2}$;(b) $5\times 10^{-3}$;(c) $1\times
10^{-3}$;(d) $5\times 10^{-4}$.

####  6.1.4 Long-time computation

In this example, we perform long-time computation with the present method and
spatial-temporal rescaling scheme [15] for simulating a large Hele-Shaw
bubble. The computation is performed in the scaled frame and then mapped back
onto the nonscaled frame. The initial shape is chosen as a nucleus

$\displaystyle(\bar{x}(\alpha,0),\bar{y}(\alpha,0))=\bar{r}(\alpha)(\cos\alpha,\sin\alpha),\quad\bar{r}(\alpha)=1.0+0.1(\sin
2\alpha+\cos 3\alpha),\quad\alpha\in[0,2\pi).$ (109)

The computational domain is chosen as $[-1.7,1.7]^{2}$. We use a $1024\times
1024$ grid and a time step of $\Delta\bar{t}=2\times 10^{-4}$. The interface
points are adaptively refined with the criteria
$\Delta\bar{s}>1.5\Delta\bar{x}$. The injection rate and the surface tension
coefficient are set as $\bar{J}=1$ and $\sigma=0.001$, respectively. The
morphologies of the curve in the nonscaled frame are shown in (scaled) time
increments of $0.2$, in Figure 6. The computation takes $5$ hours to reach the
final scaled time $\bar{T}=3$ (nonscaled time $T=203$). At the final time, the
number of the marker point on the interface is refined to $16384$, and the
enclosed area and the length of the interface are $A=\bar{R}^{2}\bar{A}=1280$
and $L=\bar{R}\bar{L}/2\pi=169$.

Figure 6: Numerical results of long-time computation of the Hele-Shaw flow:
(a) morphology histories of the curve; (b) and (c) time evolution of the curve
length and enclosed area; (d) iteration numbers of the GMRES method.

###  6.2 The Stefan problem

####  6.2.1 Grid refinement analysis

In this example, we consider a benchmark grid refinement test for the Stefan
problem. In the beginning, a solid seed is put into an undercooling liquid
domain $\mathcal{B}=[-2,2]^{2}$. The initial shape is a four-fold flower,
which is given by

$(x(\alpha,0),y(\alpha,0))=r(\alpha)(\cos\alpha,\sin\alpha),\quad
r(\alpha)=0.1+0.02\cos 4\alpha,\quad\alpha\in[0,2\pi).$ (110)

The initial undercooling is chosen as $St=-0.5$. Isotropic surface tension and
kinematic coefficients are chosen as $\varepsilon_{C}=\varepsilon_{V}=2\times
10^{-3}$. We take the time step as $\tau=0.001$ and successively refine the
grid from $64$ to $512$. The morphologies of the liquid-solid interface in
time increments of $0.05$ to the final time $T=0.8$ are shown in Figure 7. By
comparing the results with previous works by the level-set method and the
front-tracking method [16, 20], it can be observed that our method has less
grid-induced surface tension, which leads to better convergence. The method is
also more able to preserve the symmetries of the interface, even with a coarse
grid.

Figure 7: Grid refinement analysis of the Stefan problem with different grids:
(a) $64\times 64$; (b) $128\times 128$; (c) $256\times 256$; (d) $512\times
512$.

####  6.2.2 Stability test

In order to demonstrate the stability of the method, we compare the semi-
implicit scheme with an explicit Adams-Bashforth scheme. We consider the case
that $\varepsilon_{V}=0$ and $\varepsilon_{C}=0.05$, which results in a third-
order stiffness. The initial shape is chosen as a slightly perturbed circle

$(x(\alpha,0),y(\alpha,0))=r(\alpha)(\cos\alpha,\sin\alpha),\quad
r(\alpha)=1+0.02\cos 4\alpha,\quad\alpha\in[0,2\pi).$ (111)

We use a $256\times 256$ grid for the computational domain
$\mathcal{B}=[-2,2]^{2}$. The interface is discretized with $128$ points.
First, the problem is solved with an explicit scheme with time-steps
$\tau=5\times 10^{-5}$ and $\tau=2.5\times 10^{-5}$. The curves at $t=0.1$ are
plotted in Figure 8. The smaller time-step $\tau=2.5\times 10^{-5}$ is stable
for the computation, while the larger one $\tau=5\times 10^{-5}$ is unstable,
and the solution quickly blows up. Then, the problem is solved with the semi-
implicit scheme with a time-step $\tau=0.01$, and the solution at $t=0$ is
plotted in Figure 8. Compared with the explicit scheme, the semi-implicit
scheme is stable and can obtain a correct result with a much larger time step,
which makes the semi-implicit scheme more efficient.

Figure 8: Numerical results of the stability test for the Stefan problem with
different time-stepping methods: (a) the Adams-Bashforth scheme; (b) the semi-
implicit scheme.

####  6.2.3 Comparison with the solvability theory

In a dendritic growth problem, the growth rate of a dendrite can be predicted
by the solvability theory [9, 63]. We compare the numerical result with the
solvability result to assess the accuracy of the method. In this example, the
initial seed is chosen as a circle with a radius of $0.1$. The parameters in
the Gibbs-Thomson relation are chosen as $\varepsilon_{V}=0$,
$\varepsilon_{C}(\alpha)=0.001[1+0.4(1-\cos(4\alpha))]$ where $\alpha$ is the
angle between the normal to the interface and the $x$-axis. The computation
domain is chosen as $[-6,6]$. We set $\tau=0.001$ for the computation. The
morphologies of the liquid-solid interface and the time-evolution of the tip
velocity from $t=0$ to $t=2.2$ are plotted in Figure 9. The tip velocity
finally converges and agrees with the solvability result of $1.7$.

Figure 9: Numerical results of the dendritic growth problem and comparison
with the solvability theory: (a) morphology histories of the interface; (b)
time evolution of the tip velocity.

####  6.2.4 Anisotropic dendritic growth

initial seed is a circle with a radius $0.05$ at the center of the
computational domain $[-4,4]^{2}$. undercooling $St=-0.65$ four-fold
anisotropy $\varepsilon_{V}=0.002$,
$\varepsilon_{C}(\alpha)=0.002(8/3\sin^{4}(2(\alpha-\alpha_{0})))$ grid
$512\times 512$ with a time step $\tau=0.0002$.

(a) Numerical results with $\alpha_{0}=0$.

(b) Numerical results with $\alpha_{0}=\pi/4$.

Figure 10: Interface morphologies and the temperature field of the dendritic
growth problem with four-fold anisotropy. Snapshots are taken at
$t=0,0.02,0.06$, and $0.1$.

The initial seed is a circle with a radius $0.05$ at the center of the
computational domain $\mathcal{B}=[-2,2]^{2}$. The undercooling parameter
varies from $St=-0.55$ to $St=-0.65$. The parameters in the Gibbs-Thomson
relation are chosen as $\varepsilon_{V}=0.002$,
$\varepsilon_{C}(\alpha)=0.002(8/3\sin^{4}(3\alpha))$, which is a six-fold
anisotropy. The six-fold anisotropy leads to snowflake-shaped dendritic
growths, which are shown in Figure 11.

(a) Numerical results with $St=-0.55$.

(b) Numerical results with $St=-0.65$.

Figure 11: Interface morphologies and the temperature field of the dendritic
growth problem with six-fold anisotropy. Snapshots are taken at
$t=0,0.02,0.06$, and $0.1$.

####  6.2.5 Dendritic growth with external flow

We consider the convection effect in the dendritic growth problem. In this
example, the Stefan problem with natural convection is solved for simulations.
Initially, the seed is a circle with a radius $0.05$ at the center of the
computational domain $[-2,2]^{2}$ surrounded by undercooled fluid with
temperature $St=-0.5$. Four-fold anisotropic surface tension
$\varepsilon_{C}(\alpha)=0.002(8/3\sin^{4}(2\alpha))$ and local kinematic
equilibrium $\varepsilon_{V}=0$ are applied. The computation is performed with
a $512\times 512$ grid and a time step $\tau=0.0002$. Inflow and outflow
boundary conditions, $\mathbf{u}=(u_{0},0)^{T}$, are applied on the left and
right boundaries, respectively. No-slip boundary condition,
$\mathbf{u}=\mathbf{0}$, is used on the top and bottom boundaries. In this
example, we ignore the buoyancy force, i.e. $\beta=0$. As a result, the
convection of the flow is only driven by boundary conditions. When $u_{0}=0$,
the problem is identical to the classical Stefan problem without natural
convection. In order to examine the convection effect on the growth pattern,
we compare numerical results obtained with different flow velocities in the
fluid phase. In Figure 12, profiles of dendritic growth, temperature, and flow
fields at $t=0.1$ are presented. Evolution histories of $x$-components of the
left and right tips are shown in Figure 13. It can be observed that convection
of the flow leads to faster growth of the left branch and slower growth of the
right branch. This effect is more evident as the flow velocity increases. Due
to the flow convection, the released latent heat is carried by the flow from
left to right and consequently leads to non-symmetric temperature distribution
in the $x$-direction.

Figure 12: Details of dendritic growth histories, temperature fields, and flow
fields with different inflow velocities at $t=0.1$ ($u_{0}=0,2,4$ and $8$ from
left to right).

Figure 13: Evolution histories of $x$-components of the left and right tips
with different flow velocities.

####  6.2.6 Dendritic growth with buoyancy-driven flow

In the final example, we consider the dendritic growth problem with buoyancy-
driven flow. The anisotropic surface tension is chosen as a rotated one,
$\varepsilon_{C}(\alpha)=0.002(8/3\sin^{4}(2(\alpha-\pi/4)))$. The no-slip
boundary condition is applied for the fluid equation on the four boundaries.
Different thermal expansion coefficients are chosen so that the flow in the
fluid phase is driven by buoyancy force. The gravity acceleration is chosen as
$g=10$. The reference temperature is chosen as the temperature of surrounding
fluid $T_{0}=-0.5$. Other parameters are chosen as the same as those in the
last example. Numerical results with increasing thermal expansion coefficients
are presented in Figure 14. Near the solid-liquid interface, the released
latent heat increases the fluid temperature and, as a result, causes fluid
density changes and the buoyancy force. Driven by the buoyancy force, the
fluid carries the heat and flows from bottom to top, which leads to non-
symmetric temperature distribution in $y$-direction. The two upper branches
are restrained from growing due to accumulated heat, while the two lower
branches grow much faster since the heat flows away.

Figure 14: Details of dendritic growth histories, temperature fields, and flow
fields with different thermal expansion coefficients at $t=0.1$
($\beta=0,10^{3},2\times 10^{3}$ and $4\times 10^{3}$ from left to right).

##  7 Discussion

This study presents a novel numerical method for solving two representative
moving interface problems using a Cartesian grid-based boundary integral
method combined with an interface evolution approach. The proposed method
offers several advantages that combine the strengths of Cartesian grid-based
solvers and boundary integral methods.

Elliptic and parabolic PDEs with irregular boundaries or interfaces are
reformulated as well-conditioned boundary integral equations and are solved
using the KFBI method. The KFBI method utilizes a Cartesian grid-based solver
for integral evaluations, eliminating the need for analytical expressions of
Green’s functions and allowing for efficient utilization of fast PDE solvers
like FFTs and geometric multigrid methods. It is important to note that the
method is not entirely kernel-free in this work. Green’s function is still
required to impose artificial boundary conditions on the rectangular domain
boundary when dealing with infinite boundary conditions. However, the KFBI
method avoids the evaluation of singular and nearly singular integrals, which
are typically challenging in quadrature-based boundary integral methods.

The interface evolution is solved using the $\theta-L$ formulation instead of
the more commonly used $x-y$ formulation. Due to the periodicity of the
interface, an accurate Fourier pseudo-spectral method is employed for spatial
approximation in the interface evolution problem. Furthermore, the use of the
$\theta-L$ formulation enables the application of a small-scale decomposition
approach to address the curvature-induced stiffness in the Hele-Shaw flow and
Stefan problem. This approach, combined with a semi-implicit scheme, allows
for efficient and stable evolution of the interface by solving the semi-
implicit scheme in the Fourier space with the FFT.

While the current work primarily focuses on moving interface problems in two
dimensions, certain improvements are needed to solve models in three
dimensions, such as solidification problems and two-phase incompressible
flows. First, elliptic PDEs with irregular boundaries and interfaces in three
dimensions require a three-dimensional version of the KFBI method [52].
Additionally, for parabolic PDEs, dimension-splitting techniques can be
employed to accelerate computation [64]. Second, the $\theta-L$ approach is
only applicable to evolutionary curves in two dimensions. To accurately track
evolutionary surfaces in three dimensions, different numerical approaches,
such as the front-tracking method or the level-set method, are needed.
Finally, addressing the stiffness induced by mean curvature in three
dimensions poses greater challenges than in two dimensions. Developing an
efficient semi-implicit time-stepping scheme is crucial for ensuring the
stable evolution of the interface.

##  Appendix A Computation of derivative jumps

Suppose $\Gamma$ is parameterized as
$\mathbf{X}(\alpha)=(x(\alpha),y(\alpha))$ where $\alpha$ is an arbitrary
parameter. Let $s$ be the arc-length parameter. Suppose $\Gamma$ is
sufficiently smooth, at least in the class $C^{2}$. The unit outward normal is
given by $\mathbf{n}=(y_{\alpha}/s_{\alpha},-x_{\alpha}/s_{\alpha})$ where
$s_{\alpha}=\sqrt{x_{\alpha}^{2}+y_{\alpha}^{2}}$. Given an interface problem
with a constant coefficient, we shall derive jump values of the solution and
its derivatives at the point $(x(\alpha),y(\alpha))\in\Gamma$.

###  A.1 The modified Helmholtz equation

Consider the interface problem of the modified Helmholtz equation

$\displaystyle\Delta u-c^{2}u$ $\displaystyle=f,\quad\text{ in
}\Omega^{+}\cup\Omega^{-},$ (112) $\displaystyle[u]$
$\displaystyle=\Phi,\quad\text{ on }\Gamma,$ (113)
$\displaystyle[\partial_{\mathbf{n}}u]$ $\displaystyle=\Psi,\quad\text{ on
}\Gamma.$ (114)

The interface condition (113) implies the zeroth-order jump value

$[u]=\Phi.$ (115)

By taking the derivative of both sides of the interface condition (113) with
respect to $\alpha$ and combining the interface condition (114), we have a
$2\times 2$ linear system

$\displaystyle x_{\alpha}[u_{x}]+y_{\alpha}[u_{y}]=\Phi_{\alpha},$ (116)
$\displaystyle y_{\alpha}[u_{x}]-x_{\alpha}[u_{y}]=s_{\alpha}\Psi.$ (117)

Solving the linear system gives the values of $[u_{x}]$ and $[u_{y}]$. By
taking the derivative of both sides of equations 116 and 117, and using
equation 112, we have a $3\times 3$ linear system for $[u_{xx}]$, $[u_{yy}]$
and $[u_{xy}]$

$\displaystyle(x_{\alpha})^{2}[u_{xx}]+(y_{\alpha})^{2}[u_{yy}]+2x_{\alpha}y_{\alpha}[u_{xy}]=\Phi_{\alpha\alpha}-x_{\alpha\alpha}[u_{x}]-y_{\alpha\alpha}[u_{y}],$
(118) $\displaystyle
x_{\alpha}y_{\alpha}[u_{xx}]-x_{\alpha}y_{\alpha}[u_{yy}]+((y_{\alpha})^{2}-(x_{\alpha})^{2})[u_{xy}]$
$\displaystyle=s_{\alpha\alpha}\Psi+s_{\alpha}\Psi_{\alpha}-y_{\alpha\alpha}[u_{x}]+x_{\alpha\alpha}[u_{y}],$
(119) $\displaystyle[u_{xx}]+[u_{yy}]=c^{2}\Phi+[f].$ (120)

After solving the linear system, the derivative jump values are obtained.

###  A.2 The modified Stokes equation

Consider the interface problem of the modified Stokes equation

$\displaystyle\Delta\mathbf{u}-c^{2}\mathbf{u}-\nabla p$
$\displaystyle=\mathbf{f},\quad\text{ in }\Omega^{+}\cup\Omega^{-},$ (121)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,\quad\text{ in
}\Omega^{+}\cup\Omega^{-},$ (122) $\displaystyle[\mathbf{u}]$
$\displaystyle=\mathbf{\Phi},\quad\text{ on }\Gamma,$ (123)
$\displaystyle[T(\mathbf{u},p)]$ $\displaystyle=\mathbf{\Psi},\quad\text{ on
}\Gamma.$ (124)

where $\mathbf{u}=(u,v)^{T}$,
$T(\mathbf{u},p)=-p\mathbf{n}+(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})\mathbf{n}$,
$\mathbf{f}=(f_{1},f_{2})^{T}$, $\mathbf{\Phi}=(\Phi_{1},\Phi_{2})^{T}$,
$\mathbf{\Psi}=(\Psi_{1},\Psi_{2})^{T}$. The zeroth-order jump values are
obtained from equation 123,

$[u]=\Phi_{1},\quad[v]=\Phi_{2}.$ (125)

Taking the derivative of both sides of the interface condition (123) with
respect to $\alpha$ and using equations 122 and 124, we obtain a $5\times 5$
system for the first-order jump values $[u_{x}]$, $[u_{y}]$, $[v_{x}]$,
$[v_{y}]$ and $[p]$,

$\displaystyle x_{\alpha}[u_{x}]+y_{\alpha}[u_{y}]=\Phi_{1,\alpha},$ (126)
$\displaystyle x_{\alpha}[v_{x}]+y_{\alpha}[v_{y}]=\Phi_{2,\alpha},$ (127)
$\displaystyle
2y_{\alpha}[u_{x}]-x_{\alpha}[u_{y}]-x_{\alpha}[v_{x}]-y_{\alpha}[p]=s_{\alpha}\Psi_{1},$
(128) $\displaystyle
y_{\alpha}[u_{y}]+y_{\alpha}[v_{x}]-2x_{\alpha}[v_{y}]+x_{\alpha}[p]=s_{\alpha}\Psi_{2},$
(129) $\displaystyle[u_{x}]+[v_{y}]=0.$ (130)

Taking the derivative of both sides of equations 126, 127, 128 and 129 with
respect to $\alpha$ and the derivatives of both sides of equation 122 with
respect to $x$ and $y$, and using equation 121, we obtain an $8\times 8$
system for the second-order jump values $[u_{xx}]$, $[u_{yy}]$, $[u_{xy}]$,
$[v_{xx}]$, $[v_{yy}]$, $[v_{xy}]$, $[p_{x}]$ and $[p_{y}]$,

$\displaystyle(x_{\alpha})^{2}[u_{xx}]+(y_{\alpha})^{2}[u_{yy}]+2x_{\alpha}y_{\alpha}[u_{xy}]=r_{1},$
(131)
$\displaystyle(x_{\alpha})^{2}[v_{xx}]+(y_{\alpha})^{2}[v_{yy}]+2x_{\alpha}y_{\alpha}[v_{xy}]=r_{2},$
(132) $\displaystyle
2x_{\alpha}y_{\alpha}[u_{xx}]-x_{\alpha}y_{\alpha}[u_{yy}]+(2(y_{\alpha})^{2}-(x_{\alpha})^{2})[u_{xy}]$
$\displaystyle-(x_{\alpha})^{2}[v_{xx}]-x_{\alpha}y_{\alpha}[v_{xy}]-x_{\alpha}y_{\alpha}[p_{x}]-(y_{\alpha})^{2}[p_{y}]=r_{3},$
(133)
$\displaystyle(y_{\alpha})^{2}[u_{yy}]+x_{\alpha}y_{\alpha}[u_{xy}]+x_{\alpha}y_{\alpha}[v_{xx}]-2x_{\alpha}y_{\alpha}[v_{yy}]$
$\displaystyle+((y_{\alpha})^{2}-2(x_{\alpha})^{2})[v_{xy}]+(x_{\alpha})^{2}[p_{x}]+x_{\alpha}y_{\alpha}[p_{y}]=r_{4},$
(134) $\displaystyle[u_{xx}]+[u_{yy}]-[p_{x}]=c^{2}\Phi_{1}+[f_{1}],$ (135)
$\displaystyle[v_{xx}]+[v_{yy}]-[p_{y}]=c^{2}\Phi_{2}+[f_{2}],$ (136)
$\displaystyle[u_{xx}]+[v_{xy}]=0,$ (137) $\displaystyle[u_{xy}]+[v_{yy}]=0.$
(138)

where $r_{i},i=1,2,\cdots,4$ are given by

$\displaystyle
r_{1}=\Phi_{1,\alpha\alpha}-x_{\alpha\alpha}[u_{x}]-y_{\alpha\alpha}[u_{y}],$
(139) $\displaystyle
r_{2}=\Phi_{2,\alpha\alpha}-x_{\alpha\alpha}[v_{x}]-y_{\alpha\alpha}[v_{y}],$
(140) $\displaystyle
r_{3}=s_{\alpha\alpha}\Psi_{1}+s_{\alpha}\Psi_{1,\alpha}-2y_{\alpha\alpha}[u_{x}]+x_{\alpha\alpha}[u_{y}]+x_{\alpha\alpha}[v_{x}]+y_{\alpha\alpha}[p],$
(141) $\displaystyle
r_{4}=s_{\alpha\alpha}\Psi_{2}+s_{\alpha}\Psi_{2,\alpha}-y_{\alpha\alpha}[u_{y}]-y_{\alpha\alpha}[v_{x}]+2x_{\alpha\alpha}[v_{y}]-x_{\alpha\alpha}[p].$
(142)

By solving the three linear systems, the derivative jump values of $u,v,p$ can
be obtained.

## References

  * [1] H.-D. Cao, B.-L. Chen, X.-P. Zhu, Recent developments on the Hamilton’s Ricci Flow, Surveys in Differential Geometry 12 (1) (2007) 47–112.  doi:10.4310/sdg.2007.v12.n1.a3.
  * [2] X. Chen, Generation and propagation of interfaces for reaction-diffusion equations, Journal of Differential Equations 96 (1) (1992) 116–141.  doi:10.1016/0022-0396(92)90146-E.
URL https://linkinghub.elsevier.com/retrieve/pii/002203969290146E

  * [3] G. BELLETTINI, M. PAOLINI, Anisotropic motion by mean curvature in the context of Finsler geometry, Hokkaido Mathematical Journal 25 (3) (1996) 537–566.  doi:10.14492/hokmj/1351516749.
URL https://projecteuclid.org/journals/hokkaido-mathematical-
journal/volume-25/issue-3/Anisotropic-motion-by-mean-curvature-in-the-context-
of-Finsler/10.14492/hokmj/1351516749.full

  * [4] J. M. Hyman, Numerical methods for tracking interfaces, Physica D: Nonlinear Phenomena 12 (1) (1984) 396–407.  doi:https://doi.org/10.1016/0167-2789(84)90544-X.
URL https://www.sciencedirect.com/science/article/pii/016727898490544X

  * [5] J. Glimm, D. Marchesin, O. McBryan, Statistical fluid dynamics: Unstable fingers, Communications in Mathematical Physics 74 (1) (1980) 1–13.  doi:10.1007/BF01197574.
URL http://link.springer.com/10.1007/BF01197574

  * [6] J. Glimm, D. Marchesin, O. McBryan, Subgrid resolution of fluid discontinuities, II, Journal of Computational Physics 37 (3) (1980) 336–354.  doi:10.1016/0021-9991(80)90041-8.
URL https://linkinghub.elsevier.com/retrieve/pii/0021999180900418

  * [7] W. W. Mullins, Two-dimensional motion of idealized grain boundaries, Journal of Applied Physics 27 (8) (1956) 900–904.  doi:10.1063/1.1722511.
URL http://aip.scitation.org/doi/10.1063/1.1722511

  * [8] J. A. Sethian, J. Straint, Crystal growth and dendritic solidification, Journal of Computational Physics 98 (2) (1992) 231–253.  doi:10.1016/0021-9991(92)90140-T.
  * [9] D. I. Meiron, Selection of steady states in the two-dimensional symmetric model of dendritic growth, Tech. Rep. 4 (1986).  doi:10.1103/PhysRevA.33.2704.
  * [10] W. W. Mullins, R. F. Sekerka, Morphological Stability of a Particle Growing by Diffusion or Heat Flow, Journal of Applied Physics 34 (2) (1963) 323–329.  doi:10.1063/1.1702607.
URL https://pubs.aip.org/aip/jap/article/34/2/323-329/163811

  * [11] Y. Li, D. Jeong, J.-i. Choi, S. Lee, J. Kim, Fast local image inpainting based on the Allen–Cahn model, Digital Signal Processing 37 (1) (2015) 65–74.  doi:10.1016/j.dsp.2014.11.006.
URL
http://dx.doi.org/10.1016/j.dsp.2014.11.006https://linkinghub.elsevier.com/retrieve/pii/S1051200414003418

  * [12] M. Beneš, V. Chalupecký, K. Mikula, Geometrical image segmentation by the Allen-Cahn equation, Applied Numerical Mathematics 51 (2-3) (2004) 187–205.  doi:10.1016/j.apnum.2004.05.001.
  * [13] A. DEGREGORIA, L. SCHWARTZ, A boundary-integral method for two-phase displacement in hele-shaw cells, Dynamics of Curved Fronts (1988) 201–218doi:10.1016/B978-0-08-092523-3.50022-8.
  * [14] T. Y. Hou, J. S. Lowengrub, M. J. Shelley, Removing the stiffness from interfacial flows with surface tension, Journal of Computational Physics 114 (2) (1994) 312–338.  doi:10.1006/jcph.1994.1170.
  * [15] S. Li, J. S. Lowengrub, P. H. Leo, A rescaling scheme with application to the long-time simulation of viscous fingering in a Hele-Shaw cell, Journal of Computational Physics 225 (1) (2007) 554–567.  doi:10.1016/j.jcp.2006.12.023.
  * [16] D. Juric, G. Tryggvason, A front-tracking method for dendritic solidification, Journal of Computational Physics 123 (1) (1996) 127–148.  doi:10.1006/jcph.1996.0011.
  * [17] J. Strain, A boundary integral approach to unstable solidification, Journal of Computational Physics 85 (2) (1989) 342–389.  doi:10.1016/0021-9991(89)90155-1.
URL
https://www.sciencedirect.com/science/article/pii/0021999189901551https://linkinghub.elsevier.com/retrieve/pii/0021999189901551

  * [18] P. Zhao, J. C. Heinrich, D. R. Poirier, Fixed mesh front-tracking methodology for finite element simulations, International Journal for Numerical Methods in Engineering 61 (6) (2004) 928–948.  doi:https://doi.org/10.1002/nme.1098.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1098

  * [19] A. Schmidt, Computation of Three Dimensional Dendrites with Finite Elements, Journal of Computational Physics 125 (2) (1996) 293–312.  doi:10.1006/jcph.1996.0095.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999196900959

  * [20] S. Chen, B. Merriman, S. Osher, P. Smereka, A simple level set method for solving stefan problems, Tech. Rep. 1 (1997).  doi:10.1006/jcph.1997.5721.
  * [21] Y.-T. Kim, N. Goldenfeld, J. Dantzig, Computation of dendritic microstructures using a level set method, Physical Review E 62 (2) (2000) 2471–2474.  doi:10.1103/PhysRevE.62.2471.
URL https://link.aps.org/doi/10.1103/PhysRevE.62.2471

  * [22] F. Gibou, R. Fedkiw, R. Caflisch, S. Osher, A Level Set Approach for the Numerical Simulation of Dendritic Growth, Tech. Rep. 1-3 (2003).  doi:10.1023/A:1025399807998.
  * [23] L. Boledi, B. Terschanski, S. Elgeti, J. Kowalski, A level-set based space-time finite element approach to the modelling of solidification and melting processes, Journal of Computational Physics 457 (2022) 111047.  arXiv:2105.09286, doi:10.1016/j.jcp.2022.111047.
URL https://doi.org/10.1016/j.jcp.2022.111047

  * [24] A. Limare, S. Popinet, C. Josserand, Z. Xue, A. Ghigo, A hybrid level-set / embedded boundary method applied to solidification-melt problems, Journal of Computational Physics 474 (2023) 111829.  arXiv:2202.08300, doi:10.1016/j.jcp.2022.111829.
URL
http://arxiv.org/abs/2202.08300https://linkinghub.elsevier.com/retrieve/pii/S0021999122008920

  * [25] R. Scardovelli, S. Zaleski, DIRECT NUMERICAL SIMULATION OF FREE-SURFACE AND INTERFACIAL FLOW, Annual Review of Fluid Mechanics 31 (1) (1999) 567–603.  doi:10.1146/annurev.fluid.31.1.567.
URL https://www.annualreviews.org/doi/10.1146/annurev.fluid.31.1.567

  * [26] C. Hirt, B. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1) (1981) 201–225.  doi:10.1016/0021-9991(81)90145-5.
URL
https://www.sciencedirect.com/science/article/pii/0021999181901455https://linkinghub.elsevier.com/retrieve/pii/0021999181901455

  * [27] A. Karma, W. J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 57 (4) (1998) 4323–4349.  doi:10.1103/PhysRevE.57.4323.
  * [28] H. Wang, R. Li, T. Tang, Efficient computation of dendritic growth with r-adaptive finite element methods, Journal of Computational Physics 227 (12) (2008) 5984–6000.  doi:10.1016/j.jcp.2008.02.016.
  * [29] X. Hu, R. Li, T. Tang, A multi-mesh adaptive finite element approximation to phase field models, Communications in Computational Physics 5 (5) (2009) 1012–1029.
  * [30] J. Zhu, X. Chen, T. Y. Hou, An efficient boundary integral method for the Mullins-Sekerka problem, Journal of Computational Physics 127 (2) (1996) 246–267.  doi:10.1006/jcph.1996.0173.
  * [31] V. Cristini, J. Lowengrub, Three-dimensional crystal growth—I: linear analysis and self-similar evolution, Journal of Crystal Growth 240 (1-2) (2002) 267–276.  doi:10.1016/S0022-0248(02)00831-X.
URL https://linkinghub.elsevier.com/retrieve/pii/S002202480200831X

  * [32] V. Cristini, J. Lowengrub, Three-dimensional crystal growth - II: Nonlinear simulation and control of the Mullins-Sekerka instability, Journal of Crystal Growth 266 (4) (2004) 552–567.  doi:10.1016/j.jcrysgro.2004.02.115.
  * [33] C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics 25 (3) (1977) 220–252.  doi:10.1016/0021-9991(77)90100-0.
URL https://linkinghub.elsevier.com/retrieve/pii/0021999177901000

  * [34] C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 479–517.  doi:10.1017/S0962492902000077.
URL
https://www.cambridge.org/core/product/identifier/S0962492902000077/type/journal_article

  * [35] K. Taira, T. Colonius, The immersed boundary method: A projection approach, Journal of Computational Physics 225 (2) (2007) 2118–2137.  doi:10.1016/j.jcp.2007.03.005.
  * [36] R. J. Leveque, Z. Li, Immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM Journal on Numerical Analysis 31 (4) (1994) 1019–1044.  doi:10.1137/0731054.
URL https://doi.org/10.1137/0731054

  * [37] Z. Li, Immersed interface methods for moving interface problems, Numerical Algorithms 14 (4) (1997) 269–293.  doi:10.1023/A:1019173215885.
URL https://doi.org/10.1023/A:1019173215885

  * [38] R. J. LeVeque, Z. Li, Immersed Interface Methods for Stokes Flow with Elastic Boundaries or Surface Tension, SIAM Journal on Scientific Computing 18 (3) (1997) 709–735.  doi:10.1137/S1064827595282532.
URL http://epubs.siam.org/doi/10.1137/S1064827595282532

  * [39] Z. Li, M. C. Lai, The immersed interface method for the navier-stokes equations with singular forces, Journal of Computational Physics 171 (2) (2001) 822–842.  doi:10.1006/jcph.2001.6813.
  * [40] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method), Journal of Computational Physics 152 (2) (1999) 457–492.  doi:10.1006/jcph.1999.6236.
  * [41] R. P. Fedkiw, T. Aslam, S. Xu, The Ghost Fluid Method for Deflagration and Detonation Discontinuities, Journal of Computational Physics 154 (2) (1999) 393–427.  doi:10.1006/jcph.1999.6320.
  * [42] X. D. Liu, R. P. Fedkiw, M. Kang, A Boundary Condition Capturing Method for Poisson’s Equation on Irregular Domains, Journal of Computational Physics 160 (1) (2000) 151–178.  doi:10.1006/jcph.2000.6444.
  * [43] D. Q. Nguyen, R. P. Fedkiw, M. Kang, A boundary condition capturing method for incompressible flame discontinuities, Journal of Computational Physics 172 (1) (2001) 71–98.  doi:10.1006/jcph.2001.6812.
  * [44] H.-J. Jou, P. Leo, J. Lowengrub, Microstructural evolution in inhomogeneous elastic media, Journal of Computational Physics 131 (1997) 109–148.  doi:10.1006/jcph.1996.5581.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999196955813

  * [45] T. Y. Hou, Z. Shi, Removing the stiffness of elastic force from the immersed boundary method for the 2D Stokes equations, Journal of Computational Physics 227 (21) (2008) 9138–9169.  arXiv:0801.2398, doi:10.1016/j.jcp.2008.03.002.
  * [46] T. Y. Hou, Z. Shi, An efficient semi-implicit immersed boundary method for the Navier-Stokes equations, Journal of Computational Physics 227 (20) (2008) 8968–8991.  doi:10.1016/j.jcp.2008.07.005.
  * [47] V. Cristini, J. Lowengrub, Q. Nie, Nonlinear simulation of tumor growth, Journal of Mathematical Biology 46 (3) (2003) 191–224.  doi:10.1007/s00285-002-0174-6.
URL http://link.springer.com/10.1007/s00285-002-0174-6

  * [48] S. Li, J. S. Lowengrub, P. H. Leo, V. Cristini, Nonlinear theory of self-similar crystal growth and melting, Journal of Crystal Growth 267 (3-4) (2004) 703–713.  doi:10.1016/j.jcrysgro.2004.04.002.
URL https://linkinghub.elsevier.com/retrieve/pii/S0022024804004014

  * [49] S. Li, J. S. Lowengrub, P. H. Leo, V. Cristini, Nonlinear stability analysis of self-similar crystal growth: control of the Mullins–Sekerka instability, Journal of Crystal Growth 277 (1-4) (2005) 578–592.  doi:10.1016/j.jcrysgro.2004.12.042.
URL https://linkinghub.elsevier.com/retrieve/pii/S0022024804020196

  * [50] S. Li, J. S. Lowengrub, P. H. Leo, Nonlinear morphological control of growing crystals, Physica D: Nonlinear Phenomena 208 (3-4) (2005) 209–219.  doi:10.1016/j.physd.2005.06.021.
URL https://linkinghub.elsevier.com/retrieve/pii/S0167278905002691

  * [51] W. Ying, C. S. Henriquez, A kernel-free boundary integral method for elliptic boundary value problems, Journal of Computational Physics 227 (2) (2007) 1046–1074.  doi:10.1016/j.jcp.2007.08.021.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999107003774

  * [52] W. Ying, W. C. Wang, A kernel-free boundary integral method for implicitly defined surfaces, Journal of Computational Physics 252 (2013) 606–624.  doi:10.1016/j.jcp.2013.06.019.
URL http://dx.doi.org/10.1016/j.jcp.2013.06.019

  * [53] Y. Xie, W. Ying, A fourth-order kernel-free boundary integral method for implicitly defined surfaces in three space dimensions, Journal of Computational Physics 415 (2020) 109526.  doi:10.1016/j.jcp.2020.109526.
URL https://doi.org/10.1016/j.jcp.2020.109526

  * [54] W. Ying, W. C. Wang, A kernel-free boundary integral method for variable coefficients elliptic pdes, Communications in Computational Physics 15 (4) (2014) 1108–1140.  doi:10.4208/cicp.170313.071113s.
  * [55] R. Peyret, Spectral Methods for Incompressible Viscous Flow, Vol. 148 of Applied Mathematical Sciences, Springer New York, New York, NY, 2002, pp. XII, 434.  doi:10.1007/978-1-4757-6557-1.
URL http://link.springer.com/10.1007/978-1-4757-6557-1

  * [56] G. C. Hsiao, W. L. Wendland, Boundary Integral Equations, Springer International Publishing, Cham, 2021, pp. 25–94.  doi:10.1007/978-3-030-71127-6_2.
URL https://doi.org/10.1007/978-3-030-71127-6_2

  * [57] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd Edition, Society for Industrial and Applied Mathematics, 2003.  doi:10.1137/1.9780898718003.
  * [58] T. Beale, A. Layton, On the accuracy of finite difference methods for elliptic problems with interfaces, Communications in Applied Mathematics and Computational Science 1 (1) (2006) 91–119.  doi:10.2140/camcos.2006.1.91.
URL http://www.physics.nyu.edu/$∼$fh417/ensemble-sampler-with-affine-
invariance.pdfhttp://msp.org/camcos/2006/1-1/p05.xhtml

  * [59] A. Brandt, N. Dinar, Multigrid Solutions to Elliptic Flow Problems, in: S. V. PARTER (Ed.), Numerical Methods for Partial Differential Equations, Academic Press, 1979, pp. 53–147.  doi:https://doi.org/10.1016/B978-0-12-546050-7.50008-3.
URL https://www.sciencedirect.com/science/article/pii/B9780125460507500083

  * [60] A. Brandt, I. Yavneh, On multigrid solution of high-reynolds incompressible entering flows, Journal of Computational Physics 101 (1) (1992) 151–164.  doi:https://doi.org/10.1016/0021-9991(92)90049-5.
URL https://www.sciencedirect.com/science/article/pii/0021999192900495

  * [61] J. Prüss, G. Simonett, Moving interfaces and quasilinear parabolic evolution equations, Vol. 105, Springer, 2016.
  * [62] K. Sakakibara, Y. Miyatake, A fully discrete curve-shortening polygonal evolution law for moving boundary problems, Journal of Computational Physics 424 (2021) 109857.  doi:10.1016/j.jcp.2020.109857.
URL https://doi.org/10.1016/j.jcp.2020.109857

  * [63] H. Chen, C. Min, F. Gibou, A numerical scheme for the stefan problem on adaptive cartesian grids with supralinear convergence rate, Journal of Computational Physics 228 (16) (2009) 5803–5818.  doi:10.1016/j.jcp.2009.04.044.
  * [64] H. Zhou, W. Ying, A Dimension Splitting Method for Time Dependent PDEs on Irregular Domains, Journal of Scientific Computing 94 (1) (2023) 20.  doi:10.1007/s10915-022-02066-5.
URL
https://doi.org/10.1007/s10915-022-02066-5%****␣main.bbl␣Line␣525␣****https://link.springer.com/10.1007/s10915-022-02066-5