A variational method for multiphase area-preserving interface motions

We first fix the meaning of the terms “area” and “length” used in the sequel. Working in $m$-dimensional space, the word “area” shall mean the $m$-dimensional Lebesgue measure of the region corresponding to a phase, i.e., in 3 dimensions it regards the volume. The word “length”, on the other hand, shall always refer to the $(m-1)$-dimensional Lebesgue measure of the boundary of a phase region, thus in 3 dimensions corresponding to the surface area of the interface.

Mean curvature flow is related to systems whose energy depends on the length of their surface, e.g., soap films. We will explain the basic equations for curves in the two-dimensional plane, since the derivations in higher dimensions are similar but lack transparency. Accordingly, let us consider a smooth Jordan curve $\gamma$ contained in a subregion $\Omega$ as in figure 1. Let the curve be parametrized:

and in such a way that the enclosed region $P$ is on the left side of the curve. Then the length of the curve is given by

The gradient flow of the above energy can be found from its first variation. That is, for any smooth closed curve $\varphi(s),s\in[a,b]$ we compute

where $\kappa$ is the curvature, $\boldsymbol{n}$ is the unit outer normal to the curve $\gamma$ at a given point and where we integrate with respect to arc length $l$:

In the general $m$-dimensional case one obtains the same formula as (1), where $\kappa$ means the trace of second fundamental form divided by $m$. Hence, the fastest shortening of the curve occurs for the flow with normal velocity equal to minus (a multiple of) the mean curvature.

Next we consider the curve-shortening flow under the constraint of area preservation. The area functional for region $P$ reads

and its first variation is

Analogously, in any dimension the variation turns out to be the normal vector at each point of the hypersurface.

Following the construction in [24], we introduce a time-dependent Lagrange multiplier $\lambda(t)$ for the area constraint and express the velocity of the area-constrained mean curvature flow by

A precise expression for the the multiplier $\lambda$ can be obtained in the following standard way. From the fact that the area is preserved one has

Hence it follows that the integral of $-\kappa+\lambda$ over the curve $\gamma$ vanishes at each time. This yields

That is, the Lagrange multiplier expresses the average mean curvature along the interface.

We briefly derive the velocity of interfaces moving by area-preserving mean curvature flow in the multiphase two-dimensional setting. We assume that the number of phases $k$ is finite and that the interfaces between different pairs of phases form a finite collection of arcs. The boundary $\gamma_{i}$ of phase region $P_{i}$ can then be written as the union of the interfaces from all other phase regions:

Here $\gamma_{ij}$ denotes the interface between phases $P_{i}$ and $P_{j}$, and the index $l$ expresses the fact that each interface may have several disjoint parts. In the following, however, we omit the decomposition using index $l$ since it has no influence on the computations.

Each interface $\gamma_{ij}^{l}$ is considered as an oriented curve, which has the region $P_{i}$ on its left side. For any point interior to the interface $\gamma_{ij}$ we define the normal $\boldsymbol{n}$ as the unit vector pointing into the phase with larger index, i.e., $\boldsymbol{n}$ is the outer normal to $P_{i}$ if $i<j$. We remark that some of the curves $\gamma_{ij}$ may be empty.

The energy of the system, considering arbitrary surface tension $\tau_{ij}$ for $\gamma_{ij}$, is

This value is to be minimized under the condition of constant areas, that is,

Introducing Lagrange multipliers $\lambda_{i}$, $i=1,\dots,k-1$, for each of the constraints, the constrained variation reads

where $\varphi^{ij},\,i,j=1,\dots,k$, denotes a smooth perturbation within $\gamma_{ij}$.

From this it follows that the magnitude of normal velocity $v^{ij}$ for interface $\gamma_{ij}$ ($i<j$) in the direction of $\boldsymbol{n}$ will be

Here $\delta_{jk}$ denotes the Kronecker delta, which arises from the redundance of the area constraint for phase $P_{k}$ (it follows from the constraints on the other phases and the fixation of the domain).

Explicit representations of the Lagrange multipliers can be obtained as follows. Let us denote by $L_{ij}$ the length of interface $\gamma_{ij}$, by $L_{i}$ the length of the boundary of phase region $P_{i}$, and by $\kappa_{ij}$ the average mean curvature through interface $\gamma_{ij}$ times the tension $\tau_{ij}$:

Then the preservation of phase area $P_{i}$ gives the condition

Since the above holds for $i=1,\dots,k-1$, we obtain a system of linear equations for $\lambda_{1},\dots,\lambda_{k-1}$:

