A control volume finite element method for three‐dimensional three‐phase flows

A novel control volume finite element method with adaptive anisotropic unstructured meshes is presented for three‐dimensional three‐phase flows with interfacial tension. The numerical framework consists of a mixed control volume and finite element formulation with a new P1DG‐P2 elements (linear discontinuous velocity between elements and quadratic continuous pressure between elements). A “volume of fluid” type method is used for the interface capturing, which is based on compressive control volume advection and second‐order finite element methods. A force‐balanced continuum surface force model is employed for the interfacial tension on unstructured meshes. The interfacial tension coefficient decomposition method is also used to deal with interfacial tension pairings between different phases. Numerical examples of benchmark tests and the dynamics of three‐dimensional three‐phase rising bubble, and droplet impact are presented. The results are compared with the analytical solutions and previously published experimental data, demonstrating the capability of the present method.

In the past, numerous methods have been proposed and used to simulate two-phase flows on a fixed mesh, 4 such as marker-and-cell, 5 volume-of-fluid (VOF), 4,6,7 front-tracking, 8 level set, 9,10 phase field, 11 and moment-of-fluid 12 methods. In addition, particle 13,14 methods have also been developed. An alternative is to consider the use of dynamically adaptive mesh methods, where the mesh resolution can vary in time in response to the evolving solution fields. There are some examples of the use of adaptive mesh refinement for structured meshes with VOF, 15 hybrid-level set/front-tracking, 16 and moment-of-fluid 17 methods. Unstructured meshes are very attractive when dealing with complex geometries in engineering applications and there are examples of adaptive unstructured meshes with the interface tracking method, 18 level set method, 19 phase field method, 20 and algebraic VOF method. 21 Despite the fact that there are a number of numerical studies on two-phase flows, research on three-phase flows (gas-liquid-liquid) is still limited, partly due to the complexity of interactions between phases and numerical treatment of triple junction points. The moment-of-fluid method has been developed in dealing with more than two materials. 22,23 A few numerical examples with interfacial tension can be found in the literature for the level set, 24,25 VOF, 26 moment-of-fluid, 27 phase field, 28,29 smooth particle hydrodynamics (SPH), 30,31 and moving mesh 32 methods. When considering interfacial tension effects, most of the interface capturing methods use the continuum surface force (CSF) approach 33 in both level set and VOF methods. However, far less attention has been paid to level set and VOF methods with fully unstructured meshes for three-dimensional (3D) three-phase flows, even without mesh adaptivity.
In this article, we build upon our novel interface-capturing method 21 and balanced-force interfacial tension algorithm, 34 from two-dimensional (2D) two-phase flow modeling with interfacial tension and adaptive unstructured meshes 35,36 and 2D three-phase flow modeling without interfacial tension 37 to 3D three-phase flow modeling with interfacial tension. The main novelty of the present article is the coupling of interface capturing and interfacial tension on adaptive unstructured meshes for 3D three-phase flow problems, whereas most previous studies focused on two-phase flows. An adaptive unstructured mesh modeling framework is developed here, which can modify and adapt unstructured meshes to better represent the underlying physics of multiphase problems and reduce computational effort without sacrificing accuracy. The numerical framework consists of a mixed control-volume and finite-element formulation, a VOF-type method for the interface-capturing based on compressive control volume advection and second-order finite elements, and a force-balanced algorithm for the interfacial tension that minimizes spurious velocities often found in such simulations. The interfacial tension coefficient decomposition method has been employed to deal with tension pairings between different phases via a compositional approach. Numerical examples of some benchmark tests and the dynamics of 3D three-phase rising bubble and droplet impact are presented to demonstrate the capability of this method. To the best of our knowledge, this is the first 3D control-volume and finite-element method adaptive-mesh simulation of 3D three-phase rising bubble and droplet impact dynamics.
The remainder of this article is organized as follows. Description of the model and numerical methods is given in section 2. 2D and 3D numerical results are presented and compared with analytical solutions or recently published experimental data in section 3. Finally, some concluding remarks and future work are given in section 4.

