Validation of HiG-Flow Software for Simulating Two-Phase Flows with a 3D Geometric Volume of Fluid Algorithm

1. Introduction

Interfacial multiphase flows play a crucial role in a wide range of industrial fields, offering opportunities for enhanced process efficiency, product quality, and environmental sustainability. Understanding and manipulating the dynamics of multiple immiscible fluids in contact at their interfaces are key challenges that require accurate numerical computations. These systems exhibit complex behaviors due to the interplay between bulk fluid coupling and surface effects, necessitating advanced modeling and simulation techniques [1]. The ability to accurately predict and control interfacial phenomena is of paramount importance in industrial processes such as oil and gas extraction [2], chemical reactions [3], pharmaceutical manufacturing [4], food processing [5], and microfluidics [6], among others.

Interfacial multiphase flows exhibit complex and dynamic behavior, requiring sophisticated numerical techniques for their accurate calculation and analysis. Over the years, various numerical methods have been developed to tackle the challenges associated with capturing and tracking fluid interfaces [7].

One widely employed approach is the front-tracking method [8,9], which explicitly tracks the fluid interface using a moving mesh. This method accurately captures the detailed dynamics of the interface and enables precise representation of surface effects. An example of a front-tracking method is the Arbitrary Lagrangian-Eulerian (ALE) method [10]. These methods are very good in providing an accurate description of the contact surface but often come with high computational costs.

Another approach is the fixed-mesh method, which represents the fluid interface implicitly using additional fields. The Volume of Fluid (VOF) method [11], Phase-Field method [12], and Level-Set method [13] are notable examples. These methods use a fixed mesh and provide easier coupling formulations and smaller algebraic problems. However, they may sacrifice some accuracy in the description of the interface compared to front-tracking methods. In the VOF method, a mass conservation equation is solved to track the movement of the interface by computing the volume fraction of each fluid. This method is effective in capturing large-scale interface dynamics and is widely used in various industrial applications. The Phase-Field method introduces an additional scalar field to represent the interface, evolving according to a diffusion equation. This approach offers a smooth transition across the interface and is particularly useful for simulating complex interfacial phenomena such as phase separation and morphological changes. The Level-Set method evolves a scalar function called the level set function, which tracks the interface as it evolves in time. This method provides accurate interface capturing and is suitable for simulating topological changes, such as merging and breakup of droplets. Furthermore, the Immersed Boundary Method (IBM) [14] represents the interface as a boundary immersed in a fixed grid. It uses a force-based approach to model the interaction between the fluid and the boundary, enabling the simulation of complex geometries and multiphase flows involving solid boundaries.

Continued advancements in numerical methods hold the potential for further enhancing our understanding and control of interfacial multiphase flows. Specifically, the modeling and simulation of three-dimensional (3D) interfacial multiphase flows is a critical area of research that has garnered significant attention in recent years [1,15,16]. The geometric complexity and dynamic nature of the fluid interfaces demand advanced modeling techniques that can accurately resolve and track the evolving interfaces in three-dimensional space. Researchers are continuously exploring new approaches to improve conservation of mass, capture fluid interfaces, and compute interfacial surface tension. Liovic et al. [15] developed a volume tracking method for accurate calculations of interfaces in 3D geometries. The interface geometries are approximated as piecewise planar and fully multidimensional fluxes are used in a single unsplit step to advect cell volumes. The definition of these fluxes is based on backward-trajectory remapping, ensuring reliable and precise tracking of the interfaces. To assess the performance of the 3D volume tracking scheme, Liovic et al. [15] conducted 3D shearing flow and deformation field tests. Adelsberger et al. [16] conducted 3D simulations of incompressible two-phase flow to investigate the behavior of rising droplets. Three different flow solvers, namely DROPS [17], NaSt3D [18] and OpenFOAM [19] were utilized, each employing distinct numerical techniques. DROPS employed a finite element spatial discretization, an implicit theta scheme for time discretization, and a level set method for interface capturing. NaSt3D utilized a finite difference spatial discretization, a second-order Adams-Bashforth scheme for time discretization, and a level set method for interface capturing. Lastly, OpenFOAM employed a finite volume spatial discretization, an implicit Euler scheme for time discretization, and a VOF method for interface capturing. Metivet et al. [1] presented a numerical framework for simulating multiphase flows using the level-set method. The discretization of the fluid equations is done using a finite element method. For the spatial discretization the method can use finite elements of high order. This choice ensures accuracy and stability in representing the spatial domain. As for the time discretization, the implicit Backward Differentiation Formulae (BDF) of arbitrary order is used. This approach allows for precise time stepping and facilitates the accurate simulation of dynamic phenomena in interfacial multiphase flows. This ongoing research aims to enhance the fidelity and predictive capabilities of computational models, enabling more accurate simulations of 3D interfacial multiphase flows, which generally suffer from volume conservation, interface sharpness, boundedness, and shape preservation.

In this study, we present a geometric formulation of the VOF method, known as PLIC (piecewise linear interface construction), specifically designed for simulating 3D interfacial multiphase flows. The PLIC-VOF method is implemented within the HiG-Flow (hierarchical grid) framework, a recently developed Computational Fluid Dynamics (CFD) software. HiG-Flow utilizes finite differences on hierarchical grids and enables the simulation of a wide range of flow types, including Newtonian, Generalized-Newtonian, and viscoelastic flows in any dimension. The modular design of the software allows for easy implementation of new techniques and methods, making it flexible and adaptable to evolving research needs. Furthermore, HiG-Flow incorporates numerical stabilization techniques to ensure accurate and stable simulations, even in the presence of complex flow phenomena, such as 3D multiphase flows. The user-friendly interface allows researchers to choose specific simulation features, such as the dimension, modules to be used, and numerical techniques, providing a customizable and efficient simulation environment.

