Addressing Troubles with Double Bubbles: Convergence and Stability at Multi-Bubble Junctions

by   Christopher Batty, et al.

In this report we discuss and propose a correction to a convergence and stability issue occurring in the work of Da et al.[2015], in which they proposed a numerical model to simulate soap bubbles.



There are no comments yet.


page 1

page 2


A discontinuous Galerkin pressure correction scheme for the incompressible Navier-Stokes equations: stability and convergence

A discontinuous Galerkin pressure correction numerical method for solvin...

Sensitivity Analysis of Domain Decomposition method for 4D Variational Data Assimilation (DD-4DVAR DA)

We prove consistence, convergence and stability of the Domain Decomposit...

Deterministic Sparse Sublinear FFT with Improved Numerical Stability

In this paper we extend the deterministic sublinear FFT algorithm for fa...

Modified Double DQN: addressing stability

Inspired by double q learning algorithm, the double DQN algorithm was or...

Numerical analysis of a model of two phase compressible fluid flow

We consider a model of a binary mixture of two immiscible compressible f...

Addestramento con Dataset Sbilanciati

English. The following document pursues the objective of comparing some ...

Stacking and stability

Stacking is a general approach for combining multiple models toward grea...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

In this report we discuss and propose a correction to a convergence and stability issue occurring in the work of Da et al. [2015], in which they proposed a numerical model to simulate soap bubbles.

In the original implementation of their work, convergence of the geometry towards equilibrium surfaces did not behave as expected. Soap foam should converge to a configuration described by Plateau’s laws [1873], where the dihedral angles of faces incident to Plateau borders (triple junction edges) should be 120 degrees. However, the existing method of Da et al. converges to a steady state that fails to satisfy this condition: the angles differ noticeably from 120 degrees. For example, in a simple double-bubble test case, the dihedral angles of the pair of faces incident to the outer air region are approximately 135 degree, as seen in Figure 1-middle.

We traced this issue to a flaw in the original implementation, which applied a special treatment in computing the appropriate vertex area at the triple-junction that maintains stability of surface tension forces. Unfortunately, this treatment also sacrificed the correct convergence behavior of soap foam. This undesired behavior was first observed by Ishida et al. [2017]. In the following we briefly describe the problem and introduce a new method to compute the vertex area, which simultaneously maintains stability and yields the correct Plateau border angles.

2. Integral and pointwise mean curvatures

Figure 2. (a)

: A surface is intersected with a plane containing the normal vector. In the continuous setting, we can compute the mean curvature of a point on this surface by integrating through all the angles (dark gray circle);

(b): In the discrete setting, we can use the cotan-formula to compute the mean curvature from all nearby triangles, normalizing by the vertex (Voronoi) area.

To understand the problem, we begin with a discussion of mean curvature. In the continuous setting (Figure 2a), the mean curvature at a point can be interpreted as an integral of the signed curvature on the normal section curve over all directions (represented by the dark gray circle in Figure 2a).

In the discrete setting of a numerical simulation, only finitely many triangles are used to approximate the surface. The pointwise discrete mean curvature of the surface can be computed for each (volumetric) region with the following formula [Sullivan, 2008] for vertex ,

where is the edge curvature of the -th edge connected to vertex , and is normalized by the vertex area (Figure 2b). The edge curvature is given by  [Cohen-Steiner and Morvan, 2003]. The normalization by area converts integrated curvature into the pointwise curvature which is then applied as a force by Da et al.

For two volumetric regions sharing an interface the absolute pointwise mean curvatures are trivially the same; the two sides of the surface are geometric complements and therefore differ only in sign, leading to a natural force balance at equilibrium. However, when three or more regions are present, a vertex may lie on a triple (or higher-order) junction, and computing discrete mean curvature separately in each region at the junction becomes problematic. More specifically, since the discrete mean curvature is proportional to the inverse of the vertex area, the choice of how to define the vertex area per region has a critical effect, as we describe below.

3. A problematic case

Figure 3. (a): A triple junction comprised of triangles. On each vertex the mean curvatures per region are computed, with the subdivided small disks indicating the set of regions involved. (b): The same triple junction but with the triangles between region (II) and (III) stretched. In this case, under the naïve discretization (only) the discrete mean curvatures in region (II) and (III) of the vertices on the stretched triangles are affected (marked with dark gray), yet the triple junction geometry has not changed.

We will consider the triple-junction in Figure 3

