1 Introduction↩︎

Bubbly flows, which consist of gas bubbles that rise through a liquid, occur in a wide range of natural and industrial systems. They play a central role in chemical and biochemical reactors, metallurgical processes, wastewater treatment, and many other applications. In bubble swarms, rising bubbles influence each other’s trajectories, generate turbulence, induce macroscopic circulation patterns in the surrounding liquid and affect mass, momentum and energy exchange between phases. A detailed understanding of bubble motion, deformation, and interactions within swarms is essential to design efficient reactors, improve scale‑up methodologies, and advance predictive simulation tools.

Obtaining such insights experimentally remains extremely challenging due to optical occlusion and dense multiphase interactions, especially at higher gas fractions. For this reason, significant emphasis has been placed on numerical approaches, in particular Direct Numerical Simulations (DNS), with modern interface‑resolving methods such as Volume‑of‑Fluid (VoF) [1]–[4] and Front‑Tracking (FT) [5]–[7]. The FT method explicitly tracks deformable bubble interfaces and avoids the (unphysical) numerical coalescence present in VoF schemes. FT has been demonstrated to reproduce complex bubble interaction behavior, for instance wobbling bubbles and path instabilities [8], [9], bubble–bubble interactions and swarm statistics [10]–[13], turbulence generation [14], and providing physically consistent drag forces [15] and mass transfer rates [7], [16] over a wide range of gas fractions.

While Direct Numerical Simulations (DNS) provide an exceptional level of detail, their applicability is limited to small scales due to computational cost. Larger-scale simulations therefore rely on Euler–Lagrange (EL) models, which describe the liquid phase through detailed closures but approximate the bubble phase using time-averaged correlations for drag, lift and virtual mass derived from DNS or experimental data. This averaging collapses rich information on bubble-to-bubble variability, history effects, transient deformation, and local interactions into effective coefficients, eliminating variance in rise velocity, morphology, and trajectories.

This loss of information motivates data-driven approaches that aim to retain the fidelity of DNS while enabling scalable modeling. However, a critical bottleneck in this direction is the lack of high-quality, openly available datasets that capture transient, three-dimensional bubble dynamics at the level of individual bubbles within swarms. Although DNS produces such data, it is rarely disseminated in reusable form due to its size and complexity. To address this gap, we introduce a compact, high-fidelity dataset, BubbleSH, of deformable bubbles rising in swarms, derived from interface-resolving FT DNS. The data set provides a structured, time-resolved representation of bubble motion and deformation, while remaining lightweight enough for downstream analysis and learning.

As such, BubbleSH provides a unique opportunity for developing machine learning models for dynamical systems that combine n-body interaction with deformable surfaces. Interacting-particle benchmarks are commonly used for the development and evaluation of geometric deep learning models. These datasets typically consist of multiple objects whose dynamics are coupled through pairwise and collective interactions [17], [18]. In contrast, BubbleSH includes deformable objects that evolve in a turbulent flow. Here, the learning problem includes both coupled trajectories and time-varying surface geometry. Therefore, BubbleSH combines challenges from two established classes of benchmarks: first, point-particle datasets in which objects are rigid and do not carry shape information, and second, continuum simulation benchmarks, where deformation is represented by discretizing the domain into particles or a mesh (e.g., GNS [19], MeshGraphNets [20] and LagrangeBench [21]). BubbleSH is further unique, as the spherical nature of the bubbles lends for effective non-discretized representation of the shape of the surfaces as a truncated spherical harmonic expansion. Finally, the sensitive nature of bubbly flows due to the turbulent properties of the fluid with which they interact motivates the use of probabilistic modeling [22]. These properties make BubbleSH a natural benchmark for data-driven emulators.

In this paper, we introduce a new dataset of deformable bubbles that rise in swarms, derived from high‑resolution FT DNS of air–water systems. The trajectories capture the time resolved evolution of bubble swarms under periodic boundary conditions, where each bubble is parameterized by its positions, velocity, and the coefficients of the spherical harmonic orbital functions. Alongside the dataset, we provide analysis and benchmarking tools for downstream model development. For each simulation configuration, we characterize the resulting dynamics through statistics of bubble motion, shape evolution, and inter-bubble interactions. These statistics are then used to define evaluation metrics for learned emulators, combining trajectory and shape-level errors with distributional comparisons that reflect the stochastic nature of bubble-swarm dynamics. Additionally, as an initial benchmark, we have trained a probabilistic generative forecasting model, positioning BubbleSH as a benchmark for probabilistic geometric modeling of deformable multiphase dynamics.

