Efficient Implementation of an Accurate Algebraic Scheme for Sharp Interface Advection in Multiphase Flows

Introduction

Two-phase flows, a term commonly used in computational modeling, pertain to issues involving two fluid phases. On the other hand, multiphase flows encompass a broader spectrum, extending to particle-laden flows. Despite analytical studies like Plateau on such flows dating back to the 19th century, the scope of analytical work is often significantly limited, even for relatively straightforward problems [1]. Experimental observations for practical applications pose challenges due to the difficulty of adapting experimental techniques to multiphase flows [2]. Consequently, there arises a necessity for the development of numerical techniques that are not only accurate and cost-effective but also guarantee physical consistency. These methods should adeptly capture the interface, seamlessly integrating it with the momentum conservation equation. Unfortunately, developing such approaches proves to be complex due to the inherent challenges associated with multiphase flows. These difficulties encompass a range of issues, including but not restricted to Bestion: (1) Upholding the conservation of mass, momentum, and kinetic energy, (2) Modeling variations in properties across interfaces, particularly significant shifts in density and viscosity, (3) Handling intricate topologies and the differentiation of scales, (4) Ensuring stability in simulating multiphase flows, and (5) Precisely incorporating surface tension forces while maintaining accuracy [3-5].

Figure 1 presents an overview of various categories within two- phase modeling, emphasizing prominent methods that currently garner significant attention in the community. While two-fluid models prove effective for addressing simple problems, they often fall short when handling realistic scenarios [2]. In-depth analyses of two-phase flow models are extensively covered in the works of Ishii and Hibiki and Prosperetti & Tryggvason [2,6]. For readers intrigued by Marker-and-Cell (MAC) techniques, comprehensive discussions can be found in the works of McKee et al. [7]. Similarly, those interested in Front-Tracking (FT) methods can refer to the comprehensive discussions provided by Tryggvason et al. [8,9]. Additionally, the realm of diffuse- interface approaches encompasses Phase Field (PF), Constrained Interpolation Profile (CIP), and Conservative Level-Set (CLS) methods. In contrast, sharp-interface approaches involve different classes of Interface-Capturing and Interface-Tracking methods. It is important to highlight certain abbreviations in Figure 1, such as SPH, LBM, LS, and CLSVOF, which represent the Smoothed-Particle Hydrodynamics, Lattice Boltzmann Method, Level-Set and Coupled Level-Set and Volume-of-Fluid [10-13], respectively.

Figure 1: Categorization of Numerical Techniques for Handling Two-Phase Flows. The Definitions of Abbreviations can be Found in the Main Text

This study primarily centers on one-fluid models, methodologies for capturing interfaces, and Volume-of-Fluid (VOF) techniques [4]. In the VOF approach, the interface delineation occurs implicitly through the allocation of volume fractions to represent the presence of a specific fluid within computational cells. In this context, advection is accomplished by redistributing the fluid's content among neighboring cells, facilitated by its movement across the faces of the computational mesh. Since the introduction of the original VOF techniques recorded in the literature numerous varied VOF approaches have been conceptualized and created [14]. According to Figure 1, there exist two overarching classifications within this method: geometric methodologies involve a detailed reconstruction of the interface using data related to volume fractions [12,13]. On the other hand, algebraic methodologies avoid explicit reconstruction efforts. Algebraic schemes, renowned for their relative simplicity of implementation, enhanced efficiency, and adaptability to unstructured meshes, operate without the constraints of structured meshes. Nevertheless, their foundation relies on heuristic considerations, rendering them less accurate in comparison to their geometric counterparts [15]. Conversely, geometric VOF schemes necessitate intricate geometric operations, resulting in a more cumbersome implementation process and slower execution [5]. It is worth noting that the pursuit of geometric VOF methods tailored for unstructured meshes constitutes a dynamic and ongoing area of research [16-22]. Illustrated in Figure 1, algebraic VOF methods are categorized into two groups: 1) THINC and 2) compressive.

