† Corresponding author. E-mail:

Project supported by the National Natural Science Foundation of China (Grant Nos. 10972010, 11028206, 11371069, 11372052, 11402029, and 11472060), the Science and Technology Development Foundation of China Academy of Engineering Physics (CAEP), China (Grant No. 2014B0201030), and the Defense Industrial Technology Development Program of China (Grant No. B1520132012).

Motivated by inconveniences of present hybrid methods, a gradient-augmented hybrid interface capturing method (GAHM) is presented for incompressible two-phase flow. A front tracking method (FTM) is used as the skeleton of the GAHM for low mass loss and resources. Smooth eulerian level set values are calculated from the FTM interface, and are used for a local interface reconstruction. The reconstruction avoids marker particle redistribution and enables an automatic treatment of interfacial topology change. The cubic Hermit interpolation is employed in all steps of the GAHM to capture subgrid structures within a single spacial cell. The performance of the GAHM is carefully evaluated in a benchmark test. Results show significant improvements of mass loss, clear subgrid structures, highly accurate derivatives (normals and curvatures) and low cost. The GAHM is further coupled with an incompressible multiphase flow solver, Super CE/SE, for more complex and practical applications. The updated solver is evaluated through comparison with an early droplet research.

Accurate and efficient simulation of interface evolution plays an important role for a wide range of computational fluid dynamics (CFD) problems, especially for multiphase flow. It is even more crucial for problems involving derivatives such as normals or curvatures. For example, the accurate calculation of the interfacial tension, which is highly sensitive to curvatures, is essential for the simulation of droplets or bubbles. Various numerical methods have been proposed to describe the evolution of the interface. These methods are further divided into eulerian capturing methods (ECMs) (e.g. the level set method^{[1]} (LSM)), which implicitly capture the interface, and lagrangian tracking methods (LTMs) (e.g. the front tracking method^{[2]} (FTM)), which explicitly track the interface.

ECMs capture smooth interface with high accuracy and automatically handle interface topology change such as merging or breaking. Nevertheless, ECMs suffer inconveniences in the following aspects: losses of subgrid structures, inaccurate approximations of normals and curvatures, and large computational stencils required by high order numerical methods. The LTM provides sharp interface with high mass conservation. However, the application to three-dimensional (3D) problems is difficult. Main inconveniences include the inefficient interface reconstruction and complicated calculation of interface topology change.

In order to combine advantages of different methods, popular and effective hybrid methods, such as the hybrid particles level set method^{[3]} (HPLSM), have been proposed. However, they are based on ECMs and inherit common inconveniences of ECMs. Moreover, key procedures of these schemes include dissipative interpolation on fixed grids. This kind of interpolation may smooth out important sharp details.

Another possible strategy for improving interface resolution is to employ gradient information. Gradient-augmented methods (GAM), including the compact scheme^{[4]} and the space time conservation element and solution element (CE/SE) method,^{[5]} show that gradients can significantly improve efficiency, accuracy and compactness. The compactness also provides a simple treatment of adaptive mesh refinement and boundary conditions. GAMs achieve success in a wide range of computational fluid dynamics problems.^{[6–12]} As far as we know, few papers^{[13,14]} report applications of this strategy for capturing interface.

Motivated by inconveniences of hybrid methods, a gradient-augmented hybrid method (GAHM) is presented for capturing interfaces in an incompressible two-phase flow. The method uses a modified FTM as its skeleton and combines the LSM and cubic Hermit interpolation (CHI).^{[14]} Gradients are coupled in all steps of the GAHM to maintain high resolution within a single spatial cell. Smooth eulerian level set values are calculated from the interface. They are used for physical simulation as done in the LSM and also for a local eulerian interface reconstruction. The reconstruction avoids marker particles redistribution and enables an automatic treatment of interfacial topology change.

In the next section, we will briefly present the numerical approach used for simulating an incompressible two-phase flow. Then we describe the GAHM in Section 3. Finally, we demonstrate properties of the GAHM including accuracy, convergence and cost.

Consider two incompressible fluids, fluid 1 with density *ρ*_{1} and viscosity *μ*_{1} and fluid 2 with density *ρ*_{2} and viscosity *μ*_{2}. Let *σ* denote the surface tension coefficient. The dimensionless governing equations for two-phase incompressible flow with interfacial forces are

*u*

*ρ*,

*p*,

*μ*,

*g*

*n*

*κ*denote the velocity, density, pressure, viscosity, gravity, outward unit normal vector, and curvature, respectively;

*δ*(

*x*

*x*

_{Γ}) is a delta function that is zero everywhere except at the interface where

*x*

*x*

_{Γ};

*Re*=

*ρ*

_{1}UL/

*μ*

_{1},

*Fr*=

*U*

^{2}/

*gL*, and

*We*=

*ρ*

_{1}

*U*

^{2}

*L*/

*σ*denote the Reynolds, Froude, and Weber numbers, respectively. To prevent numerical instability, discontinuous variables

