A hybrid particle volume-of-fluid method for curvature estimation in multiphase flows Note: Published in \@journal
   https://doi.org/10.1016/j.ijmultiphaseflow.2020.103209 Journal: International Journal of Multiphase Flow

We consider the numerical model describing two-component incompressible flows with surface tension available in the open-source solver Basilisk [1, 29]. The system consists of the Navier-Stokes equations for the mixture velocity $\bm{u}$ and pressure $p$

and the advection equation for the volume fraction $\alpha$

with density $\rho=(1-\alpha)\rho_{1}+\alpha\rho_{2}$, dynamic viscosity $\mu=(1-\alpha)\mu_{1}+\alpha\mu_{2}$, gravitational acceleration $\bm{g}$ and constant material parameters $\rho_{1}$, $\rho_{2}$, $\mu_{1}$ and $\mu_{2}$. The surface tension force is defined as $\bm{f}_{\sigma}=\sigma\kappa\nabla\alpha$ with the surface tension coefficient $\sigma$ and interface curvature $\kappa$. A finite volume discretization is based on Chorin’s projection method for the pressure coupling [12] and the Bell-Colella-Glaz scheme [8] for convective fluxes. The advection equation is solved using the volume-of-fluid method PLIC with piecewise linear reconstruction [43] where the normals are computed using the mixed Youngs-centered scheme which is a combination of Youngs’ scheme and the height functions. The approximation of the surface tension force is well-balanced [18] (i.e. the surface tension force is balanced by the pressure gradient if the curvature is uniform) and requires face-centered values of the curvature, which are computed as the average over the neighboring cells. We refer to [4] for further details about the algorithm.

To compute the cell-centered interface curvature, we use the proposed algorithm described in the following sections. Our implementation of the curvature estimator is available online¹¹ 1 Implementation with Basilisk: https://cselab.github.io/aphros/basilisk_partstr.zip and can be directly used in Basilisk for simulations of multiphase flows with surface tension. We also provide a visual web-based demonstration of the method²² 2 Visual demonstration: https://cselab.github.io/aphros/curv.html and a reference implementation in Python³³ 3 Implementation in Python: https://cselab.github.io/aphros/curv.py.

This section introduces a particle method for estimating the curvature of an interface represented by a set of line segments, which are given by the piecewise linear reconstruction of the interface as described in the next section. The particles are constrained to circular arcs and equilibrate on the interface due to attraction forces.

Consider a set of $N_{L}$ line segments $[\bm{a}_{l},\bm{b}_{l}],\;l=1,\dots,N_{L}$ with endpoints $\bm{a}_{l},\bm{b}_{l}\in\mathbb{R}^{2}$. One is distinguished as the target line segment $[\bm{a}^{*},\bm{b}^{*}]$ at which the curvature is to be estimated. The union of all line segments is denoted as

Given an odd number $N$, we introduce a string of $N$ particles $\bm{x}_{i}\in\mathbb{R}^{2},\;i=1,\dots,N$ and denote the index of the central particle as $c=(N+1)/2$. The particles are constrained to circular arcs, which leads to parametrization of their positions

where $\bm{p}$ is the origin, $\phi$ is the rotation angle, $\theta$ is the bending angle, $h_{p}={H_{p}h}/(N-1)$ is the distance between neighboring particles and $\bm{e}(\psi)=\cos{\psi}\,\bm{e}_{x}+\sin{\psi}\,\bm{e}_{y}$. Parameter $H_{p}$ defines the length of the string relative to the mesh step $h$. Values of $N$ and $H_{p}$ are chosen in Section 2.7. The curvature of the circular arc is related to the bending angle as

We define the force attracting a particle at position $\bm{x}$ to the nearest point on the interface

where $\eta\in[0,\,1]$ is a relaxation factor. Figure 1 illustrates the parametrisation of positions and computation of forces.

Further derivations use vector notation of the form $\bm{X}=[\bm{x}_{i},\;i=1,\dots,N]$, where $\bm{x}_{i}$ is component $i$ of $\bm{X}$. The scalar product is defined as $\bm{X}\cdot\bm{Y}=\sum\limits_{i=1}^{N}\bm{x}_{i}\cdot\bm{y}_{i}$. In this notation, the positions and forces combine to

The particles evolve until equilibration according to the following iterative procedure. Initially, the particles are arranged along the target line segment $[\bm{a}^{*},\bm{b}^{*}]$. The central particle is placed at the segment center $\bm{p}^{0}=(\bm{a}^{*}+\bm{b}^{*})/2$, the rotation angle $\phi^{0}$ is given by vector $\bm{b}^{*}-\bm{a}^{*}$, and the bending angle is zero $\theta^{0}=0$. Iteration $m+1$ consists of three steps, each correcting one parameter:

We estimate the interface curvature from a discrete volume fraction field defined on a uniform Cartesian mesh. Our method uses a piecewise linear reconstruction of the interface, which is available as part of the PLIC technique for solving the advection equation. Following [7], we compute the interface normals with the mixed Youngs-centered method on a $3\times 3\times 3$ stencil in 3D (or $3\times 3$ in 2D). To define the orientation, we assume that the normals have the direction of anti-gradient of the volume fraction. Then, the interface is reconstructed in each cell independently by a polygon (or line segment) cutting the cell into two parts according to the estimated normal and the given volume fraction [34].

In two dimensions, the interface is represented as a set of line segments. To estimate the curvature in one interfacial cell, we collect the line segments from a $5\times 5$ stencil centered at the target cell. The stencil of this size includes all line segments that can be reached by the particles if the string length is $H_{p}\leq 4$. We apply the algorithm from Section 2.2 to these line segments and compute the curvature.

In three dimensions, the reconstructed interface is a set of planar convex polygons. In this case, we compute the mean curvature as the average over $N_{s}$ cross sections and thus reduce the problem to the two-dimensional case. The algorithm to estimate the mean curvature in a cell consists of the following steps:

The iterative algorithm in Section 2.2 is guaranteed to converge for a sufficiently small relaxation factor as shown in Section A4. Here we demonstrate the convergence on the test case of estimating the curvature of a sphere introduced in Section 3.1. The convergence tolerance is set to $\varepsilon_{p}=10^{-5}$ which is three orders of magnitude smaller than the typical curvature error. Figure 3 shows the number of iterations required to satisfy the convergence criterion (10) depending on the relaxation factor. We choose the relaxation factor of $\eta=0.5$ for which the equilibration takes 10-20 iterations as seen from the convergence history in Figure 3. All further computations are performed with $\eta=0.5$, $\varepsilon_{p}=10^{-5}$ and $m_{\text{max}}=20$.

We compare our method with the generalized height-function method for curvature estimation implemented in Basilisk⁴⁴ 4 We use the version of Basilisk available as of August 21, 2019  [1, 29]. The algorithm attempts to compute the curvature in each interfacial cell with a series of techniques depending on the resolution. In three dimensions they can be outlined as follows

• 1.

  Evaluate the height function along a preferred coordinate plane chosen based on the interface normal. The heights are obtained by summation over columns that span up to 13 cells. If consistent heights are available on the $3\times 3$ stencil, compute the curvature from them with finite differences.

• 2.

  Collect heights from mixed directions in the $3\times 3\times 3$ stencil. If found six consistent heights, compute the curvature from a paraboloid fitted to them.

• 3.

  If the curvature is already defined in one of the neighboring cells in the $3\times 3\times 3$ stencil, copy the curvature from the neighbor.

• 4.

  Collect centroids of interfacial cells in the $3\times 3\times 3$ stencil. If found six centroids, compute the curvature from a paraboloid fitted to them.

• 5.

  If all techniques fail, set the curvature to zero.

We refer to the documentation of Basilisk⁵⁵ 5 Curvature estimation in Basilisk: http://basilisk.fr/src/curvature.h for further details of the algorithm. As clearly seen from the above description, the generalized height-function method, while providing an algorithm applicable at all resolutions, involves four different sources of curvature: heights, parabolic fitting to heights from mixed directions, values of curvature from neighboring cells and parabolic fit to centroids of the piecewise linear interface. Such procedure is complex for implementation and lacks robustness (e.g. depends on the presence of neighboring interfacial cells for parabolic fitting).

As shown in Section 3, the error of our method asymptotically approaches a constant. To achieve second-order convergence, we follow the idea of the generalized height-function method and switch to standard heights at high resolutions. Our combined method therefore consists of two steps

• 1.

  Evaluate the height function along a preferred coordinate plane chosen based on the interface normal. If consistent heights are available on a $3\times 3$ stencil, compute the curvature from them with finite differences.

• 2.

  Otherwise, compute the curvature using the proposed particle method as described in Section 2.3.

Overall, the logic of the algorithm has greatly simplified compared to the generalized height function method. At the same time, as we demonstrate in the following sections, the combined method provides better accuracy at low resolutions and second-order convergence for well-resolved interfaces.

