A structure-preserving finite element method for the multi-phase Mullins-Sekerka problem with triple junctions

We consider a sharp interface formulation for the multi-phase Mullins-Sekerka flow. The flow is characterized by a network of curves evolving such that the total surface energy of the curves is reduced, while the areas of the enclosed phases are conserved. Making use of a variational formulation, we introduce a fully discrete finite element method. Our discretization features a parametric approximation of the moving interfaces that is independent of the discretization used for the equations in the bulk. The scheme can be shown to be unconditionally stable and to satisfy an exact volume conservation property. Moreover, an inherent tangential velocity for the vertices on the discrete curves leads to asymptotically equidistributed vertices, meaning no remeshing is necessary in practice. Several numerical examples, including a convergence experiment for the three-phase Mullins-Sekerka flow, demonstrate the capabilities of the introduced method.


Introduction
In this paper, we consider the problem of networks of curves moving under the multi-phase Mullins-Sekerka flow, see, e.g., [14].These networks feature triple junctions, at which certain balance laws need to hold.The network depicted in Figure 1 consists of three time dependent curves Γ i (t), i = 1, 2, 3, that meet at two triple junction points T 1 (t) and T 2 (t).We assume that the curve network Γ(t) := ∪ 3 i=1 Γ i (t) lies in a domain Ω ⊂ R 2 for all times t ∈ [0, T ], and it partitions Ω into the three subdomains Ω j (t), j = 1, 2, 3.The three domains correspond to different phases in the multi-component system.The evolution of the interfaces Γ 1 (t), Γ 2 (t), and Γ 3 (t) is driven by diffusion.As in [14], given a time T > 0 and the hyperplane T Σ := {u ∈ R 3 | 3 j=1 u j = 0}, we introduce a vector of chemical potentials w : (0, T ) × Ω → T Σ which fulfills the quasi-static diffusion equation for j = 1, 2, 3 and t ∈ (0, T ), ∆w = 0 in Ω j (t) (1.1a) together with ∂ ⃗ ν Ω w = 0 on ∂Ω, (1.1b) where ⃗ ν Ω denotes the outer unit normal vector to ∂Ω.To close the system, we need boundary conditions on Γ(t) and on T 1 (t) and T 2 (t).These boundary conditions are given by a Stefan-type kinetic condition and the Gibbs-Thomson law on the moving interfaces, and Young's law at the triple junctions, see [14].The kinetic condition reads [∇w]⃗ ν = −V [χ] on Γ(t), (1.1c) where χ = (χ 1 , χ 2 , χ 3 ) T denotes the vector which consists of the characteristic functions χ j = X Ω j (t) of Ω j (t), ⃗ ν is the unit normal vector on Γ(t), and V is the velocity of Γ(t) in the direction of ⃗ ν.We write ⃗ ν = 3 i=1 ⃗ ν i X Γ i (t) and use this convention for quantities defined on Γ(t) throughout the paper.The orientation of the three normal vectors is shown in Figure 1.In addition, the quantity [q] represents the jump of q across Γ(t) in the direction of ⃗ ν defined by [q](⃗ x) := lim ε↘0 {q(⃗ x + ε⃗ ν) − q(⃗ x − ε⃗ ν)}.Furthermore, the Gibbs-Thomson equations can be written as w • [χ] = σκ on Γ(t), (1.1d) where κ denotes the curvature of Γ(t) (well-defined on the interiors of Γ i (t)), and σ = 3 i=1 σ i X Γ i (t) is a surface tension coefficient on Γ(t).Our sign convention is such that unit circles have curvature κ = −1, which is different to the one used in [14].Finally, denoting by ⃗ µ i , the outer unit co-normal to Γ i (t), we further require Young's law, which is a balance of force condition at the triple junction as follows: (1.1e) In order to be able to fulfill this condition, we require σ 1 ≤ σ 2 + σ 3 , σ 2 ≤ σ 1 + σ 3 and σ 3 ≤ σ 1 +σ 2 .It can be shown that solutions to (1.1) reduce the weighted length 3 i=1 σ i |Γ i (t) of the curve network, while conserving the areas of the subdomains Ω 1 (t) and Ω 2 (t) (and hence trivially also of Ω 3 (t)), see Section 2 for the precise details.Our aim in this paper is to introduce a numerical method that preserves these two properties on the discrete level.
Prescribing an initial condition Γ(0) = Γ 0 for the interface, we altogether obtain the following system: (1.2) The system (1.2) at present is written for the setup from Figure 1, i.e. a network of three curves, meeting at two triple junctions and partitioning Ω into three phases.We will later generalize this to an arbitrary network of curves.The simplest case is given by a single closed curve that partitions the domain into two phases.Then we obtain the classical two phase Mullins-Sekerka problem, see [18] and the references given below.Indeed, let (w, {Γ(t)} 0≤t≤T ), with w = (w 1 , w 2 ) T , be a solution to the corresponding problem (1.2) with σ = 1, and let ⃗ ν point into Ω 2 (t), the interior domain of Γ(t).Then we have [χ] = (−1, 1) T on Γ(t), and it holds that w = w 2 − w 1 is a solution to the system The multi-phase Mullins-Sekerka problem (1.2) arises naturally as the sharp interface limit of a nondegenerate multi-component Cahn-Hilliard equation.Let Σ := {u ∈ R 3 | 3 j=1 u j = 1} and let T Σ be the family of all functions w : Ω → R 3 such that Im (w) ⊂ T Σ.Let ψ : R 3 → R be a potential whose restriction to Σ has exactly three distinct and strict global minima, say p i ∈ Σ, i = 1, 2, 3, with ψ(p 1 ) = ψ(p 2 ) = ψ(p 3 ).Let F : R 3 → R 3 be the projection of ∇ψ onto T Σ.According to [14, §2], the system (1.2) is derived as the limit with ε → 0 of a chemical system consisting of three species governed by the vector-valued Cahn-Hilliard equation whose form reads as: where u 0 : Ω → Σ denotes the initial distribution of each component and u : Ω × [0, T ] → Σ and w : Ω × [0, T ] → T Σ indicate the concentration and the chemical potential of each component in time, respectively.A distributional solution concept to (1.2) was proposed, and its existence was established via an implicit time discretization and under the assumption that no interfacial energy is lost in the limit in the time discretization (see [14,Definition 4.1,Theorem 5.8]).See [28] for a related work which treated the case without triple junctions and with a driving force.
Compared to the multi-phase Mullins-Sekerka problem, the binary case, namely the twophase case, has been well studied so far.For classical solutions, Chen et al. [18] showed the existence of a classical solution to the Mullins-Sekerka problem local-in-time in the twodimensional case, whereas Escher and Simonett [22] gave a similar result in the general dimensional case.When it comes to the notion of weak solutions, Luckhaus and Sturzenhecker [33] established the existence of weak solutions to (1.2) in a distributional sense.Therein, the weak solution was obtained as a limit of a sequence of time discrete approximate solutions under the no mass loss assumption.The time implicit scheme is the basis of the approach in [14].After that, Röger [40] removed the technical assumption of no mass loss in the case when the Dirichlet-Neumann boundary condition is imposed by using geometric measure theory.Recently, researches which treat the boundary contact case gradually appear.Garcke and Rauchecker [27] considered a stability analysis in a curved domain in R 2 via a linearization approach.Hensel and Stinson [29] proposed a varifold solution to (1.2) by starting from the energy dissipation property.For a gradient flow aspect of the Mullins-Sekerka flow, see e.g., [42, §3.2].
The numerical scheme that we propose in this paper is based on the BGN method, a parametric finite element method that allows the variational treatment of triple junctions and was first introduced by Barrett, Garcke, and Nürnberg in [6,7].We also refer to the review article [9] for more details on the BGN method, including in the context of the standard Mullins-Sekerka problem (1.3).Alternative front-tracking methods for geometric flows of curve networks have been considered in, e.g., [15,43,35,38,39].
Let us briefly review numerical methods being available in the literature for the Mullins-Sekerka problem and for its diffuse interface model, the multi-component Cahn-Hilliard equation (1.4).To the best of our knowledge, there are presently no sharp interface methods for the numerical approximation of the multi-phase Mullins-Sekerka problem.For the boundary integral method for the two-phase case, we refer the reader to [10,11,16,34,44].A level set formulation of moving boundaries together with the finite difference method was proposed in [17].For an implementation of the method of fundamental solutions for the Mullins-Sekerka problem in 2D, see [23].For the parametric finite element method in general dimensions, see [8,37].Numerical analysis of the scalar Cahn-Hilliard equation is dealt with in the works [3,4,13,21].Feng and Prohl [25] proposed a mixed fully discrete finite element method for the Cahn-Hilliard equation and provided error estimates between the solution of the Mullins-Sekerka problem and the approximate solution of the Cahn-Hilliard equation which are computed by their scheme.The established error bounds yielded a convergence result in [26].Aside from the sharp interface model, the Cahn-Hilliard equation for the multi-component case has been computed in several works, see [5,12,24,31,36].The multi-component Cahn-Hilliard equation on surfaces has recently been considered in [32].This paper is organized as follows.In the first part, we focus on the three-phase case, as outlined in the introduction.In Section 2, we show a curve-shortening property and an area-preserving property of strong solutions to the system.In Section 3, we introduce a weak formulation of the system, which in Section 4 will then be approximated with the help of an unfitted parametric finite element method.The scheme, which is linear, can be shown to be unconditionally stable.Section 5 is devoted to discussing solution methods for the linear systems that arise at each time level.In Section 6, we adapt our approximation to allow for an exact area-preservation property on the fully discrete level, leading to a nonlinear scheme.Section 7 is devoted to generalizations of the problem formulation and our numerical approximation to the general multi-phase case.Finally, we will show several results of numerical computation in Section 8, including convergence experiments for a constructed solution in the three-phase case.