THINC schemes, as proposed by Xiao et al. encompass a collection of recently developed techniques within algebraic VOF methods [23]. These methods assume a hyperbolic-tangent profile for the volume fraction within the computational cell containing the interface and obtain the fluxes algebraically. As asserted by Xie et al. these approaches exhibit accuracy levels similar to those of geometric schemes but with a lower computational cost [21]. As a result, THINC schemes are increasingly garnering attention in current research. Additionally, unlike compressive schemes, these methods eliminate necessity for artificial compression and are not dependent on the Courant number at the local cell level [24]. The second group is also known as Interface Compression (IC) methods, which are high- resolution schemes. These formulations integrate an artificial compression factor into the advection equation for volume fraction, leading to an unconventional diffusion coefficient with a negative value. This term functions by constricting the volume fraction distribution perpendicular to the fluid interface, thereby mitigating interface dispersal resulting from numerical diffusion. Consequently, it supports the limitations and preservation of the proportion of phases [25]. The compressive impact varies based on the mesh resolution as well as intrinsic IC coefficient, which represents a parameter linked to the artificial compression term. Various investigations have been carried out to examine the grid resolution and the coefficient of the interaction, as outlined in prior studies [15,26-28].

Deshpande et al. identified the occurrence of parasitic currents induced by the implementation of the IC technique in flows primarily governed by the effects of interfacial tension [15,29]. Parasitic flows manifest as resilient abstract vortices in close proximity to the interface, arising autonomously without external influences due to inaccuracies in calculating interface curvature or disparities between surface tension and pressure gradient forces. Typically, these issues manifest in simulations of static bubbles and droplets within fluids characterized by high-density ratios [30]. Hoang et al. investigated how the IC coefficient influences various factors such as highest velocity attained, the thickness of the VOF interface, and the occurrence of parasitic currents [31]. They convincingly demonstrated that the optimal condition for preventing parasitic flows along with numerical diffusion corresponds to an IC coefficient of 1. Furthermore, their research emphasized the importance of cell dimensions in determining the optimal conditions for IC coefficient.

In this investigation, our aim is to utilize a precise algebraic approach for resolving the interface advection equation within the context of two-phase flow. This method is commonly known as the Multidimensional Universal Limiter for Explicit Solution (MULES) algorithm [15]. Our approach involves implementing this method through Python coding, as it is currently accessible within the OpenFOAM software. While this open-source software provides users with convenient utilization of this method, it has yet to be made available as an in-house code. The rest of the article is organized as follows: Firstly, the problem statement and modeling approach are introduced in sections 2.1 to 2.3. Subsequently, section 2.4 provides details on the numerical solution procedure, including the discretization of equations. The verification follows in section 3, where results from various schemes, such as first-order upwind second-order central van Leer flux limiter, and geometric VOF, are compared. Section 4.1 explores the impact of changing the IC coefficient on the method's accuracy [9,32-34]. In section f4.2, the value of surface tension is determined using the results of interface capturing with two methods, MULES and van Leer flux limiter (as verification), and its accuracy is assessed. Sections f3 and f4 include criteria such as Interface Advection Error (IAE), Mass Conservation Error (MCE), Elapsed Time (ET), Dimensionless Pressure Difference (DPD), and Capillary number (Ca). The article concludes with a summary of major findings in section f5.

Physical and Mathematical Modelling

Problem Statement

In Figure 2, a bubble with a diameter D0=0.3 m is positioned at coordinates (xc = 0.5 m, yc=0.75 m) within a square domain measuring 1 m in length (l=1 m). A two-dimensional velocity field, featuring a periodicity of T = 2 s, is applied to induce a two-dimensional (2D) vortex. Within this dynamic velocity field, the droplet undergoes a reciprocating rotational motion throughout a 2-second interval [35]. The volume fraction of the droplet is denoted as α and is defined as 1, while its immediate surroundings exhibit a volume fraction of 0. Also, the interface of the droplet is precisely characterized by a volume fraction of 0.5.

Figure 2: The Computational Domain; Regions Colored in Green and Blue are the Surroundings and Droplet, Respectively.

Governing Equations

In multiphase flow modeling using the VOF method as the one- fluid formulation [9], the volume fraction advection equation is solved according to equation (1), along with the continuity and momentum equations [12]. Typically, this equation is solved for the secondary phase. Since the sum of the volume fractions of all phases is equal to 1, the volume fraction of the primary phase fluid can be obtained by a simple subtraction operation after solving the equation for the volume fraction of the secondary phase, as shown in equation (2).