2 Related work↩︎

Existing bubble datasets are concentrated primarily in the image domain. Early work such as BubGAN introduced large collections of synthetic labeled bubbly-flow images for benchmarking detection and segmentation methods [23], while subsequent studies used experimental or semi-synthetic image corpora to develop deep-learning pipelines for bubble detection, mask extraction, and tracking in two-phase flows and boiling systems [24], [25]. More recent work has moved toward 3D bubble reconstruction and tracking: Hessenkemper et al. introduced a semi-artificial dataset for 3D detection and tracking of deformable bubbles in swarms, and Yu et al. released 3DBubbles, an experimental dataset containing reconstructed 3D bubble structures together with corresponding 2D projections [26], [27]. These datasets are oriented primarily toward image analysis, reconstruction, or tracking.

A second line of related work concerns the use of spherical harmonics for compact representation of 3D surfaces. They have been used for surface parametrization and geometric modeling of closed shapes [28]–[31]. In bubbly-flow research, they have been applied to characterize oscillation modes and interface distortions [32]–[36]. These works demonstrate that SH coefficients provide a compact and physically meaningful description of bubble morphology, particularly through lower-order modes associated with large-scale deformations. However, prior work has primarily used SH as a tool for analysis or static reconstruction of individual shapes. Their use as the underlying state representation in time-resolved datasets remains limited. Recent efforts such as 3DBubbles [27] move toward data-driven applications, but remain focused on experimental reconstruction rather than the dynamics of interacting bubble swarms derived from DNS.

3 Data Generation↩︎

The simulations underlying the present dataset are performed with a three‑dimensional FT method for incompressible two‑phase flow, which explicitly resolves deformable gas–liquid interfaces.

3.1 Front-Tracking Simulations↩︎

The flow field is governed by the incompressible Navier–Stokes equations written in a one‑fluid formulation, in which both phases are described on a single Eulerian grid and the interface effects enter as localized source terms. The governing equations read

\[\rho(\mathbf{x})\left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\cdot\nabla \mathbf{u}\right) = -\nabla p + \nabla\cdot\left[\mu(\mathbf{x})\left(\nabla \mathbf{u} + \nabla \mathbf{u}^T\right)\right] + \rho(\mathbf{x})\mathbf{g} + \mathbf{f}_\sigma, \;\;\;\nabla\cdot\mathbf{u} = 0\]

where \(\mathbf{u}\) is the velocity field, \(p\) the pressure, \(\rho(\mathbf{x})\) and \(\mu(\mathbf{x})\) the spatially varying density and viscosity, \(\mathbf{g}\) gravity, and \(\mathbf{f}_\sigma\) the surface tension force acting at the gas–liquid interface. The equations are discretized on a staggered Cartesian grid using a finite‑volume approach. A velocity projection method with iterative pressure correction is employed to enforce incompressibility. All linear systems arising from the discretization are solved with an incomplete‑Cholesky conjugate‑gradient (ICCG) solver, providing robustness and efficiency for the large three‑dimensional problems considered.

The gas–liquid interface (front) is explicitly represented (tracked) by a moving triangular surface mesh, whose nodes are advected with the local fluid velocity. Velocities at the nodes are obtained by piecewise spline‑based interpolation from the Eulerian grid, and time integration of the interface position is performed using a fourth‑order Runge–Kutta (RK4) scheme. This explicit tracking yields an accurate description of interface deformation dynamics and curvature, which is essential for reliable surface tension calculations.

To maintain a high‑quality surface representation during long transient simulations, the interface is continuously remeshed. The primary operation is smoothing, which redistributes marker points so that triangles remain nearly equilateral and the accuracy of curvature and surface tension forces is preserved. Mesh quality is further controlled through local topological operations, such as node addition, node removal, and edge swapping, based on local edge length and curvature criteria. Any local change that affects the enclosed volume is immediately corrected by volume‑preserving node repositioning, thereby ensuring both local and global conservation of bubble volume throughout the simulation.