Governing equations
A three-phase flow modeling framework has been developed based on the multi-component modeling approach with information on interfaces embedded into the continuity equations. In multi-component flows, a number of components exist in one or more phases (one phase is assumed here but the method is easily generalized to an arbitrary number of phases or fluids). Let i be the volume fraction of component i, where i = 1, 2, … , N c , and N c denotes the number of components. The density and dynamic viscosity of component i are i and i , respectively. A constraint on the system is For each fluid component i, the conservation of mass may be defined as and the equations of motion of an incompressible Newtonian fluid may be written as where t is the time, u is velocity vector, p is the pressure, the bulk density is = is the gravitational acceleration vector, and F is the interfacial tension force. In the present study, we focus on interfacial flows with three components, that is, N c = 3.

2.2
Numerical methods

Computational grid
The numerical framework consists of a control-volume and finite-element formulation and also a discontinuous/continuous finite-element pair. The domain is discretized into triangular or tetrahedral elements and, in this work, they are P 1 DG-P 2 elements (linear discontinuous velocity between elements and quadratic continuous pressure between elements). 21 Figure 1 shows the locations of the degrees of freedom for the this element pair and the boundaries of the control volumes in two dimensions.

Temporal discretization
A new time discretization scheme is employed here. When high-order discretization is sought, the method is based on traditional Crank-Nicolson time stepping. This method is often used because it has the simplicity of a two-level time stepping method, is unconditionally stable, and second-order accurate. However, for interface-capturing applications, the time discretization scheme is based on the explicit forward Euler time stepping method. This introduces negative dissipation and is thus a compressive scheme, which helps maintain sharp interfaces. The use of time steps of the order of the grid Courant number and above can result in numerical oscillations and unphysical solutions. For this reason, an adaptive parameter is introduced in Reference 21, in which the forward Euler time stepping method is obtained for = 0, the Crank-Nicolson method is obtained for = 0.5, and the backward Euler method is obtained for = 1. The parameter on each control volume face is chosen using a total variation diminishing criterion. 21

Spatial discretization for the global continuity and momentum equations
In order to discretize the above governing Equations (2) and (3), a finite element representation for u and p is assumed, expressed in terms of their finite element basis functions Q j and P j , respectively, as Here,  u and  p are the total number of degrees of freedom for the velocity and pressure representations. Q j = Q j I, where I is the identity matrix and Q j is the basis function.
Testing the global continuity equation with a quadratic continuous Galerkin basis function P i and applying integration by parts once, the discrete form of the global continuity equation can be obtained as where n is the current time level, n is the outward-pointing unit normal vector to the surface Γ of the domain V, and subscript "bc" means the value on the boundary. Testing Equation (3) with a linear discontinuous Galerkin (DG) basis function Q i and applying integration by parts twice over each element e with the time stepping method, the discrete form of the momentum equation can be obtained as ) dΓ ) dΓ where n is the outward-pointing unit normal vector to the surface Γ E of the element V E , Γ bc is a boundary with prescribed pressure, u nb means the value in the neighboring element along the face Γ E , ∈ {0, 1} is the implicitness parameter, and Δt is the time step size. The upwind DG method is used for the advection terms, where u in is the upwind velocity calculated from the neighboring element or boundary and the subscript ( * ) represents the latest value during the iteration in one time step.
For the viscous terms = (∇u + ∇ T u), we use a high-order linear scheme, which results in a compact stencil with an element coupling only to its surrounding elements. In order to evaluate the viscous stress tensor on the boundary of element Γ E , we integrate over the volume of two neighboring elements in order to calculate the derivatives on the element face between the two elements. For example, the derivative in x for the x component of u is obtained as in which Γ E1∩E2 is the shared face between element 1 (E1) and element 2 (E2), u is the x component of velocity u, and u nb is the value of u in the neighboring element along the face Γ E1∩E2 . It is worth noting that this is not only validated for discontinuous elements but also can be applied for continuous elements.