where α is the volume fraction, and the subscripts p and q refer to the primary (surroundings) and secondary (droplet) phases, respectively. Also, tensorial quantities are indicated by bold symbols. For example, U is the velocity vector field created by 2D vortex.

The accurate calculation of the volume fraction is crucial because it determines the curvature of the interface, surface tension, and pressure gradients. However, in multiphase flow problems involving phases with high density ratios, even small errors in the calculation and distribution of the volume fraction can cause significant and severe changes in the effective properties. Additionally, since the phase interface is confined to only a few computational cells, the simulation's dependency on grid size increases within this range when using the VOF method. Rusche [36] introduced a convective term to equation (1), which originates from using the velocity field derived from the weighted average of the velocities of all phases. This weighted average is defined according to equation (3) under the assumption that the influence of each phase's speed on interface changes is directly related to the volume fraction of that phase [37].

The third term in this equation does not have a specific physical meaning but is used to compress the phase interface during the numerical solution of the volume fraction conservation equation. Notably, this term equals zero at the lower and upper limits of the fluid volume fraction, i.e., 0 and 1, and only operates at the interface. Consequently, it more accurately tracks the interface and reduces numerical errors resulting from numerical diffusion. Now that the problem of the numerical solution of the volume fraction equation has been solved, the surface tension force can be calculated based on the volume fraction gradient. This force arises due to the excess pressure gradient at the interface between phases [39]. Brackbill et al. proposed using the Continuum Surface Force (CSF) method to model it [40]. In this method, the surface tension force per unit volume of the fluid element is calculated using the following equation:

Therefore, there is no need to solve the continuity and momentum equations to obtain velocity field.

Dimensionless Parameters

The objective parameters described in the previous sections can be alternatively expressed as [5,17,43,44]:

Where Ar, Eo, Mo, and are Archimedes, Eotvos, and Morton numbers, respectively. Furthermore, equation (14) is derived using the analytical solution for pressure difference in a two- dimensional domain. Therefore, the value of the DPD should be equal to 1 in the ideal mode, making it a reliable criterion for evaluating different schemes in section f4.2 [41].

Numerical Solution

The schematic of computational grid is shown in Figure 3. The quantities of pressure and volume fraction are stored at the centers of the cells, while the velocity components are stored at their faces. Interpolation is used to obtain the desired quantity whenever necessary.

Figure 3: The Computational Grid; Cells Colored in Orange and Red are the Main and Ghost Cells, Respectively

Now it is time to perform the volume integration on both sides of equation (5), which is as follows [9]:

In equations above, the subscripts P, N, U, UU, and correspond to the current cell, neighbor cell, first upwind cell, second upwind cell, and first downwind cell, respectively. To enhance understanding, a row of these cells is illustrated in Figure 4.

Figure 4: The Arrangement of Important Computing Cells Including Current (P), Neighbor (N), First Upwind (U), Second Upwind (UU), and First Downwind (D) Cells

Figure 5: The Schematic for Calculating the Local Extrema of the Volume Fraction; Current and Neighbor Cells

To adhere to the lower and upper limits of the volume fraction and to minimize numerical error, it is essential to adjust equations (37) and (38) as follows:

Figure 6: The Schematic for Calculating the Inflows and Outflows of Anti-Diffusive Fluxes from each Cell Face, Along with their Cumulative Sum for Each Cell

Figure 7: The Schematic for Calculating the Mean Radius of Curvature at the Interface and the Associated Surface Tension

It is important to note that these calculations are straightforward for the other faces as well. A key aspect of calculating surface tension is its dependence on volume fraction. Smoothing this quantity through filtering may enhance accuracy. This filtering method can also effectively impact critical parameters in this study [5,9]. As shown in Table 1, this process is implemented in multiple steps.

Table 1: The Process of Smoothing the Volume Fraction, Along with its Corresponding Steps and Equations

Verification