Estimation of the interface curvature in three dimensions depends on three parameters: the number of particles per string $N$, the number of cross sections $N_{s}$ and the string length $H_{p}$. We set the parameters to $N=7$, $N_{s}=3$ and $H_{p}=4$. Then, we vary each parameter independently and examine their influence by estimating the curvature of a sphere at various resolutions following the test case in Section 3.1. The figures present the median error over 100 samples for the center from a uniform distribution over the octant of the cell, i.e. each coordinate is sampled from the range $[0,h/2]$. As seen from Figure 4, the number of particles has a minor influence on the result. Nevertheless, we observe that $N=3$ provides a two times larger error for bubbles at resolutions about one cell per radius. The influence of $N_{s}$ in Figure 5 is also small which is expected for a sphere. However, more complex shapes such as those observed during the bubble coalescence in Section 4.3, require at least $N_{s}=3$ as shown in Figure 7. Figure 6 shows a stronger influence of the string length $H_{p}$. Increasing the value from $H_{p}=2$ to $H_{p}=4$ reduces the error by a factor of ten.

We compare the computational cost of our method to that of the generalized height-function method on the test case with a Taylor-Green vortex introduced in Section 4.1. The formulation is modified to consider 890 spherical bubbles at resolution $R/h=4$ on a mesh of $256^{3}$ cells. The initial positions of bubbles are distributed uniformly in the domain. The computations have been performed in parallel with 512 cores on Piz Daint supercomputer, where each compute node is equipped with a 12-core CPU Intel® Xeon® E5-2690 v3. Table 1 reports the runtime for both methods. About 5% of the cells contain interface fragments and require the curvature estimation. The time spent on the curvature estimation with the present method amounts to 125 ms per step or 31.1% of total time, while GHF takes 33 ms per step or 10.2% of total time. Therefore, our method is about four times more computationally expensive than GHF.

The cost of our method depends on the number of iterations and the number of particles per cell, which is controlled by parameters $N$, $N_{s}$ and $\varepsilon_{p}$. Estimation of curvature for a given volume fraction field consists of two parts: extraction of line segments and iterations. Extraction of the line segments from neighboring cells takes about 40% of the time and scales linearly with $N_{s}$. Iterations for equilibration of particles take about 60% of the time and scale linearly with $N_{s}\,N$.

The volume fraction field represents a single sphere (circle) of radius $R$. We vary the number of cells per radius $R/h$ and consider 100 samples for the center from a uniform distribution over the octant (quadrant) of the cell, i.e. sampling each coordinate from the range $[0,h/2]$. We compute the relative curvature error in $L_{2}$ and $L_{\infty}$ norms

where $I$ is the indices of cells containing the interface (i.e. cells $i$ for which $0<\alpha_{i}<1$) and $\kappa_{\text{ex}}$ is the exact curvature (i.e. $\kappa_{\text{ex}}=2/R$ for sphere and $\kappa_{\text{ex}}=1/R$ for circle).

Figure 10 shows the error in comparison to GHF. At low resolutions, our method is more accurate in terms of the $L_{2}$-error up to eight cells per radius in 3D (and four cells in 2D) and at resolutions below two cells per radius the error is a factor of ten smaller. The method of particles alone as described in Section 2.3 does not converge at high resolutions, and the maximum error saturates at about 10%. The combined method from Section 2.6 switches to height functions at high resolutions and therefore converges with second order. Figures 8-9 show examples of final configuration of particles at resolutions below two cells per radius.

We apply the model described in Section 2.1 for a spherical (circular) droplet in equilibrium. The initial velocity is zero and the volume fraction field represents a single sphere (circle) of radius $R$. We assume that both components have the same density and viscosity such that $\rho_{1}=\rho_{2}=\rho$ and $\mu_{1}=\mu_{2}=\mu$. This leaves one dimensionless parameter, the Laplace number characterizing the ratio between the surface tension, inertial and viscous forces

which we set to $\text{La}=1200$. We solve the problem on a mesh of $128^{3}$ cells (or $128^{2}$ cells in 2D) placing the droplet center in the corner and imposing the symmetry conditions on the adjacent boundaries. The other boundaries are free-slip walls. We vary the number of cells per radius $R/h$ and advance the solution until time $T=\rho(2R)^{2}/\mu$.

In the exact solution, the velocity remains zero and the pressure experiences a jump at the interface given by the Laplace pressure $p_{L}=\sigma\kappa_{\text{ex}}$. The numerical solutions develop spurious currents due to an imbalance between the pressure gradient and surface tension. We compute the magnitude of the spurious velocity in the $x$-direction

the corresponding Weber number

and the pressure jump

where $C$ is the indices of all cells.