Projection method
The discretized form of the global mass balance (Equation (5)) and the momentum (Equation (6)) equations is solved using a pressure projection method. This effectively eliminates the unknown velocity and solves a system of equations for pressure or pressure correction. The discretized momentum and continuity equations, at time level n + 1, can be written in matrix form, respectively, as where u n+1 and p n+1 are the FEM solution fields for velocity and pressure, respectively, and s n+1 u and s n+1 p are discretized sources.
For generality, we have introduced a matrix A, containing advection and diffusion contributions, which may be distributed and thus cannot be easily inverted. This allows the method to be applied to inertia-dominated or viscous-dominated flows without modification. Since the velocity is discontinuous between elements, the mass matrix M u = ∫ V E Q i Q j dV is block-diagonal and thus can be easily inverted, each block being local to an element. The solution method proceeds by first solving for an intermediate velocity u n+1 * using a guessed pressure p n+1 * . On the first iteration within a time step, one may use p n+1 * = p n . Equation (9) becomes The matrix equation for velocity to be satisfied is Subtracting these two equations, the velocity correction equation is obtained Multiplying this equation by B T M −1 u (see Equation (10)) and using the global continuity equation to eliminate out u n+1 , one obtains the pressure correction equation This equation is solved for p n+1 and the velocity is corrected using Equation (13). This pressure equation is typically solved using the implementation of the GMRES Krylov subspace solver 38 in the PETSc framework. 39 The Boomer algebraic multigrid preconditioner from the HYPRE library 40 is used to accelerate the process. The velocity solutions are calculated using GMRES 38 with symmetric successive overrelaxation preconditioning.
It is worth mentioning that in the first time step after the mesh adaption, the interpolated velocity field on the new mesh is projected to a continuity-satisfied space, otherwise it will not satisfy the divergence free condition. Thus, the following calculation steps are performed: the interpolated velocity after a mesh adaption is placed into the right-hand side (RHS) of the pressure equation (Equation (14)); this is then used to produce a pressure correction that is placed in the RHS of the velocity correction equation (Equation (13)); and these velocity and pressure corrections are then added to the interpolated velocity and pressures and we commence time stepping from these values. If this modification is not made, then after a mesh adaption, there is typically a spike in pressure and velocity magnitude in isolated regions. This can result in poor accuracy in these regions and this typically leads to the simulations becoming unstable.

Interface capturing method
The interface capturing method is based on a compressive advection method, which uses a mathematically rigorous nonlinear Petrov-Galerkin method that attempts to keep interfaces between components sharp. The mass conservation for each component (Equation (2)) is solved using a control volume formulation, involving a high-order accurate finite element projection to obtain fluxes on the control volume boundaries, where these fluxes are subject to flux-limiting using a normalized variable diagram approach 41 to obtain bounded and compressive solutions for the interface. More details can be found in the work of Pavlidis et al. 21 In this study, the interface capturing method has been further extended to solve three-phase flows in 3D and coupled with interfacial tensions. In contrast to the interface reconstruction in geometrical VOF methods, the interface is represented as the isosurface of the volume fraction at the value of 0.5 in our algebraic VOF interface capturing method here.

Interfacial tension model
The interfacial tension force F in Equation (3) is obtained commonly in two-phase flows using the CSF method 33 as where is the interfacial tension coefficient, = ∇ ⋅ñ is the interfacial curvature,ñ is the interface unit normal, and is the Dirac delta function. However, the above definition has to be further developed for three-phase flows as the fluid flow might be affected by more than two interfaces at the same time near triple junction points. In addition, the interfacial tension coefficient might be different for the interface between different phases. Here, the interfacial tension coefficient decomposition method is employed to deal with tension pairings between different phases using a compositional approach. 24,31,42 The physical interfacial tension coefficients ij between phase i and phase j are decomposed into phase-specific interfacial tension coefficients as ij = i + j , such that for three-phase flows, once the physical interfacial tension coefficients ij are known, the phase-specific interfacial tension coefficients can be obtained as In this approach, the resulting interfacial tension force can be rewritten as sum of each component force: Here, we use i = |∇ i | andñ i = ∇ i |∇ i| to reformulate the CSF based on the component volume fraction: A balanced-force algorithm for the CSF model on unstructured meshes is used here. In contrast to the previous work on two-phase flows, 34 the unit normal and the curvature have to be calculated three times here for each phase. The diffused interface approach based on the original volume fraction i is employed in order to estimate accurately the unit normal and the curvature i on unstructured meshes. More details can be found in Reference 34.