Mathematical properties
In this section, we shall present two important properties of strong solutions to (1.2).Proposition 2.1 (Curve shortening property of strong solutions).Assume that (w, {Γ(t) } 0≤t≤T ) is a solution to (1.2).Then it holds that where we have defined the weighted length Here dL 2 and dH 1 refer to integration with respect to the Lebesgue measure and the onedimensional Hausdorff measure in R 2 , respectively.
Proof.Since the curvature is the first variation of the curves' length in the normal direction, we deduce, on recalling (1.1d) and (1.1c), that d dt (2.1) Now for each 1 ≤ j ≤ 3 it follows from integration by parts, (1.1a) and (1.1b) that Proof.We first deduce from the motion law (1.1c) that (2.3) Here, we note that ∇w 1 does not jump over Γ 1 (t) to derive the last equality.Meanwhile, integration by parts with (1.1a) and (1.1b) shows 0 = Combining (2.3) and (2.4) gives the assertion.The cases for j = 2, 3 can be shown in the same manner.
by discrete analogues of (3.1), (3.4) and (3.5).For each m ≥ 0 and 1 where Set ⃗ q m i,j := ⃗ X m i (q i,j ).Let T m be a sequence of triangulations of Ω and let S m be the associated scalar-and vector-valued finite element spaces, namely Let V h i be the set of all piecewise continuous functions on I which are affine on each subinterval [q i,j−1 , q i,j ], and let π h i : C(I) → V h i be the associated standard interpolation operators for 1 ≤ i ≤ 3. Similarly, V h i denotes the set of all vector valued functions such that each element belongs to We define the normal vector to each edge of Γ m i by where For two piecewise continuous functions on I, which may jump across the points q i,j (1 ≤ j ≤ N i ), we define the mass lumped inner product ⟨u, v⟩ h where u(q − i,j ) := lim [q i,j−1 ,q i,j ]∋y→q i,j u(y) and u(q + i,j ) := lim [q i,j ,q i,j+1 ]∋y→q i,j u(y) for each 1 ≤ i ≤ 3. We extend these definitions to vector-and tensor-valued functions.Moreover, we define ⟨u, v⟩ where here and throughout, the notation • (h) means an expression with or without the superscript h.
We make an assumption on the discrete vertex normals ⃗ ω m , following [6, Assumption A] and [9, Assumption 108], which will guarantee well-posedness of the system of linear equations: Here, for ⃗ ξ ∈ V h i and φ ∈ S m , we use the slight abuses of notation ⃗ ξ, φ We remark that the first condition basically means that each of the three curves has at least one nonzero inner vertex normal, something that can only be violated in very pathological cases.The proof of Theorem 4.1 shows that it is actually sufficient to require this for just two out of the three curves, but for simplicity we prefer to state the stronger assumption.The second condition in Assumption 1, on the other hand, is a very mild constraint on the interaction between bulk and interface meshes.In fact, it can only be violated if all the vectors in the set are linearly dependent, which happens, for example, if the three curves are straight lines that lie on top of each other. Given ) such that the following conditions hold: We stress that (4.2) encodes two different schemes: One that uses mass lumping in the two bulk-surface terms in (4.2a) and (4.2b), and one that uses true integration in both.The interpolation operator π h i in (4.2a) is superfluous in the former case, but necessary for the stability proof of the latter.Writing (4.2) as above allows for a compact presentation.Observe that the implementation of the scheme with mass-lumping is far easier, since there bulk finite element functions only need to be evaluated at the vertices of the curve network.We refer to [8] for more details.
Proof.Since (4.2a), (4.2b) and (4.2c) is a linear system with the same number of unknowns and equations, existence follows from uniqueness.To show the latter, it is sufficient to prove that only the zero solution solves the homogeneous system.Hence let ) ) Thus, we see that i=1 C i = 0. We deduce from (4.3a) that and so the second condition in Assumption 1 yields that ⃗ . Hence, we obtain so that the first condition in Assumption 1 implies that κ σ,i = 0 for i = 1, 2, 3, i.e. κ σ = 0. Since C 1 = C 2 = C 3 , they must all be zero, and so W = 0 also follows.Hence, we have shown the existence of a unique solution (W m+1 , ⃗ X m+1 , κ m+1 σ Before proving the stability of our scheme, we recall the following lemma from [9, Lemma 57] without the proof. where |Γ h | and | ⃗ X (Γ h )| are the lengths of Γ h and ⃗ X (Γ h ), respectively.
) be a solution to (4.2).Then the following estimate is satisfied: where we recall that Proof.We choose φ = W m+1 in (4.2a) and ξ We compute from (4.6) and (4.7), on recalling Lemma 4.2, that This proves the desired result.
Remark 4.4.Similar stability results for the discretization of other gradient flows of curve networks can be found in [7,6].For the two-phase Mullins-Sekerka problem, related stability results were derived in [8,37].

