Machine Learning Model for Gas–Liquid Interface Reconstruction in CFD Numerical Simulations ^†

Abstract

The volume of fluid (VoF) method is widely used in multiphase flow simulations to track and locate the interface between two immiscible fluids. The relative volume fraction in each cell is used to recover the interface properties (i.e., normal, location, and curvature). Accurate computation of the local interface curvature is essential for evaluation of the surface tension force at the interface. However, this interface reconstruction step is a major bottleneck of the VoF method due to its high computational cost and low accuracy on unstructured grids. Recent attempts to apply data-driven approaches to this problem have outperformed conventional methods in many test cases. However, these machine learning-based methods are restricted to computations on structured grids. In this work, we propose a machine learning-enhanced VoF method based on graph neural networks (GNNs) to accelerate interface reconstruction on general unstructured meshes. We first develop a methodology for generating a synthetic dataset based on paraboloid surfaces discretized on unstructured meshes to obtain a dataset akin to the configurations encountered in industrial settings. We then train an optimized GNN architecture on this dataset. Our approach is validated using analytical solutions and comparisons with conventional methods in the OpenFOAM framework on a canonical test. We present promising results for the efficiency of GNN-based approaches for interface reconstruction in multiphase flow simulations in the industrial context.

1. Introduction

Gas–liquid interfacial flows are found in a wide range of industrial applications, including cooling systems for electrical engines, chemical reactor modeling, and pore-scale flows in porous media. Among the numerous numerical methods developed to simulate such gas–liquid flows, an important issue is tracking the motion of the free interface. Two classes of methods have emerged in recent decades: front-tracking methods and interface-capturing methods. The latter are more suitable for industrial applications, as they can more easily handle topological changes in the free interface (e.g., breakup or coalescence) in complex flows. The volume of fluid (VoF) method [1,2] is an interface-capturing method that retrieves interface properties (i.e., normal, location, curvature) from the volume fraction of each fluid. However, a bottleneck in the VoF method is the cost of accurately computing the local interface curvature, which is essential for evaluating the surface tension force at the interface. Additionally, while the VoF method is accurate on regular structured grids, it introduces unphysical fluid motions (Harvie et al. [3]) when applied to unstructured tetrahedral grids.

Several attempts to address these shortcomings based on machine learning algorithms have been made recently. Qi et al. [4] proposed a two-layer neural network to predict the interface curvature on a 2D Cartesian grid. The model was trained on a synthetic dataset to predict the curvature for a given volume fraction in neighboring cells. This model was then implemented in a CFD solver and evaluated on test cases with moving interfaces. The prediction was accurate enough to consider this approach a valid solution. Patel et al. [5] extended this approach to 3D Cartesian geometry using a two-layer neural network. They trained their model on a synthetic dataset generated from a spherical surface and evaluated its performance on analytical and CFD-based test cases. Their approach outperformed the convolution method, and even matched the accuracy of state-of-the-art VoF methods such as the height function method. Svyetlichnyy [6] also trained a neural network to predict the interface properties for 2D and 3D surface reconstruction. Their model showed good performance in predicting normals, but not in predicting curvature and surface locations. While these results are promising, they are limited to Cartesian grids due to the fixed size input restriction of the machine learning algorithms used by the authors.

Unfortunately, most industrial problems involving multiphase flow simulations are discretized on unstructured grids, as this enables efficient handling of complex geometries with irregular shapes. An overview of the algorithms used to generate such grids is provided in Owen [7]. As described in Katz and Sankaran [8], persistent challenges remain due to the poor quality of the generated meshes. Finally, due in part to the costs associated with interface reconstruction, the application of VoF methods is generally restricted to simulations within relatively small physical domains. For complex chemical processes that require simulating flows with a large number of bubbles, as described in [9,10], the VoF method is not feasible for the entire simulation. Due to its accuracy, the VoF method is sometimes employed in conjunction with other methods that are better suited for large-scale simulations. In multiscale coupled strategies, it typically helps to accurately determine local parameters that model the impact of bubbles (treated as particles) on the overall flow dynamics.

Table 1. Summary of issues for the two interface reconstruction approaches.

	Structured Grid	Unstructured Grid	Computational Cost
VoF	Accurate	Unstable	Expensive
Neural Network	Accurate	Not available	Improved

presents a comparative analysis of various interface reconstruction approaches, focusing on their accuracy, stability, and computational cost. Accuracy is assessed based on the ability to compute exact solutions, while stability is evaluated by the algorithms’ performance irrespective of mesh quality. The computational cost is considered in terms of the total expense incurred during each timestep calculation in the Volume of Fluid (VoF) method.

In this work, we propose a graph neural network (GNN)-based algorithm designed to enhance the accuracy of interface reconstruction on unstructured meshes. To reach this objective, we also propose a novel data generation process that is able to generate surfaces in three dimensions that are more representative of CFD simulations, such as deformed grid cells intersecting highly curved surfaces. We subsequently train the GNN on our generated dataset and integrate the resulting interface reconstruction model in OpenFoam.

In Section 2, we provide a foundation for the volume of fluid (VoF) method. Next, in Section 3, we present a new GNN architecture for recovering the interface properties of each mesh element from the neighboring mesh geometry and volume fraction, generalizing the work presented in Patel et al. [5] to unstructured meshes. In Section 4, we describe how to generate a training dataset that enables the model to make accurate predictions on a wide range of unseen test data that are representative of industrial applications. Section 5 details model training and hyperparameter optimization. We then validate our methodology in Section 6 using several synthetic test cases and a canonical test simulated with the OpenFOAM framework, on which we benchmark our method to the conventional interface reconstruction method.

This paper is an extended version of our paper published in [11]

2. Volume of Fluid Method