For a given configuration, the solution to this system gives the Lagrange multipliers. For example, in the case of three phases, as in figure 2, one has the velocities:

Here $\alpha=L_{13}L_{23}+L_{23}L_{12}+L_{13}L_{12}$ and $\kappa_{a}^{P_{i}}$, $i=1,2,3$, represent the average mean curvatures along the whole boundary of each phase, weighted by the surface tension.

When phases $P_{1}$ and $P_{2}$ are separated, i.e., when $L_{12}=0$, the above formulae reduce to the form of the two-phase flow from section 2.1.1. The extension of the above calculations to higher space dimensions is natural, as we use only the notions of interfaces and their oriented normals.

We describe the BMO algorithm for the case when only two phases are present. This algorithm works in any space dimension. Given an initial interface $\gamma$, we take $P$ to be the region enclosed by this interface (and possibly by the boundary of the region $\Omega$ where the motion is considered) and define its characteristic function $\chi$ as

The new interface after a time $\Delta t$ will be the boundary of the $\tfrac{1}{2}$-level set of the solution to the heat equation at time $\Delta t$ with initial datum $\chi$. The two-phase algorithm thus reads as follows:

1. 1.

   Given a region $P$, set $\chi$ to be its characteristic function.

2. 2.

   Solve the heat equation with initial condition $\chi$:

   	$\displaystyle u_{t}(t,x)$	$\displaystyle=$	$\displaystyle\Delta u(t,x)\qquad\;\,\text{for}\;(t,x)\in(0,\Delta t]\times\Omega,$
   	$\displaystyle\frac{\partial u}{\partial n}(t,x)$	$\displaystyle=$	$\displaystyle 0\qquad\qquad\quad\;\,\text{on}\;(0,\Delta t]\times\partial\Omega,$
   	$\displaystyle u(0,x)$	$\displaystyle=$	$\displaystyle\chi(x)\qquad\qquad\text{in}\;\Omega.$