The FT code simulates multiple bubbles in a single, fully periodic domain, effectively representing an infinite bubble swarm without wall effects. As bubbles approach, the intervening liquid film is naturally resolved and drains dynamically under hydrodynamic forces. Unlike front‑capturing methods such as Volume‑of‑Fluid or Level‑Set techniques [4], the FT method does not allow automatic interface merging; coalescence is therefore prevented by design, enabling bubble-swarm interactions to be studied without artificial merging. If numerical operations cause interfaces to intersect, these events are detected explicitly, and corrected by locally conservative node repositioning, again preserving bubble volume.

Surface tension forces are computed directly from the triangular mesh geometry and mapped onto the Eulerian grid as localized force densities. The exact interface position also allows an analytical determination of the local phase fraction in each computational cell, from which the spatially varying density and viscosity fields are constructed using appropriate averaging in mixture cells. This tight coupling between the Lagrangian interface description and the Eulerian flow solver is a defining feature of the FT approach and underpins its accuracy for deformable multiphase flows.

At prescribed output intervals, typically of order 1E-4 s, the code records detailed bubble dynamics, including the full surface mesh of each bubble, centroid positions, translational velocities, and additional geometric and hydrodynamic quantities. These high-fidelity outputs form the basis of the dataset presented in this work.

3.2 Spherical Harmonics Representation↩︎

The raw FT output contains a detailed triangulated interface for every bubble at every recorded time step. A single bubble mesh typically contains 1.1E4 to 1.5E4 nodes (2.2E4–3.0E4 triangular cells), distributed over the surface, depending on the amount of deformation. Although this mesh-based representation is geometrically detailed, meshes are storage-intensive and inconvenient for downstream statistical analysis and tasks, which typically require compact, fixed-dimensional descriptors. In order to reduce the footprint with several orders of magnitude, the triangulated bubble surface meshes are projected onto a basis of spherical harmonics.

Spherical harmonics provide a smooth and orthonormal basis on the sphere, and are well suited to bubble interfaces, which are typically near-spherical but may exhibit moderate deformations such as ellipsoidal or wobbling disk-like shapes. The bubble surface is described in a local spherical coordinate system centered at the bubble centroid. Let \((\theta,\phi)\) denote spherical coordinates, where \(\theta \in [0,2\pi]\) is the azimuthal angle and \(\phi \in [0,\pi]\) is the polar angle, and let \(r(\theta,\phi)\) be the radial distance from the centroid to the interface in direction \((\theta,\phi)\) [27]. The surface is approximated by a truncated spherical harmonic expansion, \[r(\theta,\phi) \approx \sum_{\ell=0}^{L}\sum_{m=-\ell}^{\ell} c_{\ell m}\, Y_{\ell}^{m}(\theta,\phi),\] where \(Y_{\ell}^{m}\) are spherical harmonic basis functions given by Equation 1 and \(c_{\ell m}\) are the associated coefficients. \[\label{shbasis} Y_{\ell}^{m}(\theta,\phi) = \sqrt{\frac{(2\ell+1)(\ell-m)!}{4\pi(\ell+m)!}} \,P_{\ell}^{m}(\cos\phi)\,e^{im\theta},\tag{1}\] where \(P_{\ell}^{m}\) denotes the associated Legendre polynomial [27]. For a truncation order \(L\), the representation yields \((L+1)^2\) coefficients per bubble. Lower-order modes capture the dominant large-scale deformations of the interface, while higher-order modes encode progressively finer surface structure. The representation assumes that the bubble is star-like with respect to its centroid, so that each ray emanating from the centroid intersects the interface only once. This excludes extreme cases such as toroidal or highly concave interfaces, but is appropriate for the class of deformable bubbles considered in the present simulations.

The coefficients are computed independently for every bubble at every recorded time step. The vertices of each triangulated interface are first expressed relative to the bubble centroid and converted from Cartesian to spherical coordinates, producing sampled values of \(r(\theta,\phi)\) over the interface. A truncated spherical harmonic expansion is then fitted to these samples in a least-squares sense, and the resulting coefficient vector is stored as the shape descriptor of the bubble. This representation preserves the essential geometry of deformable bubble interfaces while reducing the dimensionality of the original mesh-based simulation.