*ρ*and

*μ*are smoothed as

*ρ*=

*ρ*

_{2}/

*ρ*

_{1}+ (1 −

*ρ*

_{2}/

*ρ*

_{1})

*H*(

*φ*) and

*μ*=

*μ*

_{2}/

*μ*

_{1}+ (1 −

*μ*

_{2}/

*μ*

_{1}) with

*H*(

*φ*) denoting a heaviside function of

*ϕ*, which is defined by a level set function as follows:

*Ω*

_{1}is the region bounded by the interface

*Γ*and

*Ω*

_{2}is the region outside

*Ω*

_{1}. The total computational domain is the union of

*Ω*

_{1},

*Ω*

_{2}, and

*Γ*.

*d*denotes the shortest distance from point

*x*

*Γ*.

*φ*is governed by the level set equation and captured with the HPLSM

^{[3]}or GAHM. Normal vector

*n*

*κ*are derived as

^{[15]}The interface force is reformulated as a volume force

*κ*(

*φ*)

*φ*

*δ*(

*φ*)/

*We*with

*δ*(

*φ*) = d

*H*(

*φ*)/d

*φ*and

*κ*(

*φ*) = (

*φ*/|

*φ*|).

*H*(

*φ*) is written as

Equation (^{[16]} The main idea is to first predict an intermediate velocity * u** by splitting Eq. (

*u*

*u*

^{n + 1}= 0 and take divergence of Eq. (

*n*

*p*= 0. Then, pressures are submitted into Eq. (

As shown in Fig. ^{[17]} (RKM).

In the GAHM, the lagrangian interface is composed of points and cubic curves. A curve is governed by a CHI equation and stored as a line segment. The two end points of the line segments are called markers. Unlike the FTM, the GAHM also evolves level set value *φ* and its gradient *ψ*^{[18]} Secondly, eulerian level set values are projected from the interface for physical calculations. Finally, the interface is locally reconstructed with eulerian level set values.

In the GAHM, the interface is evolved by solving the lagrangian level set equation

*φ*and its gradient

*ψ*

*x*

*t*=

*u*

^{[18]}and the RKM.

^{[17]}Velocities and their gradients are approximated with the CHI

^{[14]}in eulerian grid cells. Readers can consult Ref. [14] on the CHI.

Eulerian level set values are calculated from the tracked interface. Figure *PQ* in Fig. *PQ* is divided into three line segments *PP*_{0}, *P*_{0}*P*_{1}, and *P*_{1}*Q*. Points *P*_{0} and *P*_{1} are yielded using the CHI in line segment *PQ*. *Q*_{0} and *Q*_{1} are points of tri-section of line segment *PQ*. Then, level set values are calculated from line segments *PP*_{0}, *P*_{0}*P*_{1}, and *P*_{1}*Q*. Take line segment *P*_{0}*P*_{1} for example. For a grid node *M*, we project vector *P*_{0}*M**P*_{0}*P*_{1}, yielding vector *P*_{0}*P*_{M} and point *P*_{M}. Local level set values near *P*_{0}*P*_{1} is determined by

if (*P*_{M} is in line segment *P*_{0}*P*_{1})

*n*

_{ave}is the average of normals at points

*P*

_{0}and

*P*

_{1}. Fuction

*S*(

*x*) is defined as

*M*for example, we draw a circle of radius 1.5 h centered at

*M*, here

*h*is the eulerian spatial step. The circle covers five interfacial points (

*P*

_{1}∼

*P*

_{5}) and 6 line segments (

*P*

_{0}

*P*

_{1}∼

*P*

_{5}

*P*

_{6}). Six level set values are yielded. The final level set value is the one with minimum absolute value (

*φ*(

*M*) =

*φ*(

*M*,

*P*

_{2}

*P*

_{3})) Geometric information (such as

*n*

*κ*) is approximated with eulerian level set values as done in the LSM. An obvious observation is that all derived level set values exactly follow the definition of the level set function. Moreover, level set values are independent of eulerian datum. Hence, the GAHM is capable of performing the adaptive simulations and calculations on unstructured grids.

As the interface deforms, the distances of neighboring markers are non-uniform. This may cause extra resources, where markers are crowded, or loss of the subgrid structures, where markers are insufficient. Numerical errors may also lead to oscillations on the tracked interface. This problem is hard to handle by the FTM, particularly difficult in 3D applications. Du *et al.*^{[2]} presented three methods: grid free tracking (GF), grid-based tracking (GB), and locally grid-based tracking (LGB). The GF reconstructs an interface with information about markers. It is accurate, but expensive and inefficient. The GB reconstructs an interface with crossings of the interface with the fixed grid edges. It is robust but dissipative. The LGB combines GF with GB. It uses the GF for interface propagation, and the GB for local interface reconstruction in bifurcating regions. The concept of the LGB is borrowed and combined with the CHI. We call the modified LGB the gradient-augmented LGB (GALGB).

Figure *M* and *N*). In step 2, a new line segment is generated by connecting boundary markers of the bifurcating region (line segment *PQ*). Points of *n*-section on the new line segment are used as intermediate points for step 3 (points *I*_{0} and *I*_{1}). The *n* is the number of crossing eulerian cells of the new line segment. The *φ* and *ψ* on intermediate points are calculated using the CHI in eulerian cells. In step 3, the position of new markers *x*_{nm} and corresponding variables are determined by

A difficult case^{[19]} is used to examine the capability of capturing subgrid structures. The initial interface is a circle of radius 0.15 centered at (0.5, 0.75) in a unit square with grid size 100 × 100. The prescribed velocity field is defined as

^{[3]}obtained by the LSM on a grid of size 64 × 64 completely disappears. In order to analyze the accuracy, the velocity field is reversed at

*t*= 8 when the interface reaches a maximal deformation. Then the circle is recovered at

*t*= 16. In this test, the reconstruction of the interface is an essential part of the scheme.

Figures *t* = 0, *t* = 4, *t* = 8, and *t* = 16, respectively. Thin spiral filament is well captured on a scale of the grid spacing by the GAHM, while serious mass losses are observed at both the head and tail of the filament for the HPLSM. According to Fig.

To further investigate grid dependence and accuracy, the values of considered grid cell spacing *l* are 1/50, 1/100, 1/200, and 1/400, respectively. Table *t* = 16 and particle numbers at *t* = 0 and 8 respectively. Errors by the GAHM are indicated to be close to second-order convergence of the lost area. Particle number reveals that the GAHM saves over 90% of memories, which are occupied by particles in the HPLSM. The growth rate of *n*_{p} (the number of marker particles at *t* = 8) with respect to *l* in the GAHM is smaller than that in the HPLSM. We can easily prove it by deriving the relationship between *n*_{p} and *l*. The reason is that marker particles in the GAHM describe a curve while particles in the HPLSM describe a local domain. This feature can help save memories and increase efficiency. Note that the lost area for HPLSM with grid cell spacing 1/50 is significantly large. It implies that either the grid resolution or distributed particles are insufficient.

Table *t* = 16. The GAHM guarantees approximately third order accurate level set function and second order accurate derivatives. While the HPLSM is only first order accurate for all information.

A symmetrical and laminar 2D droplet falling case^{[20]} is employed to evaluate the performance of the GAHM for incompressible multiphase flows. The computational domain is a rectangle with 100 × 200 uniform grid cells. The boundary condition is non-slip wall. The three numbers are *Re* = 100, *Fr* = 1, and *We* = 50. Density and viscosity ratios are 1.125 and 50. The interface is captured with the HPLSM or GAHM. An initial stationary circular droplet with a radius of 1 is centered at (4.5, 16.5) inside a fluid. Driven by the gravity, the droplet falls toward the bottom and induces several vortices. During the falling, it deforms and stretches into a dumbbell shape. The handle of the final shape at the time of 105 is much thinner than the size of a grid cell and it is very difficult to capture.

Table

In Fig. ^{[20]} The droplet by our approach remains connected and symmetrical. The comparison implies that our approach is less dissipative and preserves more subgrid structures.

Based on a combination of the accurate FTM, smooth LSM and compact CHI, a gradient-augmented interface capturing method is presented for an incompressible two-phase flow. Gradient is of great help to reduce numerical error and to optimize computational stencil. Smooth and sharp level set contours are projected from the tracked interface. As a result, normals are approximated more easily and accurately. The local eulerian interface reconstruction eliminates marker redistribution and automatically handles interfacial topology changes. Test results demonstrate significant improvements on mass loss, highly accurate derivatives, clear subgrid structures, and low cost.

The GAHM is further coupled with an incompressible two-phase flow, Super CE/SE, for more practical evaluations. The updated solver maintains a sharp dumbbell-shaped interface, which ensures correct interfacial tension and physical characteristics. Results reveal that the GAHM is capable of accurately simulating the incompressible two phase flow with complex interface deformation.

**Reference**

1 | J. Comput. Phys. 79 12 |

2 | J. Comput. Phys. 213 613 |

3 | J. Comput. Phys. 176 205 |

4 | J. Comput. Phys. 24 10 |

5 | J. Comput. Phys. 119 295 |

6 | Comput. Phys. Commun. 182 1589 |

7 | Chin. Phys. Lett. 28 124705 |

8 | J. Comput. Phys. 229 7503 |

9 | Comput. Fluids 56 92 |

10 | Chin. Phys. 21 040202 |

11 | Chin. Phys. 23 020202 |

12 | |

13 | Comput. Math. Appl. 29 15 |

14 | J. Comput. Phys. 229 3802 |

15 | J. Comput. Phys. 114 334 |

16 | J. Comput. Phys. 101 146 |

17 | |

18 | J. Comput. Phys. 161 512 |

19 | J. Comput. Phys. 85 257 |

20 | Int. J. Heat Mass Transfer 49 366 |