Implementation of integral surface tension formulations in a volume of fluid framework and their applications to Marangoni flows

TL;DR

The paper addresses the challenge of accurately modeling surface tension with variable gradients in multiphase flows by extending the integral formulation to a volume-of-fluid framework. It introduces three geometry-based implementations—CLSVOF, HF, and HF2D—and integrates them into Basilisk, yielding mass- and momentum-conserving surface-tension forces without relying on a discrete delta function. Through rigorous validation on capillary-wave damping, spurious currents, rising bubbles, and Marangoni-driven phenomena, the work shows that CLSVOF is most accurate but slower, HF2D provides a practical speed–accuracy trade-off, and HF offers the fastest option with some oscillatory behavior; together, they advance robust, scalable Marangoni simulations in 2D/axisymmetric settings. The results have broad implications for simulating thermocapillary and solutal Marangoni flows in engineering and physics, with potential extension to 3D and complex surfactant scenarios.

Abstract

Accurate numerical modeling of surface tension has been a challenging aspect of multiphase flow simulations. The integral formulation for modeling surface tension forces is known to be consistent and conservative, and to be a natural choice for the simulation of flows driven by surface tension gradients along the interface. This formulation was introduced by Popinet and Zaleski [1] for a front-tracking method and was later extended to level set methods by Al-Saud et al. [2]. In this work, we extend the integral formulation to a volume of fluid (VOF) method for capturing the interface. In fact, we propose three different schemes distinguished by the way we calculate the geometric properties of the interface, namely curvature, tangent vector and surface fraction from VOF representation. We propose a coupled level set volume of fluid (CLSVOF) method in which we use a signed distance function coupled with VOF, a height function (HF) method in which we use the height functions calculated from VOF, and a height function to distance (HF2D) method in which we use a sign-distance function calculated from height functions. For validation, these methods are rigorously tested for several problems with constant as well as varying surface tension. It is found that from an accuracy standpoint, CLSVOF has the least numerical oscillations followed by HF2D and then HF. However, from a computational speed point of view, HF method is the fastest followed by HF2D and then CLSVOF. Therefore, the HF2D method is a good compromise between speed and accuracy for obtaining faster and correct results. Keywords: Multiphase flows; Surface tension modeling; Marangoni flows

Figures (15)

• Figure 1: A representation of the control volume $\Omega$ that is cut by a fluid-fluid interface at points $A$ and $B$.
• Figure 2: $(a)$ An example of a computational cell with length and width $\Delta$. In this cell, the diagonal terms of the surface tension stress tensor $\sigma^{xx}$ and $\sigma^{yy}$ are defined at the cell centers and the off-diagonal term $\sigma_{xy}$ and $\sigma_{yx}$ are defined at the vertices. $(b)$ A staggered computational cell $CDEF$ cut by the interface $AB$. Unit vectors tangent to the interface at point $A$ and $B$ are represented as $\mathbf{t}_A$ and $\mathbf{t}_B$ respectively.
• Figure 3: Representation of an interface with in a computational grid with height functions, shown with dotted arrows in $(a)$ Vertical direction $h^x$ and $(b)$ Horizontal direction $h^y$.
• Figure 4: Capillary wave damping test case for all three methods CLSVOF, HF, HF2D for a fixed Laplace number $\rm{La} = 6000$$(a)$ The evolution of the maximum perturbation amplitude. $(b)$ Convergence of the $L2$ norm of the numerical errors with the number of grid points per wavelength of the capillary wave $N/\lambda$. The code to reproduce these results is available at http://basilisk.fr/sandbox/msaini/Marangoni/capwave.c.
• Figure 5: Evolution of the RMS value of spurious velocity currents for a translating circular interface for a fixed Weber number $\rm{We} = 0.4$ and a varying Laplace number $\rm{La} = 120, 1200, 12000$ with $(a)$ the standard CSF method of Basilisk $(b)$ the CLSVOF method $(c)$ the HF2D method. The code to reproduce these results is available at http://basilisk.fr/sandbox/msaini/Marangoni/spuriousMov.c.