4 Dataset and Metrics↩︎

4.1 Data Content↩︎

r0.35

1pt

The dataset is organized as a collection of simulations of deformable bubble swarms in an air–water system rising under the influence of gravity. Each simulation corresponds to a fixed combination of bubble diameter \(d\) and gas volume fraction \(\varepsilon\), and records the time-resolved evolution of all bubbles in a cubic periodic domain. Periodic boundary conditions are applied in all three spatial directions, so bubbles exiting one side of the domain re-enter through the opposite side. See Table 1 for more simulation details.

r0.35

3.5pt

The parameter space spans three bubble diameters, 4, 5, and 6 mm, and eight gas volume fractions ranging from 5% to 40% in increments of 5%, yielding 24 parameter configurations. For each configuration, the simulation tracks a swarm of \(N=32\) bubbles over time. The domain size \(L_\mathrm{box}\) is determined by \(\varepsilon\) and bubble radius, \(L_{box}=(N \frac{4}{3} \pi r^3_{bub} \;/ \;\varepsilon)^{\frac{1}{3}}\), where \(r_{\mathrm{bub}}=d/2\). \(L_\mathrm{box}\) is in the range of 6 to 41 mm.

At each time step \(t\), and for each bubble \(i\), the dataset contains the centroid position \(\mathbf{x}_i(t) \in \mathbb{R}^3\), velocity \(\mathbf{v}_i(t) \in \mathbb{R}^3\), and shape descriptor \(\mathbf{c}_i(t) \in \mathbb{R}^{225}\)   represented by a fixed set of spherical harmonic coefficients. We use up to order \(L=14\) of the spherical harmonics basis functions, yielding \((14+1)^2=225\) coefficients per bubble. This truncation order was selected such that the mean absolute percentage error (MAPE) on the reconstructed surface area remains below \(0.1\%\), while the reconstruction error decreases rapidly for increasing \(L\) [27]. The FT simulator represents each bubble shape using \(3.9\times10^3\) floats, compressing these to 225 spherical harmonic coefficients yields a compression ratio of \(\sim 173\times\). These quantities describe both the translational motion of each bubble and its time-varying surface deformation.

Each trajectory covers between 1.5 s and 2 s of simulated time. The temporal resolution depends on the parameter configuration and is set to either 0.1 ms or 1 ms. The initial 0.2 s are discarded to eliminate startup transients, resulting in a usable duration of approximately 1.3 s–1.8 s. Table 2 summarizes the main characteristics of the dataset and the parameter ranges considered. Visualizations are provided in Appendix 7.6.

4.2 Dataset Statistics↩︎

To characterize the dynamics captured in BubbleSH, we analyze kinematic, morphological, and interaction statistics. These descriptors reveal how bubble behavior varies with diameter \(d\) and gas volume fraction \(\varepsilon\), and form the basis for the evaluation metrics introduced in Section 4.3.

Kinematics We look at bubble centroid velocity magnitude \(||\mathbf{v}(t)||=||\dot{ \mathbf{x}}(t)||\) to analyze the system’s kinetic energy, we use the rise angle \(\theta_z=\arccos(v_z(t) / ||\mathbf{v}(t)||)\), to quantify buoyancy-to-drag dominance and the tortuosity \(\tau=\sum_{t=0}^{T}{||\dot{\mathbf{x}}(t)||} \;/ \;||\mathbf{x}(T)-\mathbf{x}(0)||\), which is a geometric property describing how "twisted" a bubble trajectory is, expressed by the ratio of arc length to displacement, to investigate movement efficiency.

Morphology Shape information of the bubble interfaces is expressed using three quantities, of which surface area \(A\) is the first. The sphericity, \(\Psi=\pi^{\frac{1}{3}}(6V)^{\frac{2}{3}} / A\), is the ratio of surface area \(A\) of a sphere with the same volume \(V\) to the object’s surface area, measuring how closely the bubble’s shape resembles a perfect sphere. We also look at the amount of change of the surface over time by transforming SH coefficients \(\mathbf{c}\) back into point clouds and calculating the radial point displacements over time, while taking the mean over the displacements, to get surface velocity at each timestep, \(\langle | \dot{r}(\theta, \phi) | \rangle\).

