\magnification = \magstephalf
\hsize=14truecm
\hoffset=1truecm
\parskip=5pt
\nopagenumbers
\input amssym.def
\font\eightrm=cmr8 \font\eighti=cmti8 \font\eightbf=cmbx8
\headline={\ifnum\pageno=1 \hfill\else%
{\tenrm\ifodd\pageno\rightheadline \else
\leftheadline\fi}\fi}
\def\rightheadline{EJDE--1995/11\hfil A Numerical Scheme\hfil\folio}
\def\leftheadline{\folio\hfil P.W. Bates, X. Chen, and X. Deng
\hfil EJDE--1995/11}
\voffset=2\baselineskip
\vbox {\eightrm\noindent\baselineskip 9pt %
Electronic Journal of Differential Equations
Vol. {\eightbf 1995}(1995) No. 11, pp. 1--28.\hfill\break
ISSN 1072-6691: URL: http://ejde.math.swt.edu (147.26.103.110)\hfil\break
telnet (login: ejde), ftp, and gopher access:
ejde.math.swt.edu or ejde.math.unt.edu}
\footnote{}{\vbox{\hsize=10cm\eightrm\noindent\baselineskip 9pt %
1991 {\eighti Subject Classification:} 35R35, 65C20, 65M06, 65R20, 82C26.
\hfil\break
{\eighti Key words and phrases:} Boundary Integral, Coarsening, Free
Boundary Problem, Motion by Curvature, Ostwald Ripening.
\hfil\break
\copyright 1995 Southwest Texas State University and
University of North Texas.\hfil\break
Submitted June 14, 1995. Published August 18, 1995.\hfil\break
Supported by: NSF grant DMS 9305044 (PWB), NSF grant DMS 9404773 and
Alfred P. Sloan Research Fellowship (XC).} }
\bigskip\bigskip
\centerline{A NUMERICAL SCHEME FOR THE TWO}
\centerline{PHASE MULLINS-SEKERKA PROBLEM}
\smallskip
\centerline{Peter W. Bates, Xinfu Chen, and Xinyu Deng}
\bigskip\bigskip
{\eightrm\baselineskip=10pt \narrower
\centerline{\eightbf Abstract}
An algorithm is presented to numerically treat a free boundary problem
arising in the theory of phase transition. The problem is one in which a
collection of simple closed curves (particles) evolves in such a way that
the enclosed area remains constant while the total arclength decreases.
Material is transported between particles and within particles by
diffusion, driven by curvature which expresses the effect of surface
tension. The algorithm is based on a reformulation of the problem, using
boundary integrals, which is then discretized and cast as a semi-implicit
scheme. This scheme is implemented with a variety of configurations of
initial curves showing that convexity or even topological type may not be
preserved.
\bigskip}
\input epsf
\def\gammat{\Gamma_{\!t}}
\def\pupn{\partial u \over \partial n}
\def\pwpn{\partial W_g \over \partial n}
\def\intgamma{\int_\Gamma}
\def\intgammao{\int_{\Gamma_0}}
\def\intgammat{\int_{\gammat}}
\def\xtoinf{|x|\rightarrow\infty}
\def\ov2pi{1\over2\pi}
\def\xyeps{x_0-y\pm\epsilon n(x_0)}
\def\xyepsf{\xyeps\over{|\xyeps|^2}}
\def\xzeroyf{{x_0-y}\over{|x_0-y|^2}}
\def\xzeroyepsf{{x_0-y}\over{|\xyeps|^2}}
\def\epsxyepsf{\epsilon n(x_0)\over{|\xyeps|^2}}
\def\epsxyepsfnon{\epsilon\over{|\xyeps|^2}}
\def\xyf{{x-y}\over{|x-y|^2}}
\def\tfpmeps{t^2+{\left(f(t)\pm\epsilon\right)}^2}
\def\sqrfprim{\sqrt{1+f^{'2}(t)}}
\def\xf{x\over|x|^2}
\def\oxsqr{O\left(1/|x|^2\right)}
\def\ox{O\left(1/|x|\right)}
\def\zrzh{z_R(t+h)-z(t+h)}
\centerline{\bf 1. Introduction} \medskip
Many formulations of time-dependent, multibody free-boundary problems
involve the solution of Laplace's equation in an irregular domain at each
instant of time as the free boundary evolves. For such problems boundary
integral techniques have been used and are a natural choice since they are
adaptive in that information only at the free surface is used to predict its
motion. Thus the dimension of the problem is reduced by one. Moreover, the
discretisation of the free surface can be done to resolve areas of
high curvature, a much more difficult task for other methods. In addition,
boundary integral techniques automatically account for far field asymptotic
behaviour.
Boundary integral techniques have been used successfully to study a
variety of free surface flows. Examples include the propagation of waves on a
fluid-fluid interface [4], the Rayleigh-Taylor instability problem [3, 18],
Hele-Shaw flow [9, 15], and crystal growth [14]. McFadden, Voorhees, Boisvert
and Meiron examined a multibody free boundary problem as it applies to Ostwald
ripening [17, 21, 22]. The model is a quasistatic version of one proposed by
Mullins and Sekerka to study the morphological development of a particle [19].
In this case, the free boundary consists of the interface separating two
distinct thermodynamic phases, and the dynamics of the process is due to the
tendency of the system to minimize the total interfacial surface area by
diffusing material from smaller particles to larger ones through the surrounding
region of second-phase material. The computational results in [17, 21], based on
a boundary integral formulation, give a detailed description of the local
behavior of single particles and small collections of particles. However,
diffusion within particles is neglected and the numerical analysis employs an
explicit time-stepping procedure. Because of the explicit nature of the
scheme, to maintain stability the algorithm requires the length of the
time-step to be scaled as the cube of the distance between mesh points.
Therefore, this procedure is unable, in practice, to compute close to the
disappearance of a particle. Implementation was done on a supercomputer with
vectorized code. Even then for a four particle system over five and a half
hours of CPU time was required to compute until the smallest particle reached a
mean radius of $.3$, its initial radius being $.9$.
Our purpose is to give a numerical treatment of the two phase
Mullins-Sekerka problem in two dimensions employing an algorithm which does not
suffer from severe stability constraints on the time step. We present a
semi-implicit scheme based on a different boundary integral formulation of the
problem. Implementation is done on a workstation, performing experiments with a
variety particle configurations, including a four particle system examined in
[21]. While the single phase model seems more appropriate when considering solid
particles of solute in a saturated solvent, the two phase model may be more
realistic when considering diffusion of atomic species in binary solids.
The particular model we treat can be described as follows: Given a
simple closed curve (or finite collection of nonintersecting simple closed
curves),
$\Gamma$, in the plane, find a harmonic function in ${\Bbb R}^2\backslash\Gamma$
which on $\Gamma$ is equal to the curvature of $\Gamma$. This function is
continuous but not smooth across $\Gamma$ (unless $\Gamma$ is a collection of
circles of equal radii). Now evolve
$\Gamma$ with normal velocity equal to the jump in the normal derivative across
$\Gamma$ of the harmonic function.
The problem can also be posed in a bounded domain, $\Omega$, by
imposing homogeneous Neumann conditions on $\partial\Omega$. When $\Omega={\Bbb
R}^2$ we require a certain decay on the gradient as $\vert x\vert\to\infty$.
Pego [20] derived this model, using formal asymptotic expansions, to
describe the evolution of level sets of layered solutions to the Cahn-Hilliard
equation. In [1] this connection was rigorously proved assuming the existence
of smooth solutions to the Mullins-Sekerka flow. At the same time uniqueness
of solutions was established. In [5], the existence of weak solutions was
proved and later in [6] the existence of smooth solutions was established (see
also Constantin and Pugh [8] for results concerning a related problem).
In [5] it is shown that the evolution is curve shortening and area
preserving so one is naturally drawn to the questions of whether the flow
preserves convexity of a single particle, whether or not a single particle can
split in two, and whether or not two particles can coalesce.
In [16], Mayer proved that for the one (exterior) phase problem,
convexity could be lost. In our work, numerical evidence is presented which
suggests that this is also the case for the two phase problem. Experiments are
also presented which suggest that pinching-off and coalescence are not excluded
by this model.
In section 2 we present the model in its original form and establish
an equivalent formulation using boundary integrals. This is then used to devise
a semi-implicit algorithm for advancing the curve. This involves linearizing the
operator which maps normal velocity to curvature of the curve at the next time
step under the corresponding explicit scheme. In section 3 we discretize the
algorithm, computing geometric quantities such as tangent vectors and curvatures
using natural differencing ideas. In section 4 we present the results of
experiments performed by using the procedure developed, starting with an accuracy
check with a configuration for which we know the solution, namely, the evolution
of concentric circles. Here we use a variety of time step sizes and numbers of
mesh points to aid in our selection of these parameters for subsequent
experiments and also to estimate the convergence rate as mesh size and time step
tend to zero. However, no attempt is made here to rigorously establish
convergence of our general scheme. The next experiment shows the evolution of a
nonconvex ``rose'' curve as it quickly becomes convex and tends to a circle,
keeping the enclosed area constant. The third experiment starts with a convex
particle, a tube with circular endcaps, which at first loses its convexity and
then regains it, eventually becoming circular. We then conduct experiments to
show that a single particle can ``pinch off'' to become two and two particles
can coalesce to become one. Our final experiment uses a configuration of four
particles arranged as in an experiment reported in [21] for the one phase
problem.
Our code, written in $C$, is contained in an appendix to this paper.
After this work was completed we were made aware of several interesting results
for related problems. In [12] the authors present a way to handle integrals
with singular parts arising in fluid interface problems such as Hele-Shaw. A
detailed analysis of singularity formation is given in [2] where numerical and
analytical methods provide convincing evidence for the formation of
singularities in Hele-Shaw flow. Also, the second author, with Hou and Zhu,
have incorporated some ideas from [12] with other time saving techniques to
substantially accelerate an algorithm similar to that presented here.
\noindent{\bf Acknowledgement:} The first author would like to thank Giorgio Fusco for helpful
suggestions and discussions. We would also like to thank the referee for a
careful reading of the manuscript. Much of this work was reported in the
Master's thesis [10] of the third author.
\bigskip
\centerline{\bf 2. Equivalent formulations and algorithm }\medskip
Let $\Omega$ be a bounded and simply connected domain in ${\Bbb R}^2$
and $\Gamma_0$ be
a smooth simple closed curve (or finite collection of such) in $\Omega$.
Consider the free boundary problem of finding a function
$u(x,t), x \in \Omega, t \geq 0$, and a free boundary $\Gamma_{\!0,T} =
\cup_{0\leq t\leq T} \left(
\gammat \times \{t\}\right) $ for some $T>0$ satisfying
$$\cases{\ a)\quad \Delta u(\cdot,t) = 0&in $\Omega \backslash \gammat$, $0 \leq
t \leq T$, \cr
\ \cr
\ b)\quad {\pupn} = 0&on $\partial\Omega \times [0, T]$, \cr
\ \cr
\ c)\quad u = K&on $\gammat$, $0 \leq t \leq T$, \cr
\ \cr
\ d)\quad -[{\pupn}]_{\gammat}
= V&on $\gammat$, $0 \leq t \leq T$, \cr
\ \cr
\ e)\quad \Gamma_{\!0,T} \cap \{t = 0\} = \Gamma_0, \cr}
\eqno(2.1)'
$$
where $n$ is the unit outward normal to $\partial\Omega$ or to $\gammat$,
$[{\pupn}]_{\Gamma_{\!\!t}}$ is the sum of the
{\it outward} normal derivatives
of $u$ from each side of $\gammat$ enclosing a bounded region (which is also
equal to the jump of the normal derivative of $u$ across
$\gammat$), $K$ is the curvature, taking the sign convention that convex curves
have positive curvature, while $V$ is the normal velocity of $\gammat$ and the
normal velocity of an expanding curve enclosing the bounded region is
positive. For simplicity, we only consider the case in which $\Omega = {\Bbb
R}^2$. In this case, we replace the boundary condition (2.1b)$'$ by
$$ \nabla u(\cdot, t) = \oxsqr \qquad\qquad \hbox{as }\xtoinf$$
and call the resulting system (2.1).
We first derive an integral representation for the solution to Eqs.
(2.1).
\medskip
\noindent {\bf Lemma 2.1.} {\sl Let $\Gamma = \cup_{i=1}^M\Gamma^i$ $(M\geq 1)$,
where
$\Gamma^i$ are disjoint simple closed curves such that $\Gamma$ separates
${\Bbb R}^2$ into one unbounded and $m$ bounded regions. Let $n$ be the
outward unit normal to $\Gamma$. For each
$g
\in L^2(\Gamma)$, define
$$ W_g(x) = {\ov2pi} \int_\Gamma \ln |x-y|\>g(y)\>ds_y,\qquad x \in {\Bbb
R}^2\backslash \Gamma.$$
Then the following holds:
(1) $\Delta W_g = 0 \qquad$ in ${\Bbb R}^2 \backslash \Gamma$;
(2) $\displaystyle{-\left[{\pwpn}\right]}_\Gamma \equiv {\pwpn}^{out}
- {\pwpn}^{in} = g$ \qquad on $\Gamma$;
(3) if in addition we assume that$\>\displaystyle{\intgamma
g(y)\>ds_y\>\>=\>\>0}$, \quad then
\qquad$\displaystyle{|{\nabla W_g}| = \oxsqr}$ \quad and \quad
$\displaystyle{W_g = \ox}$ \qquad as $\xtoinf$;
(4) The mapping $g\in\{ g \in L^2(\Gamma)|\intgamma g =
0\}\longmapsto W_g$ is negative definite, i.e.
\qquad $\displaystyle{(g,W_g) \equiv \intgamma g\>W_g}\> < 0$\qquad if
$g\not\equiv0$.}
\noindent {\bf Proof:}
The first assertion of the lemma follows from the fact that
$$ \Delta_x \ln |x-y|=0 \qquad \hbox{if } x\not=y.$$
The second is a standard calculation in potential theory (see for instance [11]
section 3D or [13]).
To see (3), note that when $|x|$ is large,
$$\eqalign{\nabla W_g &= {\ov2pi}\intgamma {\xyf}\>g(y)\>ds_y\cr &=
{\ov2pi}\intgamma {\left({\xyf}-{\xf}\right)}\>g(y)\>ds_y,\cr}$$
since $\intgamma g(y)\>ds_y=0$. Also,
$$\eqalign{{\xyf}-{\xf}
&={{|x|^2-|x-y|^2}\over{|x-y|^2\cdot|x|^2}}\>x-{y\over{|x-y|^2}}\cr &=
\oxsqr.\cr}$$
Hence,
$$\nabla W_g(x)=\oxsqr.$$
Similarly,
$$\eqalign{W_g &= {\ov2pi}\>\intgamma
\left[{\ln|x-y|\>g(y)-\ln|x|\>g(y)}\right]\>ds_y\cr &= {\ov2pi}\>\intgamma
\ln\left({|x-y|\over|x|}\right)\>g(y)\>ds_y\cr &= {\ov2pi}\>\intgamma
\ln\left({1+\ox}\right)\>g(y)\>ds_y\cr &= \ox.\cr}$$
The third assertion of the lemma then follows.
Finally, set $u=W_g$, then
$$\eqalign{-\intgamma W_g(y)g(y)\>ds_y &= \intgamma u
\left[{\pupn}\right]_{\Gamma} \>ds_y\cr
&= \int\!\!\!\int_{{\Bbb R}^2} {u\Delta u}
+\int\!\!\!\int_{{\Bbb R}^2} {\nabla u\cdot \nabla u}\cr
&=
\int\!\!\!\int_{{\Bbb R}^2} {|\nabla u|^2}\>dxdy > 0,\cr}$$
where in the second equation we have used the assertion of (3).
This completes the proof of the lemma.\quad\vrule height4pt width3pt
The previous lemma specifies a zero mean value jump across $\Gamma$ of the
normal derivative of a proposed function which is harmonic on each side of
$\Gamma$ and produces such a harmonic function. Clearly, any constant can be
added to the harmonic function produced by the lemma giving another suitable
function. The next lemma shows that any harmonic function on ${\Bbb
R}^2\backslash\Gamma$ having a sufficiently regular trace on $\Gamma$ has this
trace realized by the harmonic function constructed above from the jump in the
normal derivative, up to an additive constant.
\medskip
\noindent {\bf Lemma 2.2.} {\sl Suppose that $\Gamma$ is as in Lemma 2.1.
Suppose that $u$ defined on ${\Bbb R}^2$, $f\in H^1(\Gamma)$ and $g\in
L^2(\Gamma)$ are related by
$$
\left.{\matrix{\Delta u &=0\hfill \hbox{ \quad \ \ in } {\Bbb R}^2\backslash\Gamma\hfill\cr
\vert\nabla u\vert &=O\left({{1\over \vert x\vert^2}}\right)\hfill \hbox{ as
}\vert x\vert\to\infty\hfill\cr
u &=f\hfill \hbox{ \ on }\Gamma\hfill\cr
-\left[{{\partial u\over\partial n}}\right]_{\Gamma}&= g\hfill \hbox{ \ on
}\Gamma.\hfill}}\right\}\eqno(2.2)$$
Then there exists a constant $c$ such that
$$f(x)={1\over 2\pi}\int_{\Gamma} \ln\vert x-y\vert g(y)ds_y+c\hbox{ for }
x\in\Gamma.\eqno(2.3)$$
Furthermore,}
$$\int_{\Gamma} g(y)ds_y=0.\eqno(2.4)$$
\noindent {\bf Proof:}
Define $W_g$ as in Lemma 2.1 and
$$c(x)=u-W_g.$$ Then $\Delta c=0$ in ${\Bbb R}^2\backslash\Gamma$ and $c$ is
continuous across
$\Gamma$ since both $u$ and $W_g$ are. Furthermore, $\left[{{\partial c\over
\partial n}}\right]_{\Gamma}=0$ and hence $\Delta c=0$ in ${\Bbb R}^2$.
Finally, $\nabla c= \oxsqr$ as $\vert x\vert\to\infty$ and so $c$ is a
constant. This yields (2.3). To establish
(2.4), let $B_R$ be a disk of radius $R$, then, as $R\to\infty$,
$$\int_{\Gamma} g=-\int_{\Gamma} \left[{{\partial u\over\partial n}}\right]
=\int_{\partial B_R} {\partial u\over\partial
n}-\int\!\!\!\!\int_{B_R\backslash\Gamma}
\Delta u=2\pi R\cdot O(1/R^2)\to 0.\eqno \quad\vrule height4pt width3pt$$
The above lemmas can now be combined to show the equivalence of (2.1) and the
system (2.3)--(2.4) with $f=K$ and $g=V$, the curvature and normal velocity,
respectively, of $\Gamma$.
\medskip
\noindent {\bf Theorem 2.3} {\sl If $\Gamma_{0,T}\equiv\cup_{0\leq t\leq T}
\left({\Gamma_t\times\{t\}}\right)$ is a continuous family of $C^3$ curves which
satisfy (2.3) and (2.4) with $\Gamma=\Gamma_t$ for some $g=G(x,t)$ and
$c=c(t)$, where $f=K(x,t)$ is the curvature of $\Gamma_t$ at $x$, then
$\Gamma_{0,T}$ is the interface associated with the solution to (2.1) and
$g$ is the normal velocity, $V$, of $\Gamma_t$.
Conversely, if $(u,\Gamma_{0,T})$ is a solution to (2.1) then (2.3)--(2.4)
hold for each $t\in [0,T]$ with $\Gamma=\Gamma_t, g=V$ and $f=K$.}
\noindent {\bf Proof:}
Let $(\Gamma_t, g, c)$ be a solution to (2.3)--(2.4) at each time $t\in [0,T]$
with $f=K(\cdot, t)$. Then Lemma 2.1 shows that defining $u(\cdot, t)$ on
${\Bbb R}^2\backslash\Gamma_t$ by
$$u(x,t)={1\over 2\pi} \int_{\Gamma_t} \ln\vert x-y\vert g(y,t)ds_y+c(t)$$ gives
a solution to (2.1) with the normal velocity of $\Gamma_t$ being given by $g$.
Solutions to (2.1) with initial curve $\Gamma_0$ are unique by the results in
$[1]$ and so $\Gamma_{0,T}$ is that solution.
Conversely, if $(\Gamma_{0,T}, u)$ is a solution to (2.1) then the solution is
smooth by the results in [6]. By Lemma 2.2, (2.3)--(2.4) hold for some $c$ with
$\Gamma=\Gamma_t, f=K(\cdot, t)$, and $g=V(\cdot, t)$, for each
$t\in [0,T]$.\quad\vrule height4pt width3pt
Before we present a numerical scheme for solving (2.3)--(2.4), we first solve an
inverse problem:
\medskip
\noindent {\bf Lemma 2.4} {\sl Given $f\in H^1(\Gamma)$, with $\Gamma$ as in
Lemma 2.1, there exists a unique $g\in L^2(\Gamma)$ and constant $c$ such that
(2.3) and (2.4) hold.}
\noindent {\bf Proof:}
We first prove uniqueness. Since (2.3), (2.4) is a linear system it suffices to
show that $f\equiv 0$ implies $g\equiv 0$ and $c=0$. Let $w=W_g+c$, then
$$\Delta w=u\, \hbox{ in }\, {\Bbb R}^2\backslash\Gamma$$ and on $\Gamma, w=f=0$.
Furthermore, by Lemma 2.1, $\vert\nabla w\vert=\oxsqr$
as $\vert x\vert\to\infty$. It follows that $w\equiv 0$ on
${\Bbb R}^2$ and by Lemma 2.1 (2), $g\equiv 0$ on $\Gamma$. Consequently,
$c=0$ also.
To prove existence of a solution let $u$ satisfy
$$\eqalign{\Delta u&=0 \hbox{ on } {\Bbb R}^2\backslash\Gamma,\cr
u&=f \hbox{ on } \Gamma,\cr
\vert\nabla u\vert&= \oxsqr \hbox{ as } \vert x\vert\to\infty.}$$
Since $f\in H^1, u$ exists, is unique, and $u\in H^2$. Set
$$g= -\left[{{\partial u\over\partial n}}\right]_{\Gamma}.$$ Then, $g\in
L^2(\Gamma)$ and by Lemma 2.2, $\int_{\Gamma}g=0$ and $c\equiv u-W_g$ is
constant.\quad\vrule height4pt width3pt
From the previous discussion, (2.1) can now be solved step by step as
follows:
Assume that $\Gamma_t$ is known. Calculate $f=K(\cdot,t)$ on $\Gamma_t$, then
find
$g=V(\cdot, t)$ and $c(t)$ by solving
$$\eqalign{&K(x,t) = {\ov2pi}\>\int_{\Gamma_t} \ln|x-y|V(y,t)\>ds_y + c(t),\cr
&\int_{\Gamma_t} V(y,t) \>ds_y = 0.\cr}\eqno(2.5)$$
Once we know $V$, we can advance the curve by
$$x(t+dt)=x(t)+Vndt.\eqno(2.6)$$
This is an explicit scheme. The weakness of the scheme is that the time
step has to be extremely small due to instablity. A small time step, however,
introduces more machine error. We avoid these problems through the following
semi-implicit scheme. It linearizes the mapping which uses the explicit method
to get the curvature at one time step from the normal velocity at the previous
step:
Write
$$\eqalign{&\Gamma_t =\left\{ x=x(s,t)\>\>\vert \>\>s\in S^1=R^1 \>\>
(\hbox{mod} \>\>
2\pi)\right\},\cr
&k(s,t) = K(x(s,t), t),\cr
&V(s,t) = {\partial x\over\partial
t}\cdot n.\cr}$$
Assume $x(s,t)$ is known for some t, we find $V$ by solving
$$\left\{\eqalign{&k(s,t)+BV\Delta t = {\ov2pi}\>\int_{\Gamma_t}
\ln|x(s,t)-y|V\>ds_y+c,\cr
&\int_{\Gamma_t} V(y,t) \>ds_y =
0,\cr}\right.\eqno(2.7)$$
where
$$BV={\partial K(x(s,t)+hnV)\over\partial h}\Big\vert_{h=0}$$
defines a linear operator.
We then advance $\Gamma$ according to (2.7). Note that $k(s,t)+BV\Delta t$ is
an approximation of the curvature at $t+\Delta t$.
\bigskip
\centerline{\bf 3. Discretization }\medskip
Assume that $\Gamma_t=\cup_{m=1}^M{\Gamma^m_t}$ consists of $M$ disjoint simple
closed curves and take $N$ points from each $\Gamma^m_t$, labelling them by
$z_{(m-1)N+1}$, $z_{(m-1)N+2}$, $\dots$, $z_{mN}$ listed counterclockwise. Then,
to deal with the periodicity, we define the following permutation functions:
$$\eqalign{L[i]&=\cases{\ i-1&if ${{i-1}\over N} \not= $ integer, \cr
\ i-1+N&if ${{i-1}\over N} = $ integer, \cr}\cr R[i]&=\cases{\ i+1&if
${i\over N} \not= $ integer, \cr
\ i+1-N&if ${i\over N} = $ integer, \cr}\cr} \eqno(3.1)$$
to denote the
indices of the points clockwise and counterclockwise from $z_i$, respectively.
Assume that $z_L$, $z$, $z_R$, are three counterclockwise consecutive points on
$\Gamma$, we interpolate $\Gamma$ near $z$ as a segment of the circle passing
through $z_L$, $z$, $z_R$. We use the unit tangent,
$\tau$, the outward unit normal, $n$, and curvature $K$ of the circle passing
through $z_L$, $z$, $z_R$ as the corresponding quantities of $\Gamma$(see Fig.
1). If $z_L$, $z$, $z_R$ are colinear, then the curvature is zero and the
tangent is the obvious one.
We denote
$$\eqalign{T_L&={z-z_L\over|z-z_L|}, \qquad T_R={z_R-z\over |z_R-z|},\cr
d_L&=|z-z_L|, \qquad d_R=|z_R-z|,\cr
\hbox{ and }\qquad\qquad d_{RL}&= |z_R-z_L|.}\eqno(3.2)$$ Note that if a tangent
vector has coordinates $\tau=(\tau^x,\tau^y)$ then
$n=(\tau^y,-\tau^x)$ is a unit normal. This motivates the definitions
$$N_L=(T^y_L,-T^x_L)\hbox{ and } N_R=(T^y_R,-T^x_R)$$ as approximate normals,
which are outward by our ordering of the $z_j$'s. The following is easily
demonstrated and we omit the proof:
\centerline{\epsfbox{Fig1.ps}}
\centerline{Fig. 1 Interpolation}
\medskip
\noindent{\bf Lemma 3.1.} {\sl Let $C$ be the circle passing through noncolinear
$z_L,z,z_R$ and let $\tau,n$ and $K$ be the counterclockwise unit tangent,
outward unit normal, and curvature of $C$ at $z$, respectively. Then,}
$$\eqalign{\tau&={d_RT_L+d_LT_R\over d_{RL}},\cr n&={d_RN_L+d_LN_R\over
d_{RL}},\cr K&={2T_L\cdot N_R\over d_{RL}}=-{2T_R\cdot N_L\over
d_{RL}}.}\eqno(3.3)$$
Let $\Gamma=\cup_{j=1}^{MN}{\Gamma_j}$, where $\Gamma_j$ is a segment of $\Gamma
$ containing $z_j$ in its interior. We take these segments small so that a
given function $g$ on $\Gamma$ is almost constant on
$\Gamma_j$. Define
$$W_g(z') = {\ov2pi}\int_{\Gamma} \ln|z'-z|g(z)\>ds_z\eqno(3.4)$$ and set $d_i
=
\int_{\Gamma_i} ds_i$, $W^i={1\over d_i}\int_{\Gamma_i} W_g(z')\>ds_z$, and
$g_j=g(z_j)$. Then
$$W^i \simeq \sum_{j=1}^{MN} a_{ij}g_j \quad \hbox{for}\quad i= 1,2,\dots,MN,
\eqno(3.5)$$ where
$$a_{ij} = {1\over{2\pi d_i}} \int_{\Gamma_{\!i}}\!\int_{\Gamma_{\!j}} \ln|z'-z|
\>ds_{z'}ds_z.\eqno(3.6)$$ The remaining work is for the approximation of
$a_{ij}$.
\vbox{
\epsfbox{Fig2.ps}
\centerline{Fig. 2. Approximation of $a_{ij}$}}
We are seeking second order approximation, so we can replace $\Gamma_i$
and $
\Gamma_j$ in (3.6) by line segments $\hat \Gamma_i$ and $\hat \Gamma_j$ (Fig. 2).
Define $\tau_i$ as the unit tangent at point $z_i$, as in (3.3), and take
$$\eqalign{d &= d_{ij} = |z_i-z_j|,\cr
d_i &= {1\over 2}
\left(d_{R[i]}+d_{L[i]}\right),\cr
\tau_{ji} &= \cases{\ \tau_i&if $i=j,$ \cr
\ {{z_j-z_i}\over d}&if $i\not= j,$ \cr}\cr
A &= \cos\alpha =
-\tau_i\cdot \tau_{ji},\cr
B &= \cos\beta = \tau_j\cdot \tau_{ji}.\cr}$$
We have
$$ \eqalign{ a_{ij} &= {1\over {2\pi} d_i} \int_{\Gamma_{\!i}}\!\int_{\Gamma_{\!
j}} \ln|z'-z|\>ds_{z'}ds_z\cr
&\cong {1\over {2\pi} d_i}\int_{-{d_{L[i]}\over
2}}^{d_{R[i]}\over 2}\!\int_{-{d_{L[j]}\over 2}}^{d_{R[j]}\over 2}
\ln|d+Ax+By|\>dy dx \cr
&=
\Biggl[{1\over{4{\pi}d_iAB}}[d+Ax+By]^2\left(\ln|d+Ax+By|-{3\over 2}\right)
\Biggr]\Biggl|_{y=-{d_{L[j]}\over 2}}^{y={d_{R[j]}\over 2}}
\Biggr|_{x=-{d_{L[i]}\over 2 }}^{x={d_{R[i]}\over 2}} \cr
&={1\over
{2{\pi}AB\left(d_{R[i]}+d_{L[i]}\right)}}
\Biggl\{a_1^2\left(\ln|a_1|-{3\over 2}\right)+a_2^2\left(\ln|a_2|-{3\over
2}\right)}$$
$$-a_3^2\left(\ln|a_3|-{3\over 2}\right)-a_4^2\left(\ln|a_4|-{3\over 2}\right)
\Biggr\} \equiv \hat a_{ij},\eqno(3.7)$$
where
$$ \eqalign{a_1 &= d+A{d_{R[i]}\over 2}+B{d_{R[j]}\over 2}, \cr
a_2 &=d-A{d_{L[i]}\over 2}-B{d_{L[j]}\over 2}, \cr
a_3 &= d+A{d_{R[i]}\over2}-B{d_{L[j]}\over 2}, \cr
a_4 &= d-A{d_{L[i]}\over2}+B{d_{R[j]}\over 2}. \cr}$$
A more detailed analysis shows that $|a_{ij}-\hat a_{ij}| = O(r^2)$,
where
$r=\max |z_i-z_{R[i]}|$.
Note that $W^i\simeq W_g(z_i)$ and so with $V_j=V(z_j)$ and $K_j=K(z_j)$ for
fixed $t$, system (2.5) is discretized to give the following system, where we
denote $\hat a_{ij}$ by $a_{ij}$.
$$ \eqalign{&\sum_{j=1}^{MN}a_{ij}V_j + C = K_i\quad \hbox{for} \quad i=1,\dots,
MN,\cr &\sum_{j=1}^{MN}d_jV_j = 0,\cr}\eqno(3.8)$$
where
$$ \eqalign{U &\equiv (V_1, \dots, V_{MN}, C)^T \hbox{ are unknowns,}\cr
K_i &=
{-2T_{R[i]}\cdot N_{L[i]}\over d_{RL[i]}},\cr
d_j &= {1\over 2}
\left(d_{R[j]}+d_{L[j]}\right),\cr
a_{ij} &= {1\over
{2{\pi}AB\left(d_{R[i]}+d_{L[i]}\right)}} \Biggl\{a_1^2\left(
\ln|a_1|-{3\over 2}\right)+a_2^2\left(\ln|a_2|-{3\over 2}\right) \cr
&\qquad-a_3^2\left(\ln|a_3|-{3\over 2}\right)-a_4^2\left(\ln|a_4|-{3\over
2}\right)\Biggr\}, \cr}$$
where $a_1$, ..., $a_4$ are as above with $d=d_{ij}=|z_i-z_j|$.
Hence, we find the solutions $U$, by solving (3.8) and the explicit
scheme updates the location of the interface through the formula
$$ z(t+h) = z(t)+hn(z(t))V,\eqno(3.9)$$
where $h$ is the time step and $n$ is the
outward unit normal.
In (3.8), if we let $K_i$ be the curvature evaluated at $z_i(t)$ then the
scheme is explicit and unstable unless $h$ is sufficiently small. If we evaluate
$K_i$ at
$z_i(t+h)$ given in (3.9), then the scheme is implicit and stable for larger
$h$. In this case, $K_i$ depends on $V$ in a non-linear and non-local manner. For
ease of computation, we take a semi-implicit scheme by extracting the linear
part of this dependence and ignoring the remainder. From experiments we
performed, we see that this is enough for the stability of the scheme, as well as
the accuracy (see section 4). To extract the linear part of the dependence of
curvature on velocity, we compute the first derivative of $K$ with respect to
$h$. Since
$$K=-{{2T_R\cdot N_L}\over d_{RL}}$$ then
$$ {\partial K \over \partial h} = -{{2{{\partial T_R}\over{\partial h}}\cdot N_
L}\over d_{RL}}-{{2T_R\cdot{{\partial N_L}\over{\partial h}}}\over d_{RL}}+{{2T_
R\cdot N_L {\partial d_{RL}\over\partial h}}\over d_{RL}^2}.$$ Since
$$\eqalign{T_R\cdot {\partial N_L\over\partial h} &=(T_R^x, T_R^y)\left (
{\partial N_L^x\over\partial h}, {\partial N_L^y\over\partial h}\right )\cr
&=(T_R^x, T_R^y) \left ( {\partial T_L^y\over\partial h}, -{\partial
T_L^x\over\partial h}\right )\cr &=T_R^x{\partial T_L^y\over\partial h} -
T_R^y{\partial T_L^x\over\partial h}\cr &=-\left(N_R^y{\partial
T_L^y\over\partial h} + N_R^x{\partial T_L^x\over\partial h}\right)\cr
&=-N_R\cdot {\partial T_L\over\partial h}\cr}$$ we get
$$ {\partial K \over \partial h} = -{{2{{\partial T_R}\over{\partial h}}\cdot N_
L}\over d_{RL}}+{{2N_R\cdot{{\partial T_L}\over{\partial h}}}\over d_{RL}}-{{K{
\partial d_{RL}\over\partial h}}\over d_{RL}}.\eqno(3.10)$$
From (3.2) and (3.9), We have
$$ \eqalign{{\partial T_R\over\partial h} &= {\partial \over\partial
h}\left(\zrzh\over {|\zrzh|}\right)\cr &= {{{\partial \over \partial
h}\left(\zrzh\right)}\over {|\zrzh|}}\cr &-{\zrzh\over{|\zrzh|}^3} \left[(\zrzh)
\cdot {\partial
\over \partial h} \left(\zrzh\right)\right]\cr &= {{n_RV_R-nV}\over{|\zrzh|}}\cr
&-{\zrzh\over{{|\zrzh|}^3}} \left[(\zrzh)\cdot (n_RV_R-nV)\right],\cr}$$
where we
use $n$, $n_L$ and $n_R$ to denote unit normals at $z$, $z_L$ and
$z_R$, respectively (see Fig. 1).
Sending h $\rightarrow$ 0, since $T_R = {z_R-z\over |z_R-z|}$ and
$N_R(N_R\cdot x)+T_R(T_R\cdot x) = x,$ we obtain
$$ \eqalign{{\partial T_R\over\partial h}\bigg|_{h=0} &= {1\over {|z_R-z|}} \left
\{(n_RV_R-nV)-T_R(T_R\cdot (n_RV_R-nV))\right\}\cr &= {N_R\over d_R}
\left\{N_R\cdot (n_RV_R-nV)\right\}.\cr}\eqno(3.11)$$
Similarly,
$${\partial T_L\over\partial h}\bigg|_{h=0} = {N_L\over d_L} \left\{N_L\cdot (nV-
n_LV_L)\right\}.$$
We also compute
$$\eqalign{&{\partial d_{RL}\over\partial h} = {1\over d_{RL}} \left\{(z_R(t+h)-
z_L(t+h))\cdot (n_RV_R-n_LV_L)\right\}\cr &{\partial d_{RL}\over\partial
h}\bigg|_{h=0} = {1\over d_{RL}} \left\{(z_R-z_L)
\cdot (n_RV_R-n_LV_L)\right\}.\cr}$$
Then (3.10) can be written as
$$ \eqalign{{\partial K \over \partial h}\bigg|_{h=0} &= -{2\over {d_R d_{RL}}}((
n_RV_R-nV) \cdot N_R) (N_R \cdot N_L)\cr &\quad +{2\over {d_L
d_{RL}}}((nV-n_LV_L) \cdot N_L) (N_L \cdot N_R)\cr &\quad-{K\over
d_{RL}^2}(z_R-z_L)\cdot (n_RV_R-n_LV_L)\cr &=V_L\left(-{2(N_L \cdot n_L) (N_R
\cdot N_L) \over{d_L d_{RL}}}+K{(z_R-z_L)
\cdot n_L\over d_{RL}^2}\right)\cr &\quad+V\left({2(N_R \cdot n) (N_R \cdot
N_L) \over{d_R d_{RL}}}+{2(N_L \cdot
n) (N_R \cdot N_L) \over{d_L d_{RL}}}\right)\cr &\quad+V_R\left(-{2(N_R \cdot
n_R) (N_R \cdot N_L) \over{d_R d_{RL}}}+K{(z_L-z_R) \cdot n_R\over
d_{RL}^2}\right).\cr}\eqno(3.12)$$ Hence, we have
$$K(t+h)\simeq K(t)+hBV,\eqno(3.13)$$ where
$$h B = \{b_{ij}\}\quad \hbox{for} \quad i,j=1,\dots, MN $$ and
$$ \eqalign{&b_{iL[i]} = h\left(-{2(N_L \cdot n_L) (N_R \cdot N_L) \over{d_L
d_{RL}}}+K{(d_RT_R+d_LT_L) \cdot n_L\over d_{RL}^2}\right)\cr
&b_{ii} =
h\left({2(N_R \cdot n) (N_R \cdot N_L) \over{d_R d_{RL}}}+{2(N_L \cdot
n) (N_R \cdot N_L) \over{d_R d_{RL}}}\right)\cr
&b_{iR[i]} = -h\left({2(N_R
\cdot n_R) (N_R \cdot N_L) \over{d_R d_{RL}}}+K{( d_RT_R+d_LT_L) \cdot n_R\over
d_{RL}^2}\right)\cr &\qquad\qquad \hbox{for } i=1,\dots, MN,\cr
&b_{ij}= 0
\qquad\hbox{for } j\not= i, \> R[i], \hbox{or } L[i]; \>\> i,j=1,
\dots, MN.\cr}\eqno(3.14)$$
Now (3.8) can be modified as
$$ \left\{\eqalign{&\sum_{j=1}^{MN}(a_{ij}-\beta b_{ij})V_j + C = K_i\quad
\hbox{for}\quad
i=1,\dots, MN,\cr &\sum_{j=1}^{MN}d_jV_j = 0\cr}\right.\eqno(3.15)$$ where
$\beta\in [0,1]$ is a factor indicating the extent of the implicit nature of the
scheme. Several experiments using different values for $\beta$ (not reported
here) convinced us that taking $\beta=1$ provides a stable scheme while
maintaining accuracy, therefore all our subsequent experiments take
$\beta=1$.
We may write (3.15) in matrix form
$$GU = P,\eqno(3.16)$$ where
$$G=\pmatrix{a_{ij}-b_{ij}&1\cr d&0\cr}, \quad P=\pmatrix{K\cr 0\cr}, and \quad
U=\pmatrix{V\cr C\cr}.$$
After solving (3.16) for $U$ we update $z$ using (3.9).
Our algorithm is based on the premise that each boundary component $\Gamma_t^m$
is described
by $N$ mesh points. During initialization, the $N$ points are chosen to be
equispaced in arc length of the boundary. After each iteration, new mesh points
will be generated, representing the evolution of the boundary. These new points
may not be equispaced and, unless we are careful, may concentrate at certain
locations, leading to computational errors. To avoid this problem, we shall
redistribute these newly generated points so that
they are almost equispaced.
\centerline{\epsfbox{Fig3.ps}}
\centerline{Fig. 3 Reparameterization}
Let $z_0$ be the center of the circle passing through three newly created
points: $z$, $z_L$, and $z_R$, and let $z^*$ be the mid-point of the arc
$\widehat {z_Lz_R}$. We shall take $z^*$ as the new location for the point $z$.
The outward unit normals $n$, $N_L$, and $N_R$ are given by the definition
following (3.2), while
$n_{LR}$ represents the outward unit normal at $z^*$ (refer to Fig. 3). To solve
for
$z^*$, we write $$ \eqalignno{z_0 &= z-{1\over K}n &(3.17)\cr z^* &= z_0 +{1\over
K}n_{LR} = z+{1\over K}\left(n_{LR}-n\right). &(3.18)\cr}$$
Replacing $n$ and $K$ in (3.18) using (3.3), we have
$$\eqalign{z^*&=z+{d_{RL}\over 2T_L\cdot N_R}
\left(n_{LR}-{d_RN_L+d_LN_R\over d_{RL}}\right)\cr
&=z+{d_{RL}n_{LR}-d_RN_L-d_LN_R\over 2T_L\cdot N_R}.}\eqno(3.19)$$ Since
$$d_{RL}n_{LR} = {\left(d_{RL}T_{RL}\right)}^\perp={\left(d_RT_R+d_LT_L\right)}^
\perp= d_RN_R+d_LN_L,$$ (3.19) can be written as
$$\eqalign{z^* &= z+{{d_RN_R+d_LN_L-d_RN_L-d_LN_R}\over{2T_L\cdot N_R}}\cr &=
z+{(d_R-d_L)(N_R-N_L)\over{2T_L\cdot N_R}}.\cr}\eqno(3.20)$$ If we write
$$N_R-N_L = aT_L +bT_R,$$ then we find
$$ \eqalign{a=b &={{1-N_L\cdot N_R}\over {T_L\cdot N_R}}={{1-T_L\cdot T_R}\over
{T_L\cdot N_R}}={{1-{(T_L\cdot T_R)}^2}\over {\left(1+(T_L\cdot T_R)\right)(T_L
\cdot N_R)}}\cr &= {(T_L\cdot N_R)^2\over {\left(1+(T_L\cdot
T_R)\right)(T_L\cdot N_R)}}= {{T_L
\cdot N_R}\over {1+(T_L\cdot T_R)}}\cr}\eqno(3.21)$$ and so
$$N_R-N_L = {{(T_L\cdot N_R)(T_L+T_R)}\over{1+(T_L\cdot T_R)}}.\eqno(3.22)$$
Substituting into (3.20) gives
$$ z^* = z+{{(d_R-d_L)(T_L+T_R)}\over{2\left(1+T_L\cdot T_R\right)}}.\eqno(3.23
)$$
Therefore, after obtaining $z_i=z(t+h)$ , we rearrange $z_i$ as follows before
calculating the positions predicted by the next time step. Calculate
$$\eqalign{d_{R[i]}&=\vert z_{R[i]}-z_i\vert, \quad d_{L[i]}=\vert
z_i-z_{L[i]}\vert,\cr T_{R[i]}&={z_{R[i]}-z_i\over d_{R[i]}}, \quad
T_{L[i]}={z_i-z_{L[i]}\over d_{L[i]}},\cr
\hbox{then set } \hfil z_i: &=
z_i+{{(d_{R[i]}-d_{L[i]})(T_{L[i]}+T_{R[i]})}\over{2\left(1+T_{L[i]}\cdot
T_{R[i]}\right)}}.\cr}\eqno(3.24)$$
So our numerical scheme has two steps: Given a curve discretization $\Gamma_t$,
compute a discretization of $\Gamma_{t+h}$ by using the semi-implicit algorithm
to find $V$, then redistribute the mesh points on $\Gamma_{t+h}$ and repeat the
whole process. The results of implementing this in several examples are given
below.
\bigskip
\centerline{\bf 4. Experiments }\medskip
First, to test the accuracy of our numerical method, we choose a case which we
can solve analytically, namely, the case for two concentric circles.
Take concentric circles of radii $R_1dR_1$$ and so
by using (4.5) we have
$$ t = \int_{R_1}^{R_1^0} {{R_1\ln{R_2\over R1}}\over {{1\over R_1}+{1\over R_2}
}}\>dR_1 = {1\over 2}\int_{R_1}^{R_1^0} {{ R_1^2 \ln{{A^2+R_1^2}\over
R_1^2}\sqrt{R_1^2+A^2}}\over {R_1+\sqrt{R_1^2+A^2}}}\>dR_1.\eqno(4.6)$$
Let $R_1 = Ar$ then (4.6) can be written as
$$t={A^3\over 2}\int_{R_1\over A}^{R_1^0\over A} {{r^2\sqrt{1+r^2}
\ln\left(1+{1
\over r^2}\right)}\over{r+\sqrt{1+r^2}}}\> dr.\eqno(4.7)$$ Or we may write
$r={1\over k}$ in (4.7), to obtain the curvature-time relationship
$$t= {A^3\over 2}\int_{Ak^0}^{Ak} \ln\left(1+k^2\right){\sqrt{1+k^2}\over{\left(
1+\sqrt{1+k^2}\right)k^4}}\> dk,\eqno(4.8)$$ where
$$ k^0 = {1\over R_1^0}.$$
Given curvature values $k_1$ and $k_2$ obtained from the numerical simulation of
(2.1) over time step $\Delta t$, setting $a=Ak_1$ and $b=Ak_2$,
we will compare the results with those obtained by integrating (4.8).
To calculate the integral in (4.8), we use Simpson's rule. Denote
$$ f(k) = {\sqrt{1+k^2}\over {k^4\left(1+\sqrt{1+k^2}\right)}}
\ln\left(1+k^2\right).\eqno(4.9)$$
The discretization of the integral will use step length
$h < {\Delta t\over 10}$ by setting $n=[10(b-a)/\Delta t]+1$ and
$h=(b-a)/n$. Then the integral over the interval [$a$, $b$] is
$$\int_a^b f(k)dk \simeq \sum_{i=0}^{n-1} {h\over
6}\left(f(a+ih)+4f\left(a+ih+{h\over 2}\right)+f(a+ih+h)\right).\eqno(4.10)$$
The above formula has accuracy given by
$$\eqalign{\int_\alpha^{\alpha+h} f(k)dk = &{h\over
6}\left(f(\alpha)+4f\left(\alpha+{h\over 2}\right)+f(\alpha+h)\right)\cr
&+{1\over 2880} f^{(4)}(\tilde\theta)h^5 \hbox{ for some }\ \
\tilde\theta\in(\alpha,\alpha+h).\cr}$$
Clearly, $f^{(4)}(k)$ is bounded for $k\geq 1$. Hence, $t$ can be computed to an
accuracy of order $O(h^4)$ by
$$t = {A^3\over 2}\Biggl[\sum_{i=0}^{n-1} {h\over
6}\left(f(a+ih)+4f\left(a+ih+{h\over
2}\right)+f(a+ih+h)\right)\Biggl].\eqno(4.11)$$.
We apply our algorithm, (3.16), taking concentric circles of radii $1$
and $3$ and terminate when the radius of the small circle is approximately $0.1$.
We display the curvature-time relationship of the ``exact'' solution (obtained
using (4.11)) against that produced using (4.16) for $N=8$, $16$,
$32$, $64$, $128$ as $\Delta t=0.01$, $0.002$, and $0.0004$
respectively in Fig. 4.
\smallskip
\centerline{\epsfbox{outcmp8b1.ps} \epsfbox{outcmp16b1.ps}}\smallskip
\centerline{\epsfbox{outcmp32b1.ps} \epsfbox{outcmp64b1.ps}}\smallskip
\centerline{\epsfbox{outcmp128b1.ps}}\smallskip
\centerline{Fig. 4 Accuracy Check for for Concentric Circles}
Fig.~4 reveals that our simulation becomes more accurate as the time step
decreases and as more mesh points are used, which is to be expected. To make a
detailed comparison, first we fixed a time (0.3 in the figures below) and for
each value of $\Delta t$ estimated the ``spatial error'', $se$, as the
difference between the curvature of the small circle when 128 mesh points were
used and that when $N(=8, 16, 32, 64)$ points were used.
Then assuming the spatial error satisfies the relationship
$$ se = c{\left({1\over N}\right)}^\alpha,\eqno(4.12)$$
for two choices,
$N_1$ and $N_2$, we have
$$ \alpha={{\ln\left({se_1\over se_2}\right)}\over {\ln\left({N_2\over
N_1}\right)}}.\eqno(4.13)$$
This allows us to estimate the spatial convergence
rate $\alpha$ according to (4.13):
\settabs\+\indent&(n1/n2)/$\Delta t$\qquad&2.275846725508301\quad&2.2758
46725508301\quad&2.275846725508301\quad&\cr
\vbox{\hrule width4.8in height1pt}
\+&(N1/N2)/$\Delta t$&0.01&0.002&0.0004\cr
\vbox{\hrule width4.8in height1pt}
\+&8/16&2.173 &2.167 &2.165\cr
\+&16/32&2.121 &2.119 &2.119\cr
\+&32/64&2.341 &2.340 &2.340\cr
\vbox{\hrule width4.8in height1pt}
\centerline{Table 1 Spatial Convergence Rate Estimation}
Similarly, compared with the data in the case of $\Delta
t=0.0004$, we calculated the relative errors of the other cases
and estimate the temporal
convergence rate:
$$\vbox{
\settabs\+\indent&($n/(\Delta t_1/\Delta t_2)$\qquad&\cr
\vbox{\hrule width2.5in height1pt}
\+&$N/(\Delta t_1/\Delta t_2)$&0.01/0.002\cr
\vbox{\hrule width2.5in height1pt}
\+&8&1.091\cr
\+&16&1.091\cr
\+&32&1.091\cr
\+&64&1.091\cr
\+&128&1.092\cr
\vbox{\hrule width2.5in height1pt} }$$
\centerline{Table 2 Temporal Convergence Rate Estimation}
Tables 1 and 2 suggest that for the algorithm given by
(4.16), the spatial and temperal convergence rates are of
orders around 2 and 1, respectively.
Now we report the results of implementing the algorithm with a variety of
initial curves $\Gamma_0$. Since a component of the interface accelerates as it
shrinks to a point, we make a dynamical choice of time step $\Delta t$ according
to the maximum interfacial velocity
$$\Delta t = \Delta t_0 *\min\left\{{1\over V_{max}}, 1\right\},$$
where $\Delta
t_0$ is the initial setting of the time step and
$V_{max}= \max(|V_i|)$, for $i=1$, $\dots$, $MN$ at each time $t$.
We first consider the relaxation of a single body as it becomes circular.
The initial shape is given parametrically by $x=(2+0.5\sin3\theta)\cos\theta$,
$y=(2+0.5\sin3\theta)\sin\theta$. $N=132$ mesh points are distributed around the
boundary. The interface shape and curvature as a function of arc length at four
subsequent
times are displayed in Fig.~5, where the arc length is measured from the point
where the interface intersects the positive $x$ axis.
\settabs 2\columns
\+\epsfbox{outinrose132.ps} &
\epsfbox{outkrose132c0.ps}\cr \medskip
\+\epsfbox{outfrose132c1.ps} &
\epsfbox{outkrose132c1.ps}\cr \medskip
\+\epsfbox{outfrose132c2.ps} &
\epsfbox{outkrose132c2.ps}\cr \medskip
\+\epsfbox{outfrose132c3.ps} &
\epsfbox{outkrose132c3.ps}\cr\smallskip
\centerline{Fig. 5 Interface Evolution for Rose Curve}
Notice that the initially nonconvex body becomes convex fairly rapidly and then
remains convex with the interfacial velocity decreasing as the boundary becomes
closer to being circular. It was proved in [5] that a curve sufficiently
close to being circular converges to a circle. It would be natural to conjecture
that convexity is preserved by this evolution, however, the next example
suggests this is not the case.
In [16] Mayer proved that convexity could be lost when the evolution was
governed by the one phase flow. The initial curve in [16] is basically that of
our next experiment.
We consider a shape given by a long thin tube with semicircular end caps. The
straight part has length of 16, while the radius of the circular part is 0.125.
We take $N=120$ mesh points around the boundary. Shown in Fig. 6 are the
interface shape and curvature as a function of arc length at five subsequent
times, where the arc length is measured from the point (8. 0163, -0.1239) (the
body center is located at (0, 0)).
\+\epsfbox{outgintube120.ps}&
\epsfbox{outktube120c0.ps}\cr \smallskip
\+\epsfbox{outftube120c1.ps}&
\epsfbox{outktube120c1.ps}\cr \smallskip
\+\epsfbox{outftube120c2.ps} &
\epsfbox{outktube120c2.ps}\cr \smallskip
\+\epsfbox{outftube120c3.ps} &
\epsfbox{outktube120c3.ps}\cr \smallskip
\+\epsfbox{outftube120c4.ps}&
\epsfbox{outktube120c4.ps}\cr \smallskip
\centerline{Fig. 6 Interface Evolution for a Thin Tube}\medskip
Convexity is lost immediately as the ends become bulbous but eventually
convexity is regained and the curve tends to a circle. A careful analysis shows
that even though convexity is lost, the initial motion widens the tube in the
center but it widens more quickly near the endcaps. This raises the question of
whether or not a curve can ever pinch off, increasing the number of components.
Our next simulation considers an initial curve which is a dumbbell formed by
taking two ellipses joined by a narrow tube. If the ellipses are small then the
curve does not form self-intersections. However, for large ellipses, such as
shown in Figure 7, self intersections are formed. Of course, the model ceases
to be valid as soon as singularities form but our simulation suggests that the
model allows pinching-off. Even though geometric singularities occur, it can
be shown that no singularities form in the integral when one part of the
evolving curve becomes tangent to another part. However, for computational
ease $10^{-8}$ is added to the terms $|a_1| - |a_4|$ in (3.7), thereby avoiding
possible problems arising due to discretization. Figure 7 superimposes the
initial curve and its state at two later times (smooth curves). It is the
nonlocal nature of the problem that allows the curvature to drive a curve to
self-intersect. It is well known that the explicit motion by curvature in the
plane does not allow self intersections and also preserves convexity. For the
case at hand it is interesting to observe the undulations in the connecting
tube prior to the occurance of self intersection. One is led to speculate
whether simultaneous multiple self intersections could occur in the connecting
tube.
\epsfbox{d10d2.ps}
\centerline{Fig. 7 Pinch-Off}
\epsfbox{Fig82.ps}
\centerline{Fig. 8 Coalescence}
The reverse change in topology, coalescence of two particles, also seems
possible as shown by the simulation starting with two equal ellipses having
vertices close together with colinear minor axes. The symmetry ensures that
neither particle gains area at the expense of the other, so one expects the
particles to evolve to become equal circles, which is a stationary
configuration. However, the particles are so close that these circles could not
exist without overlapping unless their centers are further apart than the
centers of the original ellipses. In fact the centers do move apart but not far
enough (see Figure 8).
Figures 9 and 10 show the morphological evolution of four particles in an
initial configuration considered in [21], but here we use our algorithm for the
two phase flows, while in [21] the evolution was according to single phase
diffusion. The initial radii, from left to right, are 1, 0.9, 0.8 and 0.8. The
centers of the particles are located at (-2.5, 0), (0, 0), (2, 1.2) and (2,
-1.2), respectively.
$N=64$ mesh points are taken equally spaced around the boundary of each
particle. Originally, the center of the overall mass of the particles is at
(0.019, 0). In Fig. 9, we use a small diamond to indicate the center of mass
as the particles evolve.
\+\epsfbox{outgin4c64_1.ps}&
\epsfbox{outf4c64_1c1.ps}\cr
\+\epsfbox{outf4c64_1c2.ps}&
\epsfbox{outf4c64_1c3.ps}\cr
\+\epsfbox{outf4c64_1c4.ps}&
\epsfbox{outf4c64_1c5.ps}\cr
\centerline{Fig. 9 Evolution of Four Particles}\smallskip
Notice that the leftmost (largest) particle grows slowly with time, maintaining
approximately circular symmetry, while the middle particle expands more quickly
accompanied by the shrinking of the smaller particles. This is because there is
a strong local diffusional interaction between the middle particle and nearby
small particles. With a dynamic change in the time step, our algorithm is able
to follow to the disappearance of the smallest particles. The removal of
vanishing particles is done by deleting particles whose
average curvature is greater than 300. It is interesting to note the motion
of the center of mass during the evolution. For the single phase problem, one
can show that the centroid remains fixed.
Fig.~10 gives a clearer description of shape distortions before the smaller
particles disappear. It is clear that the middle particle recedes slightly on
its left side, material being transported to the leftmost particle, which has a
corresponding advance on its right side. By the time the smallest particles
disappear, the middle particle has become largest and so this is the only
survivor (see Figure 9).
\bigskip
\centerline{\bf Appendix }\medskip
This contains two programs. The first, called {\tt calcinit.c} constructs
initial curves and after compiling, execution is done by typing {\tt calcinit
file1 file2}, where the user gives names to {\tt file1} and {\tt file2}. The
second, called {\tt hs.c} executes our algorithm. It is executed with the
command {\tt hs file1 file3 file4 file5 file6}. In constructing the curve the
user is asked to choose the number of curves, their type from a menu (e.g.
ellipse, tube, etc.), and then dimensions (e.g. semiminor and semimajor axes)
and center locations are required.
When executing {\tt hs}, the desired values of $N, \Delta t$ and the number of
time steps are requested.
The output is placed in {\tt file4} and can be viewed using gnuplot or some other
data plotting program.
\epsfbox{outf4c64_1cal.ps}
\centerline{Fig. 10 Evolution of Four Particles}
\bigskip
\centerline{\bf References }\medskip
\item{[1]} Alikakos, N.D., Bates, P.W. and Chen, X., ``Convergence of the
Cahn-Hilliard equation to the Hele-Shaw Model,'' {\it Arch. Rat. Mech. Anal.}
{\bf 128} (1994), 165--205.
\item{[2]} Almgren, R., Singularity formation in Hele-Shaw bubbles, Preprint
May 1995.
\item{[3]} Baker, G.R., Meiron, D.I., and Orszag, S.A., ``Boundary integral
methods for axisymmetric and three-dimensional Rayleigh-Taylor instability
problems,'' {\it Physica D} {\bf 12} (1984), 19--31.
\item{[4]} Baker, G.R., ``Generalized vortex methods for free-surface flows,''
in Meyer, R.E. (ed.), Waves on fluid interfaces, Academic Press, New York,
(1983), 53--81.
\item{[5]} Chen, X., ``The Hele-Shaw problem and area-preserving, curve
shortening motions,'' {\it Arch. Rat. Mech. Anal.} {\bf 123} (1993), 117--151.
\item{[6]} Chen, X., Hong, J., Yi, F., ``Existence, uniqueness, and regularity
of classical solutions of the Mullins-Sekerka problem,'' preprint.
\item{[7]} Chen, X., Hou, T.Y., and Zhu, J., Preprint, 1995.
\item{[8]} Constantin, P., and Pugh, M., ``Global solutions for small data to
the Hele-Shaw problem,'' {\it Nonlinearity} {\bf 6} (1993), 393--415.
\item{[9]} DeGregoria, A.J., and Schwatz, L.W., ``A boundary-integral method
for two-phase displacement in Hele-Shaw cells,'' {\it J. Fluid Mech.} {\bf 164}
(1986), 383--400.
\item{[10]} Deng, X., ``A numerical analysis approach for the Hele-Shaw
problem,'' M.S. Thesis, Brigham Young University, April 1994.
\item{[11]} Folland, G. B., Introduction to partial differential equations,
Princeton University Press, Princeton, NJ, 1976.
\item{[12]} Hou, T.Y., Lowengrub, J.S., and Shelley, M.J., Removing the
stiffness from interfacial flows with surface tension, {\it J. Comp. Phys.}
{\bf 114} (1994), 312--338.
\item{[13]} Kellogg, O.D., Foundations of potential theory, Springer-Verlag,
New York, 1929.
\item{[14]} Kessler, D.A., Kopik, J., and Levine, H., ``Numerical simulation of
two- dimensional snowflake growth,'' {\it Phys. Rev. A} {\bf 30} (1984),
2820--2823.
\item{[15]} Kessler, D.A., and Levine, H., ``Theory of the Saffman-Taylor
`finger' pattern. I \& II,'' {\it Phys. Rev. A} {\bf 33} (1986), 2621--2639.
\item{[16]} Mayer, U.F., ``One-sided Mullins-Sekerka flow does not preserve
convexity,'' {\it Electr. J. Diff. Eqns.} {\bf 1993} (1993), No. 8, 1--7.
\item{[17]} McFadden, G.B., Voorhees, P.W., Boisvert, R.F., and Meiron, D.I.,
``A boundary integral method for the simulation of two-dimensional particle
coarsening,'' {\it J. Sci. Comp.} {\bf 1} (1986), 117--143.
\item{[18]} Menikoff, R., and Zemach, C., ``Methods for numerical conformal
mapping,'' {\it J. Comp. Phys.} {\bf 36} (1980), 366--410.
\item{[19]} Mullins, W.W. and Sekerka, R.F., ``Morphological stability of a
particle growing by diffusion and heat flow,'' {\it J. Appl. Phys.} {\bf 34}
(1963), 323--329.
\item{[20]} Pego, R., ``Front migration in the nonlinear Cahn-Hilliard
equation,'' {\it Proc. Roy. Soc. London A} {\bf 422} (1989), 261--278.
\item{[21]} Voorhees, P.W., McFadden, G.B., and Boisvert, R.F., ``Numerical
simulation of morphological development during Ostwald ripening,'' {\it Acta
Meta.} {\bf 36} (1988), 207--222.
\item{[22]} Voorhees, P.W., ``Ostwald ripening of two-phase mixtures,'' {\it
Annu. Rev. Mater. Sci.} {\bf 22} (1992), 197--215.
\bigskip
\vbox{Peter W. Bates \hfil\break
Mathematics, Brigham Young University \hfil\break
Provo, UT 84602 \hfil\break
E-mail address: peter@math.byu.edu}
\bigskip
\vbox{Xinfu Chen \hfil\break
Mathematics, University of Pittsburgh \hfil\break
Pittsburgh, PA 15260 \hfil\break
E-mail address: xinfu$+$@pitt.edu}
\bigskip
\vbox{Xinyu Deng \hfil\break
Mathematics, Brigham Young University \hfil\break
Provo, UT 84602 \hfil\break
E-mail address: cindy@math.byu.edu}
\vfil\eject
\centerline{\bf Erratum} \medskip
{\bf September 26, 1995}.
I did not write Theorem 2.3 correctly so I enclose the following new
version with the necessary changes in the first few lines of the
proof.\hfil\break
Best, Peter Bates.
\bigskip
\noindent {\bf Theorem 2.3} {\sl If $\Gamma_{0,T}\equiv\cup_{0\leq
t\leq T}\left({\Gamma_t\times\{t\}}\right)$ is a continuous family of
$C^3$ curves which satisfy (2.3) and (2.4) with $\Gamma=\Gamma_t$ for
$g=V(x,t)$, the normal velocity of $\Gamma_t$, and some
$c=c(t)$, where $f=K(x,t)$ is the curvature of $\Gamma_t$ at $x$, then
$\Gamma_{0,T}$ is the interface associated with the solution to (2.1).
Conversely, if $(u,\Gamma_{0,T})$ is a solution to (2.1) then (2.3)--(2.4)
hold for each $t\in [0,T]$ with $\Gamma=\Gamma_t, g=V$ and $f=K$.}
\noindent {\bf Proof:}
Let $(\Gamma_t, g, c)$ be a solution to (2.3)--(2.4) at each time
$t\in [0,T]$ with $f=K(\cdot, t)$ and suppose that $g = V$, the normal
velocity of $\Gamma_t$.
Then Lemma 2.1 shows that defining $u(\cdot, t)$ on
${\Bbb R}^2\backslash\Gamma_t$ by
$$u(x,t)={1\over 2\pi} \int_{\Gamma_t} \ln\vert x-y\vert g(y,t)ds_y+c(t)$$
gives a solution to (2.1) with $-\left[{{\partial c\over
\partial n}}\right]_{\Gamma_t}$ being given by $V$.
\bye