There is a wide variety of numerical methods for solving multiphase flow problems. Most of these are based on the resolution of the Navier–Stokes (NS) equations based on the momentum and mass conservation (Equation (1)), which can be written as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial \rho \mathbf{U}}{\partial t}}+\nabla ·\left(\rho \mathbf{U}\mathbf{U}\right)=-\nabla p+\rho \mathbf{g}+\nabla ·(\mu \nabla \mathbf{U})+\rho \mathbf{S}+{\mathbf{F}}_{\mathsf{\sigma }}\hfill \\ {\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathbf{U}\right)=0\hfill \end{array}\right.$$

(1)

where $\rho$, $\mathbf{U}$, and p are the density, velocity, and pressure unknowns, $\mu$ and $\mathbf{g}$ are the dynamic viscosity and gravity, and $\mathbf{S}$ is an extra force acting on the fluid (e.g., due to a moving object). The interface between the two fluids induces a surface tension force ${\mathbf{F}}_{\mathsf{\sigma }}$, modeled as follows:

$${\mathbf{F}}_{\mathsf{\sigma }}={\sigma }_s\kappa \mathbf{n}\delta \left(x-x_s\right)$$

(2)

where $x_s$ represents the coordinates of the interface, ${\sigma }_s$, $\kappa$, and $\mathbf{n}$ are the surface tension coefficient, mean curvature of the interface, and normal direction to the interface, and $\delta$ is the Kronecker $\delta$, which ensures that surface tension is only considered for points on the interface.

The volume of fluid (VoF) method allows for the resolution of two-phase flow problems through a single-field formulation of the NS equations. Each fluid is represented in the VoF method with the characteristic function $\chi \left(x\right)$, defined as $\chi \left(x\right)=1$ if x is in the fluid and $\chi \left(x\right)=0$ otherwise. Let us consider two fluids related to a primary phase A and secondary phase B partitioning the simulation domain $\mathsf{\Omega }$. The VoF method is a Eulerian–Eulerian approach that is capable of accurately capturing the interface from the volume fraction of the primary phase in each cell, denoted by $\alpha ={\displaystyle \frac{{\int }_{\tau }{\chi }_A}{V_{\tau }}}$, where $\tau$ is a cell with volume $V_{\tau }$ of a mesh discretization ${\mathsf{\Omega }}_h$ of $\mathsf{\Omega }$. Figure 1a illustrates a bubble defined by a characteristic function ${\chi }_A$ in blue and the discretization on a 2D tetrahedral mesh with cells colored depending on the value of the volume fraction $\alpha$. Physical quantities such as fluid density, viscosity, and velocity are expressed as characteristic function-weighted sums. For instance, letting ${\rho }_A$ and ${\rho }_B$ be the respective densities of the primary and secondary phases, the fluid density is defined as $\rho ={\chi }_A{\rho }_A+(1-{\chi }_A){\rho }_B$.

The evolution of the interface implicitly defined by $\alpha \left(x\right)=0$ is tracked with the additional advection equation shown below as Equation (3).

$${\displaystyle \frac{\partial \alpha }{\partial t}}+\nabla ·\alpha \mathbf{U}=0$$

(3)

In the VoF method, the main NS equations in (1) are solved with the PIMPLE algorithm, a combination of the PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) methods (more details can be found in Devolder et al. [12]). The advection equation in (3) is solved with a numerical scheme, either MULES (Deshpande et al. [13]) or IsoAdvector (Roenby et al. [14]). MULES is designed to ensure boundedness and stability in the face of sharp gradients such as those found at fluid interfaces. It employs a flux limiter that adjusts the interface to prevent the overshoots and undershoots that can lead to non-physical values in volume fraction fields. On the other hand, the IsoAdvector method focuses on achieving a more accurate representation of the interface geometry by using an isosurface-based advection approach. This method tends to provide better interface sharpness and reduced numerical diffusion compared to traditional VoF methods.

The interface reconstruction step is carried out with the piecewise linear interface construction (PLIC) method; the interface is then uniquely determined in each cell by the interface normal $\mathbf{n}$ and the interface center M as illustrated in Figure 2. The interface normal $\mathbf{n}$ is usually calculated as the gradient of the volume fraction $\alpha$:

$$\mathbf{n}=-{\displaystyle \frac{\nabla \alpha }{\parallel \nabla \alpha \parallel }}$$

(4)

while the interface center M is determined such that the interface truncates the cell with a volume that matches the volume fraction.

Various procedures have then been developed for the estimation of the interface curvature $\kappa$ required to compute the surface tension force ${\mathbf{F}}_{\mathsf{\sigma }}$. $\kappa$ is usually defined as

$$\kappa =\mathbf{\nabla }·\mathbf{n}.$$

(5)

Currently, the convolution (CV) (Williams et al. [15], Williams et al. [16]), height function (HF) (Lörstad et al. [2]) and radial distance function (RDF) (Cummins et al. [17]) methods are the most standard and accurate approaches.

The objective of our work is to replace computation of the interface reconstruction and interface curvature with a machine learning algorithm that directly predicts these surface properties from the volume fraction in the neighboring cells and the coordinates of these cells.

3. Graph Neural Network Model

Interface reconstruction in the VOF method consists of geometric algorithms that cut each diphasic cell into two monophasic subcells and provide the center, normal, area, and mean curvature of the subcells’ interface. To tackle the computational complexities inherent in these geometric algorithms, recent efforts have been made to develop deep learning-based surrogate models. These methodologies typically involve the creation of a dataset where input features (i.e., geometric configurations defined by the local mesh and values of the volume fraction) are paired with corresponding output labels (i.e., interface normal, curvature, and position). Subsequently, a neural network model architecture is designed and the model is trained and validated to ensure its accuracy and efficiency. The objective is to produce a surrogate model that can rival the performance of the original geometric algorithms, offering a more computationally efficient alternative. To date, the methods put forth have predominantly relied on multilayer perceptrons (MLP) and convolutional neural networks (CNN); however, these architectures face limitations when confronted with irregularly-organized data such as unstructured meshes. In this section, we present an innovative approach utilizing graph neural networks (GNN). The rationale behind this choice lies in the ability to represent unstructured meshes as graphs.