Mesh adaptivity algorithm
As the spatiotemporal dynamics of three-phase flows are highly complex, the dynamically adaptive mesh method has the advantage to resolve small interfacial features during a typical simulation. Mesh adaptivity can also reduce the computational effort, as finer mesh is placed around interfaces during their development, while coarser mesh is used away from these regions, as necessary.
The present model adapts the mesh to the solution without sacrificing the integrity of the boundary (geometry) or internal boundaries (regions) of the domain. It circumvents the complexities of boundary-conforming Delaunay methods by operating on the existing mesh. The error measure employed to adapt the mesh is based on the solution of the phase volume fraction i and provides a directional measure. The objective is to obtain a mesh that has a uniform interpolation error in any direction. This is accomplished with use of a metric, which is related with the Hessian of the solution field. 43 Appropriate scaling of the metric enables the resolution of multiscale phenomena as encountered in multiphase flows. The resulting metric is used to calculate element size and shape. The mesh optimization method is based on a series of mesh connectivity and node position searches of the landscape, defining mesh quality that is gauged by a functional. The mesh modification thus fits the solution field(s) in an optimal manner. The anisotropic mesh adaptivity technique developed by Pain et al. 43 is used here. In this article, the pressure and volume fraction projection uses consistent interpolation, which sets the node value on the new mesh based on its physical position on the old mesh, while the velocity projection uses a conservative Galerkin interpolation technique. 44 The computational cost of the metric construction and the mesh optimization procedure are approximately equal to the time spent on solving the governing equations in one time step. The frequency of mesh adaption is based on our experience and we intend to use the Courant number × frequency < 2 for choosing the frequency of refinement by considering the computational efficiency and accuracy. As we normally adapt the mesh every 10 or 20 time steps, the overhead and extra computational cost introduced by the adaptive procedure would be around 10% of the total CPU time.

Spiral-shearing flow
Before we apply our method to coupled fluid flow problems, we first test our control volume-based interface-capturing method for a pure advection problem, where the fluid interfaces move under a prescribed velocity field. The single vortex problem, 7 which is widely used as a benchmark test for two-phase flows, is extended here for three-phase flows. Two semicircles (radius 0.15) are initially placed on the left and right of the center point at (0.50, 0.75) in a unit square computational domain. The velocity field is defined by the stream function as where u and v are the horizontal and vertical components of the velocity field, respectively, and T corresponds to the total time in the simulation. The initial interface shape is deformed by the velocity field and then returned to its initial state at t = T, where T = 4 is used in the simulation. In order to avoid the influence of time step size on the results, Δt = 2.5 × 10 −4 is used in the computations with minimum mesh size h min = 1∕512 and the Courant number is 0.025. Figure 2 shows the interface shape alongside the mesh during the simulation. It can be seen that the semicircles are stretched from t = 0 under the specified velocity field, reaching its maximum deformation at t = 2. At this stage, they are spiral in shape with elongated filament, which is very thin at the tail. Due to the shearing low, the left semicircle moves to the inner part of the spiral. The fluid interfaces have been efficiently captured in the computation by using the adaptive unstructured mesh, which provides fine resolution equivalent to that of a 512 × 512 uniform mesh ("equivalent" means two meshes have the same minimum grid size in the present article). After t = 2, the velocity field is reversed and the interface shape is returned to its initial shape.
In order to demonstrate the width of the calculated interface, the volume fraction distribution for three phases near the largely deformed interface at t = 2 is shown in Figure 3. In addition, the enlarged areas in the vicinity of the tip are also presented on the bottom of Figure 3 to show the relation between the adaptive mesh and the interface. It can be seen that sharp interfaces are captured by the present compressive VOF method, where the interfaces normally span one or two elements. It is worth noting that it is not easy to present three-phase flow results compared with two-phase flows results; thus, we will only show the region for each phase F I G U R E 2 Volume fraction fields (left) and associated mesh (right) for the single-vortex shearing flow test case (Equation (19)). Three time levels are shown at T = 0 (top), T = 2 (middle), and T = 4 (bottom). The adaptive mesh with minimum length h min = 1∕512 provides a fine resolution equivalent to that of a 512 × 512 uniform mesh [Colour figure can be viewed at wileyonlinelibrary.com] (volume fraction from 0.5 to 1.0) for all the results hereinafter, rather than present three different volume fraction distributions for each phase.
The comparison of the final deformation of the interface against the analytical solution is shown in Figure 4, where the simulation result for the initial circle which is not artificially cut in half as carried out in Reference 34 is also presented here for comparison. It can be seen that the three-phase result is in good agreement with the analytical solution and is similar to the performance of two-phase flows. The normalized deviation of the outer fluid interface position after one rotation can be calculated as where N is the number of interpolation points along the interface and L is a reference length scale, which is the radius here. It is found that the errors for the simulations with minimum mesh size h min = 1∕512 are 0.28% and 0.51% for two-phase and three-phase flows, respectively. Figure 5 shows the time history of the total number of elements during the simulations, N adaptive , which has been normalized by the number of elements on an equivalent 512 × 512 fixed mesh, N adaptive . It can be seen that the ratio increases gradually before reaching its maximum value when the interfaces are at their maximum deformation and then decreases gradually to its initial state. It also shows that the maximum number of total elements for the adaptive mesh simulation is less than 4% of total degrees of freedom for its equivalent fixed mesh. In order to assess the mesh quality, the maximum aspect ratio and maximum angle within the elements during the simulation are shown in Figure 6. It can be seen that the unstructured mesh is anisotropic with large angle and aspect ratio; however, good results can still be obtained. This test demonstrates the power of the adaptive mesh approach, which can refine the mesh in the vicinity of the interface or an area of interest, and reduce computational effort without sacrificing accuracy.