In this section, the verification of the MULES algorithm using a spherical droplet in a 2D vortex (Nx = Ny=32 and t=0.005s), as detailed in section f2.1, is examined. The results, including IAE, MCE, and ET, are compared with those obtained from the first-order upwind, second-order central, van Leer flux limiter, and geometric VOF methods. First, it is important to note that if the MULES limiter, λm, in Equation (22) is set to 0, the results of the MULES solver should exactly match those of the first- order upwind solver. Table 2 presents three parameters including interface advection error, mass conservation error, and elapsed time (in seconds) for both the MULES and first-order upwind methods.

Table 2: The Comparison Between the two Interface Capturing Strategies, MULES and First-Order Upwind, Based on Three Criteria: IAE, MCE, and ET

As observed in Table 2, for two different values of IC in the MULES scheme, the results for the specified parameters are consistent with those of the first-order upwind method, confirming the accuracy of the MULES solver. The slight difference in elapsed time (approximately 40 seconds) between the two methods is due to the more complex coding of the MULES solver and its use of different functions, which require more time to achieve the same results as the first-order upwind method. According to Amani [5], the MULES method has a lower computational cost than the geometric VOF method, but the geometric VOF method offers higher accuracy. It is also important to ensure that the MULES method achieves good convergence and stability, and performs more accurately than other methods such as first-order upwind, second-order central, and van Leer flux limiter. To demonstrate this, the results for these methods are provided in Table 3. Additionally, Figure 8a-e displays the volume fraction contours for various schemes at three selected times (0,T / 2,T).

Table 3: The Comparison between the Different Interface Capturing Schemes, Based on Three Criteria: IAE, MCE, and ET

As shown in Table 3, the computational cost of the MULES method is higher than that of the other methods (approximately 40 seconds). This is expected due to the more complex calculations and formulations involved in MULES. While the accuracy of the MULES method is lower than that of the geometric method, it is on par with, or even better than, the other methods. This is confirmed by comparing the percentages of IAE and MCE. For example, the IAE in the first-order upwind method is nearly twice that of the MULES, second-order central, and van Leer flux limiter methods. Consequently, this method is likely to experience significant numerical diffusion, leading to inaccurate interface tracking (see Figure 8c), whereas the other methods achieve more consistent interface tracking. Additionally, the MCE for the first-order upwind, van Leer flux limiter, and geometric VOF methods is zero, ensuring high stability and preventing interface noise due to mass differences. In contrast, the MCE for the second-order central difference method is 10.68%, which is 13 times and 1.89 times higher than the MULES method with two different IC values of 0.1 and 0.5, respectively. This higher error indicates that the second-order central difference method is prone to fluctuations at the phase interface, resulting in lower stability (see Figure 8d). Therefore, both the MULES (see Figure 8a-b) (Multimedia available online) and van Leer flux limiter (see Figure 8e) schemes are stable numerical approaches that accurately track the phase interface.

Figure 8: The volume fraction contours at t=0 (left), t=T/2 (center), and t=T (right) for a) MULES with IC=0.1 (Multimedia available online)

2. MULES with IC=0.5 (Multimedia available online)
3. first-order upwind
4. second-order central
5. van Leer flux limiter methods. The interface is depicted with a black iso-line of α=0.5

Results and Discussion

Effect of Interface Compression (IC)

As mentioned in section 2.4, the IC value is typically fixed between 0 and 2. However, Berberovic et al. and Deshpande et al. expanded the range and examined this constant value up to 4 [15,48]. If the IC value is zero, no compression will occur at the interface. A value between 0 and 1 indicates moderate compression, while values greater than 1 indicate enhanced compression. In this section, IC values between 0 and 2 are considered to evaluate their impact on interface dynamics and mass conservation within the droplet (secondary phase). This investigation focuses on a spherical droplet within a 2D vortex problem, using a grid resolution of Nx = Ny = 32 and a time step of t=0.005s, as detailed in section 2.1. Figure 9a-i illustrates the volume fraction contours derived from the MULES method with varying IC values at the final time, T.

Figure 9: The volume fraction contours at t=T for MULES method with a) IC=0, b) IC=0.25, c) IC=0.5, d) IC=0.75, e) IC=1, f) IC=1.25, g) IC=1.5, h) IC=1.75, and i) IC=2. The interface is depicted with a black iso-line of α=0.5.