3.1. GNN Architecture Design

The work by Zhou et al. [18] provides a comprehensive review of how GNNs can be used for unstructured meshes. A succinct summary of the fundamental principles underlying these GNN methods is presented in Appendix A.

3.1.1. Data Modeling

The objective of the proposed GNN approach is to predict the interface normal, curvature, center, and area for each interface cell (i.e., where $\alpha \in [0,1]$) from a graph containing the information from all neighboring tetrahedral cells. Therefore, each interface cell is associated with its own graph and interface properties, representing a single data point in the training process. In each graph, nodes correspond to neighboring cells, while edges represent cell connections through a shared face. The dimension of the graph is defined by all tetrahedral elements sharing at least one vertex with the current cell. Node features consist of the cell vertices coordinates and volume fraction. Figure 3b sketches the reconstructed surface properties for a three-dimensional tetrahedral cell, while Figure 3c shows the graph for a given interface cell in two dimensions.

In the training algorithm, all coordinates are normalized by $L=\max (L_x,L_y,L_z)$, where $L_x$, $L_y$ and $L_z$ are the maximum dimensions in the $x,y$, and z directions over all mesh cells, ensuring that each tetrahedral cell belongs to the unitary cube. The corresponding dimensionless curvature H is $H=H^{\prime }\ast L$, where $H^{\prime }$ is the original curvature and the corresponding dimensionless area of the reconstructed surface S is $S=S^{\prime }/L^2$, where $S^{\prime }$ is the original area. The coordinates of the interface center M are expressed in the barycentric system as $M=p_1{P}_1+p_2{P}_2+p_3{P}_3+p_4{P}_4$, where $P_1,P_2,P_3,P_4$ are the coordinates of each vertex of the tetrahedral cell (see Figure 3b) and $p_1+p_2+p_3+p_4=1$.

3.1.2. Neural Network Architecture

The complete architecture is illustrated in Figure 4d. It consists of three interconnected submodels, designated as NN-1, NN-2, and NN-3. The inputs to these submodels are depicted by the respective arrows directed toward each model, which outline the flow of information within the architecture. The predictive process unfolds as follows:

• NN-1 determines the interface normal and curvature.
• NN-2 identifies the center of the plane within the VOF approximation.
• NN-3 computes the area of the plane intersecting the target element.

The predictive sequence is inherently hierarchical: NN-2 relies on the outputs generated by NN-1; in turn, NN-3 depends on the results provided by NN-2. The entire architecture shown in Figure 4d is simultaneously trained by backpropagation during the training process.

NN-1 (Figure 4a) shows a graph neural network that predicts both the normal and the curvature from the input graph associated with a given interface cell. The node features are the coordinates of the cell vertices and the volume fraction $\alpha$ of the cells. More specifically, at each node we have the three-dimensional coordinates for the four vertices of a tetrahedron $\{{\mathbf{P}}_v\in {\mathbb{R}}^{12},\forall v\in \mathcal{V}\}$. Here, the coordinates of the vertices are flattened in a manner such that ${\mathbf{P}}_v=(P_1,P_2,...,P_4)=(x_1,y_1,z_1,x_2,y_2,\dots ,z_4)$, where $1,2,3,4$ are the number of the vertices of the tetrahedron. We then concatenate ${\mathbf{P}}_v$ with the volume fraction of the same node as ${\mathbf{x}}_v=\left(\alpha ,{\mathbf{P}}_v\right)\in {\mathbb{R}}^{13}$. Therefore, an input feature of NN-1 for a graph having n nodes can be denoted as ${\mathbf{X}}_1\in {\mathbb{R}}^{13\times n}$. The first two layers of NN-1 are SAGEConv layers that process the input graph and return a new graph with updated features. Node features are further processed by two linear layers. Skip connections are employed to prevent the vanishing gradient problem. The following global mean pool layer returns the graph-averaged node features. A latent representation of the initial graph is then recovered, then the model splits into two branches of linear layers for prediction of the normal $\mathbf{n}\in {\mathbb{R}}^3$ and curvature H.

NN-2 (Figure 4b) shows a simple MLP network that predicts the interface center M from the local information available in the current interface cell, including the normal $\mathbf{n}$ predicted by NN-1. We concatenate the volume fraction $\alpha$, coordinates of the vertices ${\mathbf{P}}_{cell}\in {\mathbb{R}}^{3\times 4}$, and normal $\mathbf{n}$, which is fed in the form of the dot product with ${\mathbf{P}}_{cell}$, that is, $\mathbf{n}·{\mathbf{P}}_{cell}\in {\mathbb{R}}^4$. This helps the neural network to discriminate vertices that are above or below the interface. Finally, a concatenated form ${\mathbf{X}}_2=\left(\alpha ,\mathbf{n}·{\mathbf{P}}_{cell},{\widehat{\mathbf{P}}}_{cell}\right)\in {\mathbb{R}}^{17}$ is used as input, where ${\widehat{\mathbf{P}}}_{cell}=(P_{cell1},P_{cell2},...,P_{cell4})\in {\mathbb{R}}^{12}$ is the flattened ${\mathbf{P}}_{cell}$ and $P_{cell1},P_{cell2},...$ are the coordinates of the vertices. At the output, the barycentric coordinates $p_{j=1,2,3,4}\in {\mathbb{R}}^4$ are recovered, then a final $L1$ normalization is applied on the output to enforce ${\sum }_{j=1}^4{p}_j=1$.