Solution of the linear system
In this section, we discuss solution methods for the systems of linear equations arising from (4.2) at each time level.To this end, we make use of ideas from [6,8].Here the crucial idea is to avoid having to work with the trial and test spaces V (Γ m ) and S m Σ directly, and rather employ a technique that is similar to a standard treatment of periodic boundary conditions for ODEs and PDEs.In particular, following [6, (2.44)] we introduce the orthogonal projections Σ .On letting 1 = (1, 1, 1) T , it is easy to see that for W ∈ S m it holds that ) be the unique solution to (4.2) whose existence has been proven in Theorem 4.1.Let N := 3 i=1 (N i + 1) be the sum of the vertices on each individual curve, and let K m Ω be the number of vertices of the mesh T m inside Ω.From now on, as no confusion can arise, we identify (W m+1 , κ m+1 σ , δ ⃗ X m+1 ) with their vectors of coefficients with respect to the bases {Ψ m i } 1≤i≤K m Ω and {{Φ i,j } 1≤j≤N i } 3 i=1 of the unconstrained spaces S m and V (Γ m ).In addition, we let the Euclidean space equivalent of P, and similarly for the equivalent Q : Then the solution to (4.2) can be written as (QW m+1 , κ m+1 σ , ⃗ X m + P δ ⃗ X m+1 ) for any solution of the linear system where for each 1 ≤ c ≤ 3. The advantage of the system (5.1) over a naive implementation of (4.2) is that complications due to nonstandard finite element spaces are completely avoided.A disadvantage is, however, that the system (5.1) is highly singular, in that due to the presence of the projections the dimension of its kernel is larger than the dimension of the scalar bulk finite element space S m .This makes it difficult to solve (5.1) in practice.A more practical formulation can be obtained by eliminating one of the components of W m+1 completely.In particular, on recalling that W m+1 • 1 = 0, we can reduce the unknown variables where I M denotes the identity matrix of size M for M ∈ N.
Then the solution to (4.2) can be written as ) for any solution of the reduced linear system where In contrast to (5.1), the kernel of (5.3) is small.In fact, it has dimension 8 due to the fact that P has a kernel of dimension 8. Hence, iterative solution methods, combined with good preconditioners, work very well to solve (5.3) in practice.For our numerical results in Section 8, below, we employ a GMRes iterative solver with least squares solution of the block matrix in (5.3) without P as preconditioners.
6 Obtaining a fully discrete area conservation property Although the linear scheme (4.2) introduced in Section 4 can be shown to be unconditionally stable, recall Theorem 4.3, in general the areas occupied by the discrete approximations of the three phases will not be conserved.In this section, we state how to modify the previously introduced scheme (4.2) in such a way, that it satisfies both of the structure defining properties from Section 2. To this end, we follow the discussion in [37, §3] in order to obtain an exact area preservation property on the fully discrete level.We remark that our approach hinges on ideas first presented for area-conserving geometric flows for closed curves in [30,2].See also [1, §3.2] for related work in the context of the surface diffusion flow for curve networks with triple junctions.Let us define families of polygonal curves {Γ h i (t)} t≥0 , i = 1, 2, 3, that are parameterized by the time variable.In particular, for each 0 ≤ m ≤ M , t ∈ [t m , t m+1 ] and 1 ≤ i ≤ 3, we define the polygonal curve Γ h i (t) by Precisely speaking, the vertices of Γ h i (t) are defined as follows: ⃗ q h i,j (t) := for 0 ≤ j ≤ N i , while we write each edge of Γ h i (t) as e h i,j (t) := [⃗ q h i,j−1 (t), ⃗ q h i,j (t)] for 1 ≤ i ≤ 3 and 1 ≤ j ≤ N i .Lemma 6.1.For each m ≥ 0 and (i, j, k) ∈ Λ, it holds that Proof.The desired result for i = 1, 2 is shown in [1, Lemma 3.1], and the result for i = 3 follows analogously on noting that ∂Ω is fixed.Now the weighted vertex normal vector ⃗ ω is defined through the following formula: Consequently, we obtain a nonlinear system with the aid of ⃗ ω We can now prove the area preserving property of each domain surrounded by the polygonal curve in the discrete level.Theorem 6.2 (Area preserving property for the discrete scheme).Let m ≥ 0 and let ) be a solution of (6.3).Then, for each 1 Proof.We will argue similarly to the proof of [37,Theorem 3.3].
Σ in (6.3a), we obtain from (6.2) and Lemma 6.1 that This yields the desired result for k = 1.The other cases can be treated analogously.
) be a solution of (6.3).Then, it holds that Proof.The proof is analogous to the proof of Theorem 4.3 once we replace ⃗ ω m by ⃗ ω m+ 1 2 .