Figure 9 demonstrates that as the IC value increases, the interface thickness becomes more significant, leading to more accurate interface tracking. For instance, in Figure 9a, when the IC value is zero, there is no phase within the droplet, and the maximum volume fraction inside the droplet is 0.9. However, as the IC value increases, the maximum volume fraction inside the droplet reaches 1, and the phase inside the droplet is considered in the simulation. As the IC value increases, the compressibility of the interface also increases, leading to the volume fraction of 1 inside the droplet advancing towards the interface. However, if the IC value exceeds 1, non-physical instabilities may form during the simulation, disturbing the mass conservation of the droplet. Therefore, it is recommended to keep the IC value between 0 and 1. To substantiate this claim, the results of IAE and MCE are presented in Table 4 and Figure 10. According to the data in Table 4 and Figure 10, the IAE decreases from 4.8% to 3.95% as the IC value increases from 0 to 2. This decline continues steeply until the IC reaches 1.4, where the IAE reaches 3.8%, after which this error fluctuates around 4%. Therefore, it can be concluded that increasing the IC value initially has a favorable effect on interface tracking, but this effect diminishes once the IC surpasses the critical value of 1.4. Conversely, the MCE rises from 0% to 23.19% as the IC increases from 0 to 2. This steep increase in MCE is due to parasitic currents at the interface and instability in the numerical solution. By comparing both errors in Figure 10, it is evident that IAE is almost independent of IC changes, while IC significantly affects the droplet mass conservation. Additionally, the time results reported in Table 4 indicate that changes in IC have no effect on the elapsed time for each simulation. The ET for all simulations remains consistently around 150 seconds.

Table 4: The Comparison Between the Different Values of IC in MULES Scheme, Based on Three Criteria: IAE, MCE, and ET.

Effect of Smoothing the Volume Fraction

In this section, surface tension is calculated using the relationships introduced in section f2.4. The impact of volume fraction smoothing and noise reduction on this quantity is examined using the filters outlined in Table 1 (C1 to C4). This investigation focuses on a stationary spherical droplet problem (x_c = y_c=0.5 m and u = v = 0 m/s), using a grid resolution of N_x = N_y =32 and a time step of t=0.00125s. Table 5 and Figure 11a-b present various parameters for both the MULES and van Leer flux limiter methods, including interface advection error, dimensionless pressure difference, capillary number, and elapsed time. Additionally, Figure 12a-e displays the volume fraction contours along with parasitic current for various schemes at final time, t = 0.5 s.

Table 5: The Comparison Between the two Interface Capturing Strategies, MULES and Van Leer Flux Limiter, Based on Four Criteria: IAE, DPD, Ca, and ET. (C1 to C4 are introduced in Table 1).
Figure 11: The Results of two Interface Capturing Methods, MULES and The Van Leer Flux Limiter, include: a) DPD and Ca, and b) ET and IAE Versus Number of Filtering.

According to the data in Table 5 and Figure 11a-b, both the MULES and van Leer flux limiter methods perform similarly, indicating that volume fraction smoothing has a consistent effect on both schemes of tracking the interface (see Figure 12a-e) (Multimedia available online). A closer examination of Table 5 and Figure 11a reveals that, as the number of filtering steps increases from 0 to 4, DPD and Ca initially decrease. However, after the first filtering, the values of these two quantities increase. With one filtering step, the interface fluctuations are reduced, leading to a decrease in the maximum velocity of the parasitic flow around the droplet by almost 75%, as shown in Figure 12b (Multimedia available online) and Figure 13. Consequently, the capillary number decreases significantly from 0.02484 to its lowest value of 0.00663, representing a 73.31% reduction. A smaller capillary number, closer to zero, indicates higher stability of the solution. Moreover, with one filtering step, DPD reduces by 0.23, approaching the analytical DPD value of 1. However, with further increases in the number of filtering steps, DPD and CA increase, likely due to numerical errors from excessive smoothing of the volume fraction (see Figure 12c-e) (Multimedia available online). Consequently, DPD tends toward an approximate limit value of 1.016, experiencing increases of 0.0162 and 0.01509 compared to the analytical DPD value in the third and fourth filters, respectively. Additionally, according to Figure 13, the maximum parasitic velocity around the droplet increases by 30.92% percent from 0.0663 m/s value to 0.0868 m/s, causing an upward trend in CE from 0.00663 to 0.00868 (also refer to Figure 12c-e) (Multimedia available online). Furthermore, according to Table 5 and Figure 11b, filtering the volume fraction once increases IAE from 0.7 to 0.9. This is likely because minimizing parasitic currents and approaching the real pressure difference inside and outside the droplet makes tracking the interface more challenging, thus increasing the numerical error for both methods. However, with an increased number of filtering steps, IAE decreases from 0.9 to 0.63. As more filtering is performed, fluctuations in the volume fraction are further reduced, resulting in smoother interfaces that the solver can track with less error. It should be noted that filtering imposes additional computational cost, approximately 10 to 15 seconds for the MULES method and 2 to 14 seconds for the van Leer flux limiter method, on the simulation.