NN-3 (Figure 4c) shows another simple MLP that predicts the interface area from the available cell quantities ${\mathbf{X}}_3=\left(\mathbf{n}·{\mathbf{P}}_{cell},{\widehat{\mathbf{P}}}_{cell},p_{1,2,3,4}\right)\in {\mathbb{R}}^{19}$. A sigmoid function is used after the final layer. Due to the nondimensionalization by L of the coordinates, the maximum possible area is $\sqrt{2}/2$, which is less than 1, allowing us to skip normalization of the target area.

4. Dataset

Patel et al. [5] employed spherical surfaces with varying radii to create their synthetic dataset and train their model. We can observe their prediction of the normal on a sinusoidal surface, as shown in Figure 5. We observe that the largest errors are located on saddle point regions; it is likely that the model trained on a sphere-based dataset is biased and cannot generalize to non-spherical surfaces. Because a surface in three dimensions is locally characterized by its principal curvatures ${\kappa }_1$ and ${\kappa }_2$, we created a new dataset derived from paraboloidal surfaces, defined by the equation

$${\displaystyle \frac{x^2{\kappa }_1}{2}}+{\displaystyle \frac{y^2{\kappa }_2}{2}}=z.$$

(6)

Figure 6 illustrates an example of such a paraboloidal surface with curvatures of opposite signs. The dataset creation process was then adapted from Patel et al. [5] to account for sampling principal curvatures. The exact data generation process is detailed in Appendix A.2 and illustrated in Figure 7.

With this procedure, we generated a dataset with $N_{all}=700,000$ paraboloids. The distribution of the graph features of the dataset elements is shown in Figure 8. We observe that ${\kappa }_1$ and ${\kappa }_2$ follow a uniform distribution, as designed. Figure 8b shows the histogram of the volume fraction $\alpha$ of the cell corresponding to the first node of each graph. Even if this distribution were not controllable a priori in our generation algorithm, the figure indicates that quasi-uniformity was reached by our dataset generation process. Figure 8c displays the distributions of the rotations in each direction, as explained in step 6 of the data generation process. Figure 8d shows the normals at the center point of the rotated paraboloid, which appear to be evenly distributed, confirming the 3D uniformity of the rotation. The mean curvature, defined as $H={\displaystyle \frac{{\kappa }_1+{\kappa }_2}{2}}$, follows the expected triangular shape in the histogram shown in Figure 8e, as it results from the sum of two independent uniform distributions. The area histogram in Figure 8f exhibits a unimodal distribution with a long tail, which is not controllable a priori. As observed, most properties exhibit a uniform distribution; therefore, we expect such a dataset generation procedure to better represent general surfaces than a dataset solely based on spherical surfaces.

5. Model Training and Optimization

5.1. Training

We utilized a dataset comprising 500,000 graphs selected from a total of $N_{all}$ = 700,000 to train our model. The choice of this specific data size aligns with the considerations outlined in Section 5.3. The $\mathrm{500,000}$ graph dataset was partitioned into a 7:2:1 ratio for training, validation, and testing, respectively.

In a deep architecture such as the one shown in Figure 4, the backpropagation of gradients becomes less effective due to the depth of the entire architecture. At the same time, the errors from intermediate predictions are propagated to the subsequent models.

To address challenges associated with ineffective training due to model depth and mitigate error propagation from intermediate predictions, auxiliary loss functions are introduced. These functions incorporate ground truth labels (“wl”, with-label) as inputs to models responsible for computing the center and area. By addressing the vanishing gradient problem, these auxiliary losses help to prevent error compounding from intermediate predictions. In NN-2, the input normal can either be the predicted value from NN-1 (“wol”, without-label) or the ground truth label (“wl”, with-label), necessitating two distinct loss functions. Similarly, this principle applies to the input center in NN-3. Each loss function employs an L2-normalized mean squared error (MSE). The overall loss is then defined as follows:

$${\mathrm{Loss}}_{\mathrm{total}}={\mathrm{Loss}}_{\mathrm{normal}}+{\mathrm{Loss}}_{\mathrm{curv}}+{\mathrm{Loss}}_{\mathrm{center}}^{\mathrm{wl}}+{\mathrm{Loss}}_{\mathrm{center}}^{\mathrm{wol}}+{\mathrm{Loss}}_{\mathrm{area}}^{\mathrm{wl}}+{\mathrm{Loss}}_{\mathrm{area}}^{\mathrm{wol}}$$

where, ${\mathrm{Loss}}_{\mathrm{normal}}$, ${\mathrm{Loss}}_{\mathrm{curv}}$, ${\mathrm{Loss}}_{\mathrm{center}}$, and ${\mathrm{Loss}}_{\mathrm{area}}$ represent the loss functions associated with the normal, curvature, center, and area predictions, respectively.

The Pytorch-integrated adaptive learning rate was employed and configured within the range of $1\times {10}^{-2}$ to $1\times {10}^{-6}$. The learning rate was reduced by a factor of $0.9$ when the validation loss function failed to improve for five consecutive epochs. The batch size was set to 512. To prevent overfitting, the early stopping patience parameter was established at 80 epochs. Additionally, gradient clipping was implemented with a threshold of $1\times {10}^{-2}$, and the optimization was conducted using the Adam optimizer.

Figure 9 illustrates the learning curves for both total and individual train/validation losses during the training process. Notably, no improvement was discerned for the total validation loss beyond the 588th epoch, leading to the termination of training at the 588th epoch. It is observed that all types of losses converge without evident signs of major overfitting.