Interactions The motion and deformation of bubbles in close proximity are strongly influenced by hydrodynamic interactions, primarily arising from drag forces [15]. We quantify spatial organization using the radial distribution function \(g(r)\), which describes how bubble density varies with distance from a reference bubble centroid. To characterize velocity interactions, we compute the angle between velocity vectors for bubble pairs within a distance of \(4d\), \(\theta_{i,j}=\arccos(\mathbf{v}_i \cdot \mathbf{v}_j / ||\mathbf{v}_i|| ||\mathbf{v}_j||)\).

The empirical probability mass functions of the eight mentioned properties from dataset configurations \(d\in \{4,5\}\) and \(\varepsilon\in \{ 15,30 \}\) are plotted in Figures 3 and 4, using \(d\) as unit of length. Looking at the densities and their medians in Figure 3, represented by the dotted lines, we see notable shifts from \(\varepsilon=15\%\) to \(\varepsilon=30\%\) across all descriptors. Bubbles are slower, rise at a less vertical angle and move in more tortuous paths when increasing \(\varepsilon\) as a result of increased drag forces [37]–[39]. This difference is more apparent for smaller bubbles [15]. Examining the morphology densities on the right side of Figure 3, we note that bubbles have lower surface area, are more spherical and have a higher rate of interface deformation at higher \(\varepsilon\). These effects are again dampened at \(d=5\). In terms of surface velocity, we notice an increase in deformation for \(d=4\) which is not present at \(d=5\) when raising gas density. The interaction statistics are shown in Figure 4. A peak in bubble density at distance \(r\approx2d\) can be observed across all regimes, while the median depends inversely on \(\varepsilon\) as expected. Lastly, the pairwise angles between velocity vectors increases with \(\varepsilon\), indicating that more streamlined upwards motion is replaced by inefficient horizontal motion.

4.3 Evaluation Metrics↩︎

Measuring the performance of a model is often done using displacement errors in the field of trajectory forecasting and point-based distance metrics and volumetric measures in the scope of shape prediction. These standard metrics like Average Displacement Error, Chamfer Distance and Intersection over Union, penalize deviations from individual ground truth trajectories or shapes without regard for predicted collective physical behavior of the system. Distribution-based evaluation additionally measures whether a model captures the statistical structure of dynamically meaningful quantities, ensuring that predictions are physically consistent at the population level rather than merely geometrically close to specific realizations. We advocate that density-based metrics of relevant domain-informed quantities are crucial in a machine learning evaluation pipeline, especially for probabilistic generative models, designed to sample multiple different possible predictions. For stochastic, turbulent systems such as bubble swarms, where multiple physically valid futures exist for given initial conditions, standard metrics would conflate diversity with model error and therefore provide an incomplete picture of model quality.

As such, in addition to displacement and volumetric errors, we measure the similarity of the predicted and ground truth distributions \(p,q\) using the 1-Wasserstein distance (Earth Mover’s Distance): \[\mathcal{W}_1(p,q)= \inf_{\gamma\in\Pi(p,q)} \mathbb{E}_{(x,y) \sim \gamma} \;|x-y|\] where \(\Pi(p,q)\) denotes the set of all couplings of \(p\) and \(q\). Unlike KL divergence, \(\mathcal{W}_1\) is a proper metric and has the natural interpretation as the minimum cost of transporting one distribution to another, making it well suited to physical quantities that may be multimodal or heavy-tailed.

We apply this to the marginal distributions of the kinematic, morphological, and interaction quantities described in Section 4.2. For each quantity, we compute \(\mathcal{W}_1\) between the empirical distribution of predicted and ground truth values aggregated across all bubbles and time steps. As these quantities have different scales, raw \(\mathcal{W}_1\) distances are normalized by the interquartile range of the ground truth distributions before aggregation, yielding a dimensionless score per quantity, insensitive to heavy-tailed distributions. A scalar benchmark metric for model comparison can then be obtained by averaging across all normalized \(\mathcal{W}_1\) distances. To anchor the reported scores, a reference baseline should always be included in the evaluation.