Figure 12: The volume fraction contours at final time, t=0.5 s, for MULES (left) and van Leer flux limiter (right) schemes with varying numbers of filtering: a) zero, b) one, c) two, d) three and e) four. The interface is depicted with a black iso-line of α=0.5. (Multimedia available online)

Figure 13: The Maximum Velocity of Parasitic Flow Versus Number of Filtering.

Conclusion

In the present work, an efficient algebraic scheme known as MULES was implemented for sharp interface advection. This method was rigorously verified and compared with several other schemes, including first-order upwind, second-order central, van Leer flux limiter, and Geometric VOF. Two problems involving a spherical droplet in a 2D vortex and a stationary droplet were investigated. The verified model was then employed to examine the effects of the Interface Compression (IC) coefficient, varying from 0 to 2, representing no compression, moderate compression, and enhanced compression. Several key objective parameters, including Interface Advection Error (IAE), Mass Conservation Error (MCE), and Elapsed Time (ET), were analyzed. As the IC value increases, the thickness of the interface becomes more distinct, enhancing the accuracy of the interface tracking within the simulation. At very low IC values, the droplet lacks a discernible phase, and the volume fraction remains limited. However, as the IC value rises, the simulation effectively captures the internal phase of the droplet, leading to a more realistic representation. Interestingly, while the accuracy of interface tracking improves with increasing IC values, this enhancement comes with a caveat. When the IC value surpasses a certain threshold, the simulations begin to exhibit non-physical the droplet. These instabilities manifest as erratic behaviors in the simulation, particularly in the form of parasitic currents that arise at the interface, compromising the overall integrity of the results. While IAE showed an almost steady improvement as the IC increased, MCE displayed a heightened sensitivity to changes in IC. This disparity highlights the complex relationship between the IC value and mass conservation, underscoring the potential pitfalls of selecting an excessively high IC value. Additionally, MULES and van Leer flux limiter schemes were employed to examine the effects of volume fraction smoothing. Initially, filtering reduces the Dimensionless Pressure Difference (DPD), but beyond the first step, DPD begins to rise again, eventually approaching a value slightly above the analytical prediction. This pattern suggests that while some filtering can stabilize the system, excessive filtering introduces numerical errors that offset the initial benefits. The Capillary Number (Ca) experiences a significant drop after the first filtering step, indicating improved solution stability. However, additional filtering causes the Ca to increase again, reflecting the same trend observed with DPD. The velocity of parasitic flow around the droplet initially decreases markedly with the first filtering step, leading to reduced interface fluctuations and a more stable system. However, with further filtering, the parasitic flow velocity starts to increase again, suggesting that excessive filtering can destabilize the system. The IAE behaves differently compared to other parameters. Initially, IAE increases with one filtering step, likely due to the challenge of accurately tracking a more stable interface. However, as more filtering steps are applied, IAE decreases, indicating that the solver can better track the smoother interfaces produced by extensive filtering.

References