The paper is organised as follows: in Section 2 we present the governing equations that describe interfacial multiphase flows. Subsequently, Section 3 describes the numerical method used to discretize the governing equations using the PLIC-VOF method. In Section 4, we focus on the validation of the numerical algorithm through two 3D advection tests: the shearing flow and deformation field problems. Furthermore, the validation of the PLIC-VOF method is also conducted using a 3D multiphase benchmark problem involving the rising of a bubble. The paper ends with the conclusion in Section 5.

2. Governing Equations

This study examines the characteristics of an unsteady, laminar, isothermal and incompressible two-phase flow. The flow consists of two immiscible fluid phases, with no mass transfer occurring across their interface. The governing equations for this flow are the mass conservation equation (Equation (1)),

$$\nabla ·\mathbf{u}=0,$$

(1)

and the linear momentum balance equation (Equation (2)),

$$\rho \left[\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\left(\nabla ·\mathbf{u}\right)\right]=-\nabla p+\nabla ·\left\{\mu \left[\nabla \mathbf{u}+{(\nabla \mathbf{u})}^T\right]\right\}+\rho \mathbf{g}+\sigma \kappa \mathbf{n}\delta ,$$

(2)

where $\mathbf{u}=\left(u,v,w\right)$ is the velocity field of the flow, $\rho$ is the density, p the pressure, $\mu$ the absolute viscosity, $\mathbf{g}$ the gravitational field, $\sigma$ the interfacial tension coefficient, $\kappa$ and $\mathbf{n}$ the curvature and the normal vector to the interface, respectively, and $\delta$ is a function that takes a value of one at the interface and zero elsewhere in the domain. Equation (2) can be written in dimensionless form as:

$$\rho \left[\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\left(\nabla ·\mathbf{u}\right)\right]=-\nabla p+\frac{1}{Re}\nabla ·\left\{\mu \left[\nabla \mathbf{u}+{(\nabla \mathbf{u})}^T\right]\right\}+\rho \frac{\mathbf{g}}{g}+\frac{1}{Bo}\kappa \mathbf{n}\delta ,$$

(3)

where the non-dimensional parameters $Re=\left({\rho }_1{g}^{1/2}{L}^{3/2}\right)/{\mu }_1$ and $Bo=\left({\rho }_{1}g{L}^2\right)/\sigma$ represent, respectively, the Reynolds number, which describes the ratio between inertial and viscous forces, and the Bond number (also known as the Eötvös number), which describes the ratio between gravitational force and surface tension. The parameters g, ${\rho }_1$, L and ${\mu }_1$ represent the magnitude of the gravitational field, the characteristic density, the characteristic length scale and the characteristic viscosity, respectively.

In the VOF method, assuming the flow is incompressible, the advection of the volume fraction is modeled by the following transport equation:

$$\frac{\partial f}{\partial t}+\nabla ·\left(\mathbf{u}f\right)=0.$$

(4)

The method used to model the evolution of the interface plays a fundamental role in the accuracy of the VOF method. In the HiG-Flow software, a geometric formulation was taken to solve the transport equation (Equation (4)), where the fluid-fluid interface is reconstructed at each time step using the PLIC-VOF method. In this sense, an advection algorithm with a fractional step (split) was implemented. This technique simplifies the method and prevents the transportation of the same fluid portion twice, as happens in single-step (unsplit) methods [20].

3. Numerical Procedure

3.1. HiG-Tree and HiG-Flow software

The HiG-Tree software is responsible for the data structure, domains, approximations and interpolations, solvers of linear and non-linear systems, used in numerical simulations in the HiG-Flow software.

The HiG-Tree data structure is a hierarchical Cartesian mesh representation which can be constituted of elements with varying sizes. This structure is based on the recursive spatial subdivisions (refinements) of its elements, enabling their representation in an m-tree (see Figure 1). In the HiG-Tree data structure, arbitrary mesh refinement processes can be applied at any level of the tree, generating an unstructured mesh. Classic examples of this type of mesh are the data structures Quadtree (two-dimensional) and Octree (three-dimensional). Although the illustration in Figure 1 is in two dimensions, the framework is designed to be generic and can be extended to any number of dimensions.

In the module responsible for domain and subdomain partitioning, the generated domains are composed of blocks that can be complex to simulate problems in channels with intricate geometries. Each block is discretized using the HiG-Tree structure, which enables appropriate spatial refinement for the simulated problems. These domains are partitioned using the Zoltan-Trilinos library [36,37], which ensures a good load distribution among processes during execution. This module also allows for the storage and enumeration of properties represented in the cells and facets of the mesh.

The module for approximations and interpolations handles the estimation of properties within cells and facets. It allows for the selection of the polynomial approximation order utilized by the moving least squares (MLS) method, as described by Sousa et al. [32]. In terms of meshes composed of elements with different levels of refinement, numerical methods based on approximations of finite differences in Cartesian uniform meshes become limited, since they will require spatial interpolations in unknown points of the stencil, as illustrated in Figure 2. The accuracy of this type of method heavily depends on the geometric characteristics of the mesh. The MLS method avoids this geometric dependence, using a cloud of neighboring points to the unknown to apply the interpolation.