In addition to density-based evaluation, we provide improvements on common point-based metrics measuring errors on positional and shape data. The Average/Final Displacement Error, ADE \(=\langle ||\mathbf{x}_i(t)-\hat{\mathbf{x}}_i(t)|| \rangle\), FDE \(=\langle ||\mathbf{x}_i(T)-\hat{\mathbf{x}}_i(T)|| \rangle\), depend on both prediction horizon and spatial scale, which makes comparison across different time scales and units of measurement harder. We propose normalizing the mean displacement errors by ground truth arc length, as Relative Average/Final Displacement Error, for a dimensionless and more interpretable position-based error metric. Defined as R-ADE \(=\langle ||\mathbf{x}_i(t)-\hat{\mathbf{x}}_i(t)|| \;/ \;\sum_{t=0}^T || \dot{\mathbf{x}}_i(t)|| \rangle\), it measures the mean ratio of displacement error to distance traveled. We extend this to evaluating deformation errors over time, as the Relative Average Chamfer Distance (R-ACD), which is normalized by the total change in \(r(\theta, \phi)\) over time, as R-ACD \(=\langle \;CD(r_i(t),\hat{r}_i(t)) \;/ \;\sum_{t=0}^T \langle |\dot{r}_i(t)|\rangle \;\rangle\) with \(CD(A,B)=\frac{1}{|A|}\sum_{a\in A}min_{b\in B}||a-b||+\frac{1}{|B|}\sum_{b\in B}min_{a\in A}||a-b||\) for bubble point clouds \(A\) and \(B\) [40].

5 Evaluation↩︎

BubbleSH targets a gap in multiphase flow modeling, DNS resolve individual bubble dynamics with high fidelity but are computationally expensive, while practical large-scale Euler-Lagrange models collapse this rich per-bubble information into averaged correlations. Data-driven emulators trained on BubbleSH allow for bridging this gap, learning to reproduce coherent bubble dynamics at a fraction of the computational cost. To demonstrate this, along with our proposed metrics, we have trained a probabilistic generative forecasting model on BubbleSH that jointly predicts bubble trajectories and shape evolution.

5.1 Model↩︎

Concretely, we implemented a data-driven surrogate based on STFlow [41], a generative model using the conditional flow matching framework [42], [43]. We frame trajectory simulation as learning the joint probability distribution over future states conditioned on observed initial frames. This formulation reflects the sensitivity of bubble-swarm dynamics to turbulent interactions and local perturbations, where multiple physically valid futures may exist for a given initial conditions. Our model constructs an informed prior \(p_0\), where the observed part is kept fixed, while the unobserved future is initialized using a stochastic random-walk of all 32 bubbles over velocities \(\mathbf{v} \in \mathbb{R}^{(T-T_c) \times 32 \times 3}\) and changes in spherical harmonic coefficients \(\dot{\mathbf{c}} \in \mathbb{R}^{(T-T_c) \times 32 \times (L+1)^2}\) for observed trajectory length \(T_c\). Positions and velocities are represented in the cylindrical coordinate system to align with the translational and vertical-axis rotational symmetry of the dynamics. Only velocity, acceleration and relative position information is used to represent the centroid trajectories, making the emulator equivariant to translation.

The architecture consists of alternating spatial and temporal processing layers. Our permutation-equivariant spatial message-passing graph neural network layer is based on Neural Message Passing [44] and captures multi-bubble interactions by aggregating messages along edges defined by inter-bubble distances, operating independently at each time step. The temporal layer, implemented as a 1D UNet, models the dynamics of individual trajectories across time. Both layers share high-dimensional latent node-level embeddings as part of a spatiotemporal graph, which are updated after every layer. Separate decoders for \(\mathbf{v}\) and each SH order within \(\dot{\mathbf{c}}\) produce the outputs that denoise the prior \(p_0\) into targets. Additional model details can be found in Appendix 7.1, 7.2 and 7.3.

5.2 Experiments↩︎

We evaluate the model on the presented dataset and study its performance using both the point displacement and distributional-based metrics discussed in Section 4.3.

Data The model is trained and evaluated on a large subset of BubbleSH spanning multiple diameters and gas fractions. All trajectories are sampled at a temporal resolution of \(10^{-3}\, \mathrm{s}\), providing a consistent setting across configurations. Fixed length trajectory windows of 60 ms and \(T=30\) frames are extracted using a sliding-window approach. For each window, the first 20 ms (\(T_c=10\)) are used as the observed conditioning prefix, and the remaining portion are treated as the prediction target. The dataset is split into training, validation, and test sets with proportions of \(80\%\), \(5\%\), and \(15\%\), respectively, containing \(2094\) samples.