In terms of the center and area, the with-label loss yields a lower loss function compared to the without-label scenario. This outcome aligns with expectations, as in the absence of label input the input variables are predicted variables from another neural network, potentially introducing errors. The total and individual loss functions at the 508th epoch for both training and validation are summarized in

Table 2. Model performance at the 508th epoch in training (top) and on the test set (bottom).

Type of Loss	Training	Validation	MSE
${\mathrm{Loss}}_{\mathrm{total}}$	$7.50\times {10}^{-3}$	$9.07\times {10}^{-3}$	$9.30\times {10}^{-3}$
${\mathrm{Loss}}_{\mathrm{normal}}$	$2.30\times {10}^{-3}$	$2.39\times {10}^{-3}$	$2.41\times {10}^{-3}$
${\mathrm{Loss}}_{\mathrm{curv}}$	$4.13\times {10}^{-3}$	$5.26\times {10}^{-3}$	$5.40\times {10}^{-3}$
${\mathrm{Loss}}_{\mathrm{center}}^{\mathrm{wl}}$	$1.10\times {10}^{-4}$	$2.01\times {10}^{-4}$	$1.84\times {10}^{-4}$
${\mathrm{Loss}}_{\mathrm{center}}^{\mathrm{wol}}$	$6.66\times {10}^{-4}$	$9.76\times {10}^{-4}$	$1.05\times {10}^{-3}$
${\mathrm{Loss}}_{\mathrm{area}}^{\mathrm{wl}}$	$1.03\times {10}^{-4}$	$5.53\times {10}^{-5}$	$5.51\times {10}^{-5}$
${\mathrm{Loss}}_{\mathrm{area}}^{\mathrm{wol}}$	$1.99\times {10}^{-4}$	$1.95\times {10}^{-4}$	$2.03\times {10}^{-4}$
Predicted Variable		Median
$n_{error}=\left|\right|{\mathbf{n}}_{label}-{\mathbf{n}}_{pred}\left|\right|$		$5.61\times {10}^{-2}$
$H_{error}=\left|H_{label}-H_{pred}\right|$		$4.16\times {10}^{-2}$
$M_{error}=\left|\right|M_{label}-M_{pred}\left|\right|$		$1.98\times {10}^{-2}$
$A_{error}=\left|A_{label}-A_{pred}\right|$		$7.57\times {10}^{-3}$

. Furthermore, the table includes predictions on the test dataset generated by the trained model. All types of losses on the test dataset remain within the same range as those observed during training and validation. This outcome substantiates the model’s generalization capability beyond the training and validation datasets.

To conduct a detailed analysis of the results, the errors in the $L_1$ norm for each predicted property outlined in Table 2 (right) are depicted in Figure 10. The error histograms display a unimodal distribution centered around zero. Large errors (errors $\ge 0.2$) occur significantly less frequently than those near zero. These large errors typically arise in configurations where the volume fractions $\alpha$ are close to 0 or 1. In such cases, the center of the interface is positioned very near to one of the nodes that define the corners of the grid cells. This proximity can lead to inaccuracies in the interface representation, affecting the results. The median of the errors is also presented in Table 2.

In addition, Figure 11 displays the $R^2$ score between the predictions and their corresponding labels in the test dataset for the scalar variables, namely, the curvature and the area. The scattered points are color-coded based on the volume fraction $\alpha$ of the cell where the prediction was made. The $R^2$ value exceeds $0.9$ for both variables. Interestingly, the value of $\alpha$ does not appear to exhibit a clear correlation with the error; high and low $\alpha$ values are distributed uniformly in the figure, irrespective of whether the prediction underestimates or overestimates the label.

5.2. Hyperparameter Optimization

The complex architecture employed in our study is characterized by numerous hyperparameters, and optimizing these parameters is crucial for effective training; to facilitate this process, we utilized an automated tuning tool called Optuna Akiba et al. [19], which is a free library based on Bayesian optimization. The hyperparameter study was conducted independently for the three neural networks illustrated in Figure 12.

In Optuna, optimal hyperparameters are determined through a series of “trials’, each representing a training iteration with a specific set of hyperparameters chosen from suggested values, as detailed in

Table 3. Hyperparameter study for NN-1.

Hyperparameters	Suggested Values	Best Parameter
$s_1$	$\{5,6,7,8\}$	6
$s_2$	$\{s_1,s_1+1,\dots ,8\}$	6
$l_{2\left(i\right)}$	$\{5,6,\dots ,s_2\}(l_{2\left(i\right)}\ge l_{2(i+1)})$	$l_{2\left(1\right)}=6$
$l_{3\left(j\right)}$	$\{5,6,\dots ,l_{2\left(n_2\right)}\}(l_{3\left(j\right)}\ge l_{3(j+1)})$	$l_{3\left(1\right)}=6$
$l_{4\left(k\right)}$	$\{5,6,\dots ,l_{3\left(n_3\right)}\}(l_{4\left(k\right)}\ge l_{4(k+1)})$	$l_{4\left(k\right)}=6(\mathrm{for}k=1,2,3)$
$n_1$	$\{0,1,2,3,4\}$	1
$n_2$	$\{0,1\}$	1
$n_3$	$\{0,1,2,3,4\}$	1
$n_4$	$\{1,2,3,4\}$	3
$n_b$	$\{8,9,10,11,12,13\}$	9

,

Table 4. Hyperparameter study for NN-2.

Hyperparameters	Suggested Values	Best Parameter
$l_1$	$\{7,8,9,10,11\}$	11
$l_{2\left(i\right)}(i=1,2,\dots ,n_2)$	$\{4,5,\dots ,l_1\}(l_{2\left(i\right)}\ge l_{2(i+1)})$	$l_{2\left(1\right)}=5,l_{2\left(2\right)}=5$
$n_1$	$\{1,2,3\}$	3
$n_2$	$\{0,1,2\}$	2
$n_b$	$\{8,9,10,11,12\}$	8

