Comparison of Surface Tension Models for the Volume of Fluid Method

Abstract

With the increasing use of Computational Fluid Dynamics to investigate multiphase flow scenarios, modelling surface tension effects has been a topic of active research. A well known associated problem is the generation of spurious velocities (or currents), arising due to inaccuracies in calculations of the surface tension force. These spurious currents cause nonphysical flows which can adversely affect the predictive capability of these simulations. In this paper, we implement the Continuum Surface Force (CSF), Smoothed CSF and Sharp Surface Force (SSF) models in OpenFOAM. The models were validated for various multiphase flow scenarios for Capillary numbers of 10${}^{-3}$–10. All the surface tension models provide reasonable agreement with benchmarking data for rising bubble simulations. Both CSF and SSF models successfully predicted the capillary rise between two parallel plates, but Smoothed CSF could not provide reliable results. The evolution of spurious current were studied for millimetre-sized stationary bubbles. The results shows that SSF and CSF models generate the least and most spurious currents, respectively. We also show that maximum time step, mesh resolution and the under-relaxation factor used in the simulations affect the magnitude of spurious currents.

1. Introduction

For a comprehensive understanding of flow physics in multiphase systems, which is ubiquitous in both nature and technological processes, consideration of surface tension is important. For instance, the break down of a fluid jet into droplets is used to form droplets in inkjets [1] and lab-on-chip devices [2] while the thinning and breakdown dynamics of non-Newtonian fluid filaments is critical in its application in jetting [3,4]. Flow scenarios such as underground water flows [5], oil recovery [6] and paper-based microfluidics [7] are examples of flow through porous media where dominance of surface tension may produce a capillary rise. The detachment diameter of the bubble [8,9] and shape of rising bubble [10] during bubble evolution in champagne, boiling and electrochemical gas evolution is dependent on surface tension, as is the droplet size produced during atomisation of fuels [11], spraying [12,13] and growth of a bubble in confined geometries [14]. The effect of surface tension is also important in events such as nucleation of bubbles [15,16] and droplets [17].

Due to the importance and complex nature of multiphase flows, numerical simulations, especially computational fluid dynamics (CFD), are commonly used to study and understand these processes. The CFD strategies used to model multiphase flows can broadly be divided into four categories: Euler–Euler (EE), Euler–Lagrange (EL), interface tracking and capturing methods. The EE approach assumes that phases are interpenetrating, which is efficient when modelling large-scale industrial processes [18,19], while EL tracks the dispersed phases individually, which can be computationally expensive [20,21]. As both EE and EL approaches do not resolve the complete interactions between the phases, they require so called “closure laws” (see [18,19,20,21]). Interface tracking methods, such as the moving mesh method, use a separate boundary-fitted moving mesh for each phase [22]. Although interface tracking methods are quite accurate, they are typically used to model bubble or droplets with mild-moderate deformations [22,23] but to handle complex interface deformations these methods require a global or local re-meshing [24]. Interface capturing methods use a fixed grid with functions to capture the interface such as the Volume Of Fluid method (VOF) [25], level-set [26] and diffuse interface methods [27]. Other methods available in the literature employ a hybrid interface tracking-capturing approach, such as the immersed boundary [28] and front tracking method [29]. Due to its ability to conserve mass (both level-set [30] and phase-field [31] models have difficulties in conserving mass), robustness and ability to produce reasonably sharp interfaces VOF is very popular in multiphase simulations [32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57] and implemented in both commercial (ANSYS Fluent^® and Flow-3D^®) and open source (OpenFOAM^®) CFD packages.

Due to the popularity of open source CFD packages, this paper predominantly delves into the VOF formulation and reported development in interFoam, which is the VOF-based solver available in OpenFOAM^®. In the VOF method, a scalar function representing the volume fraction of phases in the computational cells is advected. The advection of the volume fraction equation is done using specific discretisation schemes, such as the interface compression method [58], to prevent the excessive smearing of the interface thickness. Apart from interface compression method, recent work has explored reconstruction of interface using techniques such as the isoAdvector method [59,60] and piece wise linear interface calculation (PLIC) algorithm [61]. Although the VOF approach in theory produces a sharp interface, the “real” VOF, which is implemented in solvers such as interFoam, produces a non-sharp interface, which stretches over a few computational cells. This non-sharp nature of the volume fraction leads to errors in the calculated curvature which generates spurious currents that is quantified in the work by Harvie et al. [62], appearing as vortices around the interface (see [63,64]). The presence of these spurious currents introduces non-physical velocities near the interface, which can increase the interfacial mass transfer while modelling condensation [32] and evaporation [57] scenarios and adversely effects the accuracy of simulations. In the literature, spurious currents in VOF methods can be reduced by the following approaches:

• force balance, which is achieved by discretising the surface tension and pressure forces at the same location [65];
• accurate calculation of the curvature (see

  Table 1. An overview of improved curvature calculations and surface tension models developed for VOF method.

  Publication	Remarks
  Brackbill et al. [66]	Introduced the Continuum Surface Force (CSF) and density scaled CSF models. These methods are very common due to its relatively straightforward implementation in a VOF framework.
  Ubbink [67]	Proposed using a smoothed volume fraction to calculate curvature (referred to as “Smoothed CSF” in this paper). Using a smoothed volume fraction to compute the curvature instead of non-smoothed volume fraction in CSF model reduced spurious current up to one order of magnitude [56]. This method has been used in modelling condensing bubbles [32] and droplets in microfluidic devices [56]. A similar smoothening of ${\alpha }_1$ was proposed by Heyns and Oxtoby [68].
  Sussman and Puckett [69]	Developed a fully coupled level-set VOF (CLSVOF) method which combines the mass conservativeness of VOF method with smoothness of the level-set function to reduce spurious currents. The CLSVOF method has been used to applications such as splashing droplets [45], flow through microfluidic devices [46], wave breaking [47] and droplet evaporation [43]. Another variant of coupled level set approach is the simple coupled level-set VOF (S-CLSVOF) proposed by Albadawi et al. [70].
  Raessi et al. [71]	Proposed a method to calculate $\kappa$ based on advected normals. The spurious currents were lower than CSF (within the same order) while modelling cases such as stationary bubble, rising bubble and Rayleigh–Taylor instability [71].
  Renardy and Renardy [72]	Introduced parabolic reconstruction of surface tension (PROST) algorithm which uses a least-squares fit of the interface to a quadratic surface. The spurious current produced by the algorithm is lower by two orders of magnitude compared to CSF [72]. The model was used to simulate droplet deformation including breakup [48,72].
  Cifani et al. [61]	Implemented piecewise linear interface calculation (PLIC) algorithm (proposed by Youngs [73]) to reconstruct the interface in interFoam and managed to reduce spurious currents.
  Pilliod and Puckett [74]	Developed an efficient least squares volume-of-fluid interface reconstruction algorithm (ELIVRA) which reconstructs the interface using a least square method to fit the interface to a linear surface.
  Popinet [75]	Proposed calculating curvature using height functions. Use of height functions have reduced spurious current (slightly outperformed PROST algorithm [75]) and has been shown to model flow in microchannels [49], rising bubble [34,44] and other multiphase flows [50].
  Raeini et al. [76]	Introduced a sharp surface force formulation to calculate the capillary force, which is then filtered to reduce the spurious currents (known as FSF model). Neglecting the filtering terms in the FSF model provides a sharp surface formulation of surface tension known as SSF, which is described in [42]. The SSF has been reported to be reduce the spurious currents by two to three orders in comparison to CSF [42]. The FSF model has been reported to provide periodic bursts in the velocity fields but lower spurious current than SSF [42]. The approach has been used to model bubbles in microfluidic devices [51] and flow through porous media [52].
  Denner et al. [77]	Proposed the use of artificial viscosity model, which applies artificial shear stress in the tangential direction to interface, to dampen the spurious currents. The model has been used to model rising bubble and capillary instability of a water jet [77].
  Lafaurie et al. [78]	Proposed an alternative to the CSF model, known as the Continuum Surface Stress (CSS) model, determines surface tension as divergence of stress tensor without relying on complex curvature calculations. Due to imbalance in the surface tension and pressure, CSS model can also produce spurious currents [35] which has reported to be in the same order as CSF [72]. CSS model has been used to model static droplets and rising bubbles [35], but it does not provide reliable results for falling films [41].

  ); and