3. 3.

   Update $\chi$ as the $\tfrac{1}{2}$-level set of $u(\Delta t,\cdot)$:

   $$\chi(x)=\left\{\begin{array}[]{l}1\qquad\text{if}\;\;u(\Delta t,x)>\tfrac{1}{2},\\ 0\qquad\text{if}\;\;u(\Delta t,x)\leq\tfrac{1}{2}.\end{array}\right.$$

   The evolved interface is now the boundary of the set $\{x\in\Omega;\;\chi(x)=1\}$.

4. 4.

   Go back to step 2 to proceed with the computation for the next time step.

It has been rigorously shown, in a general framework including topological changes, that this algorithm converges to motion by mean curvature as $\Delta t\to 0$ [25, 26, 27].

Here we remark that the Neumann boundary condition in the diffusion step guarantees that the interface will touch the boundary of $\Omega$ with right angle. Other boundary conditions, such as Dirichlet conditions pinning the interface at the boundary, may be used according to necessity.

We next address the case of multiple phases. The idea of sharpening separately diffused characteristic functions (one for each phase) was introduced in [22]. The algorithm is as follows.

1. 1.

   Given regions $P_{i},i=1,\dots,k$, set $\chi_{i}$ to be the characteristic function of $P_{i}$.

2. 2.

   For $i=1,\dots,k$, obtain $u_{i}(\Delta t,x)$ by solving the heat equation with initial condition $\chi_{i}$ up to time $\Delta t$.

3. 3.

   Update $\chi_{j}$ as the characteristic function of the set where $u_{j}$ has the largest value amongst all solutions $u_{i}$:

   $$\chi_{j}(x)=\left\{\begin{array}[]{l}1\qquad\text{if}\;\;u_{j}(\Delta t,x)\geq u_{i}(\Delta t,x)\;\;\forall i\neq j,\\ 0\qquad\text{otherwise}.\end{array}\right.$$

   The new interfaces are the boundaries of the sets $\{x\in\Omega;\;\chi_{i}(x)=1\}$.

4. 4.

   Go back to step 2 to proceed with the computation for the next time step.

The above algorithm can be reformulated to obtain an equivalent algorithm using a single vector-valued heat equation. This is essential for implementing the area constraint and for dealing with more general motions.

We prepare $k$ reference unit vectors $\boldsymbol{p}_{i},i=1,\dots,k$, of dimension $k-1$, each corresponding to a phase $P_{i}$. They are defined as the vectors pointing from the centroid of a standard $k$-simplex to its vertices (cf. figure A.12 in the appendix). Hence, there are infinitely many possible $k$-tuples but the relative distributions of the vectors are identical. See A for a simple way to construct these vectors.

Using the vectors $\boldsymbol{p}_{i}$, the multiphase algorithm can be written as follows:

1. 1.

   Given regions $P_{i},i=1,\dots,k$, set

   $$\boldsymbol{u}_{0}(x)=\boldsymbol{p}_{i}\qquad\text{for}\;\;x\in P_{i}.$$

2. 2.

   Solve the vector-valued heat equation with initial condition $\boldsymbol{u}_{0}$:

   	$\displaystyle\boldsymbol{u}_{t}(t,x)$	$\displaystyle=$	$\displaystyle\Delta\boldsymbol{u}(t,x)\qquad\;\;\text{for}\;(t,x)\in(0,\Delta t]\times\Omega,$		(3)
   	$\displaystyle\frac{\partial\boldsymbol{u}}{\partial n}(t,x)$	$\displaystyle=$	$\displaystyle 0\qquad\qquad\quad\;\;\;\text{on}\;(0,\Delta t]\times\partial\Omega,$
   	$\displaystyle\boldsymbol{u}(0,x)$	$\displaystyle=$	$\displaystyle\boldsymbol{u}_{0}(x)\qquad\qquad\text{in}\;\Omega.$

3. 3.

   Update $\boldsymbol{u}_{0}$ by identifying the reference vector which is closest to the solution $\boldsymbol{u}(\Delta t,x)$:

   $$\boldsymbol{u}_{0}(x)=\boldsymbol{p}_{j},\qquad\text{where}\quad\boldsymbol{p}_{j}\cdot\boldsymbol{u}(\Delta t,x)=\max_{i=1,\dots,k}\boldsymbol{p}_{i}\cdot\boldsymbol{u}(\Delta t,x).$$

   (4)

   This redistribution of reference vectors determines the approximate new phase regions after time $\Delta t$.

4. 4.

   Go back to step 2 to proceed with the computation for the next time step.

The equivalence of our algorithm with the original BMO can be shown by considering the functions

Since $\boldsymbol{u}$ is the solution to the linear vector-valued heat equation, we have

Moreover, since from the construction of $\boldsymbol{p}_{i}$ it holds that $\boldsymbol{p}_{i}\cdot\boldsymbol{p}_{j}=1/(1-k)$ for $i\neq j$, we can readily check the identity

where $\chi_{i}$ is the characteristic function of $i$-th phase region. It follows that $w_{i}$ is identical to the solution $u_{i}$ of the scalar heat equation from the original algorithm for each $i=1,\dots,k$. From the definition of $w_{i}$ it is immediate that

which proves the equivalence of both algorithms.

The paper [23] presents an area-preserving BMO method for 2 phases. The authors adjust the height at which the thresholding occurs in such a way that the resulting area of the level set is preserved. Such a level set is guaranteed to exist by the maximum principle and one can compute that the thresholding height has to be changed from $\tfrac{1}{2}$ to

where $\kappa_{a}$ is the average mean curvature along the interface.

However, when three or more phases are present, each interface has a different average mean curvature and the thresholding heights become different. One could try to use a different thresholding height for each phase but then the global interaction would be ignored and the phases may overlap or create vacuums, especially when they are initially touching. Therefore, we suggest a different approach incorporating the area constraint into the diffusion process. In this method, a heat equation with a nonlinear source term expressing the area preservation of level sets is solved. Subsequently, the solution is sharpened at a fixed height.

The mentioned nonlinear heat equation corresponds to the gradient flow of the functional

in the constrained set

where $A_{i}$ denotes the given area of phase $P_{i}$. For simplicity, we consider only the case when one type of interface divides one type of phase into multiple regions, like in a soap froth. When multiple phases and interfaces with nonuniform properties are present, a different form of the energy $J$ has to be adopted [28].

Let us consider the two-dimensional two-phase case as in figure 1 to understand the meaning of this gradient flow. Our argument is formal and we simplify the exposition by discretizing time and adopting the numerical approach from the next section, which is based on the method of Rothe [29].

Let the curve $\gamma$ denote the $\tfrac{1}{2}$-level set of a scalar function $u$ and call the region enclosed by this curve $P$. We assume that an initial function $u_{*}$ is given and we search for the minimizer of the integral

under the constraints

Here $h$ denotes the length of the discrete time step and $A$ is the required area.

Note that by perturbing function $u$ in regions away from the $\frac{1}{2}$-level set, one can readily deduce that the minimizer $u$ satisfies, in a weak sense , the following

Now consider a perturbation of $u$ of the form $u+\delta u$, which is allowed to affect the $\frac{1}{2}$-level set. The corresponding change in the level curve can be written in the form $\alpha(\gamma+\delta\gamma)$, $\delta\gamma=(\delta\gamma_{1},\delta\gamma_{2})$, where $\alpha$ is a constant depending on $\delta\gamma$ and whose role is to adjust the area to the correct value. Because of constraints (6) and (7), we cannot choose $\delta u$ arbitrarily. However, we can select an arbitrary $\delta\gamma$, find an appropriate constant $\alpha$, and use the corresponding $\delta u$.

We compute the variation of functional $J$ remembering that $u$ may not be smooth across the interface $\gamma$:

Here $w_{P}$ denotes the value of a function $w$ taken as a limit from inside of $P$. The symbol $w_{\Omega\setminus\bar{P}}$ has the analogous meaning.

Due to (8), on the interface we have

We next to rewrite this identity in terms of the arbitrary perturbation $\delta\gamma$. To this end, we use the fact that the phase area is preserved to obtain

which yields

Here, the symbol $\simeq$ means that the equation holds up to second order in terms of $\delta\gamma$ or $\delta u$. From this we can compute the value of $\alpha-1$ which will be needed later:

Condition (6) means

Using a Taylor expansion and (9) we obtain that, on $\gamma$, one has

Expressing $\delta u$ from (2.2.3) as

we obtain

Noting that

we arrive at the interface condition

Similar calculation can be carried out in the multiphase setting, see C for more details.

In view of the derived condition on $\gamma$ and the results regarding two-phase free boundary problems in [30], we can reformulate the variational problem in the following way: Find a minimizer $u_{\lambda}\in H^{1}(\Omega)$ of the functional

According to the results in [30], we can expect that this minimization problem is in a sense equivalent to looking for a weak solution of the problem

where $\mu$ is a Radon measure given by

for a suitable space-independent function $\lambda$. Here the symbol ${\mathcal{H}}^{1}$ means the one-dimensional Hausdorff measure. This type of problem is known as the two-phase parabolic free boundary problem [31, 32, 33].

Having deduced the partial differential equation corresponding to the constrained gradient flow, let us return to the considerations concerning the BMO method. Formal calculation shows that the BMO method using the parabolic equation (12) instead of the heat equation evolves interfaces with normal velocity equal to minus mean curvature, plus a space-independent term:

Thus, if $\lambda$ is chosen appropriately, the area-preserving mean curvature flow is realized. Since the derivation of the above relation is rather technical, we present it in B.

Minimization of (11) is much easier to implement than the constrained minimization (5) or the partial differential equation (12). In each time step of the BMO algorithm, if for some $\lambda$ we obtain that the area of $\{x;\;u_{\lambda}(x)>\tfrac{1}{2}\}$ is equal to the given value $A$, then $u_{\lambda}$ becomes a solution to our original problem. However, it is not clear how to calculate such a $\lambda$ or even if such a function exists. For the two-phase case we see that $\lambda$ should be half of the average mean curvature, but determining the Lagrange multiplier for the multiphase case is complicated. Moreover, one of the advantages of the BMO approach is that it does not require the computation of curvature values, so we want to avoid the direct calculation of $\lambda$.

Therefore, based on the results in [34] we approach (5) by using the method of penalization and consider minimization of the functional

where $\varepsilon$ is a small positive number and the penalty function is defined by

Penalization techniques of this type for stationary problems were analyzed in [34, 35, 36] and other works. The general conclusion of these treatments is that the penalized solutions converge as $\varepsilon\to 0$ to a solution of the corresponding PDE or of the original constrained minimization problem. Moreover, it was discovered that, in order to obtain the solution of the original problem, it is not necessary to carry out the limit $\varepsilon\to 0$, since the exact minimizer is obtained for some sufficiently small but positive value of $\varepsilon$. Although the analysis of the time-dependent problem has not yet been addressed, this observation is essential to justify the robustness of the penalty method, at least in the stationary case.

In closing this section we note that, except for the paper [36], the constrained multiphase (vector-type) problem is almost unexplored. We provide here a formal analysis of the convergence of the BMO algorithm for the constrained multiphase problem. Although this analysis is an important part of the present work, we defer it to C in order to make the main text less tortuous.

This section examines the behavior of our method in comparison to the standard BMO. We refer to the standard method as BMO, and to our own algorithm (which utilizes the thresholding from section 3.2) as BMO^∗.

We examine the error of our method when applied to a simple test problem. By symmetry, a circle of initial radius $r_{0}$ which is evolving by mean curvature flow remains a circle whose radius $r(t)$ satisfies the following ordinary differential equation: