Bubble Meshing: Techniques & Applications

• Bubble meshing is a physically-based surface triangulation method that leverages bubble dynamics and conformal mapping to efficiently generate Delaunay triangles.
• It integrates approaches such as adaptive front refinement, level set interface capture, and spectral filtering to handle evolving geometries and ensure numerical stability.
• The method is applied in fluid dynamics, fracture-driven growth, and analytical constructions, optimizing computation by focusing mesh resolution where needed.

Bubble meshing most specifically denotes a physically-based surface triangulation method in which physical bubbles are simulated to automatically generate mesh vertices, resulting in high-quality Delaunay triangles (Shangyu et al., 7 Aug 2025). In adjacent computational literature, related bubble-centered meshing practices also appear in interface-capturing two-fluid FEM/level set formulations, polygonal constructions for compacted bubble clusters, adaptive fracture-front meshing for growing gas bubbles, and stabilized surface discretizations for boundary element bubble dynamics (Doyeux et al., 2012). This suggests that the term is used both for a named triangulation algorithm and, more broadly, for mesh-generation or mesh-management strategies organized around bubble geometry, bubble packing, or bubble-interface evolution.

1. Scope and meanings

The literature represented here does not use “bubble meshing” in a single narrow sense. In "An Improved Physically-Based Surface Triangulation Method" (Shangyu et al., 7 Aug 2025), bubble meshing is explicitly a surface triangulation method. In "Simulation of two-fluid flows using a Finite Element/level set method. Application to bubbles and vesicle dynamics" (Doyeux et al., 2012), the central issue is not explicit bubble packing but the meshing of a moving bubble-fluid interface on a fixed, globally unstructured mesh. In "Symmetry break in the eight bubble compaction" (Bevilacqua, 2020), bubble geometry generates a polygonal tiling with flat and curved interfaces. In "Fracture-Driven Single Bubble Grows and Migration Model in Aquatic Muds" (Katsman, 26 Sep 2025), the emphasis is adaptive refinement and mesh motion at a bubble front governed by LEFM. In "GPU accelerated fast multipole boundary element method for simulation of 3D bubble dynamics in potential flow" (Gumerov et al., 2019), bubble surfaces are triangulated and then stabilized by filtering. In "Doubly-periodic array of bubbles in a Hele-Shaw cell" (Silva et al., 2010), “meshing” appears in the geometric sense of constructing periodic bubble networks by conformal mapping.

Context	Mesh treatment	Representative paper
Surface triangulation	Bubble packing induces planar Delaunay triangles	(Shangyu et al., 7 Aug 2025)
Two-fluid interface simulation	Fixed unstructured mesh with immersed interface	(Doyeux et al., 2012)
Bubble-cluster geometry	Dual tessellation and polyhedral interfaces	(Bevilacqua, 2020)
Fracture-driven growth	Front-focused adaptive 1D/2D/3D refinement	(Katsman, 26 Sep 2025)
Potential-flow dynamics	Triangular surface mesh with spectral filtering	(Gumerov et al., 2019)
Hele-Shaw arrays	Conformal maps for periodic bubble networks	(Silva et al., 2010)

A useful way to organize the subject is therefore by mechanism: bubble packing for triangulation, immersed-interface treatment for moving boundaries, adaptive refinement near bubble fronts, and analytic or geometric bubble-network construction. This broader organization is an inference from the papers rather than a formal taxonomy stated in any single source.

2. Physically-based surface triangulation

The named method “bubble meshing” is a physically-based surface triangulation procedure in which the target surface is first flattened onto a plane by a conformal mapping $f : M \rightarrow \mathbb{C}$, circles are packed and relaxed in the planar domain, planar Delaunay triangulation is generated from bubble centers, and the resulting mesh is inverse-mapped back to the original surface by barycentric coordinates. The inverse mapping is described through

$$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$

The improved method employs conformal mapping to simplify surface bubble packing by flattening the surface onto a plane, avoids direct bubble movement on the surface, optimizes bubble quantity control, and separates that control from the relaxation process; the reported consequence is a computation-time reduction of over $70\%$. The enhanced method enables efficient triangulation of disk topology surfaces and supports local size control, curvature adaptation, and re-meshing of discrete surfaces. Reported examples include a sample surface remesh improving from $104\text{s}$ to $18\text{s}$, a flat plate with hole triangulated in $0.495\text{s}$, minimum angles above $30^\circ$, and more than $80\%$ of triangles with minimum angles above $45^\circ$ (Shangyu et al., 7 Aug 2025).

Algorithmically, the paper’s central change is boundary-focused overlap control. Instead of repeatedly adjusting all bubbles in each iteration, the method only monitors boundary bubbles for overlap, removes bubbles when the overlap exceeds a threshold, and invokes relaxation once after excessive overlaps have been removed. The physical relaxation stage is written as