Figures 11-12 show the values of $\text{We}_{\text{max}}$ and the relative pressure jump $\Delta p/p_{L}$ at $t=T$ depending on the resolution. The evolution of $\text{We}_{\text{max}}$ for two selected resolutions is shown in Figure 13. We observe that the solutions from GHF convergence with time to a zero spurious flow in most cases. This demonstrates the existence of a volume fraction field for which the method of height functions provides a uniform curvature field. In other cases, such as $R/h=2$ and $26.9$ in 3D, the spurious flow does not converge to zero. Meanwhile, our method of particles more accurately predicts the pressure jump at low resolutions below $R/h=4$ cells in 3D (and 3.36 cells in 2D) and provides lower magnitude of the spurious flow in cases such as $R/h=2$ in 3D shown in Figure 13.

We note, however, that the equilibration of a static droplet is largely of theoretical interest as the equilibrium shape is symmetric and aligned with the mesh directions. Such conditions are incompatible with advection and therefore rarely observed in practically relevant simulations. The spurious flow occurs in a more realistic scenario such as the translating droplet case discussed in Section 3.3.

We extend the previous case by adding a uniform initial velocity field $\bm{u}(\bm{x},0)=\bm{U}=U\bm{l}/|\bm{l}|$ where $\bm{l}=(1,\;0.8,\;0.6)$ in 3D and $\bm{l}=(1,\;0.8)$ in 2D. The additional parameter of the problem is the Weber number

which we set to $\text{We}=0.1$ while keeping $\text{La}=1200$. We solve the problem in a periodic domain on a mesh of $128^{3}$ cells (or $128^{2}$ cells in 2D). We vary the number of cells per radius $R/h$ and advance the solution until time $T=2R/U$.

The magnitude of the spurious flow is computed relative to the initial velocity as the maximum over all cells

and definitions of $\text{We}_{\text{max}}$ and $\Delta p$ follow (18) and (19). Figures 14-15 show $\text{We}_{\text{max}}$ and the relative pressure jump $\Delta p/p_{L}$ depending on the resolution. The quantities are averaged over $t\in[T/2,T]$ and the evolution of $\text{We}_{\text{max}}$ for two selected resolutions is shown in Figure 16. Our method provides lower magnitudes of the spurious flow than GHF and more accurate values of the pressure jump at resolutions below $R/h=2.82$ cells in 3D (and $R/h=2$ cells in 2D), and at higher resolutions the error is comparable to GHF. We note that more accurate estimates of curvature in Section 3.1 are observed at similar resolutions.

The Taylor-Green vortex is a classical benchmark for the capabilities of flow solvers to simulate single-phase turbulent flows [40]. Here we extend the formulation by adding a gaseous phase. The problem is solved in a periodic domain $[0,2\pi]^{3}$ with the initial velocity

and a single bubble of radius $R=0.1$ placed at $(2,2,2)$. Parameters of the problem are the Reynolds number $\text{Re}=\rho_{1}/\mu_{1}$ and the Weber number $\text{We}=2\rho_{1}R/\sigma$. Here we choose $\text{Re}=800$ and $\text{We}=2$. The density and viscosity ratios are set to $\rho_{2}/\rho_{1}=0.01$ and $\mu_{2}/\mu_{1}=0.01$.

Figure 17 shows the trajectory of the bubble at various resolutions in comparison to GHF. Both methods converge to the same solution with the mesh refinement. However, at lower resolutions $R/h\leq 3.06$ our method provides qualitatively more accurate trajectories, while in GHF the trajectory is dominated by the spurious flow. Snapshots of the vorticity magnitude and the bubble shape are shown for both methods at $R/h=3.06$ in Figure 18

To evaluate our method on non-spherical shapes, we consider the Taylor-Green vortex with a droplet initially elongated in the $z$-direction. We use the same initial conditions (22) for the velocity field and define the volume fraction to describe a single elliptical droplet with $R=0.1$

The density and viscosity ratios are set to $\rho_{2}/\rho_{1}=10$ and $\mu_{2}/\mu_{1}=10$. Other parameters of the problem are kept the same: the Reynolds number $\text{Re}=\rho_{1}/\mu_{1}=800$ and the Weber number $\text{We}=2\rho_{1}R/\sigma=2$.

The deformation of the droplet is characterized by the gyration tensor