To obtain an approximation of the second derivative of variable u in the x direction at point c, a standard second-order finite difference scheme can be used. The approximation is given by the following expression:

$$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{1}{\Delta x^2}(u_l-2u_c+u_r).$$

(5)

Note that the point $u_l$ on the stencil in Figure 2 does not coincide with any point on the mesh and must be approximated by some interpolation technique on the mesh of unknowns in the neighborhood of $u_c$. In the HiG-Tree software the following interpolation is used:

$$u_l=\sum _{k\in {\mathcal{I}}_l}{w}_k^l{u}_k,$$

(6)

where ${\mathcal{I}}_l=i_k,k=1,...,N_l$ is the set of indices for the unknowns that are in the neighborhood of $u_c$ and the weights $w_k^l$ are calculated using the MLS method. The number of neighbors $N_l$ depends on how many points are needed to maintain the order of accuracy of the global approximation.

Lastly, the module designed for solving linear and non-linear systems consists of two powerful libraries that utilize parallel computation techniques for efficient solution. These libraries are specifically tailored for handling both shared and distributed memory architectures, making use of high-performance routines. These libraries are HYPRE (High Performance Preconditioners) [38] and PETSc (Portable, Extensible Toolkit for Scientific Computation) [39]. These libraries offer significant advantages in the development of parallel code, and their performance has been extensively optimized and tested in various computational fluid dynamics codes.

The HiG-Flow software is a versatile tool for numerical simulation of fluid flows in various configurations. It enables the simulation of single-phase flows of Newtonian, generalized Newtonian, and viscoelastic fluids, as well as two-phase flows of Newtonian fluids using the volume of fluid method (VOF) for interface representation. The software is designed with a modular structure, allowing for easy implementation of new techniques and methods. It provides flexibility to users, who can choose the dimensionality of the simulation, select the desired modules (e.g., monophasic, Newtonian, generalized Newtonian, viscoelastic, biphasic), and specify numerical techniques through a data input file. This user-friendly approach facilitates customization and adaptation to specific simulation requirements.

In the HiG-Flow software, the most traditional methods for temporal advancement were implemented. The explicit methods implemented are: explicit Euler, second-order TVD (total variation diminishing) Runge-Kutta (or modified Euler method), and third-order TVD Runge-Kutta. The implicit methods implemented in the system are: implicit Euler, Crank-Nicholson, and second-order Backward Differentiation Formula (BDF). To obtain the pressure and velocity fields, the incremental (second order in time) and non-incremental (first order in time) projection methods are implemented.

3.2. Discretization procedure in HiG-Flow

To obtain the pressure fields and flow velocity through the resolution of the Navier-Stokes Equations (1) and (3), the incremental and non-incremental projection methods proposed by Chorin [29], in which the pressure and velocity fields are decoupled, are implemented in HiG-Flow. As stated by Lages [30], in terms of the error arising from the pressure-velocity decoupling, it has been observed that the non-incremental version of the method exhibits first-order accuracy in time, whereas the incremental version achieves second-order accuracy in time. In this work, the incremental projection method was used to solve the equations. Therefore, discretizing Equations (1) and (3) by the explicit Euler method we obtain:

$$\begin{array}{ccc}\hfill \nabla ·{\mathbf{u}}^{n+1}& =& 0,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \frac{{\rho }^{n+1}{\mathbf{u}}^{n+1}-{\rho }^n{\mathbf{u}}^n}{\Delta t}& =& -{\rho }^n({\mathbf{u}}^n·\nabla {\mathbf{u}}^n)-\nabla p^{n+1}+\frac{1}{Re}\nabla ·\left\{\mu \left[\nabla {\mathbf{u}}^n+{(\nabla {\mathbf{u}}^n)}^T\right]\right\}+{\mathbf{F}}^n,\hfill \end{array}$$

(8)

where

$${\mathbf{F}}^n={\rho }^n\frac{\mathbf{g}}{g}+\frac{1}{Bo}{\left(\kappa \mathbf{n}\delta \right)}^n.$$

(9)

The superscripts n and $n+1$ refer to the time level and $\Delta t$ represents the time step. The elements with superscript $n+1$ are unknowns and elements with superscript n are known.

In the incremental projection method, initially using $p^n$ instead of $p^{n+1}$, an intermediate velocity is calculated, ${\mathbf{u}}^{*}$, as follows:

$$\begin{array}{ccc}\hfill \frac{{\rho }^{*}{\mathbf{u}}^{*}-{\rho }^n{\mathbf{u}}^n}{\Delta t}& =& -{\rho }^n({\mathbf{u}}^n·\nabla {\mathbf{u}}^n)-\nabla p^n+\frac{1}{Re}\nabla ·\left\{\mu \left[\nabla {\mathbf{u}}^{*}+{(\nabla {\mathbf{u}}^{*})}^T\right]\right\}+{\mathbf{F}}^n,\hfill \end{array}$$

(10)

Subsequently, the difference between the pressure gradients $p^n$ and $p^{n+1}$ is calculated to correct the intermediate velocity and, lastly, the final velocity ${\mathbf{u}}^{n+1}$ is obtained:

$$\frac{{\rho }^{n+1}{\mathbf{u}}^{n+1}-{\rho }^n{\mathbf{u}}^{*}}{\Delta t}=-\nabla (p^{n+1}-p^n)=-\nabla \Delta p,$$

(11)

or rather,

$${\rho }^{n+1}{\mathbf{u}}^{n+1}={\rho }^n{\mathbf{u}}^{*}-\Delta t\nabla \Delta p.$$

(12)

To calculate the pressure, the divergent operator is applied on both sides of Equation (12) and the incompressibility hypothesis, Equation (7), is invoked, to give the Poisson equation for pressure:

$$\nabla ·\left(\frac{1}{{\rho }^n}\nabla \Delta p\right)=\frac{1}{\Delta t}\nabla ·{\mathbf{u}}^{*}.$$

(13)

Since the phenomena treated in this work involve two-phase flows, that is, they involve two fluids with different physical properties, it follows that $\rho$ and $\mu$ are not constant throughout the domain. Thus, the resolution of the Poisson equation for the pressure and momentum equation presents greater difficulties and requires special care. To remediate this issue, in this work, the BiCGS (BiConjugate Gradient Squared) iterative method [40] with Block Jacobi preconditioner [41] was used, considering the convergence tolerance on the residual norm of $1\times {10}^{-7}$.

3.3. Interface reconstruction

In the existing literature, various methods have been proposed in conjunction with the VOF technique to handle the geometric reconstruction of the interface. These methods aim to locally approximate the interface within each interfacial cell by utilizing the volume fraction as a fundamental piece of information. Traditional examples of methods for the reconstruction of interfaces are the SLIC (Simple Line Interface Calculation) [21], Hirt–Nichols-VOF (original VOF method) [22] and PLIC (Piecewise-Linear Interface Calculation) [23]. A comparison of these methods was performed by Rudman [24], where the PLIC method showed more consistent results. Therefore, the PLIC method was implemented in the HiG-Flow solver to represent the interface.

In the PLIC-VOF method, assuming that the volume fraction and the normal vector are known in each interfacial cell, the interface is approached by a plane defined as follows:

$$\mathbf{n}·\mathbf{x}=\alpha ,$$

(14)

where $\mathbf{n}$ is the normal vector over the interface, $\mathbf{x}$ is the position of a point on the plane and $\alpha$ is the smallest distance from the plane to the origin of the computational cell. The distance $\alpha$ must be adjusted so that the volume fraction below the reconstructed plane (say, $\overline{f}$) is equal to f, i.e.,

$$\overline{f}(\alpha ,\mathbf{n},\Delta )-f=0.$$

(15)

Given a cell with dimension $\Delta =(\Delta x,\Delta y,\Delta z)$ and the normal vector, $\mathbf{n}=(n_x,n_y,n_z)$, the fluid volume below the interface within the interfacial cell is defined as follows:

$$\begin{array}{ccc}\hfill \overline{V}& =& \frac{1}{6|n_x\left|\right|n_y\left|\right|n_z|}[{\alpha }^3\hfill \\ \hfill & & -R(\alpha -|n_x|\Delta x)·(\alpha -|n_x{|\Delta x)}^3\hfill \\ \hfill & & -R(\alpha -|n_y|\Delta y)·(\alpha -|n_y{|\Delta y)}^3\hfill \\ \hfill & & -R(\alpha -|n_z|\Delta z)·(\alpha -|n_z{|\Delta z)}^3\hfill \\ \hfill & & +R(\alpha -|n_x|\Delta x-|n_y|\Delta y)·(\alpha -|n_x|\Delta x-|n_y{|\Delta y)}^3\hfill \\ \hfill & & +R(\alpha -|n_x|\Delta x-|n_z|\Delta z)·(\alpha -|n_x|\Delta x-|n_z{|\Delta z)}^3\hfill \\ \hfill & & +R(\alpha -|n_y|\Delta y-|n_z|\Delta z)·(\alpha -|n_y|\Delta y-|n_z{|\Delta z)}^3],\hfill \end{array}$$

(16)

where

$$R\left(a\right)=\left\{\begin{array}{c}0,ifa\le 0,\hfill \\ 1,ifa>0.\hfill \end{array}\right.$$

(17)

Equation (16) can be simplified based on the intersection of the plane with the cell, allowing to estimate the distance from the plane to the origin of the cell and the volume of fluid below the plane. To calculate the normal vector, the second-order accurate method described by Mehta et al. [25] was used. In this method, the normal vector is defined as the gradient of the volume fraction (i.e., $\mathbf{n}=\nabla f/\left|\right|\nabla f\left|\right|$), where the gradient is approximated using a $3\times 3\times 3$ stencil around each interfacial cell.

3.4. Interface advection

Following the same strategy that Puckett et al. [26] adopted in the two-dimensional case, in order to maintain the conservation of f, the incompressibility condition was taken into account (i.e., $\nabla ·\mathbf{u}=0$), resulting in the following modification to the volume fraction transport equation (Equation (4)):

$$\frac{\partial f}{\partial t}+\nabla ·\left(\mathbf{u}f\right)=f\nabla ·\mathbf{u},$$

(18)

or rather,

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}+\nabla .\left(uf\right)& =& f\nabla .u,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}+\nabla .\left(vf\right)& =& f\nabla .v,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill \frac{\partial f}{\partial t}+\nabla .\left(wf\right)& =& f\nabla .w.\hfill \end{array}$$

(21)