$m_i \ddot{x}_i + c_i \dot{x}_i = f_i.$

This makes the method less dependent on repeated global updates and less sensitive to overlap-threshold tuning.

The significance of the conformal-map formulation is that computationally expensive surface operations are replaced by planar operations. The paper explicitly lists avoided costs such as geodesic distance calculation, surface tangency or attraction, and dynamic boundary adaptation on the original curved surface. Within the papers considered here, this is the clearest instance in which “bubble meshing” refers to a direct mesh-generation algorithm rather than to a mesh-handling strategy for bubble simulations.

3. Fixed-mesh interface treatment for bubble flows

A distinct line of work treats bubble meshing as the representation of a moving bubble-fluid interface on a fixed mesh. In the FEM/level set framework for two-fluid flow and vesicle dynamics, the interface $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$0 is represented implicitly by a level set function $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$1, with $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$2 in $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$3, $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$4 in $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$5, and $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$6 on $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$7. The interface evolves by

$$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$8

The mesh is globally unstructured and not fitted to the moving interface; instead, the interface is immersed and tracked via the level set on a fixed mesh. There is no explicit remeshing around the interface. Fluid properties and interfacial forces are smeared over a few elements near the interface, with $$u = \sum_{i=1}^3 \lambda_i(u)\,u_i, \qquad \Psi(u) = \sum_{i=1}^3 \lambda_i(u)\,r_i.$$9 typically taken as $70\%$0, where $70\%$1 is the mesh size near the interface. Smoothed Heaviside and delta functions are used to interpolate discontinuous material properties and replace surface integrals by volume integrals near the interface. High-order polynomial approximation is used to increase simulation accuracy, and reinitialization by fast marching is used to maintain the signed-distance property $70\%$2 (Doyeux et al., 2012).

For bubble dynamics, the governing two-fluid equations are the Navier–Stokes equations with level-set-dependent density and viscosity, together with a surface-tension contribution approximated by

$70\%$3

For vesicles, the model is enriched by a Lagrange multiplier enforcing inextensibility,

$70\%$4

The discrete multiplier $70\%$5 is only assembled in elements intersected by the interface, because it is weighted by a smoothed delta concentrated in a narrow band.

The method’s principal advantage is that no remeshing is required as the interface moves, while arbitrary topology changes such as merging or splitting are handled naturally by the level set function. Its principal limitation is equally explicit: the method does not employ adaptive mesh refinement or local mesh adaptation. Interface diffusion, numerical thickness over a few mesh cells, sensitivity in highly deformed cases, and perimeter loss remain central issues. Within the broader topic of bubble meshing, this paper establishes that high-order approximation on a fixed mesh is an alternative to interface-fitted remeshing.

4. Adaptive and stabilized bubble-surface meshes

Another group of methods keeps the mesh explicitly attached to the bubble boundary and modifies it as the geometry evolves. In the fracture-driven single-bubble growth model for muddy aquatic sediments, the domain begins with a penny-shaped crack, and $70\%$6 mesh elements of maximum size $70\%$7 are placed along the initial curved bubble front. These edge elements seed refined $70\%$8 triangular elements on the bubble surface and local symmetry plane, and the $70\%$9 layer then serves as the basis for $104\text{s}$0 tetrahedral elements in the bulk. Coarsening proceeds away from the front with a growth rate of $104\text{s}$1. Mesh movement is coupled directly to LEFM: local crack increments are computed along the crack front, mesh nodes near the front are moved by those increments, Laplace smoothing maintains quality after node movement, and full remeshing is performed only if mesh quality falls below an acceptable threshold. The local update criterion is

$104\text{s}$2

and crack propagation occurs when $104\text{s}$3 (Katsman, 26 Sep 2025).

This strategy differs sharply from the fixed-mesh level set approach. Here, fidelity is concentrated at the bubble surface and particularly at the advancing fracture front, and aggressive coarsening controls the total element count. The paper explicitly states that no further mesh refinement is needed during simulation after validation against analytical solutions. A plausible implication is that this formulation is designed for problems where the most important gradients are localized at the bubble front and cannot be represented adequately by an interface smeared over several cells.

A separate stabilization problem arises in $104\text{s}$4 BEM for bubble dynamics in potential flow, where only the bubble surfaces are discretized. There the surface meshes are triangular, and a new smoothing technique using a surface filter is introduced. For bubbles topologically equivalent to a sphere, Cartesian coordinates are expanded in spherical harmonics,

$104\text{s}$5

and then truncated at bandwidth $104\text{s}$6. The filtering operator is

$104\text{s}$7