Metrics Following Section 4.3, for centroid position accuracy, we report the Relative ADE and FDE metrics. For assessment of the bubble deformations, we report the Intersection over Union and Relative Average Chamfer Distance. To assess whether the model reproduces physically meaningful dynamics, we also report the normalized 1-Wasserstein distance between predicted and ground-truth marginal distributions of the kinematic, morphological, and interaction quantities.

3pt

We compare four configurations: (i) a baseline where the generated trajectory is sampled from the random walk prior \(p_0\) itself; (ii) the presented model using \(L=5,(5+1)^2=36\) SH coefficients; (iii) a variant where the bubbles are modeled as being conditionally independent by removing the spatial interaction layers and (iiii) the model using \(L=3,(3+1)^2=16\) SH coefficients.

R0.35 4pt

The results in Table 3 demonstrate that the full generative model with \(L=5\) achieves the strongest overall performance, improving over the random walk prior across all metrics. In Figure 5 we plot the PMF’s calculated from predictions and targets of the \(L=5\) experiment and highlight their overlap, from which we qualitatively confirm that the learned dynamics are aligned with the data. The distributional Wasserstein scores show substantial improvement over the prior, indicating that the model captures realistic behavior at the population level rather than memorizing individual trajectories. The ablation with conditional independence performs surprisingly similar to the full model, indicating that the model does not yet learn to capture spatial inter-bubble interactions well. The \(L=3\) model produces the lowest displacement errors but higher morphology dissimilarity, confirming that truncating the spherical harmonics too aggressively limits the model’s ability to reproduce accurate deformations. Beyond predictive accuracy, we compare the speedup in simulation time by using the generative model in Table 4. The trained emulator achieves four to five orders-of-magnitude speedup over the FT simulator. This speedup, combined with the physically grounded evaluation framework, positions BubbleSH as practical stepping stone towards scalable data-driven multiphase flow modeling that retains per-bubble fidelity.

6 Conclusion↩︎

We presented BubbleSH, a dataset of transient, three-dimensional bubble swarm dynamics derived from high-fidelity interface-resolving direct numerical simulation. The dataset pairs time-resolved centroid trajectories with continuous spherical harmonic shape descriptors across multiple bubble sizes and gas volume fractions. BubbleSH sits between rigid particle dynamics benchmarks and mesh-based simulation benchmarks; the bubbles interact collectively like point particles but carry evolving shape states. As demonstrated by our generative baseline, models trained on BubbleSH can reproduce coherent bubble dynamics at four to five orders of magnitude lower computational cost than the underlying DNS, offering a practical path toward predictive swarm-scale modeling that preserves individual bubble variability. The accompanying evaluation framework provides both point-wise and distribution-based metrics for assessing trajectory accuracy and physical consistency. Together, these provide a foundation for developing and comparing data-driven multiphase flow emulators that must jointly handle complex many-body interactions and high-resolution deformable geometry.

Limitations and Future Work By simulating in a fully periodic domain, the representation of large-scale inhomogeneities, wall effects, and long-range plume dynamics is restricted. The simulations are idealized air-water systems without breakup, coalescence, surfactants, or contaminants, and therefore do not capture physicochemical effects that may be important in industrial applications. The current parameter space covers strongly deformable bubbles in water at selected \(\varepsilon\) values, corresponding roughly to Reynolds numbers of 100-1000 and Eötvös numbers of 2-5, but excludes nearly spherical, weakly deformed ellipsoidal, and spherical-cap regimes. Although spherical harmonics provide a compact shape representation, truncation can exhibit high-frequency Gibbs-type oscillations that may affect local geometric quantities. While our initial generative benchmark highlights the use of BubbleSH, the model remains limited in its treatment of symmetry and long-horizon stability. Future work could incorporate the rotational structure of spherical harmonics and study more stable rollout methods beyond fixed-window forecasting.

7 Technical appendices and supplementary material↩︎

7.1 Model Training and Inference Procedures↩︎

The training and inference algorithms, inspired from [41] are denoted below in Algorithm 6 and Algorithm 7.