Generalization to multi-component systems
In order to simplify the presentation, in the previous sections we concentrated on the simple three phase situation depicted in Figure 1.However, it is not difficult to generalize our introduced finite element approximations to the general multi-phase case.We present the details in this section, following closely the description of the general curve network used in [1, §2].

Problem setting
For later use, we let N ≤M := {1, • • • , M } for each M ∈ N. First, let us introduce counters which will be used frequently later on.Given a curve network Γ(t), I C ≥ 1 denotes the number of curves which are included in Γ(t).Thus, we have Γ(t) = ∪ I C i=1 Γ i (t), where each Γ i (t) is either open or closed.Let I P ≥ 2 be the number of phases, i.e. the not necessarily connected subdomains of Ω which are separated by Γ(t).This means that Ω\Γ(t) = I P j=1 Ω j (t).Finally, each endpoint of an open curve included in Γ(t) is part of a triple junction.We write the number of triple junctions as I T ≥ 0.
Assumption 2 (Triple junctions).Every curve Γ i (t), i ∈ N ≤I C , must not self-intersect, and is allowed to intersect other curves only at its boundary ∂Γ i (t).
Assumption 3 (Phase separation).The curve network Γ(t) is equipped with a matrix O : {−1, 0, 1} I P ×I C that encodes the orientations of the phase boundaries.In particular, each row contains nonzero entries only for the curves that make up the boundary of the corresponding phase, with the sign specifying the orientation needed for the curves normal to make it point outwards of the phase.I.e. for p ∈ N ≤I P and i ∈ N ≤I C we have For every i ∈ N ≤I C , there exists a unique pair In this situation, we say that Ω p 1 (t) is a neighbour to Ω p 2 (t).
Clearly, given the matrix O, the boundary of Ω p (t) can be characterized by where we have assumed that Ω I P (t) is the only phase with contact to the external boundary.
Remark 7.1 (Examples).We note that for the three-phase problem shown in Figure 1, we have Examples for more complicated networks, and their description in our general framework, can be found in the numerical results section, see §8.
Under these preparations, we can generalize the system (1.2) to the multi-phase case.Where no confusion can arise, we use the same notation as before, e.g., T Σ = {u ∈ R I P | I P j=1 u j = 0}.Given an initial curve network Γ 0 , our aim is to find w : Ω → T Σ and the evolution of a curve network {Γ(t)} t≥0 which satisfy where χ denotes the vector of the characteristic functions χ j = X Ω j (t) , j = 1, . . ., I P , as before.
Remark 7.2.The system (7.1)does not depend on the choice of normals ⃗ ν i , for i ∈ N ≤I C .Indeed, if we take −⃗ ν i as the unit normal vector to Γ i (t), then the sign of κ i reverses.On the other hand, the sign of the jump [χ] is also reversed.Thus, the second law of (7.1) does not change.Meanwhile, the sign of the normal velocity V i in the third condition is also reversed, balancing with the sign change of ⃗ ν i on the left hand side.