and

Table 5. Hyperparameter study for NN-3.

Hyperparameters	Suggested Values	Best Parameter
$l_1$	$\{5,6,7,8,9\}$	7
$l_{2\left(i\right)}$	$\{4,5,\dots ,l_{2\left(i\right)}\}(l_{2\left(i\right)}\ge l_{2(i+1)}),l_{2\left(0\right)}=l_1$	$l_{2(1,2,3,4)}=6$
$n_1$	$\{1,2,3,4\}$	1
$n_2$	$\{1,2,3,4\}$	4
$n_b$	$\{8,9,10,11\}$	9

. Each trial’s effectiveness is gauged by the final validation loss, which serves as its score. Hyperparameters for subsequent trials are automatically selected based on these scores using a Bayesian optimization approach. In our study, we conducted 100 trials, with each trial capped at a maximum of 50 epochs. This cap represented a practical balance between optimizing performance and efficiently managing computational resources. The tables further highlight the parameters that achieved the lowest validation loss.

For NN-1, the node numbers for the first two SAGEConv layers were set to $2^{s_1}$ and $2^{s_2}$, with $2^{s_2}$ maintained in the subsequent $n_1$ linear layers to preserve the skip connection. Another $n_2$ linear layers followed with node numbers $2^{l_{2\left(i\right)}}$, where $i=1,2,\dots ,n_2$; here, $2^{l_{2\left(i\right)}}$ could vary at each layer i under the condition $l_{2\left(i\right)}\ge l_{2(i+1)}$. After a global max pooling operation, an additional $n_3$ linear layers followed with node numbers $2^{l_{3\left(j\right)}}$ ($j=1,2,\dots ,n_3$ with $l_{3\left(j\right)}\ge l_{3(j+1)}$). Following a divergence in the latent space into normal and curvature predictions, the curvature side included another block of layers, beginning with a linear layer with node number $2^{l_{3\left(n_3\right)}}+3$ and succeeded by $n_4$ layers with node numbers $2^{l_{4\left(k\right)}}$ ($k=1,2,\dots ,n_4$ with $l_{4\left(k\right)}\ge l_{4(k+1)}$).

For NN-2, the constant node number $2^{l_1}$ was established for the first layer and the subsequent $n_1$ linear layers preceding the dropout layer. Following the dropout operation, an additional block of $n_2$ linear layers was introduced, where the node number $2^{l_{2\left(i\right)}}$ could vary at each layer i under the condition $l_{2\left(i\right)}\ge l_{2(i+1)}$.

For NN-3, the constant node number $2^{l_1}$ was assigned to the first layer and the subsequent $n_1$ linear layer, followed another block of $n_2$ linear layers, where the node number $2^{l_{2\left(i\right)}}$ could vary at each layer under the condition $l_{2\left(i\right)}\ge l_{2(i+1)}$.

The impact of batch size, denoted as $N_b=2^{n_b}$, was investigated for all three neural networks; the results indicated $n_b=9$ to be the optimal parameter for two of the networks. Consequently, a batch size of $N_b=512$ was selected for training.

Considering all the optimizations discussed above, the architecture of our model was finalized as presented in Section 3.1.2.

5.3. Data Size

The size of the training dataset significantly influences the quality of the model. Generally, a larger dataset tends to result in a better model, albeit at the expense of longer training times. From an engineering standpoint, understanding the tradeoff between model quality and training time is crucial. In this study, we investigated the impact of dataset size on training by examining the results of two generalization tests. The predictive performance comparison across all variables in the two test cases is depicted in Figure 13a,b. The use of colors in these figures aids in swiftly identifying discrepancies in errors among different scenarios. The training was conducted using 50, 100, 200, 300, 400, 500, 600, and 700 thousand graphs. In both tests, we observe the expected downward trend in error with an increase in the number of graphs. However, the reduction in losses becomes less pronounced for datasets larger than $500,000$ graphs, especially in test-No.2.

Table 6. Training times of different data sizes.

Data size (in thousand)	50	100	200	300	400	500	600	700
Time (mins.)	$9.50\times {10}^1$	$1.25\times {10}^2$	$2.79\times {10}^2$	$4.90\times {10}^2$	$7.71\times {10}^2$	$1.03\times {10}^3$	$1.23\times {10}^3$	$1.41\times {10}^3$

provides a summary of training time versus dataset size. As anticipated, larger datasets result in longer training times. Considering the tradeoff between time and accuracy, the model was ultimately trained with a dataset size of $500,000$ graphs.

6. Model Validation

6.1. Generalization Test on a Sinusoidal Surface

To evaluate the generalization of the model, we conducted predictions on a 3D sinusoidal surface, as illustrated in Figure 14. The surface was defined by $f\left(x\right)=0.1sin9xcos9y$. The surrounding space was voxelized (generating a tetrahedral mesh) with a constant discretization of $\Delta =0.05$. Graphs were then generated within each of these voxels, with each voxel containing the center point of a graph, and its neighboring voxels containing other nodes of the same graph. Therefore, a voxel can act as both a container for a center point and a neighbor point.

The maximum absolute nondimensionalized curvature of the test surface was ${\left|H\right|}_{max}=0.703$, falling within the range of the training dataset ($H_{training}\approx [-0.75,0.75]$). With the volume fraction in the range $0.01<\alpha <0.99$, each interface cell was then input to the model, while interface cells with very high or low volume fractions were filtered out, as we considered these to be outside the range of applicability for our model. Figure 15 displays the errors of the prediction for the four variables and their corresponding histograms. The errors for each variable are defined in

Table 7. Median of the error for the entire surface.