3D reversible vortex
As in Reference 27, we test the reversible vortex in 3D. A sphere with radius 0.15 and center (0.35, 0.35, 0.35), which is artificially cut in half along the vertical axis, is placed in the following flow field: where u, v, and w are the components of the velocity field in the x, y, and z directions, respectively, and t corresponds to the time in the simulation. The initial interfaces are deformed by the velocity field for 0 < t < 3∕2 and then the flow is reversed for 3∕2 < t < 3, and finally the interfaces return to their initial state at t = 3. Δt = 1.0 × 10 −3 is used in the computations with minimum mesh size h min = 1∕100 and maximum mesh size h max = 1∕4 and the Courant number is 0.1. Only fine mesh is used in the vicinity of the interface and coarse mesh is used away from the interface in order to reduce the computational effort. Figure 7 shows the interface shape alongside the mesh during the simulation, in which a reasonable agreement is  obtained for this 3D shear flow field. For comparison, we also performed a two-phase simulation in which the sphere is not artificially cut in half and we can see that similar results are obtained.
As the location of the interface is known (which is represented as the isosurface of the outer fluid volume fraction at = 0.5), the normalized deviation of the interface position after the reversible vortex can be calculated as where N is the number of interpolation points along the interface and L is a reference length scale, which is the radius here. Figure 8 shows the convergence for the computations for two-phase and three-phase cases with three different simulations with minimum mesh size h min = 1∕25, 1∕50, 1∕100. It is shown that the three-phase case has slightly bigger error than the two-phase case due to the additional calculation for the third phase. It can be seen that the present method is close to second-order accurate, which is consistent with the quadratic polynomial function used for the P 2 finite element type. It is worth mentioning that the maximum total number of element for the three-phase case with h min = 1∕100 is 61506, which is one order of magnitude less than the total number of cells for an equivalent fixed mesh (10 6 ) and can reduce the computational effort for the simulations.