7.2 Definitions of Message Passing Layer and the UNet↩︎

The spatial component of the model updates latent bubble embeddings through message passing over the bubble graph. For an edge \((i,j)\in\mathcal{E}\) at time \(t\), the message from bubble \(j\) to bubble \(i\) is computed as \[\mathbf{m}_{ij}^{t} = \varphi_{\mathbf{m}} \left( \mathbf{h}_i^{t}, \mathbf{h}_j^{t}, \mathbf{e}_{ij}^{t}, \|\mathbf{x}_i^{t}-\mathbf{x}_j^{t}\|^2, \gamma_\tau(\tau) \right),\] where \(\mathbf{h}_i^{t}\) and \(\mathbf{h}_j^{t}\) are latent node embeddings, \(\mathbf{e}_{ij}^{t}\) denotes edge features, \(\mathbf{x}_i^{t}-\mathbf{x}_j^{t}\) is the relative displacement, and \(\gamma_\tau(\tau)\) is the embedding of the flow-matching time.

Messages are aggregated over neighboring bubbles, \[\bar{\mathbf{m}}_i^{t} = \sum_{j\in\mathcal{N}(i)} \mathbf{m}_{ij}^{t},\] and the hidden state is updated with a residual MLP block, \[\mathbf{h}_i^{t\,\prime} = \mathrm{LayerNorm} \left( \mathbf{h}_i^{t} + \varphi_{\mathbf{h}} \left( \mathbf{h}_i^{t}, \bar{\mathbf{m}}_i^{t}, \gamma_\tau(\tau) \right) \right).\] Here, \(\varphi_{\mathbf{m}}\) and \(\varphi_{\mathbf{h}}\) are MLPs with SiLU activations. The layer updates latent node features using geometric pairwise information; The layer preserves translational equivariance by constructing edge features from pairwise relative displacements and distances rather than absolute coordinates.

Our temporal convolution layer follows the U-Net design of [41], [45] with the same parameters. The U-Net takes the latent node embeddings h together with a 16-dimensional sinusoidal embedding of the flow-matching time \(\tau\). Since the trajectory data also depends on physical time within the window, we add a sinusoidal frame-index embedding to \(\mathbf{h}\). The U-Net applies convolutions along the temporal axis over the T frames of each trajectory. Its output is an updated latent representation \(\mathbf{h}\), from which a two-layer \(\varphi_{\mathbf{v}}\) predicts the velocity-field update \(\Delta \mathbf{v}\).

7.3 Hyperparameters↩︎

The hyperparameters used to train and evaluate the model are included in Table ¿tbl:appendix:hyper?. We use AdamW with default parameters and a ReduceLROnPlateau scheduler based on the minimum validation loss, with a decay factor of 0.5 and patience of 60 epochs. The table lists the graph connectivity, trajectory window length, observed prefix length, and the prior-noise parameter s, which controls the additional noise in the informed prior by scaling the velocity variance estimated from the observed prefix.

7.4 Compute resources↩︎

Data Generation. Simulations were run as single CPU processes on Debian Trixie or Ubuntu 24.04 desktop/workstations with 24-core AMD CPUs. Runtime varied substantially with the gas volume fraction and simulation stability. Lower gas fractions, such as \(\varepsilon\)=5% or 10%, typically required more than 30 days to generate over 1.5 s of simulated time, while higher gas fractions generally required at least 15 days for comparable physical durations.

Model training and experiments. All model experiments were run on a single NVIDIA Tesla V100 GPU. The prior baseline required less than 30 minutes to run, while the other model configurations required approximately 10-13 hours of training depending on the spherical harmonic order and whether message passing was used. Inference was comparatively inexpensive and took only a few minutes per experiment.

7.5 Data Generation Pipeline↩︎

Multiple simulations were run concurrently on any system alongside other activities. A daemon was created to detect the simulator writing an output file containing the Eulerian fields and original bubble meshes, and would then trigger the conversion program [46]. The conversion of the bubble shapes to spherical harmonics was done ad hoc on the same machines, in parallel (up to 48 threads using OpenMP in C++, or multiprocessingPool in Python), after which the resulting coefficients were stored and the original output file was removed to save space.

7.6 Visualizations↩︎

References↩︎