Weak formulation
Similarly to Section 3, it is possible to derive a weak formulation of the multi-phase problem (7.1).On noting that the jump of the j-th characteristic function χ j across the curve Γ i is −O ji , we can compute from the third condition in (7.1) that Hence, overall we obtain the following weak formulation:

Finite element approximations
We now generalize our finite element approximation (4.2) to the multi-phase case.The necessary discrete function spaces are the obvious generalizations, for example where ρ k i ∈ {0, 1} encodes at which if its two endpoints the discrete curve Γ m c k i meets the k-th triple junction.
Given ⃗ X m ∈ V (Γ m ), we find (W m+1 , ⃗ X m+1 , κ m+1 σ ) such that the following conditions hold: (7.3c) Remark 7.3 (Linear system).The linear system of equations arising at each time level of (7.3) is given by the obvious generalization of (5.1), where the block matrix entries of (5.Finally, on using the techniques from Section 6, we can adapt the approximation (7.3) to obtain a structure preserving scheme that is unconditionally stable and that conserves the areas of the enclosed phases exactly.Given ( Γ m i = 0. (7.4b) We conclude this section by stating theoretical results for the generalized schemes.Their proofs are straightforward adaptations of the proofs of Theorems 4.1, 4.3 and 6.2.

Numerical results
We implemented the fully discrete finite element approximations (7.3) and (7.4) within the finite element toolbox ALBERTA, see [41].The arising linear systems of the form (5.3) are solved with a GMRes iterative solver, applying as preconditioner a least squares solution of the block matrix in (5.3) without the projection matrices P .For the computation of the least squares solution we employ the sparse factorization package SPQR, see [19].
For the triangulation T m of the bulk domain Ω, that is used for the bulk finite element space S m , we use an adaptive mesh that uses fine elements close to the interface Γ m and coarser elements away from it.The precise strategy is as described in [8, §5.1] and for a domain Ω = (−H, H) d and two integer parameters N c < N f results in elements with maximal diameter approximately equal to h f = 2H N f close to Γ m and elements with maximal diameter approximately equal to h c = 2H Nc far away from it.For all our computations we use H = 4.An example adaptive mesh is shown in Figure 3, below.
We stress that due to the unfitted nature of our finite element approximations, special quadrature rules need to be employed in order to assemble terms that feature both bulk and surface finite element functions.For all the computations presented in this section, we use true integration for these terms, and we refer to [8,37] for details on the practical implementation.Throughout this section we use (almost) uniform time steps, in that τ m = τ for m = 0, . . ., M − 2 and τ M −1 = T − t m−1 ≤ τ .Unless otherwise stated, we set σ = 1.

3 phases
In the next set of experiments, we investigate how a standard double bubble and a disk evolve, when one phase is made up of the left bubble, and the other phase is made up of the right bubble and the disk.With the notation from §7.1 we have I C = 4, The two bubbles of the double bubble enclose an area of about 3.139 each, while the disk has an initial radius of 5  8 , meaning it initially encloses an area of 25π 64 ≈ 1.227.During the evolution the disk vanishes, and the right bubble grows correspondingly, see Figure 3.We note that our theoretically framework does not allow for changes of topology, e.g., the vanishing of curves.Hence, in our computations we perform heuristic surgeries whenever a curve becomes too short.Here a closed curve is simply discarded, while a curve that was part of a network is removed.This will leave two triple junctions, where only two curves meet, and the involved curves can be glued together so that the simulation can continue.Repeating the simulation with a bigger initial disk gives the results in Figure 4.Here the radius is 5  4 , so that the enclosed area is 4.909.Now the disk grows at the expense of the right bubble, so that eventually two separate phases remain.With the next simulation we demonstrate that in the given setup, which of the two components of phase 1 survives is not down the initial size.In particular, we allow the initial disk to have area 3.300, so that it is bigger than the other component of the same phase: the right bubble in the double bubble.And yet, due to the perimeter of the bubble being overall cheaper than the boundary of the disk, the latter shrinks to extinction.See Figure 5.
For the next experiment we start from a nonstandard triple bubble, where we choose the right most bubble to have area 1.5, while the other two bubbles have unit area.We assign the two outer bubbles to belong to the same phase.In particular, with the notation from §7.1 we have I C = 6, 3 ) = (2, 4, 6),     We observe that during the evolution the larger bubble on the right grows at the expense of the left bubble, until the latter one vanishes completely.The remaining interfaces then evolve towards a standard double bubble with enclosed areas 1 and 2.5.See Figure 6.
In the final numerical simulation for the setting with 3 phases, we consider the evolution of two double bubbles.With the notation from §7.
The first bubble is chosen with enclosing areas 3.14 and 6.48, while the second double bubbles encloses two areas of size 3.64.In each case, the left bubble is assigned to phase 1, while the right bubbles are assigned to phase 2. In this way, the lower double bubble holds the larger portion of phase 1, while the upper double bubble holds the larger portion of phase 2. Consequently, each double bubble evolves to a single disk that contains just one phase.See Figure 7.

4 phases
In the next set of experiments, we investigate simulations for a nonstandard triple bubble, with one of the bubbles making a phase with a separate disk.In particular, with the notation from §7.The three bubbles of the triple bubble enclose an area of unity each, while the disk has an initial radius of 1  2 , meaning it initially encloses an area of π 4 ≈ 0.785.During the evolution the disk vanishes, and the right bubble grows correspondingly, see Figure 8. Repeating the experiment with initial data where the disk is a unit disk leads to the evolution in Figure 9, where the disk now expands and survives, at the expense of the right bubble in the triple bubble.

3 phases
As an example for non-equal surface energy densities for the various curves, we repeat the simulation in Figure 3, but now weigh curves 1 and 3 in the double bubble with σ 1 = σ 3 = 2, while keeping the other two densities at unity.This now means that in contrast to Figure 3, it makes energetically more sense to increase the size of the single bubble, while shrinking the bubble that is surrounded by the more expensive interfaces.See Figure 3 for the observed evolution.

4 phases
Similarly, if we make the interface of the single circular bubble in the initial data in Figure 9 more expensive, it will no longer grow but shrink to a point.Setting the weight for the curve to σ 7 = 2 leads to the evolution seen in Figure 11.−1/R 2 (t) on Γ 2 (t) and w 3 − w 1 = 1/R 3 (t) on Γ 3 (t).Thus, w satisfies the second condition in (7.1).We move on the motion law.A direct calculation shows R 3 (t) on Γ 3 (t). (A.2) Hence, the third condition of (7.1) is valid by (A.1) and (A.2).Therefore, w = (w 1 , w 2 , w 3 ) T given by (8.1) is an exact solution of (7.1) with σ = 1.

Figure 1 :
Figure 1: Three open curves with two triple junctions.

Theorem 4 . 1 (
Existence and uniqueness).Let Assumption 1 hold and let m ≥ 0. Then there exists a unique solution
1 we have I C = 3, I P = 3, I T = 0 and O =

Figure 2 :
Figure 2: A network of curves which consists of three concentric circles. h

Figure 3 :
Figure 3: The solution at times t = 0, 0.4, 0.8, 1, and a plot of the discrete energy over time.Below we show the adaptive bulk mesh at times t = 0 and t = 1.

Figure 4 :
Figure 4: The solution at times t = 0, 4, 4.5, 6, and a plot of the discrete energy over time.

Figure 5 :
Figure 5: The solution at times t = 0, 4, 5, 7, and a plot of the discrete energy over time.

Figure 6 :
Figure 6: The solution at times t = 0, 0.4, 0.6, 1, and a plot of the discrete energy over time.

Figure 8 :
Figure 8: The solution at times t = 0, 0.4, 0.7, 1, and a plot of the discrete energy over time.

Figure 9 :
Figure 9: The solution at times t = 0, 0.4, 0.8, 1, and a plot of the discrete energy over time.

Figure 10 :
Figure 10: The solution at times t = 0, 1, 2, 3, and a plot of the discrete energy over time.

Figure 11 :
Figure 11: The solution at times t = 0, 1, 1.8, 3, and a plot of the discrete energy over time.