Floating lens
In order to check that the present method produces the correct equilibrium contact angles with triple junction points, we follow References 24 and 28 to investigate the floating lens problem, where a droplet of one liquid spreads between two other immiscible liquids, eventually forming a lens shape at steady state (see Figure 9(A)). A circular droplet (fluid 2) is initialized between the interface of fluid 1 and fluid 3 with zero velocity as shown in Figure 9(A). The droplet (radius 0.15) is placed at the center point (0.5, 0.5) in a unit square computational domain, whereas no-slip boundary condition is applied at the walls. The fluids have the same density and viscosity ( 1 = 2 = 3 = 1 and 1 = 2 = 3 = 1), and Re = L * V * = 60 and We 2 = 2 L * V 2 * 2 = 60, where L * = 60 and V * = 1 denote characteristic scales of length and velocity and 2 is the phase-specific interfacial tension for phase 2. We vary the value of the phase-specific interfacial tension coefficient of phases 1 and 3 in order to consider three different cases when increasing the interfacial tension with We 1 = We 3 = 108, 60, 36.
In Figure 9(B-D), the equilibrium droplet shapes are shown for the three cases for We 2 = 60 and We 1 = We 3 = 108, 60, 36, respectively. It can also be seen that the fully unstructured meshes with minimum length h min = 1∕256 and Courant number of 0.25 are adapted in order to capture the evolution of the interface shapes. The equilibrium contact angle is known when specifying the interfacial tension coefficients, and there is an analytical solution for the lens area and its length (the distance between two triple junction points). 24,28 Thus, we compare the length predicted by the numerical simulations and compare it against analytical solutions and a previous study in Table 1; this shows a good accuracy of the present method.
A convergence study with three different fixed meshes and an adaptive mesh has also been carried out for case 2 (We 1 = We 3 = 60). Figure 10 shows the comparison between the analytical solution and numerical simulations for the lens' length TA B L E 1 A comparison of the equilibrium distance between triple junctions obtained from analytical solutions and past (Reference 28) and present numerical simulations  during the interface evolution. It can be seen that the result is converged and all the numerical errors for different meshes are less than 2% when compared with the analytical solution. It is worth noting that the error for the adaptive mesh is slightly larger when compared with the fixed mesh simulations, which might be due to the fully unstructured mesh. We note, however, that the adaptive simulations require fewer degrees of freedom (the number of elements is less than 1300) as the mesh is only refined near the interface. In order to demonstrate the efficiency and the benefits of complexity introduced by adaptive scheme, Figure 11 shows the CPU wall time running in a single core for the simulations with different meshes used to generate Figure 10. As expected, the CPU time increases when refining the mesh. It follows the h −3 slope, which is consistent with the theory of order of 3 as we double the resolution in each direction in 2D and also double the total time step when keeping the constant Courant number. It is worth noting that the CPU time for adaptive mesh is less than all three fixed mesh simulations, being only 0.55% of its CPU time on an equivalent fixed mesh.

Drop levitation
In this section, we consider a droplet leaving an interface under surface tension forces. 24,31 We followed Reference 24 and the same computational setup for floating lens in Section 3.3 is used here, where a circular droplet (fluid 2) is initialized between the interface of fluid 1 and fluid 3 with zero velocity as shown in Figure 9A. Here, we consider when wetting occurs by varying the interfacial tension as 13 (interfacial tension between top and bottom fluids) and 23 (interfacial tension between droplet and top fluid) are equal to unity, whereas 12 (interfacial tension between droplet and bottom fluid) are set to 10. Figure 12 shows the predicted results for the evolution of the 2D droplet due to surface tension forces. The present method is able to deal with large deformation of the interface, which is similar to the results in Reference 24. From t = 0 to t = 0.3, the area for the 2D droplet for fluid 2 has lost only 0.11% during the simulation.
We also extend this case in a fully 3D simulation. Figure 13 shows the predicted evolution of the 3D droplet at the initial state, during leaving and about to leave the interface. At t = 0.25, it can be seen that there is an attachment of the  Figure 14 Variables red color fluid and the blue color fluid at the center and it will eventually leave the bottom fluid during the rising process. From t = 0 to t = 0.25, the volume for the 3D droplet for fluid 2 has lost 0.25% during the simulation.

Air bubble crossing an oil-water interface
After validating the numerical predictions against analytical solutions involving an equilibrium three-phase problem and a clinical test case, in this section, we study 3D three-phase flows of an air bubble crossing a totally stratified oil-water interface. The experimental setup and results have been reported in Reference 26 and a spherical cap bubble case (computational parameters are presented in Table 2) is considered here. The Bond number is Bo = 1 gd 2 ∕ 12 = 29.5 and the Archimedes number is Ar = 1 (gd 3 ) 1/2 ∕ 1 = 42.6 for this case, where d denotes the bubble diameter. For this simulation, we use a relative large computational domain 6d × 12d × 6d to avoid boundary effects; on the other hand, the no-slip condition is applied; the center of bubble is initialized at y = 2d, and the oil-water interface is located at y = 5d in the vertical direction. The minimum length h min = 0.02d is set in the computation and the maximum Courant number is about 0.5 during the whole simulation. Figure 14 shows the evolution of the three-phase flow of an air bubble crossing an oil-water interface. 2D profiles ( Figure 14B) in a central plane associated with the adaptive unstructured mesh ( Figure 14C) and 3D side-view results ( Figure 14D) are presented, along with experimental results in Reference 26 ( Figure 14A) for comparison. It can be seen that the initial spherical cap bubble deforms to a slightly dimpled shape when reaching the oil-water interface. There is a thin film of water above the bubble (shown in the second column of Figure 14B) as it penetrates through the interface. The thin film ruptures (as shown in the third and fourth columns of Figure 14B) and the bubble recovers its spherical cap shape, albeit with a smaller radius in the oil phase. The predicted results are in agreement with the experimental observations at the early stage and the deviation in the late stage (as shown in the fifth and last columns of Figure 14B) is mainly due to the difficulty to simulate the rupture of the thin film. The adaptive meshes on a 2D plane are also shown in Figure 14C, which only refine the resolution in the vicinity of the interfaces and demonstrate the efficiency of the present method. It is worth noting that although the predicted results in Reference 26 are in better agreement with the experimental observations for this case, they studied this problem using a 2D axisymmetric computations with a finer grid, whereas a fully 3D simulation is carried out here.