For toroidal bubbles, a $104\text{s}$8 Fourier basis replaces spherical harmonics. Filtering is applied twice per time step, and additional stabilization is obtained by moving surface mesh points with only the normal component of velocity. Reported behavior is explicit: without shape filtering, all tested meshes destabilized within $104\text{s}$9–$18\text{s}$0 time steps; with the filter, simulations ran indefinitely without instability. The method was used for bubble clusters with thousands of bubbles, including $18\text{s}$1 bubbles and more than $18\text{s}$2 million mesh vertices (Gumerov et al., 2019).

Taken together, these two papers show two complementary approaches to bubble-surface meshing. One uses front-focused refinement, node motion, and occasional remeshing; the other uses fixed-connectivity surface meshes stabilized by spectral filtering. Both reject the assumption that accurate bubble simulation necessarily requires repeated global remeshing.

5. Geometric and analytical construction of bubble networks

Bubble meshing also appears in geometric constructions of bubble aggregates and periodic bubble arrays. In the eight-bubble compaction problem, the spatial arrangement of seven peripheral bubbles around a central bubble is dictated by maximization of the minimum mutual distance between centroids on the sphere, i.e. a Tammes-type maximin criterion,

$18\text{s}$3

From this arrangement, the central sphere is tiled by a dual tessellation whose faces are $18\text{s}$4 triangle, $18\text{s}$5 quadrilaterals, and $18\text{s}$6 pentagons. Peripheral bubbles are then moved radially inwards under volume conservation. The contact polygons are flat planar interfaces, while the remaining free surfaces are curved spherical caps. Mechanical balance is imposed by edge-force equilibrium,

$18\text{s}$7

together with Laplace–Young relations for free surfaces and volume conservation. The paper reports anisotropy in the field of surface tensions, with circumferentially oriented bubble-bubble interfaces carrying larger tensions than radially oriented ones, and proposes this anisotropy as a mechanical cue for symmetry breaking (Bevilacqua, 2020).

The Hele-Shaw literature develops a more analytical construction. Exact solutions for doubly-periodic arrays of steadily moving bubbles are obtained by conformal mapping when surface tension is neglected. Symmetry assumptions reduce the flow domain to a simply connected region, and the physical map is written as

$18\text{s}$8

The mappings $18\text{s}$9 and $0.495\text{s}$0 are built from generalized Schwarz–Christoffel formulae with slits encoding bubble locations. The paper explicitly formulates a principle of meshing: adding or removing slits in the parameter domain corresponds to adding or removing rows or columns of bubbles, enabling systematic construction or decimation of periodic bubble networks, including multi-file arrays and doubly-periodic solutions in unbounded cells (Silva et al., 2010).

These works broaden the topic beyond numerical discretization in the narrow finite-element sense. Here, “meshing” concerns explicit geometric decomposition of space into bubble-related cells, interfaces, and periodic units. The common feature is not a single data structure but a constructive rule linking bubble arrangement to a discretized geometry.

6. Limitations, misconceptions, and related terminology

A recurrent misconception is to equate bubble meshing with interface-fitted remeshing. The FEM/level set formulation for two-fluid flow does not fit the mesh to the bubble boundary at all; the interface is immersed in a fixed, globally unstructured mesh, and accuracy is improved by high-order approximation rather than by local mesh adaptation (Doyeux et al., 2012). Conversely, the fracture-driven mud model is explicitly front-fitted and adaptively graded near the bubble front (Katsman, 26 Sep 2025). The literature therefore supports no single assumption about whether bubble meshing is fixed-grid, front-tracking, or remeshing-intensive.

A second misconception is that bubble meshing is purely geometric. In the eight-bubble compaction study, the geometric construction is coupled to mechanical equilibrium and volume conservation (Bevilacqua, 2020). In the LEFM-based growth model, mesh updates are directed by stress intensity factors and crack-growth criteria (Katsman, 26 Sep 2025). In the BEM study, mesh stabilization is inseparable from the numerical treatment of potential-flow dynamics (Gumerov et al., 2019). Even the improved triangulation method is “physically-based,” not merely combinatorial (Shangyu et al., 7 Aug 2025).

A third point of clarification concerns terminology. “Bubble meshing” should not be confused with “bubble-enriched” finite elements. In "Bubbles Enriched Quadratic Finite Element Method for the 3D-Elliptic Obstacle Problem" (Gaddam et al., 2016), the bubble function

$0.495\text{s}$1

is an element-wise enrichment on tetrahedra for an obstacle problem. It vanishes on $0.495\text{s}$2, is positive in the interior, and supports optimal convergence and a posteriori estimation, but it is not a bubble-interface meshing strategy.

The current arXiv record suggests a field organized around several stable motifs: conformal parameterization for efficient surface triangulation, immersed representations for moving interfaces, graded front refinement for fracture-coupled bubble growth, spectral filtering for long-time surface-mesh stability, and dual or conformal constructions for periodic bubble networks. The shared objective is consistent across these variants: to control geometry resolution, numerical stability, and computational cost in problems where bubbles are themselves the dominant geometric objects.