Predicted Variable	Median	Color Scale
$n_{error}=\left|\right|{\mathbf{n}}_{label}-{\mathbf{n}}_{pred}\left|\right|$	$6.05\times {10}^{-2}$	$[min\left(n_{error}\right),max\left(0.3n_{error}\right)]$
$H_{error}=\left|H_{label}-H_{pred}\right|$	$4.32\times {10}^{-2}$	$[min\left(H_{error}\right),max\left(0.2H_{error}\right)]$
$M_{error}=\left|\right|M_{label}-M_{pred}\left|\right|$	$1.84\times {10}^{-2}$	$[min\left(M_{error}\right),max\left(0.1M_{error}\right)]$
$A_{error}=\left|A_{label}-A_{pred}\right|$	$6.84\times {10}^{-3}$	$[min\left(A_{error}\right),max\left(0.2A_{error}\right)]$

. It is observed that errors do not exhibit a clear correlation with the geometry for any variable. The histograms display a unimodal shape, and no abnormal behavior is observed. Table 7 also presents the median error values.

The test surface utilized in this study exhibited a nonuniform distribution about $\alpha$, as depicted in Figure 16. A higher concentration of samples is observed at $\alpha$ values close to 0 and 1, referred to as “marginal $\alpha$”. This uneven distribution may impact the prediction results. In Figure 16, the prediction is compared to the label and colored by $\alpha$. Notably, cases with marginal $\alpha$ values exhibit larger errors. For instances with small marginal $\alpha$, the model tends to predict values smaller than the label, whereas it tends to predict larger values for instances with large marginal $\alpha$. This influence of $\alpha$ was not observed in the test case during training, where $\alpha$ was quasi-uniform. This suggests limited capacity of the model for marginal $\alpha$, and should be taken into consideration when implementing the model in a real CFD solver and applying it to real-world use cases.

6.2. Validation Simulating a Two-Phase Flow Model with OpenFOAM

To validate our approach, we implemented our methodology in OpenFOAM (“Open-source Field Operation And Manipulation”), the C++ toolbox for the development of customized numerical solvers in computational fluid dynamics (CFD). We extended the IsoAdvector solver of OpenFOAM described in Roenby et al. [14] with the implemented geometrical interface reconstruction method, denoted IsoAlpha. To compute the geometrical features of interfaces, we introduced a new reconstruction method, denoted IsoAlphaML, based on the GNN model described in Section 3.

We benchmarked the two methods using a 3D study case inspired by the original 2D rising bubble study case published in 2009 by Hysing et al. [20]. This study case involves simulating a two-phase flow within a cylindrical domain. The simulation domain is discretized with a tetrahedral mesh using the open-source GMSH tool presented in Geuzaine and Remacle [21]. Three meshes, $Mes{h}_1$, $Mes{h}_2$, and $Mes{h}_3$, are built. The original mesh $Mes{h}_1$ of size ${\Delta }_1=0.1$ is refined twice to obtain $Mes{h}_2$ of size ${\Delta }_2={\displaystyle \frac{{\Delta }_1}{2}}=0.05$ and $Mes{h}_3$ of size ${\Delta }_3={\displaystyle \frac{{\Delta }_2}{2}}=0.025$. The primary goal of this study is to assess the impact of the mesh discretization size on the finite volume schemes employed across various models. Additionally, we seek to illustrate the performance of our method with different mesh sizes. Notably, as the grid cell size diminishes, the relative curvature of interfaces tends to decrease to lower values after normalization to the reference root cell of unit size. This significant reduction may trigger the emergence of unexpected phenomena.

The two phases are modeled as Fluid1 and Fluid2 with their specific dynamical properties (density and viscosity). Initially, the entire domain is occupied by Fluid1 except for the inside of a spherical bubble occupied by Fluid2. We impose boundary conditions as described in Figure 17.

We denote the height of the domain $\mathsf{\Omega }$ by L and the radius of the initial spherical bubble by $r_0$. We performed different simulations with a constant advected velocity $U_0$ along the z-axis imposed within the domain and with a surface tension coefficient denoted $\sigma$. These simulations can be classified and characterized by the scaling factors ${\displaystyle \frac{L}{r_0}}=2$, ${\displaystyle \frac{L}{U_g}}$, where $U_g=\sqrt{gL}$ is the gravitational velocity, the Reynolds number $Re={\displaystyle \frac{{\rho }_1{U}_{0}L}{{\mu }_1}}$ describes the ratio of inertial to viscous effects, and the Eötvös number $E_0={\displaystyle \frac{\Delta \rho g{L}^2}{\sigma }}={\displaystyle \frac{\Delta \rho {U_g}^{2}L}{\sigma }}$ describes the ratio of gravitational forces to surface tension effects.

The surface tension force plays a crucial role in these models, making its correct discretization necessary for achieving accurate simulation results. As seen in Section 2, computation of the surface tension force based on the continuous surface force model (2) relies on approximation of the mean curvature $M={\displaystyle \frac{{\kappa }_1+{\kappa }_2}{2}}$. We compared our IsoAlphaML approach, where M is predicted with our GNN-based model, to the IsoAlpha method, where M is computed by discretizing Equations (5) with a finite volume scheme, as explained in Section 2.

The study case was designed in three progressive test cases:

• Test-1: The Stationary Droplet test case, in which Fluid1 and Fluid2 have the same dynamic properties and $U_0=0$. Nothing should happen. This test case aims to detect residual streams that may appear due to errors in the curvature computation.
• Test-2: The Translating Droplet test case, in which the two fluids have the same properties. The imposed velocity is different from 0. The bubble should be translated without deformation.
• Test-3: The Original 3D Hysing test case, in which the properties of the fluids are different. Fluid1 represents the liquid phase, while Fluid2 represents the gas phase. No velocity is imposed, and $U_0=0$. The bubble is supposed to rise with the gravity effects. The bubble is deformed under the effect of surface tension at the gas–liquid interface, which depends on the interface curvature.