Water drop impact onto a deep oil bath
In this section, we further investigate 3D three-phase flows by considering the dynamics of a water droplet impacting upon a deep oil bath. The experimental measurements have recently been reported in Reference 45, and we follow the same setup and parameters in that experimental study (shown in   Table 3. A, The experimental results in Reference 45 obtained for the same parameters as in Table 3 The Weber number (We = 1 V 2 d∕ 12 ) based on the water drop is 21. A large computational domain 6d × 9d × 6d, with the air-oil interface at y = 6d, is used to avoid boundary effects, whereas no-slip condition is applied at boundaries, and the minimum length h min = 0.02d is set in the computation. During the simulation, the maximum Courant number is about 0.5. Figure 15 shows the evolution of the three-phase flow of a water drop impacting upon a deep oil bath. 3D view results ( Figure 15B) and 2D profiles ( Figure 15C) in a central plane associated with the adaptive unstructured mesh ( Figure 15D) are presented, along with experimental results in Reference 45 ( Figure 14A) for comparison. It can be seen that the water drop opens up a crater and sets the surrounding oil liquid into motion during impact. The drop becomes flat and spreads below the crater interface. Afterward, the crater collapses upward due to buoyancy effects, and the water drop decelerates and deforms mainly due to the interfacial tension force. There is some discrepancy between the simulation and experimental measurements, which is mainly due to the mesh not being sufficiently fine to resolve triple region between the air, water, and oil.  Figure 16 shows the time history of the droplet volume and its bottom position during the simulation. It can be seen that the droplet volume has lost less than 2.5% during the simulation and the droplet velocity is higher during the initial impact and gradually decreases afterward.
In order to demonstrate the efficiency of the adaptive unstructured mesh simulation, Figure 17 shows the time history of the total number of elements during the simulations, N adaptive , which has been normalized by the number of elements on an equivalent fixed mesh, N fixed . It can be seen that N adaptive of the adaptive mesh gradually increases after the droplet impact in order to resolve the complex interfacial structures. The mesh is finest during the crater-generation process, and the number of elements gradually decreases during droplet spreading and reaches a nearly steady value after the droplet penetrates the air-oil interface. It can be also seen that the maximum number of elements used in the adaptive simulation is less than 0.5% of the fixed mesh case with equivalent mesh resolution; this reduces significantly the computational efforts without sacrificing accuracy.

CONCLUSIONS
In this article, a novel control volume finite element method with adaptive anisotropic unstructured meshes has been presented for 3D three-phase flows with interfacial tension. This method can modify and adapt unstructured meshes to better represent the underlying physics of interfacial problems and reduce computational effort without sacrificing accuracy. The numerical framework consists of a mixed control volume and finite element formulation, a VOF-type method for the interface capturing based on a compressive control volume advection method and second-order finite element methods, and a force-balanced algorithm for the interfacial tension implementation. The interfacial tension coefficient decomposition method has been employed to deal with tension pairings between different phases using a compositional approach. The numerical framework has been validated with several benchmark problems for interface advection, interfacial tension pairings for a floating lens, and three-phase rising bubble and droplet impact problems. The proposed method has the ability to handle a range of complex three-phase flow problems in 3D, involving different phase-phase interactions and triple junctions, in an efficient way. All simulations here were performed in serial using a single processor (Intel Xeon E5-2697v3 2.6 GHz), and typical runtimes are a couple of hours for 2D and approximately 2 weeks for 3D. Future work will include high-performance computing for large scale problems and the inclusion of surfactant effects.