\section{2.\ Numerical method} \subsection{2.1.\ Governing equations} We consider the isothermal motion of two incompressible, Newtonian fluid phases and we adopt a whole-domain formulation, which considers the two-fluids system as a single-fluid system with space dependent physical properties (density and viscosity)~\cite{Brackbill1992,Ferrari2013,Ferrari2014,Scardovelli1999}. A single set of Navier--Stokes equations is solved with a unique pressure and velocity field that is shared among the phases. The total mass conservation for the mixture is expressed by the divergence free condition on the velocity: \begin{equation} \nabla \cdot \mathbf{u} = 0 \tag{1} \end{equation} while the momentum conservation is expressed as \begin{equation} \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot (2\mu \mathbf{D}) + \rho \mathbf{g} + \sigma \kappa \mathbf{n}\delta_\Gamma \tag{2} \end{equation} In the above equations, $\mathbf{u}$ is the fluid velocity, $\rho$ the density, $\mu$ the viscosity, $p$ the pressure, $\mathbf{g}$ the gravitational acceleration and $\mathbf{D} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the rate-of-strain tensor. The last term of Eq.~(2) represents the effect of the Laplace pressure, where $\sigma$ is the surface tension, $\kappa$ is the curvature of the interface, $\mathbf{n}$ is the unit normal vector to the interface and $\delta_\Gamma$ is the Dirac function centered on the surface of discontinuity $\Gamma$. The physical properties are defined for the mixture in terms of an indicator function, $\alpha(\mathbf{x})$, that is one in fluid 1 and zero elsewhere, \begin{equation} \begin{aligned} \rho(\mathbf{x}) &= \alpha(\mathbf{x}) \rho_1 + \bigl(1 - \alpha(\mathbf{x})\bigr)\rho_2, \\ \mu(\mathbf{x}) &= \alpha(\mathbf{x}) \mu_1 + \bigl(1 - \alpha(\mathbf{x})\bigr)\mu_2. \end{aligned} \tag{3} \end{equation} \subsection{2.2.\ Volume of fluid method} In the VOF method, the indicator function used to identify each phase is defined as a volume fraction; it represents the ratio of the cell volume occupied by fluid 1 and is $\alpha=1$ if the computational cell is completely filled with the primary phase, $\alpha=0$ if the cell is filled with the secondary phase and it assumes values between zero and one in the interfacial cells. Therefore, the interface is smeared over the cells in which $0<\alpha<1$. The volume fraction is transported as a passive scalar by the flow field with the velocity resulting from the solution of Eqs.~(1) and (2) and, hence, the interface location is updated by solving an advection equation of the form: \begin{equation} \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{u}) + \nabla \cdot \bigl[\alpha(1-\alpha)\mathbf{u}_r\bigr] = 0 \tag{4} \end{equation} The last term in the above equation is an artificial compression term used to limit the numerical diffusion and $\mathbf{u}_r$ is a suitable compression velocity. The factor $\alpha(1-\alpha)$ guarantees the compression only in the interface region, while the boundedness of the volume fraction, and thus the mass conservation, is guaranteed by the MULES (Multidimensional Universal Limiter with Explicit Solution) algorithm implemented in the interFoam solver and used here to discretize Eq.~(4). The artificial compression velocity is given by: \begin{equation} \mathbf{u}_r = \min\bigl[C_f |\mathbf{u}|, \, (|\mathbf{u}|)_{\max}\bigr]\cdot \mathbf{n} \tag{5} \end{equation} The scalar product with $\mathbf{n}$ ensures a compression perpendicular to the interface, whereas the magnitude of the compression is controlled by the coefficient $C_f$ (a value of $C_f=1$ is set for the simulations presented in this paper). The surface force (last term of Eq.~(2)) is calculated using the Continuum Surface Force (CSF) method~\cite{Brackbill1992}, which consists in replacing the force acting at the surface of discontinuity by a continuum force. In the CSF method the Dirac function $\delta_\Gamma$ is approximated (in the sense of distributions) by $|\nabla \alpha|$ and the interface unit normal vector is calculated from the volume fraction as $\mathbf{n}=\nabla \alpha/|\nabla \alpha|$. Therefore, the volume force can be written as: \begin{equation} \mathbf{f}_v = \sigma \kappa \mathbf{n}\delta_\Gamma = \sigma \kappa \nabla \alpha \tag{6} \end{equation} and the curvature is evaluated directly from the gradient of the volume fraction: \begin{equation} \kappa = -\nabla \cdot \mathbf{n} = -\nabla \cdot \left(\frac{\nabla \alpha}{|\nabla \alpha|}\right). \tag{7} \end{equation} \subsection{2.3.\ Coupled VOF and LS for uniform meshes} The VOF approach used to evaluate the curvature and the unit normal vector is known to have poor accuracy. The volume fraction is a step-function that changes sharply across the interface, and standard finite-difference approximations do not converge when applied to highly discontinuous functions~\cite{Khodaparast2015}. [...] Therefore, in order to well explain the limits of the S-CLSVOF algorithm in the case of non-uniform meshes, it is hereafter reviewed in detail. The first step for the construction of a signed distance function is to assign an initial value to the level set using the available volume fraction field after the solution of the advection equation: \begin{equation} \phi_0 = (2\alpha -1)\xi \tag{8} \end{equation} where $\xi=0.75 \Delta x$. [...] The signed distance function to the interface can be recovered by solving the following Eikonal equation: \begin{equation} \left\{ \begin{aligned} |\nabla \phi| &=1,\\ \mathrm{sgn}(\phi)&=\mathrm{sgn}(\phi_0). \end{aligned} \right. \tag{9} \end{equation} If the initial level set $\phi_0$ is already reasonably close, Eq.~(9) can be replaced by: \begin{equation} \left\{ \begin{aligned} \phi_\tau + \mathrm{sgn}(\phi_0)\bigl(|\nabla \phi|-1\bigr)&=0,\\ \phi(\mathbf{x},0)&=\phi_0(\mathbf{x}). \end{aligned} \right. \tag{10} \end{equation} The level set is numerically updated as: \begin{equation} \phi^n = \phi^{n-1} + \Delta \tau\, \mathrm{sgn}(\phi_0)\Bigl[1 - |\nabla \phi^{n-1}|\Bigr]. \tag{11} \end{equation} A proper number of iterations: \begin{equation} n_{\max} = \frac{\epsilon}{\Delta \tau} \tag{12} \end{equation} with $\epsilon=1.5\Delta x$. After the reinitialization, the interface normal vector is computed as $\mathbf{n}=\nabla \phi/|\nabla \phi|$. The volumetric surface force is calculated with the CSF method: \begin{equation} \mathbf{f}_v = \sigma \kappa(\phi)\mathbf{n}(\phi)\delta(\phi) = \sigma \kappa \nabla H, \tag{13} \end{equation} where $\delta(\phi)$ is a smooth Dirac function: \begin{equation} \delta(\phi)= \begin{cases} 0, & |\phi|>\epsilon,\\ \frac{1}{2\epsilon}\Bigl(1+\cos\bigl(\frac{\pi \phi}{\epsilon}\bigr)\Bigr), & |\phi|\le \epsilon. \end{cases} \tag{14} \end{equation} The smooth Heaviside function is: \begin{equation} H(\phi)= \begin{cases} 0, & \phi<-\epsilon,\\ \frac{1}{2}\Bigl[1+\frac{\phi}{\epsilon}+\frac{1}{\pi}\sin\bigl(\frac{\pi \phi}{\epsilon}\bigr)\Bigr], & |\phi|\le \epsilon,\\ 1, & \phi>\epsilon. \end{cases} \tag{15} \end{equation} and we have used \[ \frac{dH}{dx_i} = \frac{dH}{d\phi}\,\frac{d\phi}{dx_i} = \delta(\phi)\,\frac{d\phi}{dx_i}. \] The physical properties can be defined either using Eq.~(3) or the Heaviside function $H$. As also reported by Albadawi et al.~\cite{Albadawi2013}, both methodologies give similar results and the VOF formulation is used here. Note that the method holds in both 2D and 3D domains as long as the mesh is uniform, i.e., squares in 2D and cubes in 3D. \subsection{2.4.\ Extension to non-regular and non-Cartesian grids: flexCLV} The formulation described above is limited to regular Cartesian grids and the extension to non-uniform and unstructured meshes is not straightforward, especially in the case of elements with high aspect ratio or in the presence of highly distorted cells, which is a typical situation when simulating confined two-phase flows at small capillary numbers. The main problems lie on both the initial value assigned to the level set (Eq.~(8)) and the solution of the reinitialization equation (Eq.~(10)), for which parameters depending on the local grid size must be defined ($\xi$, $\Delta \tau$ and $\epsilon$). A special attention must be paid to the definition of the thickness of the interface, $\epsilon$, for which the correct mesh size, in a direction orthogonal to the interface regardless of the mesh topology, must be selected. Values of $\epsilon$ smaller than the desired thickness would cause $\phi$ to vary too steeply (with a consequent drop of accuracy in evaluating $\mathbf{n}$ and $\kappa$ as derivatives of $\phi$), while an $\epsilon$ too large would make $\phi$ too diffused and with $\mathbf{f}_v$ acting on too many cells on both sides of the interface (with a consequent loss of the interface sharpness). Unfortunately, simple combinations of the local mesh size along the three principal axis ($\Delta x$, $\Delta y$ and $\Delta z$) hardly provide the desired thickness of the interface on non-uniform or unstructured meshes, as illustrated below. In Fig.~2 we show some typical issues encountered with the level set reinitialization that may arise when simple combinations of the local mesh size along the three principal directions are taken. We consider the reinitialization of the level set function using Eq.~(11) on a circle of diameter $d_b=0.5$ in a square domain $2d_b \times 2d_b$. The grid is uniform along the y-direction and non-uniform along the x-direction with a maximum aspect ratio of 20. The thickness of the interface is set to $\epsilon=2.5 \Delta \lambda$, which should provide an interface spread over 5 cells. The easiest solution to define a local mesh size, needed by the definitions of $\xi$, $\epsilon$ and $\Delta \tau$, is to define a representative value of the mesh size $\Delta$ (we simply call it $\Delta$ to distinguish between the uniform size, $\Delta x$, defined before) by taking an average of the width and height (and depth in 3D) of the non-uniform cell size (see for example Ningegowda and Premachandran (2014)). However, in the presence of elements with high aspect ratios, this leads to a spreading of the interface towards the smallest cell size, as is visible in Fig.~2(a) and (d), where the interface in the region with a finer grid spreads over about 10 cells, instead of the desired 5 cells. In bubble channel flow simulations for example, the mesh is usually non-uniform close to the channel walls with aspect ratios that can reach values up to 30 (in 2D we would have $\Delta x \approx 30 \Delta y$, see Section~4.2) causing a non-physical spreading of the interface over the thickness of the liquid film. The same considerations are valid if using the minimum grid size over the domain: the interface would be of the desired thickness where the mesh is finer but very sharp where the mesh is coarser, causing difficulties in the evaluation of the level set gradient. This is evident in Fig.~2(b) and (e), where the interface thickness is of 5 cells where the mesh elements are smaller but only 1 cell where the mesh is coarser. In Abadie (2013), the local grid size is defined based on the local minimum size ($\Delta_i=\min(\Delta x_i, \Delta y_i)$ in 2D); however we found this choice not suitable for our case since the approach would yield a very sharp change in $\phi$ in the cells where the direction of the minimum size is not aligned with the interface normal vector (see Fig.~2(c) and (f)). This emphasizes an important aspect: in order to maintain the interface at the desired thickness, $\epsilon$ has to be defined not only according to the local mesh size, but also according to the direction of the unit normal vector $\mathbf{n}$. Building upon these considerations, the flexCLV algorithm proposed in this work is formulated as follows: \begin{enumerate} \item After the solution of the VOF advection equation (Eq.~(4)), the isosurface $\alpha=0.5$ is constructed from the values of the volume fraction. The isosurface is computed using a regularized marching tetrahedra algorithm that combines marching tetrahedra and vertex clustering to generate isosurfaces which are topologically consistent with the data and contain a number of points appropriate to the sampling resolution~\cite{Treece1999}. \item In order to save computational time, only a band of cells near the interface is flagged for the reconstruction of the level set field. Cells are selected if the magnitude of the interface normal vector is higher than a certain threshold. Since at this stage only the volume fraction is available, the normal vector is calculated from $\alpha$. Moreover, to ensure that the flagged zone is larger than the desired thickness of the interface, $\alpha$ is smoothed two times before calculating its gradient by recursively interpolating the volume fraction from the cell center to the cell face (a simple Laplacian smoother is used, see for example Ferrari and Lunati (2013)). \item Instead of initializing $\phi_0$ through Eq.~(8) and solving the reinitialization equation (Eq.~(10)), the level set in the flagged cells is computed as the minimum distance between the center of the specific cell and all the points on the isosurface. If the isosurface is defined by the points $\mathbf{c}_1,\ldots,\mathbf{c}_M$ and we denote with $\mathbf{x}_p$ the center of a generic cell $P$ in the flagged zone, the signed distance function, $\phi(\mathbf{x}_p,t)$ is calculated as: \begin{equation} \phi(\mathbf{x}_p,t)= \begin{cases} \min\limits_j \|\mathbf{x}_p-\mathbf{c}_j\|, & \text{if } \alpha(\mathbf{x}_p,t)>0.5 \\ -\min\limits_j \|\mathbf{x}_p-\mathbf{c}_j\|, & \text{if } \alpha(\mathbf{x}_p,t)<0.5 \\ 0, & \text{if } \alpha(\mathbf{x}_p,t)=0.5 \end{cases} \tag{16} \end{equation} This guarantees that $\phi$ well represents the actual distance from the interface in all the cells where this is needed, i.e.\ where $\mathbf{f}_v$ has to be calculated, so that no reinitialization is needed (Eq.~(10)). This step also avoids the calculation of the grid-dependent parameters defined in Section~2.3 ($\xi$, $\Delta \tau$ and $\epsilon$). \item A characteristic local mesh size, $\Delta$, is defined based on the direction of the unit normal vector at the interface: \begin{equation} \Delta(\mathbf{x}_p,t)= \begin{cases} \dfrac{\sum_{k=1}^{N_f}\bigl|\mathbf{n}(\mathbf{x}_p,t)\cdot\mathbf{S}_k\bigr|\,d_k}{\sum_{k=1}^{N_f}\bigl|\mathbf{n}(\mathbf{x}_p,t)\cdot\mathbf{S}_k\bigr|}, & \text{if } \|\mathbf{n}(\mathbf{x}_p,t)\|\neq0 \end{cases} \tag{17} \end{equation} where the sum is performed over the $k$ faces that bound cell $P$, $N_f$ is the number of faces, $\mathbf{S}_k$ is the face unit normal vector, $d_k$ is the distance between $\mathbf{x}_p$ and the center of the neighboring cell $N$ that shares the face $k$ and $\mathbf{n}(\mathbf{x}_p,t)=\nabla \alpha/|\nabla \alpha|$ is the interface unit normal vector at cell $P$ (see Fig.~3). Eq.~(17) evaluates $\Delta$ as the average distance between $\mathbf{x}_p$ and all the adjacent cell centers, weighted by the scalar product between the interface unit normal vector and the cell face unit normal vector, i.e.\ the cosine of the angle between $\mathbf{n}$ and $\mathbf{S}_k$. This yields an estimation of the mesh size in a direction orthogonal to the interface, as desired. The condition ``if $\|\mathbf{n}(\mathbf{x}_p,t)\|\neq0$'' guarantees that $\Delta$ is calculated only for the flagged cells. Note that in the particular case of an orthogonal non-uniform grid (as that in Fig.~2 for example), the cell face unit normal vector and the vector connecting two adjacent cell centers (whose magnitude is $d_k$) are parallel, whereas in the case of non-orthogonal grids (as that in Fig.~3) they are, in general, not parallel and only the component of $d_k$ perpendicular to the face contributes to the evaluation of $\Delta$. In the case of a uniform orthogonal grid, Eq.~(17) returns a constant mesh size for every angle between $\mathbf{n}$ and $\mathbf{S}_k$. The thickness of the interface is then defined as $\epsilon=1.5\Delta$ (as in the cases illustrated in Figs.~2 and 4). The value of $\epsilon=1.5\Delta$ is the traditional value in level set method~\cite{Sussman2000,Sussman1994}, which guarantees a good compromise between interface sharpness and smoothness of the physical properties at the interface. Therefore, this methodology guarantees that the interface is always of the desired thickness (measured in terms of number of mesh cells) independently of the global mesh architecture, which provides an accurate surface tension force representation. \item The Heaviside function $H$ is calculated using Eq.~(15). In the regions of the mesh far from the interface, where $\|\mathbf{n}\|=0$, $H$ is simply set to zero or to one depending on which fluid is present. \end{enumerate} The accuracy of the proposed reconstruction algorithm is shown in Fig.~4 for a grid with a uniform size along the y-direction and a non-uniform size along the x-direction, with an aspect ratio between the smallest and the largest element of 20, analogous to the case depicted in Fig.~2. Fig.~4(a) and (c) show the resulting Heaviside field after the reconstruction, from which it is evident that the interface is correctly represented with the desired thickness both in regions where the mesh is finer and coarser. Fig.~4(b) shows $\Delta$ calculated with Eq.~(17) in the flagged region around the interface: depending on the direction of the interface normal vector, the correct combination of $\Delta x$ and $\Delta y$ is selected. The characteristic mesh size given by Eq.~(17) can be, in principle, used to compute the grid-dependent parameters ($\xi$, $\Delta \tau$ and $\epsilon$) needed to solve the reinitialization equation defined in Section~2.3 for the S-CLSVOF. The level set can be initialized using Eq.~(8), by replacing $\Delta x$ with the local computed value of $\Delta$, in those cells where the interface normal vector is defined and a guess value of the level set can be assigned in cells where the normal vector is zero (a large number for example). However, problems arise in the connection region between the interfacial and non-interfacial cells. Since the new values of the level set in the reinitialization step depend on its gradient, that is on neighboring cells, the guess value of the level set assigned where $\|\mathbf{n}\|=0$ may cause very large gradients in those regions and would deteriorate the solution in the following iterations. On the other hand, the flexCLV algorithm described above, although more computationally expensive than the S-CLSVOF algorithm due to the reconstruction of $\phi$ from the $\alpha=0.5$ isosurface, allows to obtain a level set function that does not depend on any adjustable parameter and gives the desired thickness of the interface even in the case of very complex 2D or 3D meshes, including unstructured meshes (see later Section~3.2). \subsection{2.5.\ Accuracy of the level set reconstruction algorithm} Below, we measure the accuracy of the signed distance function reconstruction for the algorithm described in Section~2.4. The level set is reconstructed on a circular interface of diameter $d_b=0.5$, centered in a square domain $2d_b \times 2d_b$. Two grid architectures are used: a regular mesh of different sizes with $\Delta x = \Delta y = \{d_b/10,\, d_b/20,\, d_b/40,\, d_b/80,\, d_b/160\}$ where the S-CLSVOF~\cite{Albadawi2013} solution represents the benchmark case, and a non-uniform grid with a uniform mesh size equal to $\Delta y$ along the y-direction and non-uniform along the x-direction with an aspect ratio between the smallest and largest element of 20. Despite the different approach for calculating the level set function, the proposed flexCLV algorithm is based on similar concepts of the S-CLSVOF method and a similar order of convergence is therefore expected on Cartesian meshes. The distance function is reconstructed with the algorithm described above and the error $E_2$ with respect to the exact distance is calculated as: \begin{equation} E_2 = \sqrt{\frac{1}{N}\,\sum_{i=1}^{N}\bigl(\phi_i - d^*_i\bigr)^2} \tag{18} \end{equation} where $d^*_i$ is the exact distance from the interface at a specific cell $i$ and $N$ is the number of cells in the flagged region. The distance function error is shown as a function of the grid size in Fig.~5 for both uniform and non-uniform grids. The error obtained with the S-CLSVOF algorithm in the uniform mesh is also included as a reference solution. The new flexCLV algorithm exhibits a first order convergence rate, in agreement with the reference S-CLSVOF method~\cite{Albadawi2013}, independent of the mesh topology. The results demonstrate that the geometrical procedure introduced in Section~2.4 for the reconstruction of the level set function is robust enough to reconstruct the distance function with high accuracy on both regular and non-regular meshes. In Section~3 the methods will be validated against static and dynamic test cases. \subsection{2.6.\ Solution procedure} The flexCLV algorithm is implemented into the OpenFOAM framework by modifying the interFoam solver. The governing equations are discretized by the Finite Volume method on collocated structured or unstructured grids. The pressure-velocity coupling is handled by the Pressure Implicit Splitting of Operator (PISO) algorithm~\cite{Issa1986} and the equations are integrated with an accuracy of first order in time and second order in space. The discretization schemes, solvers and convergence criterion used in this work are shown in Table~1. The solution procedure for the flexCLV solver can be summarized as follows: \begin{enumerate} \item Define vectors and scalar field. Initialize $\alpha(\mathbf{x},0)$ and define physical properties with Eq.~(3). Initialize $\mathbf{u}(\mathbf{x},0)$, $p(\mathbf{x},0)$ and $\alpha(\mathbf{x},0)$. \item Set the time-step $\Delta t = C_0 \Delta x / |\mathbf{u}|_{\max}$, where $C_0$ is the Courant number and $|\mathbf{u}|_{\max}$ is the maximum velocity in the domain. $\Delta t$ is recalculated at every time-step to keep the Courant number below the specified value. \item Solve the advection equation for the volume fraction, Eq.~(4), and update physical properties with Eq.~(3). \item Compute the $\alpha=0.5$ isosurface. \item Calculate the distance from the $\alpha=0.5$ isosurface and compute $\Delta$ with Eq.~(17). Then, compute the new values of the Heaviside function (Eq.~(15)), curvature and surface force (Eq.~(13)). \item Solve mass and momentum conservation equations (Eqs.~(1) and (2), respectively) to obtain the new pressure and velocity using the PISO algorithm. An intermediate velocity, $\mathbf{u}^*$, is evaluated from a semi-discretized form of the momentum equation integrated over the control volume of a generic cell $P$, whose matrix form can be written as: \begin{equation} a_p \mathbf{u}^*_p + \sum_n a_n \mathbf{u}^*_n = \mathbf{r} - \nabla p + \mathbf{f}_v \tag{19} \end{equation} where $\mathbf{u}_n$ is the velocity in the neighboring cells of $P$, while $a_p$ and $a_n$ are the diagonal and off-diagonal matrix coefficients resulting from the discretization of the velocity dependent terms of Eq.~(2) and $\mathbf{r}$ includes all the source terms. By defining $H[\mathbf{u}^*] = \mathbf{r} - \sum_n a_n \mathbf{u}^*_n$, Eq.~(19) can be rewritten as: \begin{equation} \mathbf{u}^*_p = \frac{H[\mathbf{u}^*]}{a_p} + \frac{1}{a_p}\bigl(-\nabla p + \mathbf{f}_v\bigr) \tag{20} \end{equation} Eq.~(20) can be solved for the intermediate velocity $\mathbf{u}^*_p$ by treating the surface tension and the pressure gradient as an explicit source term, whereas all the other terms are treated as implicit. \item Since the intermediate velocity does not satisfy mass conservation, the continuity condition is imposed on $\mathbf{u}^*$. The discretized form of Eq.~(1) is $\sum_f \mathbf{u}_f \cdot \mathbf{S}_f=0$, where $\mathbf{S}_f$ is the face area vector and $\mathbf{u}_f$ is interpolated from Eq.~(20) on the cell faces. By imposing the continuity constraint on Eq.~(20) and rearranging the terms we obtain a Poisson equation for the pressure: \begin{equation} \sum_f \left(\left(\frac{\nabla p}{a_p}\right)_f \cdot \mathbf{S}_f \right) = \sum_f \left( \frac{H[\mathbf{u}^*]}{a_p} + \frac{\mathbf{f}_v}{a_p} \right)_f \cdot \mathbf{S}_f \tag{21} \end{equation} Special care must be paid here in the discretization of the pressure gradient and the surface force to avoid any inconsistent discrete balance between pressure and surface tension. This is accounted for in OpenFOAM by computing both the $\nabla p$ and $\nabla H$ at the cell faces using the same discrete operator. \item The pressure solution is used to correct the velocity field using Eq.~(20). Since the term $\bigl(H[\mathbf{u}^*]+\mathbf{f}_v\bigr)$ depends on the previously computed velocity $\mathbf{u}^*$ and is added explicit to Eq.~(21), iterations are needed to reduce the discretization error of the PISO splitting (3 corrections are applied in this work). \item Advance to the next time and repeat steps 2--9. \end{enumerate} \caption{Discrete representation of a circular bubble of diameter $d$ on a uniform Cartesian grid: a) volume fraction and b) level set. The interface (black line) is represented by the iso-contour $\alpha=0.5$ and $\phi=0$, respectively.} \caption{Issues on the level set reinitialization in the case of a non-uniform grid, for a circle of diameter $d_b=0.5$ in a square domain $2d_b \times 2d_b$. The grid is uniform along the y-direction and non-uniform along the x-direction with a maximum aspect ratio of 20. The color function is initialized as $\alpha=1$ in cells whose center lies inside the circle and $\alpha=0$ otherwise. The initial value of the level set is calculated using Eq.~(8) and it is then reinitialized using Eq.~(11). The Heaviside function $H$ is calculated from Eq.~(15). The value of $\Delta$, which represents the typical grid size, is calculated using three different approaches: for case (a) and (d), $\Delta$ is taken as the average of the width and height of the local non-uniform cell size ($\Delta_i=0.5\cdot(\Delta x_i+\Delta y_i)$) for $i=1:N$, where $N$ is the total number of cells in the square domain; for case (b) and (e), $\Delta$ is taken constant as the global minimum between widths and heights of every cell ($\Delta=\min(\Delta x_i,\min(\Delta y_i))$); for case (c) and (f), $\Delta$ is calculated using the local minimum in every cell ($\Delta_i=\min(\Delta x_i,\Delta y_i)$). For all the cases we set $\xi=0.75\Delta$, $\Delta \tau=0.1\Delta$ and $\epsilon=2.5\Delta$ (note that here we use 5 cells for the interface thickness just for visualization purpose, for all the other simulations presented in this work $\epsilon=1.5\Delta$). Images (a), (b) and (c) show a 2D map of the Heaviside function after the reinitialization, whereas (d), (e) and (f) show a 1D plot of the Heaviside function along the dashed line (or lines) indicated in (a), (b) and (c). When $\Delta$ is calculated using the average between $\Delta x$ and $\Delta y$ (case a and d) the interface is very diffused in the regions where the grid is finer; if $\Delta$ is calculated as the global minimum size in the domain, the interface results of the correct thickness (5 cells) where the mesh is finer but very sharp where the mesh is coarser (case b and e); finally, if $\Delta$ is calculated as the local minimum between $\Delta x$ and $\Delta y$ (case c and f) the interface results of the correct thickness only if $\Delta$ is close enough to the mesh size along the direction of the interface normal vector.} \caption{Parameters used in the flexCLV algorithm (adapted from Weller (2008)).} \caption{Level set reconstruction in the case of a non-uniform grid for a circle of diameter $d_b=0.5$ in a square domain $2d_b \times 2d_b$. The grid is uniform along the y-direction and non-uniform along the x-direction with an aspect ratio of 20. The color function is initialized as $\alpha=1$ in cells whose center lies inside the circle and $\alpha=0$ otherwise. The level set is reconstructed from the $\alpha=0.5$ isosurface only in a region around the interface to reduce computational costs with the algorithm described in Section~2.4. a) the entire $H$ field calculated with Eq.~(15) after the reconstruction of the level set, b) $\Delta$ calculated with Eq.~(17) and c) 1D plot of $H$ along the black and red dashed lines indicated in a). $\epsilon$ has been set equal to $2.5\Delta$ to improve visualization. The proposed algorithm gives the desired thickness of the interface (approximately 5 cells) in both the regions with finer and coarser meshes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)} \caption{Distance function error ($E_2$) computed around a circle of radius $d_b=0.5$ centered in a square domain $2d_b \times 2d_b$, as a function of the grid size $\Delta x$. Circles show the convergence in the case of uniform grids, while squares are in the case of non-uniform grids. The new flexCLV algorithm exhibits a first order convergence rate, in agreement with the reference S-CLSVOF method, independent of the mesh topology.} \begin{table}[htbp] \caption{Details of the discretization schemes, solvers and convergence criteria adopted for the simulations. MULES: Multidimensional universal limiter for explicit solution; TVD: Total Variation Diminishing; PbiCG: Preconditioned bi-conjugate gradient solver for asymmetric matrices; DILU: Simplified diagonal-based incomplete LU preconditioner for asymmetric matrices; PCG: Preconditioned conjugate gradient solver for symmetric matrices; DIC: Simplified diagonal-based incomplete Cholesky preconditioner for symmetric matrices; GAMG: Geometric agglomerated algebraic multigrid preconditioner. The reader is referred to the OpenFOAM user-guide for more information on the different solvers (\url{www.openfoam.com}).} \vspace{0.5em} \centering \begin{tabular}{|l|l|} \hline \multicolumn{2}{|c|}{\textbf{OpenFOAM}} \\ \hline $\partial \alpha / \partial t$ & MULES, first order explicit \\ $\partial \rho \mathbf{u} / \partial t$ & Euler, first order implicit \\ $\nabla \cdot \left(\alpha \mathbf{u}\right), \, \nabla \cdot \left[\alpha(1-\alpha)\mathbf{u}_r\right]$ & Van Leer, TVD, second order \\ $\nabla \cdot \left(\rho \mathbf{u}\mathbf{u}\right)$ & LimitedLinear, TVD, second order \\ $\nabla$ & Gauss linear (uniform/non-uniform Cartesian grid) \\ $\nabla^h,\,\nabla^2$ & Least Square (unstructured grids) \\ $p - u$ coupling & Central-difference with explicit non-orthogonal correction (if needed) \\ Solver & PISO (3 corrections)\\ & $\mathbf{u}$: PbiCG + DILU \\ & $p$: PCG + DIC or PCG + GAMG \\ Convergence Criterion & $\mathbf{u}$: Resid. $<10^{-6}$ \\ & $p$: Resid. $<10^{-7}$ \\ \hline \end{tabular} \end{table} ---Fig. 1. Discrete representation of a circular bubble of diameter d on a uniform Cartesian grid: a) volume fraction and b) level set. The interface (black line) is represented by the iso-contour $\alpha = 0.5$ and $\phi = 0$, respectively. **Chart Description:** The image contains two plots, labeled a) and b), displayed side-by-side on a uniform Cartesian grid, showing a circular feature. Both plots use a color map to represent scalar values and include a black line indicating an interface. **Chart a):** * Type: Contour plot/Heatmap on a grid. * Label: a). * Variable: $\alpha$ (alpha). * Color Bar: Positioned above the plot. Scale ranges from 0.0 (blue) to 1.0 (red). Major tick labels are 0.0, 0.3, 0.5, 0.8, 1.0. * Content: Shows the distribution of $\alpha$ on the grid. The central region is red/yellow (values near 1.0), transitioning outwards through green, cyan, and blue (values towards 0.0). * Grid: Uniform Cartesian grid. * Black Line: A circular line representing the interface. According to the caption, this corresponds to the iso-contour $\alpha = 0.5$. **Chart b):** * Type: Contour plot/Heatmap on a grid. * Label: b). * Variable: $\phi$ (phi). * Color Bar: Positioned above the plot. Scale ranges from -0.008 (blue) to 0.004 (red). Major tick labels are -0.008, -0.005, -0.002, 0.000, 0.003, 0.004. * Content: Shows the distribution of $\phi$ on the grid. The central region is red/yellow (positive values), transitioning outwards through green (near zero) to blue (negative values). * Grid: Uniform Cartesian grid. * Black Line: A circular line representing the interface. According to the caption, this corresponds to the iso-contour $\phi = 0$. **Extracted Content:** **Figure Description:** The figure contains six subplots (a-f) displaying data likely related to a circular object or distribution. Subplots a, b, and c show 2D grid data represented as heatmaps. Subplots d, e, and f show 1D profiles extracted from the 2D data. **Subplot a):** * **Type:** 2D Heatmap / Grid Plot * **Title/Label:** a) * **X-axis Label:** X * **Y-axis Label:** Y * **X-axis Range:** 1.0 to 1.6 * **Y-axis Range:** 1.2 to 1.8 * **Color Bar Label:** H * **Color Bar Scale:** 0 to 1.0 (labeled every 0.1) * **Data Representation:** Colors on a grid representing values of H. Yellow indicates high values (around 1), transitioning through orange, light blue, dark blue to low values (around 0) in the background. A roughly circular region of high H is shown. * **Annotation:** A horizontal dashed black line at Y=1.5. **Subplot b):** * **Type:** 2D Heatmap / Grid Plot * **Title/Label:** b) * **X-axis Label:** X * **Y-axis Label:** Y * **X-axis Range:** 1.0 to 1.6 * **Y-axis Range:** 1.2 to 1.8 * **Color Bar Label:** H * **Color Bar Scale:** 0 to 1.0 (labeled every 0.1) * **Data Representation:** Colors on a grid representing values of H. Yellow indicates high values (around 1). The shape appears more blocky and aligned with the grid compared to subplot a). The transition from high H to low H (dark blue background) is sharper. * **Annotation:** A horizontal dashed black line at Y=1.5. **Subplot c):** * **Type:** 2D Heatmap / Grid Plot * **Title/Label:** c) * **X-axis Label:** X * **Y-axis Label:** Y * **X-axis Range:** 1.2 to 1.9 * **Y-axis Range:** 1.2 to 1.8 * **Color Bar Label:** H * **Color Bar Scale:** 0 to 1.0 (labeled every 0.1) * **Data Representation:** Colors on a grid representing values of H. Yellow indicates high values (around 1), transitioning through orange, light blue, dark blue to low values (around 0) in the background. A roughly circular region of high H is shown, similar to a) but over a slightly different range. * **Annotation:** A horizontal dashed black line at Y=1.5. A vertical dashed red line at X=1.5. **Subplot d):** * **Type:** Line chart with data points * **Title/Label:** d) * **X-axis Label:** X * **Y-axis Label:** H * **X-axis Range:** 1.0 to 1.8 * **Y-axis Range:** 0 to 1.0 * **Data Representation:** Black circular data points connected by dashed black lines. Shows the value of H along the horizontal dashed line (Y=1.5) in subplot a). The profile shows a sharp increase from 0 to 1, stays at 1 for a range of X values, then a sharp decrease back to 0. **Subplot e):** * **Type:** Line chart with data points * **Title/Label:** e) * **X-axis Label:** X * **Y-axis Label:** H * **X-axis Range:** 1.0 to 1.8 * **Y-axis Range:** 0 to 1.0 * **Data Representation:** Black circular data points connected by dashed black lines. Shows the value of H along the horizontal dashed line (Y=1.5) in subplot b). The profile shows a sharp increase from 0 to 1, stays at 1 for a range of X values, then a sharp decrease back to 0. **Subplot f):** * **Type:** Line chart with data points * **Title/Label:** f) * **X-axis Label:** X, Y * **Y-axis Label:** H * **X-axis Range:** 1.1 to 1.9 * **Y-axis Range:** 0 to 1.0 * **Data Representation:** Two series of data points connected by dashed lines. * Black circular data points connected by dashed black lines: Represents the profile of H along the horizontal dashed line (Y=1.5) in subplot c) as a function of X. * Red circular data points connected by dashed red lines: Represents the profile of H along the vertical dashed line (X=1.5) in subplot c) as a function of Y (plotted on the same axis as X). * **Description:** Both profiles show a sharp rise from 0 to 1, a region where H is approximately 1, and then a sharp drop back to 0. The red profile appears slightly wider than the black profile, or shifted, indicating the vertical dimension of the high-H region might be slightly different or the profile line position is different. **Chart/Diagram Description:** * **Type:** Geometric diagram representing two adjacent polyhedral volumes or cells, likely part of a computational mesh. * **Main Elements:** * Two polyhedral shapes are shown, sharing a common face. * One common face between the two shapes is shaded grey. * A point labeled 'P' is located within or near one polyhedron. A grey circle marks its exact position. * A point labeled 'N' is located within or near the adjacent polyhedron. A grey circle marks its exact position. * A straight dashed line connects the point P to the point N. * A label 'f_k' is placed near the shaded face, likely representing a flux or flow across the face. * A vector labeled 'S_k' is shown originating from the shaded face and pointing towards N. It is represented by an arrow. * A vector labeled 'd_k' is shown originating from the shaded face, pointing towards N. It is represented by an arrow. The label 'd_k' also has a small line segment indicating it is perpendicular to the shaded face. * Dashed lines are used to indicate edges or lines that are hidden from view. **Textual Information:** * **Labels within the diagram:** * P * f_k * S_k * d_k * N * **Caption/Title:** * Fig. 3. Parameters used in the flexCLV algorithm (adapted from Weller (2008)). **Figure Title:** Fig. 4. Level set reconstruction in the case of a non-uniform grid for a circle of diameter $d_b = 0.5$ in a square domain $2d_b \times 2d_b$. The grid is uniform along the y-direction and non-uniform along the x-direction with an aspect ratio of 20. The color function is initialized as $\alpha = 1$ in cells whose center lies inside the circle and $\alpha = 0$ otherwise. The level set is reconstructed from the $\alpha = 0.5$ isosurface only in a region around the interface to reduce computational costs with the algorithm described in Section 2.4. a) the entire H field calculated with Eq. (15) after the reconstruction of the level set, b) $\Delta$ calculated with Eq. (17) and c) 1D plot of H along the black and red dashed lines indicated in a). $\epsilon$ has been set equal to 2.5$\Delta$ to improve visualization. The proposed algorithm gives the desired thickness of the interface (approximately 5 cells) in both the regions with finer and coarser meshes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) **Chart Descriptions:** **a) Image of Fig. 4** * **Type:** Heatmap/Pseudocolor plot on a 2D grid. * **Main Elements:** * Displays a circular region in yellow/orange representing high values, surrounded by regions of decreasing values (green, blue, dark blue). * A color bar is present on the right showing the mapping from color to values of H (from 0 to 1). * A horizontal black dashed line is at Y = 0.5. * A vertical red dashed line is at X = 0.5. * **Coordinate Axes:** * X-axis labeled "X", range approximately 0.1 to 0.9. * Y-axis labeled "Y", range approximately 0.1 to 0.9. * **Labels and Annotations:** Subfigure label "a)". Color bar labeled "H". **b) Image of Fig. 4** * **Type:** Heatmap/Pseudocolor plot on a 2D grid. * **Main Elements:** * Displays an annular region in yellow/orange/green representing non-zero values, surrounded by dark blue regions representing values close to zero. The annular region appears somewhat blocky due to the grid structure. * A color bar is present on the right showing the mapping from color to values of $\Delta$ (from 0 to 0.06). * **Coordinate Axes:** * X-axis labeled "X", range approximately 0.1 to 0.9. * Y-axis labeled "Y", range approximately 0.1 to 0.9. * **Labels and Annotations:** Subfigure label "b)". Color bar labeled "$\Delta$". **c) Image of Fig. 4** * **Type:** 1D Line plot with markers. * **Main Elements:** * Two data series are plotted: * Black dots connected by a dashed black line. * Red dots connected by a dashed red line. * Both series represent a step-like function, starting at 0, transitioning to 1, staying at 1, transitioning back to 0. * The transition region is steep. * The red line/points are slightly shifted horizontally compared to the black line/points in the transition regions. * Black dots are plotted at approximately X=0.2, 0.22, 0.24, 0.76, 0.78, 0.8. * Red dots are plotted at approximately X=0.25, 0.28, 0.31, 0.69, 0.72, 0.75. Many red dots are also at Y=0 and Y=1 for X < 0.25 and X > 0.75 respectively. * **Coordinate Axes:** * X-axis labeled "X", range from 0 to 1. * Y-axis labeled "H", range from 0 to 1. * **Labels and Annotations:** Subfigure label "c)". **Textual Information:** **Legend:** I order flexCLV (uniform grid) flexCLV (non-uniform grid) S-CLSVOFF **Y-axis Label:** $E_2$ **X-axis Label:** $\Delta x$ **Figure Caption:** Fig. 5. Distance function error ($E_2$) computed around a circle of radius $d_b = 0.5$ centered in a square domain $2d_b \times 2d_b$, as a function of the grid size $\Delta x$. Circles show the convergence in the case of uniform grids, while squares are in the case of non-uniform grids. The new flexCLV algorithm exhibits a first order convergence rate, in agreement with the reference S-CLSVOFF method, independent of the mesh topology. **Chart/Diagram Description:** * **Type:** Log-log scatter plot with connecting lines and a reference line. * **Main Elements:** * **Coordinate Axes:** * X-axis: Labeled "$\Delta x$", logarithmic scale, ranging from approximately 10⁻³ to 10⁻¹. * Y-axis: Labeled "$E_2$", logarithmic scale, ranging from approximately 10⁻⁴ to 10⁻¹. * **Data Series:** * "I order": Represented by a black dashed line. This line shows a linear relationship on the log-log scale, representing a first-order convergence rate (slope of 1). * "flexCLV (uniform grid)": Represented by a green solid line connecting green circular markers. The points generally follow the trend of the "I order" line. * "flexCLV (non-uniform grid)": Represented by a red solid line connecting red square markers. The points are very close to the "flexCLV (uniform grid)" points and follow a similar trend. * "S-CLSVOFF": Represented by a blue dashed line connecting blue square markers. The points and line are very close to the flexCLV lines, also following the trend of the "I order" line. * **Data Points:** Data points are plotted for $\Delta x$ values approximately around 0.003, 0.005, 0.01, and 0.02. For each method, the corresponding $E_2$ values are plotted, showing a decrease as $\Delta x$ decreases. * **Trend:** All three computational methods (flexCLV uniform, flexCLV non-uniform, and S-CLSVOFF) show a first-order convergence rate, as their data points lie close to the "I order" line on the log-log plot. The flexCLV results for uniform and non-uniform grids are very close to each other and to the S-CLSVOFF results.