as our motivating example. If one computes the mean curvatures at the junction separately per volumetric region using only the triangles comprising each region’s surface, then the area of a triangle on one particular branch of the junction contributes to the curvature estimate for only the two volumes that share that triangle (Figure 

3b). The curvature estimate for the third region is entirely independent of that triangle’s area; this implies that if the mesh is initially in a numerical equilibrium and that triangle is then stretched (or remeshed) in a manner that doesn’t change the angles at the triple-junction, the observed curvature values of the two incident volumes will change accordingly, but that of the third will not. That is, force balance will be lost, despite no change in the actual angles at the junction! In the limit of mesh refinement the vertex area conceptually shrinks to zero ensuring a valid discretization, but in any discrete simulation the delicate force balance is often quickly destroyed. This problem compounds from one timestep to the next as the resulting forces attempt to compensate for changing incident triangle areas by adapting the junction angles, leading to further instability. As a consequence, when simulating double bubbles, the bubbles will incorrectly accelerate and drift away (Figure 1-top).

4. A preliminary solution

In the original implementation of Da et al. [2015], the authors noticed this problem and adopted a simple solution that modifies the vertex area computation at the triple junction: the vertex area used for all the incident regions is set to be the shared edge length along the triple junction times a global constant, which is set to be a fixed fraction of the mean edge length over the entire mesh ( is used in the original implementation). In this way, the discrete mean curvatures computed at the triple junction reflect the relevant angles, but are independent of the areas of the particular triangulation.

This simple solution can indeed stabilize the simulation since the mean curvatures will not be affected by otherwise irrelevant disturbances in the local triangulation; for example, in-plane vertex-smoothing of nearby points no longer suddenly upsets the balance. Nevertheless, a global constant times the edge length cannot well approximate the actual local vertex area. In the particular example of our double-bubble, this “effective vertex area” for the exterior region at the triple-junction is larger than the actual value computed directly from the triangles, while the vertex area computed inside the bubble is smaller. As a result, the mean curvature outside the bubble is smaller and yields a larger surface tension force, which causes the dihedral angle between the bubbles to be larger than the desired equilibrium value (Figure 1-middle).

5. Our solution

# Subdivision # Faces Angular Deviation
(per bubble) (initial)
16 880
24 1,992
48 8,016
96 32,160
Table 1. Angular deviation versus resolution: the observed angular deviation from at the triple junction decreases as the resolution of the mesh increases, confirming that our method converges under refinement.

Although the remeshing process can cause unpredictable changes in the local triangulation, we identified a strategy that can always be safely applied: we evaluate the vertex areas per region, determine their average, and use that value for normalization when computing each of the per-region discrete curvatures. This technique is simple to implement, incurs minimal additional cost, and can be justified as follows.

Under this strategy, since the normalizing vertex areas used in computing the curvature for each region at the vertex are the same, when a nearby vertex in any region is disturbed, the change is immediately reflected in the computed mean curvatures for all the incident regions. This greatly improves stability, similar to the fix by Da et al. However, unlike the correction used by Da et al., the vertex areas used in our method faithfully reflect the actual areas of the surrounding triangles, rather than an arbitrary global constant. Hence convergence can be assured since valid mean curvatures are computed for the vertex in each region. In other words, the normalization of curvatures at the triple junction is now scaled appropriately by local area, consistent with how the discrete curvatures are scaled by area on other (manifold, non-junction) parts of the mesh.

In Figure 1 bottom, we demonstrate that our method produces a stable simulation of the double-bubble, and that bubbles at equilibrium have a dihedral angle much closer to the theoretical value of . In Table 1, we demonstrate that our discretization converges: as the resolution of the mesh increases, the angular deviation decreased dramatically from to , indicating that the dihedral angles in our simulation converges towards the theoretical value. The angular deviation was computed as the root-mean-square deviation compared with , across all the dihedral angles at the triple junction (consistent with Ishida et al. [2017]).


  • D. Cohen-Steiner and J. Morvan (2003) Restricted delaunay triangulations and normal cycle. In Proceedings of the nineteenth annual symposium on Computational geometry, pp. 312–321. Cited by: §2.
  • F. Da, C. Batty, C. Wojtan, and E. Grinspun (2015) Double bubbles sans toil and trouble: discrete circulation-preserving vortex sheets for soap films and foams. ACM Transactions on Graphics (TOG) 34 (4), pp. 149. Cited by: §1, §4.
  • S. Ishida, M. Yamamoto, R. Ando, and T. Hachisuka (2017) A hyperbolic geometric flow for evolving films and foams. ACM Transactions on Graphics (TOG) 36 (6), pp. 199. Cited by: §1, §5.
  • J. A. F. Plateau (1873) Statique expérimentale et théorique des liquides soumis aux seules forces moléculaires. Vol. 2, Gauthier-Villars. Cited by: §1.
  • J. M. Sullivan (2008) Curvatures of smooth and discrete surfaces. In Discrete differential geometry, pp. 175–188. Cited by: §2.