where $V=\int{\alpha\,dV}$ and $\bm{x}_{c}=\frac{1}{V}\int{\bm{x}\alpha\,dV}$. For an elliptical droplet, the principal components $\lambda_{1}<\lambda_{2}<\lambda_{3}$ of the gyration tensor relate to semi-axes $R_{1}<R_{2}<R_{3}$ of the ellipsoid as $R_{i}=\sqrt{5\lambda_{i}},\;i=1,2,3$. Figure 19 shows the trajectory of the droplet center of mass in the $z$-direction and the smallest semi-axis $R_{1}$ of the ellipsoid of gyration at various resolutions in comparison to GHF. Snapshots of the vorticity magnitude and the bubble shapes are shown in Figures 20-21 for both methods at resolutions $R/h=1.53$ and $4.08$. Both methods converge to the same solution with the mesh refinement. The oscillation of the droplet is indicated by $R_{1}(t)$ which starts at $R_{1}=0.1$ according the initial conditions and reaches the first maximum at $t=1.2$ when the droplet approaches a spherical shape shown in Figure 21. Due to inertia, the droplet further deforms to an oblate shape corresponding to a minimum of $R_{1}$ at $t=1.5$. At the lowest resolution $R/h=1.53$, our method captures the formation of the oblate shape at $t=1.8$ which is delayed compared to the finest resolution. However, with GHF the oblate shape is not captured. At the next resolution $R/h=2.04$, GHF computes the trajectory more accurately than our method.

Coalescence of bubbles and drops is an actively studied phenomenon commonly found in nature and industry [5, 35]. Here we consider coalescence of two spherical bubbles. The problem is solved in a periodic domain $[0,1]^{3}$ with zero initial velocity and two tangent bubbles of radius $R=0.15$ placed along the $x$-axis. The only parameter of the problem is the Ohnesorge number

which we set to $\text{Oh}=0.007$. The density and viscosity ratios are set to $\rho_{2}/\rho_{1}=0.01$ and $\mu_{2}/\mu_{1}=0.01$.

We refer to [35] for a detailed experimental study of bubble coalescence. The process starts with the formation of a neck connecting the bubbles which then propagates along the bubble surface. Figures 22-23 show the isosurfaces of the volume fraction with different resolutions at $t/T=0.08$ and $0.18$, where $T=\sqrt{\rho_{1}R^{3}/\sigma}$ is the capillary time.

We found that shapes from our simulations match well the experimental results (see Figures 22, 23 and 29) when a factor of 1.2 is applied to the values of time reported in [35]. We note that the simulations using the boundary integral method [35] capture the shapes reported in experiments without mentioning any such factor. We performed additional simulations using Gerris [29] (see Figure 31) and again matching shapes obtained by these simulations and those obtained experimentally required the adjustment by the factor of 1.2, in agreement with Basilisk and our own software. This discrepancy may be attributed to the use of an inviscid boundary integral method in [35] and the viscous simulations employed in Gerris, Basilisk and our own software. Furthermore, we match the results of another experimental study [39] without scaling of time as shown in Section A2 and Figures 32-33.

The evolution of the neck radius is presented in Figure 24 in comparison to the experimental data reported in [35]. We observe that the present method agrees well with experimental data while this is not the case for simulations using the generalized height function method, in particular at later times (see Figure 23). The generalized height-function method introduces spurious disturbances of the interface near the coalescence neck. Furthermore, its solution does not converge with mesh refinement. The reason for this is that the initial shape has effectively infinite curvature at the coalescence neck as refining the mesh increases the curvature resolved on the mesh. As seen from Figure 25, increasing the resolution reduces the percentage of cells where height functions are not defined and the curvature is estimated using particles. Nevertheless, more accurate estimation of curvature in such cells benefits the overall accuracy.

The following test case demonstrates the deformation and breakup of an air sheet in shear flow. The domain $[0,2]\times[0,1]\times[0,1]$ is periodic in the $x$- and $z-$directions and no-slip conditions are imposed on the boundaries in the $y$-direction. The initial velocity profile is simple shear $u_{x}=(y-0.5)U/D$ and $u_{y}=u_{z}=0$, where $D=0.05$ is the thickness of the sheet and $U=0.2$ is the velocity difference. The initial volume fraction

is shown in Figure 26. Parameters of the problem are the Reynolds number $\text{Re}=\rho_{1}UD/\mu_{1}$ and the Weber number $\text{We}=\rho_{1}U^{2}D/\sigma$. Here we choose $\text{Re}=80$ and $\text{We}=0.64$. The density and viscosity ratios are set to $\rho_{2}/\rho_{1}=0.01$ and $\mu_{2}/\mu_{1}=0.01$.

Initial deformations of the air sheet develop further in the shear flow. This leads to tearing of the interface starting in four distinct locations. Figure 27 shows the shapes computed with both methods of curvature estimation at various resolutions. Both methods converge to the same solution with mesh refinement. At lower resolutions, GHF shows more tearing of the interface, while with the present method the interface remains stable.