1. Plateau, J. (1873). Experimental and theoretical statics of liquids subject to molecular forces only.
2. Prosperetti,  A.  (2009).  Computational  Methods  forMultiphase Flow. Cambridge university press.
3. Bestion, D. (2014). The difficult challenge of a two-phase CFD modelling for all flow regimes. Nuclear Engineering and Design, 279, 116-125.
4. Mirjalili, S., Jain, S. S., & Dodd, M. (2017). Interface- capturing methods for two-phase flows: An overview and recent developments. Center for Turbulence Research Annual Research Briefs, 2017(117-135), 13.
5. Amani, E. (2022) "Numerical simulation of multiphase flows (Graduate Course)," Tehran Polytechnic University, Fall
6. Ishii, M., & Hibiki, T. (2010). Thermo-fluid dynamics of two-phase flow. Springer Science & Business Media.
7. McKee, S., Tomé, M. F., Ferreira, V. G., Cuminato, J. A., Castelo, A., Sousa, F. S., & Mangiavacchi, N. (2008). The MAC method. Computers & Fluids, 37(8), 907-930.
8. Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al- Rawahi, N., Tauber, W., ... & Jan, Y. J. (2001). A front- tracking method for the computations of multiphase flow. Journal of computational physics, 169(2), 708-759.
9. Tryggvason, G., Scardovelli, R., & Zaleski, S. (2011). Direct numerical simulations of gas–liquid multiphase flows. Cambridge university press.
10. Sussman, M., Smereka, P., & Osher, S. (1994). A level set approach for computing solutions to incompressible two- phase flow. Journal of Computational physics, 114(1), 146- 159.
11. Sussman, M., & Puckett, E. G. (2000). A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. Journal of computational physics, 162(2), 301-337.
12. Sharifi, M., & Amani, E. (2024). Electromagnetic effects on the solidification of a metallic alloy droplet impacting onto a surface. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 700, 134806.
13. Khoshnevis, A., Ahmadpour, A., & Amani, E. (2024). Double emulsion generation in shear-thinning fluids under electric field effects. International Journal of Mechanical Sciences, 281, 109556.
14. Hirt, C. W., & Nichols, B. D. (1981). Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of computational physics, 39(1), 201-225.
15. Deshpande, S. S., Anumolu, L., & Trujillo, M. F. (2012). Evaluating the performance of the two-phase flow solver interFoam. Computational science & discovery, 5(1), 014016.
16. Ahn, H. T., & Shashkov, M. (2007). Multi-material interfacereconstruction on generalized polyhedral meshes. Journal of Computational Physics, 226(2), 2096-2132.
17. Hernández, J., López, J., Gómez, P., Zanzi, C., & Faura, F. (2008). A new volume of fluid method in three dimensions—Part I: Multidimensional advection method with faceâ?ÂÃÂ?matched flux polyhedra. International Journal for Numerical Methods in Fluids, 58(8), 897-921.
18. López, J., Zanzi, C., Gómez, P., Faura, F., & Hernández, J. (2008). A new volume of fluid method in three dimensions— Part II: Piecewise-planar interface reconstruction with cubic-Bézier fit. International journal for numerical methods in fluids, 58(8), 923-944.
19. Ivey, C. B., & Moin, P. (2012). Conservative volume of fluid advection method on unstructured grids in three dimensions. Center for Turbulence Research Annual Research Briefs, 179, 192.
20. Maric, T., Marschall, H., & Bothe, D. (2013). voFoam-A geometrical Volume of Fluid algorithm on arbitrary unstructured meshes with local dynamic adaptive mesh refinement using OpenFOAM. arXiv preprint arXiv:1305.3417.
21. Xie, B., Ii, S., & Xiao, F. (2014). An efficient and accurate algebraic interface capturing method for unstructured grids in 2 and 3 dimensions: The THINC method with quadratic surface representation. International Journal for Numerical Methods in Fluids, 76(12), 1025-1042.
22. Jofre, L., Lehmkuhl, O., Castro, J., & Oliva, A. (2014). A 3-D Volume-of-Fluid advection method based on cell- vertex velocities for unstructured meshes. Computers & Fluids, 94, 14-29.
23. Xiao, F., Honma, Y., & Kono, T. (2005). A simple algebraic interface capturing scheme using hyperbolic tangent function. International journal for numerical methods in fluids, 48(9), 1023-1040.
24. Yokoi, K. (2007). Efficient implementation of THINC scheme: a simple and practical smoothed VOF algorithm. Journal of Computational Physics, 226(2), 1985-2002.
25. Weller, H. G. (2006). Anew approach to VOF–based interface capturing methods for incompressible, compressible and cavitating flow. Technical Report.
26. Piro, D. J., & Maki, K. J. (2013). An adaptive interface compression method for water entry and exit.
27. Deshpande, S. S., Trujillo, M. F., Wu, X., & Chahine, G. (2012). Computational and experimental characterization of a liquid jet plunging into a quiescent pool at shallow inclination. International Journal of Heat and Fluid Flow, 34, 1-14.
28. Kubo, S., Yoshida, H., Ichimura, T., Wijerathne, M. L. L., & Hori, M. (2021). Study on influence of prior recognition of flooding state on evacuation behavior. International Journal of Disaster Risk Reduction, 63, 102437.
29. Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., & Zanetti, G. (1994). Modelling merging and fragmentation in multiphase flows with SURFER. Journal of computational physics, 113(1), 134-147.
30. Gerlach, D., Tomar, G., Biswas, G., & Durst, F. (2006). Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows. International Journal of Heat and Mass Transfer, 49(3-4), 740-754.
31. Hoang, D. A., van Steijn, V., Portela, L. M., Kreutzer, M. T.,& Kleijn, C. R. (2013). Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method. Computers & Fluids, 86, 28-36.
32. Ferziger, J. H., Peric, M., & Street, R. L. (2019).Computational methods for fluid dynamics. springer.
33. Van Leer, B. (1974). Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme. Journal of computational physics, 14(4), 361-370.
34. Youngs, D. L. (1982). Time-Dependent Multi-Material Flow with Large Fluid Distortion.
35. Gunsing, M. (2004). Modelling of bubbly flows using volume of fluid, front tracking and discrete bubble models.
36. Rusche, H. (2002). Computational fluid dynamics of dispersed two-phase flow at high phase fractions. Ph. D. thesis, University of London.
37. Damian, S. M. (2012). Description and utilization of interFoam multiphase solver. International Center for Computational Methods in Engineering, 1-64.
38. Ubbink, O., & Issa, R. (1999). A method for capturing sharp fluid interfaces on arbitrary meshes. Journal of computational physics, 153(1), 26-50.
39. Wang, Z., Yang, J., Koo, B., & Stern, F. (2009). A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves. International Journal of Multiphase Flow, 35(3), 227-246.
40. Brackbill, J. U., Kothe, D. B., & Zemach, C. (1992). A continuum method for modeling surface tension. Journal of computational physics, 100(2), 335-354.
41. Pritchard, P. J., & Mitchell, J. W. (2016). Fox and McDonald's introduction to fluid mechanics. John Wiley & Sons.
42. Brochu, T., & Bridson, R. (2009). Robust topological operations for dynamic explicit surfaces. SIAM Journal on Scientific Computing, 31(4), 2472-2493.
43. Xiao, F., Ii, S., & Chen, C. (2011). Revisit to the THINC scheme: a simple algebraic VOF algorithm. Journal of Computational Physics, 230(19), 7086-7092.
44. Enright, D., Fedkiw, R., Ferziger, J., & Mitchell, I. (2002). A hybrid particle level set method for improved interface capturing. Journal of Computational physics, 183(1), 83- 116.
45. Zalesak, S. T. (1979). Fully multidimensional flux-corrected transport algorithms for fluids. Journal of computational physics, 31(3), 335-362.
46. Boris, J. P., & Book, D. L. (1973). Flux-corrected transport.I. SHASTA, a fluid transport algorithm that works. Journal of computational physics, 11(1), 38-69.
47. Kuzmin, D., Möller, M., & Turek, S. (2003). Multidimensional FEM-FCT schemes for arbitrary time stepping. International journal for numerical methods in fluids, 42(3), 265-295.
48. Berberovi- E., van Hinsberg, N. P., Jakirli-, S., Roisman, I. V., & Tropea, C. (2009). Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 79(3), 036306.