In

Table 8. Property of Fluid1 and 2.

	Fluid1		Fluid2		Simulation Parameters
	Density	Viscosity	Density	Viscosity	$\mathsf{\sigma }$	${\mathit{U}}_{\mathbf{0}}$	${\mathit{U}}_{\mathit{g}}$	$\mathit{Re}$	${\mathit{E}}_{\mathbf{0}}$
Test-1	1000	10	1000	10	1	0	0	0	0
Test-2	1000	10	1000	10	1	1	0	50	0
Test-3	1000	10	1	0.1	1.96	0.7	0.7	35	125

, we gather the fluids properties and simulation parameters of these three test-cases.

We first explore the distribution of the volume fraction generated by the OpenFOAM simulations to check whether it conforms to one of the datasets used to train our GNN model. Figure 18 shows the histogram of the volume fraction $\alpha$ of the diphasic cells computed during the simulation of Test-1 on the left, while the right shows the histogram of the predicted normal error $n_{error}$ on the spherical interface regarding the expected normal values, which can be easily computed. We can observe that many samples are found in the “marginal $\alpha$” range. The distribution of the $n_{error}$ is mono-modal and its peak is located near its median. These tendencies are similar to those observed in the previous generalization sinusoidal test case.

Figure 19 compares the evolution of bubble circularity over time for each test case when using the IsoAlphaML method and the IsoAlpha method. The similarity of the curves at first glance indicates that the simulation results are very similar. As expected, the bubble circularity does not change over time for Test-1 and Test-2. In Test-3, the evolution of bubble circularity using the two methods is very similar.

To better analyze the difference between the mean curvature computation of the two methods, we focus on the results of Test-1. In this test, the droplet should remain stationary. The mean curvature should also remain constant and equal to the reference value $M_{ref}={\displaystyle \frac{1}{R}}$, where R denotes the initial radius of the droplet. We have introduced two variants ${\mathit{IsoAlphaML}}_{M_{ref}}$ and ${\mathit{IsoAlpha}}_{M_{ref}}$ of the previous methods, setting the mean curvature to the reference value $M_{ref}$. The reference normal value of the interface in each diphasic voxel is also easy to evaluate using the voxel center and the stationary droplet center. We evaluate the errors in the curvature computation and use the fact that these errors impact the value of the surface tension force, which leads to the appearance of nonzero residual velocities. Then, computing the overall kinetic energy is a simple way to evaluate the accuracy of the curvature computation method. The errors regarding the reference solution of the two approaches IsoAlphaML and IsoAlpha are gathered in

Table 9. Median of the predicted/calculated value.

Variable	IsoAlphaML	IsoAlpha
$n_{error}=\left|\right|{\mathbf{n}}_{label}-{\mathbf{n}}_{pred}\left|\right|$	$5.94\times {10}^{-2}$	$1.33\times {10}^{-1}$
$M_{error}=\left|\right|M_{label}-M_{pred}\left|\right|$	$1.64\times {10}^{-2}$	$2.38\times {10}^{-2}$

. The histogram of curvature values, the normal and curvature errors are plotted in Figure 20 and  Figure 21. Figure 22 plots the time evolution of the global residual kinetic energy for the IsoAlphaML and the IsoAlpha methods along with their variants ${\mathit{IsoAlphaML}}_{M_{ref}}$ and ${\mathit{IsoAlpha}}_{M_{ref}}$. As expected, we see that the kinetic energy of the ${\mathit{IsoAlphaML}}_{M_{ref}}$ and ${\mathit{IsoAlpha}}_{M_{ref}}$ simulations remains very low. For IsoAlphaML and IsoAlpha, the kinetic energy increases with time, highlighting the effect of the curvature error on the surface tension computations. For the simulation with $Mes{h}_1$ and $Mes{h}_2$, this effect is less important for the IsoAlphaML method than for the original IsoAlpha method, while for the simulation with $Mes{h}_3$ the effect is inverted.

The errors of our GNN-based approach are reasonable. Similar to the previous generalization sinusoidal test, the results prove the capability of our method to generalize on unseen surfaces. This makes the proposed method interesting in the context of real CFD simulations, where complex unstructured meshes are required.

Notably, our approach outperforms the IsoAlpha method on both variables.

7. Conclusions and Future Perspectives

In this study, we have proposed an enhanced volume of fluid (VoF) method utilizing machine learning, specifically graph neural networks (GNNs), to expedite interface reconstruction on general unstructured meshes. Our approach involves the generation of a synthetic dataset comprising paraboloid surfaces discretized on unstructured meshes. Comprehensive data exploration was conducted to understand the characteristics of the dataset. Subsequently, we introduced a GNN-based method for computing interface reconstruction on general unstructured meshes. Extensive hyperparameter optimization was performed, focusing on the parameters of our GNN architecture and the dataset size.

We validated the efficacy of our new method through generalization tests. The results demonstrate that GNNs provide a viable alternative to conventional surface reconstruction methods, proving effective even on unstructured grids. This is in contrast to other approaches, which are limited to structured meshes. Importantly, our method showcases superior accuracy compared to traditional geometrical methods on tetrahedral meshes, making it particularly appealing in industrial contexts.

Future work will focus on efficient integration of GNN models in CFD simulators such as OpenFoam. To observe significant computational gains, extensive work is required in order to couple traditional numerical simulators, which are usually parallelized on CPU cores, with ML models, which are most efficient when evaluated using GPU accelerators. This endeavor will contribute to a more comprehensive understanding and application of our proposed GNN-enhanced VoF method in practical engineering scenarios.