Accordingly to Rider and Kothe [33], who also presented this strategy for the two-dimensional case, this "divergence correction" makes the local and global volume-filling constraints to be satisfied much more closely, since the discrete velocity divergences are not necessarily zero, but rather a function of the convergence tolerance used in obtaining the solutions of the linear system.

The calculation of the transport of the volume fraction is divided into three steps: first in the x direction, then in the y direction, and then finally, in the z direction. Thus, three consecutive updates are needed to carry f at each time step, i.e, $f^n\to f^{*}\to f^{**}\to f^{n+1}$. Thus, discretizing f implicitly on the right-hand side of Equation (19), implicitly on the right-hand side of Equation (20) and explicitly on the right-hand side of Equation (21), leads to

$$\begin{array}{ccc}\hfill f_{i,j,k}^{*}& =& \frac{f_{i,j,k}^n+\frac{\Delta t}{\Delta x}\left(F_{i-1/2,j,k}^n-F_{i+1/2,j,k}^n\right)}{1-\frac{\Delta t}{\Delta x}\left(u_{i+1/2,j,k}^n-u_{i-1/2,j,k}^n\right)},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill f_{i,j,k}^{**}& =& \frac{f_{i,j,k}^{*}+\frac{\Delta t}{\Delta y}\left(G_{i,j-1/2,k}^{*}-G_{i,j+1/2,k}^{*}\right)}{1-\frac{\Delta t}{\Delta y}\left(v_{i,j+1/2,k}^n-v_{i,j-1/2,k}^n\right)},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill f_{i,j,k}^{n+1}& =& f_{i,j,k}^{**}+\frac{\Delta t}{\Delta z}\left[f_{i,j,k}^{**}\left(w_{i,j,k+1/2}^n-w_{i,j,k-1/2}^n\right)-\left(H_{i,j,k+1/2}^{**}-H_{i,j,k-1/2}^{**}\right)\right],\hfill \end{array}$$

(24)

where $f_{i,j,k}^n$ is the volume fraction in the $(i,j,k)-\mathrm{th}$ cell, at time $t^n=n\Delta t$, u, v and w are the velocity components in each coordinate direction and F, G and H denote the fluxes of f in the facets of cell $(i,j,k)$. Suppose, without loss of generality, that the velocity $u_{i+1/2,j,k}$ is positive, then the flux volume, $F_{i+1/2,j,k}^n$, is given by:

$$F_{i+1/2,j,k}^n=\frac{u_{i+1/2,j,k}^n·{\overline{V}}_{i+1/2,j,k}}{u_{i+1/2,j,k}^n·\Delta t·\Delta y·\Delta z}=\frac{{\overline{V}}_{i+1/2,j,k}}{\Delta t·\Delta y·\Delta z},$$

(25)

where ${\overline{V}}_{i+1/2,j,k}$ is the volume of fluid that must be transported to the neighboring cell.

3.5. Surface tension force and curvature

In this work, the CSF (Continuum Surface Force) method presented by Brackbill et al. [27] was adopted to model the interfacial tension force described in the momentum equation (Equation (3)). In this model, the body force is defined by:

$$F_{\sigma }=\sigma \kappa \mathbf{n}\delta .$$

(26)

In the VOF method, Equation (26) is generally rewritten as:

$$F_{\sigma }=\sigma \kappa \frac{\nabla f}{\left|\right|\nabla f\left|\right|}\left|\right|\nabla f\left|\right|=\sigma \kappa \nabla f.$$

(27)

To obtain the curvature $\kappa$, the traditional height function method, described by Torrey et al. [28], was used. This method takes nine heights to form a stencil $7\times 3\times 3$, $3\times 7\times 3$ or $3\times 3\times 7$ around each interfacial cell, according to the orientation of the interface, where the construction of height functions is defined as follows:

$$h_{i,j}=\sum _{k-3}^{k+3}{f}_{i,j,k}{\Delta }_z,\mathrm{to}j=j^{\prime }-1,j^{\prime },j^{\prime }+1\mathrm{and}i=i^{\prime }-1,i^{\prime },i^{\prime }+1,$$

(28)

for the case where the largest component of the normal vector is in the z direction. Knowing the nine height functions, the interfacial curvature can be calculated as follows:

$$\kappa =sgn\left(\mathbf{n}\right)·\frac{h_{xx}+h_{yy}+h_{xx}{h}_y^2+h_{yy}{h}_x^2-2h_{xy}{h}_x{h}_y}{{(1+h_x^2+h_y^2)}^{3/2}},$$

(29)

where the derivatives $h_x$, $h_y$, $h_{xx}$, $h_{yy}$ and $h_{xy}$ are calculated using centered finite difference schemes.

3.6. Time-step constraint

During the geometric advection of the interface, it is necessary to consider certain constraints to ensure the proper functioning of the method. In HiG-Flow, in order to prevent the same fluid portion from being transported twice, due to the use of a fractional step advection algorithm, and to ensure that the advected fluid volume does not exceed the computational cell volume, the following condition has been imposed:

$$\Delta t<\frac{\Delta }{2\left|\mathbf{u}\right|}\equiv \Delta t_{Advection}.$$

(30)

The inclusion of the interfacial tension term in the Navier-Stokes equation demands one more restriction for the time step. According to Popinet [31], the explicit discretization of the interfacial tension term leads to a condition of temporal stability of the form:

$$\Delta t<\sqrt{\frac{({\rho }_1+{\rho }_2){\Delta }^3}{4\pi \sigma }}\equiv \Delta t_{\sigma },$$

(31)

where ${\rho }_1$ and ${\rho }_2$ represent the density of each fluid. According to the rationale provided by Popinet [31], the limitation on the time step size is justified by the need to adequately capture the propagation of the fastest capillary wave within the system.

Furthermore, the CFL (Courant–Friedrichs–Lewy) condition should be taken into account for the parabolic and hyperbolic terms on Equation (8),

$$\Delta t<Re{\Delta }^2\equiv \Delta t_{Parabolic},$$

(32)

and

$$\Delta t<\frac{\Delta }{\rho \left|\mathbf{u}\right|}\equiv \Delta t_{Hyperbolic}.$$

(33)

Therefore, the time step constraint is:

$$\Delta t=min(\Delta t_{Advection},\Delta t_{\sigma },\Delta t_{Parabolic},\Delta t_{Hyperbolic}).$$

(34)

4. Validation

In this section, we present the validation of the numerical results obtained using the HiG-Flow software by comparing them with the results reported in the scientific literature. The purpose is to assess the accuracy and reliability of the HiG-Flow software by evaluating its agreement with established findings. To assess the performance of the implemented interface reconstruction and advection methods, benchmark problems were simulated. These benchmark problems, namely the 3D shearing flow test and the 3D deforming field test, were chosen based on their relevance and availability in the literature, particularly in the work by Liovic et al. [15]. In these tests, since the velocity field is given, the sphere is not affected by the numerical dissipation of the methods used to obtain the flow velocity and pressure. In other words, the Navier-Stokes equations do not need to be solved, just the volume fraction transport equation. Therefore, it is expected that the sphere maintains its original shape at the final time of the simulation. By comparing the results obtained from the HiG-Flow simulations with the established solutions, the accuracy and fidelity of the implemented methods can be evaluated. Furthermore, the solver for the Navier-Stokes equations developed in HiG-Flow was validated using the 3D rising drop test proposed by Adelsberger et al. [16]. This test serves as a benchmark to assess the accuracy and performance of the solver in simulating the rising behavior of a three-dimensional droplet. Both quantitative and qualitative analyses of the results obtained from the simulations of these benchmark problems are presented in the following sections.

To quantitatively evaluate the performance of the advection algorithm implemented in HiG-Flow, the volume of fluid at the initial time $\left(V_I\right)$ and the volume of fluid at the final time $\left(V_F\right)$ of the simulation were compared for both the 3D shearing flow test and the 3D deforming field test. The absolute error, ${ϵ}_{L_1}$, between the volumes is defined as follows:

$${ϵ}_{L_1}=|V_I-V_F|.$$

(35)

To quantitatively evaluate the 3D rising drop test, the evolution of the bubble’s center of mass and the bubble’s rise velocity were analyzed. The evolution of the bubble’s center of mass is given by:

$$z_c={\int }_{{\Omega }_1}zdz,$$

(36)

where ${\Omega }_1$ is the region of the domain filled by the bubble. By utilizing the volume fraction f, Equation (36) can be expressed in a different form, as follows:

$$z_c=\sum _{i=1}^N{z}_i·f·\Delta x\Delta y\Delta z,$$

(37)

where N is the number of cells that contain a non-zero volume fraction.

Lastly, the bubble’s rise velocity is given by:

$$v_c=\frac{d}{dt}{z}_c.$$

(38)

4.1. 3D shearing flow test

The 3D shearing flow test is widely used in the literature to verify the accuracy of interface advection algorithms in two-phase flows [15,34,35]. In this test, a sphere is subjected to a prescribed velocity field defined by the combination of a single vortex in the $x-y$ plane and laminar pipe flow in the z direction.

To simulate this test, a sphere with a radius of 0.15 was considered, centered at the point $(0.5,0.75,0.25)$, in a domain with dimensions of $[0,1]\times [0,1]\times [0,2]$. The velocity field used to move the sphere within the domain is defined by:

$$\mathbf{u}=\left\{\begin{array}{c}u=-sin{\left(\pi x\right)}^{2}sin\left(2\pi y\right)cos\left(\frac{\pi t}{T}\right),\hfill \\ v=sin{\left(\pi y\right)}^{2}sin\left(2\pi x\right)cos\left(\frac{\pi t}{T}\right),\hfill \\ w=u_{max}{\left(1-\frac{r}{R}\right)}^{2}cos\left(\frac{\pi t}{T}\right),\hfill \end{array}\right.$$

(39)

where $T=6.0$, $r=\sqrt{{(x-0.5)}^2+{(y-0.5)}^2}$, $R=0.5$ and $u_{max}=1.0$. This velocity field generates a vortical flow, causing the initially spherical interface to take on a spiral shape. After reaching half of the total simulation time (i.e., t = 3.0), the interface reverses its trajectory and returns to its initial position at time t = 6.0, and should also reverting to its original shape.

Figure 3a shows the point of maximum interface deformation in the 3D shear flow test simulated in HiG-Flow. A mesh with a resolution of $64\times 64\times 128$ and time-step of $\Delta t=1\times {10}^{-3}$ were used in the numerical simulations. The result obtained by HiG-Flow is compared with the findings presented by Liovic et al. [15], as shown in Figure 3b, where a mesh with a resolution of $64\times 64\times 128$ was also used. The interface shape calculated by HiG-Flow demonstrates qualitative agreement with the results obtained by Liovic et al. [15].

To perform a quantitative analysis of the 3D shearing flow test, the volume fraction absolute error was calculated using Equation (35). The obtained result was then compared with the values reported in the scientific literature. This comparison is shown in Table 1, where the HiG-Flow software showed superiority, in terms of volume fraction conservation, when compared to the results presented by Liovic et al. [15], Duz et al. [34] and Lopez et al. [35].

4.2. 3D deforming field test

The 3D deforming field test is also widely used in the literature to verify interface advection methods [15,35], as the prescribed velocity field causes the interface to assume a complex geometry during its motion. To perform this test, a sphere with a radius of 0.15 was considered, centered at the point (0.35, 0.35, 0.35), within a cubic cavity with dimensions $[0,1]\times [0,1]\times [0,1]$. The velocity field used to move the sphere within the domain is defined by:

$$\mathbf{u}=\left\{\begin{array}{c}u=2sin{\left(\pi x\right)}^{2}sin\left(2\pi y\right)sin\left(2\pi z\right)cos(\pi t/T),\hfill \\ v=-sin{\left(\pi y\right)}^{2}sin\left(2\pi x\right)sin\left(2\pi z\right)cos(\pi t/T),\hfill \\ w=-sin{\left(\pi z\right)}^{2}sin\left(2\pi x\right)sin\left(2\pi y\right)cos(\pi t/T),\hfill \end{array}\right.$$

(40)

where $T=3.0$.

Similarly to the 3D shearing flow test studied in Section 4.1, the prescribed solenoidal velocity field causes the sphere to deform until reaching half of the total simulation time, and then follows a reverse trajectory, tending back to its original shape and position.

Figure 4a shows the point of maximum interface deformation in the 3D deforming field test simulated in HiG-Flow. A mesh with a resolution of ${128}^3$ and time-step of $\Delta t=1\times {10}^{-3}$ were used in the numerical simulations. The result is compared with the findings presented by Liovic et al. [15], as shown in Figure 4b, where a mesh with a resolution of ${128}^3$ was also used. The interface shape obtained by HiG-Flow exhibits qualitative agreement with the results obtained by Liovic et al. [15], with fewer regions of void volume fraction compared to their results, despite using the same level of mesh refinement.

Analogously to the previous test shown in Section 4.1, we analyze quantitatively the results obtained here for the 3D deforming field test. For that purpose, the volume fraction absolute error was calculated using Equation (35) in order to verify the ability of the volume tracking scheme to conserve volume. The results obtained by our newly-developed multiphase algorithm were compared with the results presented by Liovic et al. [15] and by Lopez et al. [35], as shown in Table 2. Across all levels of mesh refinement utilized, our numerical formulation consistently yields lower errors compared to those reported in references [15] and [35]. These results support the methodology of developing an interface advection scheme with a divergence correction approach inside the volume fraction transport equation (see right-hand side of Equation(18)), which leads to a closer satisfaction of both local and global volume-filling constraints.

4.3. 3D buoyancy-driven rising drop

To conclude the validation of the newly implemented numerical algorithm for simulating 3D multiphase flows, we now turn our attention to the problem of the buoyancy-driven rise of a bubble immersed in a Newtonian fluid. In this test case, we consider a bubble with a density lower than the fluid filling the column. As a result, the pressure gradient induced by gravity causes the bubble to move in the direction opposite to the gravitational field.

A bubble with a radius of $R=0.25$ was positioned at the point $(0.5,0.5,0.5)$ inside a tank with dimensions of $[0,1]\times [0,1]\times [0,2]$. The no-slip boundary condition was imposed on the domain walls. The simulations were performed until time reached $t=3.0$. A mesh size of $50\times 50\times 100$ elements and a time-step of $\Delta t=2.0\times {10}^{-3}$ was employed in the simulations. The two test cases proposed by Adelsberger et al. [16] were simulated. In the first case, the density and viscosity ratios between the two fluids in the tank are both equal to 10 (${\rho }_1/{\rho }_2={\mu }_1/{\mu }_2=10$). In the second case, the density ratio between the fluids is 1000 (${\rho }_1/{\rho }_2=1000$) and the viscosity ratio is 100 (${\mu }_1/{\mu }_2=100$). In test case 1, the droplet undergoes a transformation into an ellipsoidal shape. On the other hand, in test case 2, the droplet has the potential to break up due to the lower surface tension. Test case 1 can be categorized as an ellipsoidal regime, while test case 2 falls somewhere between the skirted and dimpled ellipsoidal-cap regimes. The physical and dimensionless parameters considered for simulating these two cases are presented in Table 3. These parameters were also used by Metivet et. al [1] with an high-order finite-element software, so-called FEEL.

Figure 5 illustrates a comparison of the bubble shapes at time $t=3$ simulated using HiG-Flow and OpenFOAM, for the test case 1 (ellipsoidal regime). The interface shape obtained by HiG-Flow exhibits qualitative agreement with the result obtained by Adelsberger et al. [16] using the OpenFOAM software, where the droplet reaches a stable ellipsoidal shape at steady-state regime.

For a better comparison of the actual droplet shape, Figure 6 illustrates the two-dimensional contour through the droplet’s center generated by the HiG-Flow (this work) and FEEL [1] software at time $t=3$. The contour lines of HiG-Flow (blue) and FEEL (green) exhibit a close resemblance to each other, with a small difference at the bottom of the drop.

Figure 7 shows the evolution of the bubble’s center of mass and bubble’s rise velocity for test case 1. The results obtained using the newly-developed algorithm in HiG-Flow are compared to the results computed by FEEL [1], as well as OpenFOAM, DROPS, and NaSt3D [16]. In Figure 7(a), at the initial stage of the simulation, up until approximately $t\approx 1$, the drop position is nearly identical in all simulations. Afterwards, the vertical position of the droplet in OpenFOAM and FEEL software is lower than the one obtained in HiG-Flow, DROPS and NaSt3D. Our results are in excellent agreement with the findings obtained by DROPS and NaSt3D. Figure 7(b) illustrates the rise velocity of the droplet. In the initial stages of the simulation, the vertical velocity component $\left(w\right)$ of the fluid velocity gradually increases from its initial value of zero, reaching a maximum magnitude of approximately $0.35-0.36$ at around $t\approx 0.9$. Then, the steady-state velocity of the droplet decreases, being again in agreement with the findings obtained by DROPS and NaSt3D.

Figure 8 illustrates a comparison of the bubble shapes at time $t=3$ simulated using HiG-Flow and OpenFOAM, for the test case 2 (skirted and dimpled ellipsoidal-cap regimes). Again, the interface shape obtained by HiG-Flow exhibits qualitative agreement with the result obtained by Adelsberger et al. [16] using the OpenFOAM software, where the droplet reaches a stable cap shape at steady-state regime. In the simulated time, our results do not indicate the presence of satellite droplets, contrary to the OpenFOAM results which predicts 8 satellite droplets as stated in Adelsberger et al. [16].

Figure 9 illustrates the comparison between two-dimensional contours through the droplet’s center generated by the HiG-Flow (this work) and FEEL [1] software at time $t=3$. The shape of the contour lines of the bubbles obtained by both software are similar and reach a skirted format, since the surface tension is lower and the ratio between the densities and viscosities between the fluids present in the flow is higher. The differences between the two simulations are particularly noticeable along the bottom edge of the droplet, where the curvature is most pronounced. As depicted in the contour plot, the contour line in FEEL (green) appears more extended than in HiG-Flow. Furthermore, FEEL exhibits a higher contour line at the top edge of the droplet and a smaller contour line at the central part of the bottom side, when compared to that one obtained with HiG-Flow software.

Lastly, Figure 10 shows the evolution of the bubble’s center of mass and bubble’s rise velocity for test case 2. In Figure 10(a), we show that the droplet’s position over time is similar for all the simulation time in all the flow solvers, except for the FEEL software, which can explain the differences shown in the contour plot in Figure 9. Figure 10(b) illustrates the rise velocity of the droplet. In the initial stages of the simulation, the vertical velocity component $\left(w\right)$ of the fluid velocity gradually increases from its initial value of zero, reaching a maximum magnitude of approximately $0.37$ at around $t\approx 0.54$. Then, the steady-state velocity of the droplet decreases so that a final velocity of 0.32 at $t=3$ is obtained using the HiG-Flow software. Our result is similar to the OpenFOAM, DROPS and NaSt3D software, which predict a final velocity of 0.31. The final velocity predicted by FEEL is approximately 0.28, showing the most significant difference compared to the velocities obtained from other software.

5. Conclusions

In conclusion, this study presents a newly developed numerical method for simulating two-phase flows of incompressible, immiscible, and isothermal Newtonian fluids. The geometric Volume of Fluid (VOF) method is employed to locate and reconstruct the interface during the dynamics of the two-phase flow. Specifically, the Piecewise-Linear Interface Calculation (PLIC) technique is implemented for three-dimensional (3D) interface transport. This technique approximates the interface in each interfacial cell using line segments. The methodology of developing an interface advection scheme with a divergence correction approach inside the volume fraction transport equation, leads to a closer satisfaction of both local and global volume-filling constraints. In addition, the continuum interfacial force (CSF) model is utilized to account for the balance of forces at the interface. Laslty, the Navier-Stokes equations are solved using meshless (possibly high-order) finite difference schemes from the HiG-Flow computational fluid dynamics software.

To validate the numerical implementation, several benchmark two-phase flows are simulated with the 3D PLIC-VOF HiG-Flow algorithm. Two standard advection numerical tests, namely the 3D shearing flow test and the 3D deforming field test, are conducted, and good agreement is observed between the 3D interface shapes obtained by the 3D PLIC-VOF HiG-Flow algorithm and those found in the scientific literature. Furthermore, the absolute error of the volume tracking advection calculation is found to be smaller than the values reported in the scientific literature for both tests, even if using coarser meshes in our calculations. Additionally, the 3D bubble rising problem is simulated for different fluid density and viscosity ratios, mimicking two different regimes, the ellipsoidal regime and the skirted and dimpled ellipsoidal-cap regime. The results obtained with the newly implemented algorithm are in good agreement with those obtained using the OpenFOAM, DROPS and NaSt3D software. The bubble’s rise velocity and center of mass evolution are computed with the 3D PLIC-VOF HiG-Flow algorithm and found to be consistent with these software.

Overall, the developed method demonstrates promising accuracy and reliability in simulating complex two-phase flows, making it a valuable tool for future fluid dynamics research. Particularly, its applicability extends to the study of complex fluid rheology, including viscoelastic fluids.