• choosing the appropriate time step for the solver (see [63]).

To analyse the force balance (described in [65]), Deshpande et al. [63] evaluated interFoam and showed that there is no imbalance in the surface tension and pressure forces due to inconsistent discretisation. However, the iterative process, which is used to solve pressure equation, introduces an imbalance which is related to the user defined tolerance level of the solution [63]. An overview of literature that provides an improved estimate of the interface curvature and surface tension modelling approaches is provided in Table 1. The improved representation of the interface (which aids in accurate calculation of the interface curvature) is provided bycoupled level-set VOF (CLSVOF) method, height functions and interface reconstruction algorithms (like piecewise linear interface calculation (PLIC), parabolic reconstruction of surface tension (PROST) and efficient least squares volume-of-fluid interface reconstruction (ELIVRA) algorithms), whereas the other methods discussed in Table 1 provide alternative approaches to model surface tension. To ensure that spurious currents do not grow over time, a stability condition, proposed by Brackbill et al. [66], for explicit treatment of surface tension is

$$\Delta t<\sqrt{\frac{{\rho }_{avg}{(\Delta x)}^3}{2\pi \sigma }},$$

(1)

where $\Delta x$, $\sigma$ and ${\rho }_{avg}$ are grid spacing, surface tension and average of density of both phases, respectively. As proposed by Galusinski and Vigneaux [79], a comprehensive constraint on the time step must consider the effect of both density and viscosity which can be expressed as

$$\Delta t\le \frac{1}{2}\left(C_2{\tau }_{\mu }+\sqrt{{\left(C_2{\tau }_{\mu }\right)}^2+4C_1{\tau }_{\rho }^2}\right),$$

(2)

where ${\tau }_{\mu }$ and ${\tau }_{\rho }$ are time scales which are defined as ${\mu }_{avg}\Delta x/\sigma$ and $\sqrt{{\rho }_{avg}{(\Delta x)}^3/\sigma }$, where ${\mu }_{avg}$ is the average dynamic viscosity between the phases, respectively. An evaluation of interFoam, by Deshpande et al. [63], proposed that time step should satisfy

$$\Delta t\le \max \left(C_2{\tau }_{\mu },10C_1{\tau }_{\rho }\right),$$

(3)

along with the time step constraint discussed in Equation (2). Deshpande et al. [63] also calculated the values of $C_1$ and $C_2$ for interFoam to be equal to 0.01 and 10, respectively. Further details of the numerical methods used to model surface tension is available in the recent review work by Popinet [80].

In the literature, comparison between surface tension models is often done for a specific of flow phenomenon and at times a static scenario is used to quantify the spurious currents. Some examples of flow phenomena used to compare surface tension models are rising bubbles whose diameters are in the order of few millimetres [33,34,35], translating and rotating bubbles [64], oscillating droplets or bubbles [34], stagnant bubbles or droplets [34,35,39,64], Rayleigh–Taylor instability [37,38], Taylor bubbles [64], falling films [41], droplet splashing [38,39], capillary rise [42] and bubble evolution [37,40]. These typically compare the CSF model with height functions [33,34,64], PROST [37], PLIC [42], CLSVOF and its variants [37,38,39,40,64], FSF and SSF [42], and CSS [35,41] models. Although the flow scenarios that are used to compare surface tension models are diverse, they can be broadly categorised based on the dominance of surface tension in the flow using the Capillary number ($Ca$), which is defined as the ratio of viscous to surface tension forces in the system. Flow phenomena such as capillary rise and stationary bubbles are examples of low values of $Ca$ whereas flows with larger values of $Ca$ include rising bubbles and falling films.

During processes such as gas evolution during electrochemical reactions and boiling, nucleated bubbles grow by mass transfer across the interface [15,16] or coalescence [8], but once the bubble detaches it may deform as it rises up and/or interacts with other bubbles [53]. Other complex processes, such as splashing, involve droplet spreading on a surface which is accompanied by formation of secondary smaller droplets at the rim [81]. To reliably model these processes, surface tension models must be able to accurately handle flow scenarios with both small and large capillary numbers. In literature, the work by Hoang et al. [56] implemented the Smoothed CSF approach to model the steady motion of bubbles in a straight two-dimensional channel, the formation of bubbles in two- and three-dimensional T-junctions, and the breakup of droplets in three-dimensional T-junctions. A study by Heyns and Oxtoby [68] implemented a selection of surface tension modelling approaches (e.g., the CSF, a variant of Smoothed CSF and a force-balanced higher-resolution artificial compressive formulation) to model a stationary bubble. To the best of the authors’ knowledge, a recent study by Yamamoto et al. [36] is the only one of its kind where different surface tension models (i.e., S-CLSVOF, density scaled S-CLSVOF and CSF) are compared based on a variety of processes with various capillary numbers (e.g., rising bubbles, capillary rise, capillary wave and thermocapillary flows).

In this study, we implemented three different surface tension models, namely CSF [66], Smoothed CSF [67] and SSF [76], in interFoam available in OpenFOAM 6. To investigate the capability of the surface tension models to handle various flow scenarios, we used two benchmark cases:

• two-dimensional rising bubbles (proposed by Hysing et al. [54], Klostermann et al. [55]); and
• two-dimensional capillary rise.

These two benchmark cases were selected due to their relevance in a variety of processes. To compare the spurious currents generated by the surface tension models, a stationary bubble was simulated. For practical applications, the time step constraint can substantially increase the computational time, thus the temporal development of the spurious currents with the surface tension models were also examined. Furthermore, the evolution of spurious currents with mesh resolution and under-relaxation factor used for the simulations was also investigated. In the interest of knowledge dissemination, the solvers and the test cases (implemented in OpenFOAM 6) discussed in the paper are available in the Supplementary Materials.

2. Numerical Model

2.1. Governing Equations

The VOF approach (developed by Hirt and Nichols [25]) denotes the individual phases using a scalar function called volume fraction, represented as

$${\alpha }_1(\overrightarrow{x},t)=\left\{\begin{array}{c}0(\mathrm{within}\mathrm{Phase}2\mathrm{or}\mathrm{gas})\hfill \\ 0<{\alpha }_1<1(\mathrm{at}\mathrm{the}\mathrm{interface})\hfill \\ 1(\mathrm{within}\mathrm{Phase}1\mathrm{or}\mathrm{liquid}),\hfill \end{array}\right.$$

(4)

where ${\alpha }_1$ is the volume fraction of liquid. The fluid properties such as density ($\rho$) and viscosity ($\mathsf{\mu }$) in a control volume are calculated as

$$\chi ={\chi }_1{\alpha }_1+{\chi }_2{\alpha }_2\mathrm{where}\chi \in [\rho ,\mu ],$$

(5)

where ${\chi }_1$ and ${\chi }_2$ represent the fluid property of liquid and gas phase, respectively.

Considering the fluids to be incompressible, isothermal and immiscible, the VOF approach solves a single set of continuity and momentum equation for the whole domain. The continuity equation is written as

$$\nabla ·\overrightarrow{U}=0,$$

(6)

where $\overrightarrow{U}$ is the fluid velocity. The momentum equation is

$$\frac{\partial \rho \overrightarrow{U}}{\partial t}+\nabla ·\left(\rho \overrightarrow{U}\overrightarrow{U}\right)=-\nabla p+\rho \overrightarrow{g}+\nabla ·\mu \left(\nabla \overrightarrow{U}+\nabla {\overrightarrow{U}}^T\right)+{\overrightarrow{F}}_{st},$$

(7)

where last term represents the surface tension forces, the second last term is the viscous term, $\overrightarrow{g}$ is the acceleration due to gravity and p is the pressure. Advection of the volume fraction of liquid (${\alpha }_1$) is implemented in interFoam as

$$\frac{\partial {\alpha }_1}{\partial t}+\nabla ·\left({\alpha }_1\overrightarrow{U}\right)+\nabla ·\left({\alpha }_1(1-{\alpha }_1)\overrightarrow{U_r}\right)=0,$$

(8)

where the third term is an artificial compression term used to sharpen the interface [58,61]. The artificial compression term uses a relative velocity ($\overrightarrow{U_r}$) defined as

$$\overrightarrow{U_r}=C_{\alpha }\left|\frac{\varphi }{|S_f|}\right|\overrightarrow{n_f},$$

(9)

where $\varphi$, $S_f$, $C_{\alpha }$ and $\overrightarrow{n_f}$ are the velocity flux, face surface area, an adjustable compression factor and unit normal vector to the interface, respectively. In the literature, $C_{\alpha }$ can be set between 0 and 4, where $C_{\alpha }$ equal to zero and one correspond to the case of no and moderate compression, respectively [56]. The increase in the value of $C_{\alpha }$ sharpens the interface but increases the magnitude of spurious currents (see [51,56]). To model practical flow scenarios using interFoam, the value of $C_{\alpha }$ is generally set to unity [32,63].

2.2. Surface Tension Models

This section introduces the three surface tension models: CSF, Smoothed CSF and SSF.

2.2.1. The Continuum Surface Force (CSF) Model

Proposed by Brackbill et al. [66], the CSF model provides a volumetric representation of surface tension, represented as

$$\overrightarrow{F_{st}}=\sigma \kappa \nabla {\alpha }_1,$$

(10)

where $\sigma$ is the surface tension and $\kappa$ is the curvature, defined in Equation (13). The unit normal vector at the interface is calculated as

$$\widehat{n}=\frac{\nabla {\alpha }_1}{|\nabla {\alpha }_1|+\delta },$$

(11)

where $\delta$ is a small non-zero term to ensure that the denominator does not become zero. $\delta$ is calculated as ${10}^{-8}/{\left(\frac{{\sum }_N{V}_i}{N}\right)}^{1/3}$, where N is the number of computational cells and ${\sum }_N{V}_i$ provides the sum of the volumes of individual cells (represented by i). Once $\widehat{n}$ is calculated, it is corrected to account for wall adhesion through

$$\widehat{n}={\widehat{n}}_w\cos \theta +{\widehat{t}}_w\cos \theta$$

(12)

where $\theta$ is the contact angle of the gas–liquid interface at the walls (measured in the liquid phase), and ${\widehat{n}}_w$ and ${\widehat{t}}_w$ are unit vectors that are normal and tangential to the wall, respectively [82]. The curvature of the interface is then calculated as

$$\kappa =-\nabla ·\widehat{n}.$$

(13)

2.2.2. The Smoothed CSF Model

The Smoothed CSF model (by Ubbink [67]) proposed modifying CSF by modifying the calculation of curvature of interface by using a smoothed volume fraction of liquid (${\alpha }_1$).

The smoothed volume fraction field is calculated using a smoother proposed by Lafaurie et al. [78], which has been implemented in the literature [32,56] and is represented as

$$\widetilde{{\alpha }_1}=\frac{{\sum }_{f=1}^N<{\alpha }_1{>}_{c\stackrel{}{\to }f}{S}_f}{{\sum }_{f=1}^N{S}_f},$$

(14)

where the indices c and f are the cell and face centre indices, respectively. $<{\alpha }_1{>}_{c\stackrel{}{\to }f}$ represents the interpolation of ${\alpha }_1$ from cell to face centre. The smoothening of volume fraction, done using Equation (14), is applied twice to obtain a smooth volume fraction field, which is used in Equation (15). Implementation of Equation (14) in interFoam is done using the subroutine developed in the work by [56]. Based on the smoothed volume fraction field, the unit normal to the interface is calculated as

$$\tilde{\widehat{n}}=\frac{\nabla \widetilde{{\alpha }_1}}{|\nabla \widetilde{{\alpha }_1}|+\delta },$$

(15)

which is then corrected for wall adhesion (based on Equation (12)). The curvature of the interface is then calculated as

$$\tilde{\kappa }=-\nabla ·\tilde{\widehat{n}}.$$

(16)

The surface tension can be represented using the modified curvature ($\tilde{\kappa }$ in Equation (16)), which can be represented as

$$\overrightarrow{F_{st}}=\sigma \tilde{\kappa }\nabla {\alpha }_1.$$

(17)

2.2.3. The Sharp Surface Force (SSF) Model

In the SSF model, proposed by Raeini et al. [76], smoothened and sharpened volume fraction fields are used to calculate curvature and gradient of of volume fraction.

The smoothened volume fraction (${\alpha }_s$) is calculated based on interpolating the cell-centred values of ${\alpha }_1$ to the cell faces using a three consecutive smoothening steps described using Equations (18a)–(18c)

$${\alpha }_{s1}=\mathcal{C}{〈<{\alpha }_1{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}+\left(1-\mathcal{C}\right){\alpha }_1,$$

(18a)

$${\alpha }_{s2}=\mathcal{C}{〈<{\alpha }_{s1}{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}+\left(1-\mathcal{C}\right){\alpha }_{s1},$$

(18b)

$${\alpha }_s=\mathcal{C}{〈<{\alpha }_{s2}{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}+\left(1-\mathcal{C}\right){\alpha }_{s2},$$

(18c)

where $\mathcal{C}$ is set equal to 0.5. The unit normal to the interface is then calculated as

$${\widehat{n}}_s=\frac{\nabla {\alpha }_s}{|\nabla {\alpha }_s|+\delta },$$

(19)

which is then corrected for wall adhesion (based on Equation (12)). The curvature (${\kappa }_s$) is calculated using Equation (19) as

$${\kappa }_s=-\nabla ·{\widehat{n}}_s.$$

(20)

The interface curvature is smoothed by using a three step procedure, which can be broadly summarised into Equations (21a), (21c), and (21d). The first step involves smoothening the curvature calculated in Equation (20) as

$${\kappa }_{f1}=\left(2\sqrt{{\alpha }_c(1-{\alpha }_c)}\right){\kappa }_s+\left(1-2\sqrt{{\alpha }_c(1-{\alpha }_c)}\right){\kappa }_s^{\ast }$$

(21a)

where ${\alpha }_c$ is defined as min(1,max(${\alpha }_1$,0)) and

$${\kappa }_s^{\ast }=\frac{{〈<w{\kappa }_s{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}}{{〈<w{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}},w=\sqrt{{\alpha }_c(1-{\alpha }_c)+{10}^{-3}}.$$

(21b)

The second step further smoothens the curvature (calculated in Equation (21a)) as

$${\kappa }_{f2}=\left(2\sqrt{{\alpha }_c(1-{\alpha }_c)}\right){\kappa }_s+\left(1-2\sqrt{{\alpha }_c(1-{\alpha }_c)}\right){\kappa }_{s2}^{\ast },\mathrm{where}{\kappa }_{s2}^{\ast }=\frac{{〈<w{\kappa }_{f1}{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}}{{〈<w{>}_{c\stackrel{}{\to }f}〉}_{f\stackrel{}{\to }c}}.$$

(21c)

The final step calculates the the final curvature as

$${\kappa }_{final}=\frac{<w{\kappa }_{f2}{>}_{c\stackrel{}{\to }f}}{<w{>}_{c\stackrel{}{\to }f}}.$$

(21d)

The surface tension is then given as

$$\overrightarrow{F_{st}}=\sigma {\kappa }_{final}\nabla {\alpha }_{sh},$$

(22)

where ${\alpha }_{sh}$ is a sharpened volume fraction of liquid defined in Equation (23).

$${\alpha }_{sh}=\frac{1}{1-C_{sh}}\left[\min \left(\max \left({\alpha }_1,\frac{C_{sh}}{2}\right),1-\frac{C_{sh}}{2}\right)-\frac{C_{sh}}{2}\right],$$

(23)

where $C_{sh}$ is a sharpening coefficient. A value of $C_{sh}$ = 0 reduces ${\alpha }_{sh}$ to ${\alpha }_1$, whereas $C_{sh}$ = 1 provides sharp representation of the interface (which is numerically unstable). We used $C_{sh}$ = 0.98 for static cases and $C_{sh}$ = 0.5 for dynamic cases.

3. Solver Settings

To simplify the treatment of pressure boundary condition and density change across the interface, interFoam uses $p_{rgh}$ which is defined as $p-\rho \overrightarrow{g}·\overrightarrow{x}$, where $\rho \overrightarrow{g}·\overrightarrow{x}$ is the hydrostatic component of pressure [58]. The volume fraction evolution equation (Equation (8)) is solved using the Multidimensional Universal Limiter with Explicit Solution (MULES) algorithm, which preserves the boundedness of volume fraction [61,63]. Once volume fraction is solved, the continuity equation (Equation (6)) and momentum equation (Equation (7)) are solved using the Pressure Implicit with Splitting of Operator (PISO) algorithm [83]. In PISO, a predicted velocity is updated using a pressure correction procedure to advance velocity and pressure fields in time [58,63]. To understand the implementation and solution algorithm of the governing equations (Equations (6)–(8)) in interFoam, please refer to the work by Rusche [58] or Deshpande et al. [63]. The discretisation schemes, solvers and others parameters used to solve the governing equations for all the simulations discussed in this paper are presented in

Table 2. Discretisation schemes.

Modeling Term	Keyword	Scheme	Remarks
Time derivatives	ddtSchemes	Euler	First order implicit method (see [84])
Divergence term	$\nabla ·\left(\rho \overrightarrow{U}\overrightarrow{U}\right)$ $\nabla ·\left(\overrightarrow{U}{\alpha }_1\right)$ $\nabla ·\left(\overrightarrow{U_c}{\alpha }_1(1-{\alpha }_1)\right)$	vanLeerV vanLeer interfaceCompression	Modified vanLeer for vector fields (see [84]) See [85] See [63]
Gradient term	gradSchemes	linear	Operator with ∇ (see [84])
Laplacian term	laplacianSchemes	linear corrected	Operator with ${\nabla }^2$ (see [84])
Others	snGradSchemes interpolationSchemes	corrected linear	Surface normal gradients (see [84]) Interpolates values (see [84])

,

Table 3. Solvers used for the discretised equation.

Equation	Linear Solver	Smoother/Preconditioner	Tolerance
Pressure correction equation	PCG	DIC	10${}^{-20}$ (based on [63])
Momentum equation	smoothSolver	symGaussSeidel	10${}^{-12}$
Volume fraction equation	smoothSolver	symGaussSeidel	10${}^{-12}$

and

Table 4. Other parameters used in solving the discretised equations.

Parameter	Value	Notes
nAlphaCorr	2	Number of ${\alpha }_1$ correction [55]; typically set equal to 1 or 2 for time-dependent flows [86].
nAlphaSubCycles	1	Represents the number of sub-cycles within ${\alpha }_1$ equation [84].
cAlpha ($C_{\alpha }$)	1	Used for interface compression in Equation (9).
MULESCorr	yes	Switches on semi-implicit MULES [87].
nLimiterIter	3	Number of MULES iterations over the limiter [87].
momentumPredictor	no	Controls solving of the momentum predictor; typically set to ’no’ for multiphase and low Reynolds number flows [84].
minIter	1	Minimum number of iterations used in momentum calculation.
nOuterCorrectors	1	PISO algorithm is selected by setting this parameter equal to unity (in PIMPLE algorithm) [84].
nCorrectors	3	The number of times the PISO algorithm solves the pressure and momentum equation in each step; usually set to 2 or 3 [84].
nNonOrthogonalCorrectors	0	Used when meshes are non-orthogonal [84].

, respectively. Under-relaxation factors, if set to less than unity, cause damping of the solution, which can lead to longer computational time for the solution reach to a steady state value. In flow scenarios where there is no steady state solution, using an under-relaxation factor can lead to erroneous results due to under-prediction of the flow variables. We used an under-relaxation factor in the solver equal to one for dynamic cases and 0.9 for static cases. The effect of using an under relaxation factor of one on static cases is also investigated.

4. Validation: Benchmark Test Cases

4.1. Two Dimensional Rising Bubbles

Due to the computational overhead of modelling a three-dimensional rising bubble, we model the buoyancy driven motion of a single bubble as proposed by Hysing et al. [54], Klostermann et al. [55]. The work by Hysing et al. [54] reported benchmarking data such as the bubble shape, rising velocity and circularity for two cases. These benchmarking data are produced based on numerical simulations using codes such as TP2D, FreeLIFE and MoonNMD [54]. In the work by Klostermann et al. [55], the benchmark proposed by Hysing et al. [54] was used to evaluate the VOF solver in OpenFOAM^® (i.e., interFoam) for various meshes.

The computational domain used for the simulation is a rectangle of dimensions 1 m × 2 m where the bubble of diameter 0.5 m was initialised such that the centre of the bubble is at a distance of 0.5 m from the bottom and side walls. As mesh convergence could not be achieved perfectly in previous works [36,55], we used a uniform grid 160 × 320 for the simulations, corresponding to the fine mesh used in [54]. The pressure boundary conditions used in the simulations were zero gradient on the side and bottom walls, and a Dirchlet condition (equal to zero) at the top wall. The volume fraction of fluid used a zero gradient boundary condition on all walls. The velocity boundary conditions used for the simulations were no slip on top and bottom walls, but slip condition was implemented for the side walls. The fluid properties associated with the test cases, which are abbreviated as TC1 and TC2, are tabulated in

Table 5. Physical parameters used for the rising bubble simulations (see [54]).

Cases	${\mathit{\rho }}_1$(kg/m${}^3$)	${\mathit{\rho }}_2$(kg/m${}^3$)	${\mathit{\nu }}_1$(m${}^2$/s)	${\mathit{\nu }}_2$(m${}^2$/s)	$\mathit{\sigma }$(N/m)	$\overrightarrow{\mathit{g}}$(m/s${}^2$) *	$\mathit{R}\mathit{e}$	$\mathit{E}\mathit{o}$	$\mathit{C}\mathit{a}$
TC1	1000	100	10${}^{-2}$	10${}^{-2}$	24.5	(0 −0.98 0)	35	10	0.286
TC2	1000	1	10${}^{-2}$	10${}^{-1}$	1.96	(0 −0.98 0)	35	125	3.571

* $\overrightarrow{g}$ is the reduced gravity as described in [54].

. The maximum Courant number used by the solver was set equal to 0.01 and maximum time step permitted was based on Equations (2) and (3). The test cases are distinguished based on Reynolds ($Re$), Eötvös ($Eo$) and Capillary ($Ca$) numbers, which are defined as

$$Re=\frac{U_{g}L}{{\nu }_1},Eo=\frac{{\rho }_1{U}_g^{2}L}{\sigma },Ca=\frac{Eo}{Re}$$

(24)

with L and $U_g$ being the characteristic length scale (equal to 0.5 m) and characteristic velocity (defined as $\sqrt{|\overrightarrow{g}|L}$), respectively. The bubble shape was obtained at ${\alpha }_1=0.5$ and rising velocity was calculated based on bubble volume averaged vertical component of the velocity vector [54,55]. For validation, we used the the data reported by Klostermann et al. [55] and Hysing et al. [54] (for the predictions by the FreeLIFE solver, which is referred to as ’Benchmark’ in this paper) for a uniform grid of 160 × 320.

The first test case, TC1, corresponds to the case where surface tension effects are dominant [55]. The temporal evolution of the bubble as predicted by the various surface tension models is compared in Figure 1. Due to the stronger surface tension effects, the interface deforms into an ellipsoidal bubble (see Figure 2). The bubble shape (at t = 3 s) predicted by CSF model provides a slightly better agreement to the benchmark data compared to the other surface tension models. The surface tension models also tend to underpredict the position of the bubble at t = 3 s. This underprediction could be attributed to the lower rising velocity (see Figure 3), which has also been reported in previous studies using OpenFOAM [36,54,55]. Although bubble shape and rising velocity provide an overview of the capability of the surface tension models, the quantification of the errors associated with the models was based on the maximum rising velocity ($V_{max}$) and the time at which the $V_{max}$ occurred (tabulated in

Table 6. Benchmark quantities for TC1.

Parameter	CSF	Smoothed CSF	SSF	[55]	Benchmark ([54])
$V_{max}$	0.2366	0.2375	0.2386	0.2365	0.2419
$t\left(V_{max}\right)$	0.9632	0.9491	0.9104	0.9219	0.9270

). The benchmarking data show that SSF model provides a better agreement to the data reported by Hysing et al. [54] (absolute error is less than 2%) and Klostermann et al. [55] (absolute error is nearly 1.5%) in comparison to the other models.

The other test case, TC2, corresponds to a case where the surface tension effects are lower [55]. This results in larger deformation of interface as the bubble evolves (see Figure 4) and eventually forms a skirted bubble that has thin filaments that breaks down into smaller droplets (see Figure 5). Comparing the surface tension models to the benchmark for final bubble shape shows that the models agree quite well (see Figure 5) but there is a difference between the models with respect to the prediction of the skirted part of the bubble (see Figure 5b). Figure 6 shows that the surface tension models in comparison to the benchmark data under-predicts the rise velocity. Comparing with the benchmark, the SSF model provides the closest agreement for $V_{max1}$ (absolute error is nearly 3.5% [54] and less than 0.1% for [55]) and $t\left(V_{max1}\right)$ (absolute error is nearly 3–3.5% for both [54,55]) (see

Table 7. Benchmark quantities for TC2.

Parameter	CSF	Smoothed CSF	SSF	[55]	Benchmark ([54])
$V_{max1}$	0.2434	0.2429	0.2427	0.2431	0.2514
$t\left(V_{max1}\right)$	0.7663	0.7637	0.7502	0.7250	0.7281
$V_{max2}$	0.2302	0.2290	0.2260	0.2302	0.2440
$t\left(V_{max2}\right)$	1.9721	1.9700	1.9729	1.9594	1.9844

). On the other hand, CSF model agrees with the benchmarking data for $V_{max2}$ (absolute error is nearly 5.7% [54] and 0% for [55]) and $t\left(V_{max2}\right)$ (absolute error is nearly 0.6% for both [54,55]) (see Table 7).

In the previous work by Klostermann et al. [55], the spurious currents were reported to be the reason for the error between the benchmark ([54]) and their simulations (for both TC1 and TC2). Thus, the differences in the predictions, for the rising bubble simulations, between the three surface tension models considered in this paper and their departure from the benchmark can also be attributed to spurious currents generated by these models (which is discussed below). For TC2, the larger variation between the surface tension models after the first peak in the transient evolution of the rise velocity (see Figure 6) can be attributed to the differences in the shapes of filament or satellite droplets (based on the work of Yamamoto et al. [36]). Interestingly, there are also some differences in the predictions by the CSF model (for both TC1 and TC2) and the data reported by Klostermann et al. [55], which could be attributed to the difference in the solver settings (e.g., the discretisation schemes, linear solvers and number of iterations) and/or the variations within the different versions of OpenFOAM. The influence of the discretisation schemes on the predicting the flow variables has been previously investigated in [88,89] but further investigation into the effects of other solver settings (e.g., the choice of linear solver and number of iterations) on the solution is required to quantify its effect. As OpenFOAM gets updated, some of the functionalities and/or the algorithms are modified, for example, the artificial interface compression term used in advection of ${\alpha }_1$ (defined in Equation (9)) is computed differently in the older versions of the software (see [55]). To the best knowledge of the authors, no study has reported a comparison of the performance of various versions of OpenFOAM for specific flow scenarios. These settings, especially discretisation schemes and interface compression algorithms, would effect the generation and evolution of spurious currents, which could be the potential source of the discrepancy between our simulations and the data reported in literature.

4.2. Two-Dimensional Capillary Rise

The rise of liquid through a narrow tube or between two parallel plates, which occurs as a consequence of the wetting of the walls by the liquid, is known as capillary rise. As the liquid rises, it reaches a point of equilibrium when the vertical component of the force exerted by surface tension is balanced by the gravitational force acting on the risen liquid column. This equilibrium point (for liquid rising between two vertical parallel plates) is denoted using a height ($h_b$), which can be analytically calculated as

$$h_b=\frac{2\sigma \cos \theta }{\Delta \rho |\overrightarrow{g}|a},$$

(25)

where $\Delta \rho$ is the difference between densities of liquid and gas, and a is the distance between the plates [90].

To study capillary rise, we used a rectangular domain of dimensions 1 mm × 20 mm, where a (defined in Equation (25)) is equal to 1 mm, with a uniform mesh of 20 × 400. This mesh resolution provided the most accurate prediction of capillary rise for the same computational domain while using CSF model in the previous work by Yamamoto et al. [36]. The boundary conditions for velocity field imposes a no slip boundary condition for the walls and pressure based condition (applied to both inlet and outlet) that computes inlet velocity based on the patch-face normal component of the internal-cell velocity and outflow using the zero gradient condition. The volume fraction field uses a zero gradient condition at walls (with a contact angle of 45${}^{\circ }$) and outlet, along with a Dirchlet condition (equal to one) at inlet. The boundary condition for pressure uses a Dirichlet condition (equal to zero) at inlet and outlet whereas the walls use a Neumann boundary condition. The materials properties used for the simulations are described in

Table 8. Physical parameters used for the capillary rise simulations.

${\mathit{\rho }}_1$(kg/m${}^3$)	${\mathit{\rho }}_2$(kg/m${}^3$)	${\mathit{\nu }}_1$(m${}^2$/s)	${\mathit{\nu }}_2$(m${}^2$/s)	$\mathit{\sigma }$(N/m)	$\overrightarrow{\mathit{g}}$(m/s${}^2$) *	$\mathit{Ca}$
1000	1	10${}^{-6}$	1.48 × 10${}^{-5}$	0.07	(0 −10 0)	0.0014

* This value of $\overrightarrow{g}$ is used to study capillary rise by Yamamoto et al. [36].

. The initial volume fraction of liquid in the domain is set such that the liquid–gas interface is at a height of 8 mm from the bottom surface. The maximum time step (which satisfies both Equations (2) and (3)) and maximum Courant number were set equal to 3.5 $\mathsf{\mu }$s and 0.1, respectively.

Once the interface position stabilised (see Figure 7), the capillary height $h_{b,calc}$ was calculated approximately from the volume fraction field as

$$h_{b,calc}=\frac{{\int }_S{\alpha }_{1}dS}{a},$$

(26)

where the numerator is the area occupied by the liquid in the computational domain [36]. The capillary rise height calculated from the simulations is compared to the analytically derived $h_b$ (which was determined to be 9.9 mm using Equation (25)) in

Table 9. Errors associated with the surface tension models on prediction of capillary rise.

Surface Tension Model	${\mathit{h}}_{\mathit{b},\mathit{calc}}$ (mm)	$\mathit{E}\left(\mathit{h}\right)=({\mathit{h}}_{\mathit{b},\mathit{calc}}-{\mathit{h}}_{\mathit{b}})/{\mathit{h}}_{\mathit{b}}$
CSF	9.16	−0.076
Smoothed CSF	Capillary height did not stabilise during simulations (see Figure 7)
SSF	9.26	−0.065

.

Table 9 shows that SSF model provides a better prediction of the capillary rise height compared to CSF model. A previous work by Yamamoto et al. [36] reported an error of 7.7% for a capillary rise model using the CSF model. Interestingly, the Smoothed CSF model could not provide a reliable capillary rise prediction due to the oscillation of the water column (see Figure 7). This discrepancy can be explained based on the evolution of the spurious currents ($U_{sc}$ defined in Equation (27)), which are plotted in Figure 8. The magnitude of spurious currents ($U_{sc}$) generated in the simulations was computed at each time step as

$$U_{sc}=\max \left(\right|\overrightarrow{U}\left|\right).$$

(27)

The periodic growth and decay of the spurious currents in the Smoothed CSF model (see Figure 8) results in the unrealistic motion of the interface whereas the CSF model which has much larger magnitude of spurious currents is much more periodic (see Figure 8), which reduces the net motion of the liquid–gas interface. Compared to CSF and Smoothed CSF models, the spurious current evolution in the SSF model is lowered by nearly two orders of magnitude (see Figure 8).

5. Analysis: Spurious Current

To study the spurious currents generated during the simulations, we simulated a stationary bubble where the effect of gravity was neglected. A bubble of diameter $2R$ was set at the centre of a square domain of dimensions $4R\times 4R$. The properties of the two phases and other physical parameters used for the simulations described in this section are tabulated in

Table 10. Physical parameters used for the simulations in the analysis of spurious current.

${\mathit{\rho }}_1$(kg/m${}^3$)	${\mathit{\rho }}_2$(kg/m${}^3$)	${\mathit{\nu }}_1$(m${}^2$/s)	${\mathit{\nu }}_2$(m${}^2$/s)	$\mathit{\sigma }$(N/m)	$\overrightarrow{\mathit{g}}$(m/s${}^2$)
1000	1	10${}^{-6}$	1.48 × 10${}^{-5}$	0.07	(0 0 0)

. For these simulations, the boundaries were assigned the Dirichlet condition, equal to 101325 Pa, for pressure and zero gradient condition for both ${\alpha }_1$ and $\overrightarrow{U}$. The simulations were run until an end time of 0.05 s to ensure that initial transients (if any) were eliminated with maximum time step calculated based on Equations (2) and (3) along with maximum Courant number of 0.1.

The accuracy of the surface tension models was calculated based on the following parameters: Laplace pressure, magnitude of spurious currents and mass imbalance. For a two-dimensional bubble, the Laplace pressure can be calculated using the Young–Laplace equation as

$$\Delta p_c^{\prime }=\frac{\sigma }{R}.$$

(28)

The Laplace pressure inside the bubble was calculated from the simulation as

$$\Delta p_c=\frac{{\int }_V{\alpha }_{2}pdV}{{\int }_V{\alpha }_{2}dV}-p_0,$$

(29)

where $p_0$ is the operating pressure (which was equal to 101325 Pa). The mean error associated with the Laplace pressure calculated by the various surface tension models was determined as

$$\overline{E}(\Delta p_c)=\frac{\overline{\Delta p_c}-\Delta p_c^{\prime }}{\Delta p_c^{\prime }},$$

(30)

where the overbar represents the time averaged variables.

5.1. Stagnant Bubble of Few Millimetres

In this test case, we modelled a bubble with a radius of 2.5 mm using fluid properties described in Table 10 and under-relaxation factor of 0.9. The computations were performed using a uniform structured grid. The total number of mesh elements and maximum time step (which satisfies both Equations (2) and (3)) used in the simulations are described in

Table 11. Details of mesh and the associated maximum time step calculated based on Equations (2) and (3) used for stationary bubble simulations.

Mesh	Mesh Resolution (mm${}^2$)	Total Number of Cells	$\mathit{R}/\mathit{\delta }\mathit{x}$ *	Maximum Time Step (s)
M0	0.5 × 0.5	400	5	9 × 10${}^{-5}$
M1	0.25 × 0.25	1600	10	3 × 10${}^{-5}$
M2	0.125 × 0.125	6400	20	1 × 10${}^{-5}$
M3	0.083 × 0.083	14400	30	6 × 10${}^{-6}$

* $R/\delta x$ is the ratio of the radius of the bubble and the cell size.

.

To understand how spurious currents occur with various surface tension models, $U_{sc}$ is plotted at t = 0.05 s for the grid described by M3 in Figure 9. In the surface tension models considered in this study, the spurious currents occur around the interface but their magnitudes are much larger in the bubble than outside. To quantify the spurious currents from the simulations, the magnitude of spurious currents and capillary pressure are tabulated in

Table 12. Comparison of spurious currents based on mesh and surface tension models (using an under-relaxation factor of 0.9).

Surface Tension Model	Mesh	$\overline{{\mathit{U}}_{\mathit{sc}}}$	$\mathit{Ca}=\frac{{\mathit{\rho }}_1{\mathit{\nu }}_1\overline{{\mathit{U}}_{\mathit{sc}}}}{\mathit{\sigma }}$	$\overline{\Delta {\mathit{p}}_{\mathit{c}}}$	$\overline{\mathit{E}}(\Delta {\mathit{p}}_{\mathit{c}})$	Mass Imbalance
CSF	M0	0.133	0.002	22.29	−0.20	0
	M1	0.171	0.002	23.03	−0.18	0
	M2	0.174	0.002	24.06	−0.14	0
	M3	0.189	0.003	24.77	−0.12	0
Smoothed CSF	M0	0.096	0.001	24.12	−0.14	0
	M1	0.088	0.001	25.14	−0.10	0
	M2	0.062	0.001	25.19	−0.10	0
	M3	0.049	0.001	26.09	−0.07	0
SSF	M0	0.045	0.001	23.95	−0.14	0
	M1	0.087	0.001	25.12	−0.10	0
	M2	0.036	0.001	25.88	−0.08	0
	M3	0.041	0.001	25.55	−0.09	0

. The spurious currents generated by the surface tension models tends to reduce with finer meshes for both SSF and Smoothed CSF. On the other hand, the increase in spurious current for CSF can be explained based on the dependence on the mesh size ($\Delta x$) is given by

$$C_{\Delta x}\sim \sqrt{\frac{\sigma }{\rho \Delta x}},$$

(31)

where $C_{\Delta x}$ is the magnitude of the spurious velocities (studied for CSF model [63,66]). Equation (31) indicates that smaller mesh sizes result in larger values of spurious currents for CSF model. As shown in Table 12, the Laplace pressure predicted by the surface tension models does not perfectly match $\Delta p_c^{\prime }$ but both Smoothed CSF and SSF provides a better prediction in comparison to CSF.

5.2. Effect of Time Step

The two time step constraints were from Brackbill et al. [66] (Equation (1)) and Deshpande et al. [63] (Equations (2) and (3)). To study the effect of time step constraint, the simulations used a bubble of 2.5 mm with the M3 mesh (see Table 11) and fluid properties described in Table 10 using an under-relaxation factor of 0.9. The maximum time steps ($\Delta t$) used for the simulations are 25 $\mathsf{\mu }$s (based on [66]) and 6 $\mathsf{\mu }$s (based on [63]).

The temporal evolution of $U_{sc}$ is compared for the surface tension models in Figure 10. Using the time step dictated by Deshpande et al. [63], the spurious currents generated by the CSF model are reduced by less than half in comparison to when time step constraint proposed by Brackbill et al. [66] was used. The other models show an absolute difference in the mean spurious current of nearly 7% and 6%, respectively, for the time step constraints (see

Table 13. Comparison of spurious currents for the time stepping constraints based on M3 mesh and surface tension models (while using an under-relaxation factor of 0.9).

Surface Tension Model	$\overline{{\mathit{U}}_{\mathit{sc}}}$ based on Brackbill et al. [66]	$\overline{{\mathit{U}}_{\mathit{sc}}}$ based on Deshpande et al. [63]
CSF	0.395950	0.189170
Smoothed CSF	0.052188	0.048619
SSF	0.038550	0.040984

).

5.3. Effect of Under-Relaxation Factor

To understand the effect of under-relaxation factor, we considered a case which used an under-relaxation factor of unity for modelling the stationary bubble of 2.5 mm with M3 mesh.

Table 14. Comparison of spurious currents based on mesh and surface tension models (using no under-relaxation and time step dictated by Deshpande et al. [63]).

Surface Tension Model	Mesh	$\overline{{\mathit{U}}_{\mathit{sc}}}$	$\mathit{Ca}=\frac{{\mathit{\rho }}_1{\mathit{\nu }}_1\overline{{\mathit{U}}_{\mathit{sc}}}}{\mathit{\sigma }}$	$\overline{\Delta {\mathit{p}}_{\mathit{c}}}$	$\overline{\mathit{E}}(\Delta {\mathit{p}}_{\mathit{c}})$	Mass Imbalance
CSF	M0	0.158	0.002	22.27	−0.20	0
	M1	0.279	0.004	23.09	−0.18	0
	M2	0.510	0.007	24.34	−0.13	0
	M3	0.723	0.010	24.47	−0.13	0
Smoothed CSF	M0	0.154	0.002	24.02	−0.14	0
	M1	0.122	0.002	24.97	−0.11	0
	M2	0.104	0.001	25.21	−0.10	0
	M3	0.075	0.001	26.03	−0.07	0
SSF	M0	0.042	0.001	24.07	−0.14	0
	M1	0.065	0.001	24.86	−0.11	0
	M2	0.033	0.000	26.04	−0.07	0
	M3	0.036	0.001	25.64	−0.08	0

provides a summary of the spurious current and the Laplace pressure in the bubble. Comparison of the results from under-relaxation factor of 0.9 (see Table 12) and 1 (see Table 14) shows that spurious currents generated by Smoothed CSF model is substantially larger when using a larger under-relaxation factor (nearly twice). The SSF model provides the least amount of spurious currents for both the under-relaxation factors and the CSF model generates larger spurious currents with larger mesh density (as described by Equation (31)). It is also worth pointing out that the evolution of spurious currents for the time step constraints provide marginally higher spurious currents for CSF model (0.1% using the time step constraint by Equation (1)) but the Smoothed CSF and SSF models show a spurious current reduction by nearly 10% and 11%, respectively (see

Table 15. Comparison of spurious currents for the time stepping constraints based on M3 mesh and surface tension models (while providing no under-relaxation to the flow variables).

Surface Tension Model	$\overline{{\mathit{U}}_{\mathit{sc}}}$ based on Brackbill et al. [66]	$\overline{{\mathit{U}}_{\mathit{sc}}}$ based on Deshpande et al. [63]
CSF	0.722930	0.723390
Smoothed CSF	0.082800	0.075458
SSF	0.040016	0.035903

). Based on the evolution of spurious currents based on time step constraint, the SSF model generates the least spurious current when compared to Smoothed CSF and CSF models.

6. Conclusions

In the study, we successfully implemented CSF, Smoothed CSF and SSF models in OpenFOAM and compared them based on their ability to simulate a two-dimensional stationary bubble, rising bubbles and capillary rise. The flow scenarios modelled corresponds to a variety of capillary numbers (in the order of 10${}^{-3}$, 0.1 and 1), which is relevant in various industrial processes. The numerical simulations show that:

• For a stationary bubble with a 2.5 mm radius, CSF and SSF models generate the most and least amount of spurious currents, respectively. For the finest mesh used, Smoothed CSF and SSF models reduce spurious currents by nearly one-tenth and one-twentieth of the CSF model (when no under-relaxation factor is used), respectively. When using a lower under-relaxation factor (for the finest mesh), Smoothed CSF and SSF models reduce the spurious currents by approximately one-fourth of the CSF model.
• The time step constraints proposed by Brackbill et al. [66] and Deshpande et al. [63] show that spurious currents generated by the CSF is significantly reduced while using a lower under-relaxation factor. In Smoothed CSF and SSF models, when using the same under-relaxation factor, the time step constraint slightly reduces the spurious currents by 6–7%. Interestingly, when no under-relaxation is used, the CSF model generates marginally larger (nearly 0.1%) spurious currents with the time step constraint proposed by Deshpande et al. [63], but other models show a reduction in spurious current by less than 10%.
• The Laplace pressure in the bubbles predicted by Smoothed CSF and SSF is more accurate with an error of 7–9% for the higher mesh densities than CSF model with negligible imbalance in mass of the phases.
• Although using a lower under-relaxation factor reduces the spurious currents and predicts the Laplace pressure in the stationary bubble (for all the surface tension models considered) quite reasonably, it can adversely effect the accuracy of dynamic cases such as rising bubbles by underestimating the flow variables.
• Using a higher mesh density results in larger spurious currents for CSF model but they are reduced for both Smoothed CSF and SSF models for the static case considered.
• The effect of mesh resolution was studied only for the stationary bubble in this work. For the case of rising bubbles, previous works [36,55], using the CSF model, reported challenges in achieving a mesh independent solution. Similarly, for capillary rise using the CSF model, Yamamoto et al. [36] reported an increasing error when using a finer mesh. The meshes used in this paper correspond to the finest grid (used in FreeLIFE solver) implemented by Hysing et al. [54] and the grid that provided a most accurate model for capillary rise in the work by Yamamoto et al. [36]. We expect similar effects of mesh resolution for both Smoothed CSF and SSF models for dynamic cases, as they are variants of the same formulation. The quantification of these errors will be treated in future work.
• Rising bubbles were successfully modelled using the surface tension models and validated based on the final bubble shape and rising velocities proposed by Hysing et al. [54] and Klostermann et al. [55].
• Modelling the capillary rise with SSF was shown to provide a more accurate representation than the CSF model. Interestingly, the Smoothed CSF could not reliably simulate capillary rise due to spurious currents.

Although the surface tension models considered in this study did not eliminate spurious currents entirely, the comparison provides insights into the limitations of these models. Based on the simulations done in this study, the SSF model seems to provide a versatile surface tension formulation that generates small spurious currents and provides a more accurate representation of various processes in comparison to the standard CSF model.