In this paper, we investigate phase-field interface capturing equations for two-fluid systems to probe their accuracy and computational cost. Two different schemes are considered: In the first scheme, one of the two order parameters is numerically solved based on a phase-field equation, while the other order parameter is determined through the summation relation; the summation of order parameters equals unity. In the second scheme, the two order parameters are both obtained numerically by solving their respective phase-field equations. A phase-field model based on the color-gradient (CG) method is chosen, and available lattice Boltzmann models are employed for solving the interface-capturing equations together with the hydrodynamic equation. It is shown that for the first scheme, which includes the summation relation, numerical results become asymmetrical. Also, in some cases, it results in nonphysical interfaces. In terms of computational resources, this first scheme is about 11% faster with 25% less computational memory usage than the second scheme. It is shown that only for a zero velocity domain do the two schemes lead to equal results. Also, a theoretical analysis is conducted to highlight the differences between the two approaches.

## I. INTRODUCTION

Multiphase (multi-fluid) flows are ubiquitous both in industry and in nature, including bubbles and droplets,^{1} contact lines,^{2} and flow in porous media,^{3} to name a few. Complex phenomena accompanied by different time and space scales make both analytical and experimental investigations of these systems difficult, often impossible. Numerical modeling is, therefore, an indispensable tool in investigating multiphase systems. Over the past years, different approaches have been proposed for capturing or tracking interfaces in multiphase systems, including the phase-field,^{4,5} volume-of-fluid (VOF),^{6} and level-set^{7,8} methods. The phase-field method is a diffuse interface, and order parameters (phase fields or volume fractions) are defined in such a way that they have their bulk values in the bulk fluids and change smoothly across the interfaces.^{9} The two most popular models in the phase-field method context are the Cahn–Hilliard (CH)^{4} and Allen–Cahn (AC)^{5} models. These models were first developed for binary fluids and then were extended to handle ternary^{10,11} as well as phase-changing systems.^{12,13}

In modeling multiphase flows, one common relation is that the summation of the order parameters is equal to unity, hereafter called the summation relation. To the best of our knowledge, the effect of this relation on the results and resources when being employed in a numerical solver has not been explored. In the present work, we aim to show the effects of this analytical relation by conducting several numerical tests.

An alternative to classical methods, the lattice Boltzmann method (LBM) has experienced rapid development and attracted great attention during the past three decades, especially in fluid and thermal systems.^{14} Since its emergence, different methods have been proposed for multiphase systems.^{15} In this study, a recently developed color-gradient-based phase-field model is employed as the interface-capturing equation^{16} together with a robust lattice Boltzmann (LB) model for the hydrodynamic equations.^{10,17} Two different schemes are discussed; the first scheme is that one order parameter is determined based on the phase-field equation and the other order parameter is determined based on the summation relation (hereafter, we call this scheme the “one-eq scheme”). The second scheme is that we solve two phase-field equations to determine the two order parameters (hereafter, we call this scheme the “two-eq scheme”). The rest of this paper is organized as follows: In Sec. III, the LB model used for solving the two schemes (as one-eq and two-eq schemes) together with an LB model for recovering the hydrodynamic properties is described. In Sec. IV, the two schemes are investigated and compared with each other. We first start with a complex benchmark in which the interface-capturing equations together with the hydrodynamic equation are solved. Then, the velocity field is prescribed in a way that mimics a complex flow and a benchmark is conducted based on this velocity field. Then, over several benchmarks, the complexity of the velocity field is reduced to highlight the different behavior of the one-eq and two-eq schemes. Concluding remarks are presented in Sec. V.

## II. PHASE-FIELD EQUATIONS

*et al.*

^{16}Consider a domain Ω comprised of two incompressible and immiscible fluids, named the

*red*(subscript

*r*) and the

*blue*(subscript

*b*) fluid. The continuity equation for each fluid reads as

*k*, respectively,

**is the spatial location, and**

*x**t*is the time. The order parameter $ \varphi k$ of fluid

*k*is defined as

*ρ*is the bulk density of fluid

_{k}*k*. The order parameter $ \varphi k$ is zero (unity) when the cell is empty of (filled with) fluid

*k*and changes across the interfaces. In analytical modeling, we assume that the summation of the order parameters is always equal to unity (summation relation),

^{16}the order parameters can be captured based on the following equations:

*D*is called mobility in the phase-field context.

^{18}

*ξ*is the interface thickness, and

**is the unit vector normal to the interface with the following definition:**

*n*One can solve one of the two equations of (4) and compute the other order parameter based on the summation relation (3). In this case, relation (6) reduces to $ \varphi = \varphi r$ (assume that we solve the order parameter of the *red* fluid, as such $ \varphi b = 1 \u2212 \varphi r$). However, in the numerical simulations, due to overshooting and undershooting of the order parameters, employing relation (3) requires further care. Some remedies have been proposed such as using a lattice stencil with higher discrete velocities like D2Q17,^{19} or employing higher-order discretization schemes^{20} to reduce dispersion errors. As for the AC interface-capturing equation, Ren *et al.*^{21} showed that the numerical dispersion in the AC equation is large, and they implemented a cutoff value approach to restrict the order parameters to its theoretical range for simulating Rayleigh–Taylor instability (RTI). This is expected to violate mass concervation. Hu *et al.*^{22} proposed another approach by implementing a linear combination of the non-local (weight coefficient of zero) and local (weight coefficient of unity). They showed that a weight coefficient of 0.9 minimized the order parameter overshooting or undershooting; however, the underlying physics of this approach is not clear.

## III. LATTICE BOLTZMANN MODEL

All benchmarks in this study are conducted in 2D. A D2Q9 stencil is selected for both the interface-capturing and hydrodynamic equations. Also, for the hydrodynamic equation, a multi-relaxation time (MRT) collision operator is chosen for the RTI benchmark.^{17} A linear interpolation is used to update the relaxation time across the interface, i.e., $ \tau = \varphi r \tau r + \varphi b \tau b$ with *τ* being the relaxation time. For the sake of briefness, we do not present the lattice Boltzmann equations (LBEs) employed in this study for solving the interface-capturing equations as well as the hydrodynamic ones. The models are presented in Refs. 10, 16, and 17 in detail.

Since we are only interested in the differences between the one-eq and the two-eq schemes, we used the same hydrodynamic equation for both schemes. Since there are different types of surface tension force,^{23} and in order to avoid their effects on numerical results, we simulate all benchmarks in this study with zero surface tension. It should also be noted that the interface-capturing equations are coupled to the hydrodynamic equation through the velocity field and the surface tension force is not embedded in the interface-capturing equations.^{24,25} To verify that the employed LB model is sufficiently accurate, we simulated Poiseuille flow for our two schemes and validated the results to the analytical solution, as shown in the Appendix. In the following, we investigate the differences between the one-eq scheme and the two-eq scheme through analytical and numerical analyses.

## IV. RESULTS

*ξ*and the mobility

*D*are, respectively, determined with the Cahn number (Ch) and the Péclet number (Pe) which are defined as

*L*

_{0}and

*U*

_{0}being length and velocity references, respectively. It should be mentioned that all quantities in this study are either dimensionless or in the LB units.

### A. Rayleigh–Taylor instability

The RTI is our first benchmark for comparing the two schemes. RTI is extensively utilized for evaluating multiphase fluid models.^{16,17,26,27} This instability occurs when a heavy fluid (the *red* fluid here) is sitting on top of a lighter fluid (the *blue* fluid here) in a gravitational field. The heavy fluid moves downward pushing the lighter fluid to move upward. Two layers of fluids with the initial perturbation of $ y p = 0.1 L 0 \u2009 cos \u2009 ( 2 \pi x / L 0 )$ at the interface are placed in a domain of size $ L 0 \xd7 4 L 0$. The top and bottom boundaries are considered to be walls, while the two side boundaries are treated as periodic. The surface tension is zero and $ L 0 = 256$. The dimensionless numbers governing the physics are the Reynolds number $ Re = \rho r g L 0 L 0 / \eta r = 3000$ and the Atwood number $ At = ( \rho r \u2212 \rho b ) / ( \rho r + \rho b ) = 1 / 2$, where *η* is the dynamic viscosity and *g* is the gravitational acceleration. Other parameters are set to $ \eta r / \eta b = 1 , \u2009 Ch = 5 / 256$, and $ Pe = 1000$.

Figure 1 shows the interface, where the order parameter is $ \varphi = 1 / 2$, at different dimensionless times $ t * = t / L 0 / ( g \xd7 At )$ of the instability for the two schemes. It should be noted that for the one-eq scheme, we simulate two different cases. In one case, the *red* fluid equation is solved and the order parameter of the *blue* fluid is obtained via the summation relation (3). In another case, the order parameter of the *blue* fluid is solved based on the interface-capturing equation, and the order parameter of the *red* fluid is obtained via the summation relation (3). As can be seen, the results of the two schemes are similar at early times. However, at later times, the one-eq scheme sooner results in an asymmetrical interface with respect to a vertical line in the middle of the domain. It should be noted that due to floating-point effects,^{28} numerical results are asymmetric for both the one-eq and two-eq schemes. However, due to an additional term in the one-eq scheme for the other fluid, Eq. (10), the asymmetry happens faster in the one-eq scheme as Fig. 1 suggests. By changing the interface capturing equation from the *red* fluid to the *blue* fluid, the result changes for the one-eq scheme. Thus, the results depend on the fluid choice for which the interface equation is solved. Interestingly, the one-eq scheme for the *blue* fluid shows better symmetry. Note that since the density field is computed through the linear relation (7) and the densities of the two fluids are different, slight overshooting and undershooting of the order parameter in the one-eq scheme have different effects on the density field when the fluid of choice is changed. As shown in Fig. 4(a), the order parameter (6) used for calculation of the normal vector (5) shows higher oscillation when the *red* fluid is solved for which in return leads to higher asymmetry. On the other hand, in the two-eq scheme, the results are independent of the fluid of choice since both fluids are equally solved. It should be mentioned that, during the simulation, the droplets with radii smaller than three grid units do not shrink over time for both schemes, as apparent in Fig. 1(f).

In Fig. 2, the divergence of the velocity field ( $ \u2207 \xb7 u$), calculated based on a finite difference scheme,^{29} is shown for the one-eq scheme when the *red* fluid is solved for. As can be seen, the velocity field is not divergence-free. This affects the solution of the one-eq scheme. We will conduct simulations with the same magnitude of divergence as this benchmark to show the effects of the hydrodynamic field on the interface-phase capturing equations (Secs. IV B–IV E).

The time history of the system mass, i.e., $ M = \u222b \Omega \rho d V$, is shown in Fig. 3. Note that the red and blue lines overlap in the figure. As can be seen in Fig. 3, the mass of the system is not conserved in the two-eq scheme at early stages, which is in line with the findings of Fig. 2. At later stages, the deviation from the convergence value of unity for the two-sq scheme is about 0.01%. Based on Fig. 2, the divergence of the velocity field is more pronounced around the interfaces. This non-divergence-free condition leads to local mass conservation violation, which, in turn, results in total mass violation of the system. The reason why the one-eq scheme yields perfect conservation despite the non-divergence-free field needs further investigation which is out of the scope of the current study.

To further investigate the differences between these two schemes, the time history of the maximum and minimum values of $ \varphi , \u2009 \varphi r , \u2009 \varphi b$, and $ \varphi r + \varphi b$ are shown in Figs. 4 and 5, respectively. In order to have a clear presentation, we show the difference from the analytical value of 1 in Fig. 4. As can be seen in Figs. 4(a)–4(c), in the case of the two-eq scheme, the *red* and *blue* order parameters have a larger magnitude of oscillation around their equilibrium ones, especially at early times. Since the two order parameters are not strictly linked together through the summation relation (3) in the two-eq scheme, it takes more time for the two order parameters to reach equilibrium. Figures 5(a)–5(c) show that the variations of the minimum of $ \varphi , \u2009 \varphi r$, and $ \varphi b$ are almost the same for both schemes, except for the *blue* fluid at very early stages. Figures 4(d) and 5(d) indicate that the summation of two order parameters, i.e., $ \varphi r + \varphi b$, is not always equal to unity in numerical modelings in contrast to the theoretical relation (3). Also, the maximum and minimum values of $ \varphi r + \varphi b$ reach their steady state at almost the same time of $ t * \u2248 80$. Since the velocity field is not divergence-free, the summation cannot be unity, and imposing the restriction of Eq. (3) results in asymmetry and other problems which will be highlighted in the following benchmarks.

The computational cost of the two schemes for this benchmark is tabulated in Table I. All simulations in this paper are carried out on a PC with Intel® Core™ i7–6700 HQ 2.60 GHz CPU and 16 GB RAM. As can be seen in the table, for this benchmark the one-eq scheme is only 11% faster than the two-eq scheme. Also, the one-eq scheme consumes 25% less computational memory compared to its two-eq counterpart for data storage. Note that the computational time and memory for solving the hydrodynamic equation are included in the comparisons above. By taking the summation relation into account, one interface-capturing equation is not solved, and we expect that by increasing the number of fluids (*N*), these differences in computational time and memory will vanish as the ratio $ ( N \u2212 1 ) / N$ approaches unity.

. | Simulation time (s) . | Memory (MB) . |
---|---|---|

One-eq scheme | 178 298 | 161 |

Two-eq scheme | 198 365 | 202 |

Increase (%) | 11 | 25 |

. | Simulation time (s) . | Memory (MB) . |
---|---|---|

One-eq scheme | 178 298 | 161 |

Two-eq scheme | 198 365 | 202 |

Increase (%) | 11 | 25 |

^{30}Therefore, the velocity field is not perfectly divergence-free. As such, one should pay attention to the convection term in the phase-capturing equations. When solving two equations, we can expect the same divergence in each phase, and it looks like each equation is responsible for tackling its error and reaching its equilibrium. However, by using one equation, a non-symmetry in the results happens as the convective term is not the same for the two fluids. As such, it matters which fluid is solved since the domain is not symmetric from a fluid-occupying point of view. For instance, by considering the summation relation (3) and the one-eq scheme (assume we solve the

*red*fluid), the other fluid (the

*blue*fluid) equation would be

*red*fluid equation, and this unwanted term causes the differences and early asymmetry when we change the fluid of choice in the one-eq scheme (compare the solid red line with the solid blue line in Fig. 1).

In realistic multiphase flow numerical simulations (like the one in Fig. 2), the velocity field is complicated with all possible stretching, shearing, fluid merging, and breakup, and the additional term in Eq. (10) might be problematic. As such, it is necessary to make sure that the one-eq scheme can still give us reliable results when the flow field is complicated. We will investigate this in the following four benchmarks. The velocity is prescribed and it is reversed in the middle of the simulation. Also, the interfaces do not cross the domain boundaries during the whole simulation. As such, based on the following statement, we expect that the initial interfaces match the final interfaces.

*red*fluid):

^{31,32}

*blue*fluid and the two-eq scheme. However, they are not presented here for the sake of briefness.

The difference between Eqs. (11) and (12) results in an additional term in the recovered hydrodynamic equations as pointed out by Li *et al.*^{32} In our case, the right-hand side of Eq. (11) is the additional term. They proposed a compensating term to remove this unwanted term. However, it is not clear how to evaluate the term since it contains both the order parameters and the velocity field, and it seems that a prediction-correction procedure is needed when computing the velocity. We implemented the additional term in our code and simulated the RTI, and the results were almost the same as the ones in Fig. 2. As such, we do not attribute the asymmetry and the differences between the one-eq and two-eq schemes of the results to this additional term.

Considering the statements mentioned above, in this study, we consider four benchmark tests in which we expect that the mass of the domain remains constant and the interfaces return to their initial positions at the end of the simulation period. In the first test, which is devised for the first time (Benchmark I), the velocity field is not divergence-free, i.e., $ \u2207 \xb7 u \u2260 0$ (with the same order of magnitude as observed in Fig. 2). In the second test (Benchmark II), the velocity field is set to be divergence-free, with the added complexity of source and sink terms (in this study, they appear in terms of singularity in the velocity field). In the third test (Benchmark III), the velocity field is divergence-free and there is no source/sink term. Finally, in the fourth benchmark (Benchmark IV), the velocity field is zero. We consider two initial conditions for these four benchmarks: Either the initial interface is filled with the *red* fluid or *blue* fluid. For the one-eq scheme, we only solve the interface-capturing equation for the *red* fluid, and the order parameter of the *blue* fluid is obtained via $ \varphi b = 1 \u2212 \varphi r$. It should be noted that solving for the *red* fluid and swapping the fluids inside and outside the interface, it is analogous to solving for both the *red* and *blue* fluid when keeping the fluid inside the interface *red*. For the two-eq scheme, we compute the order parameters of both fluids by solving two interface-capturing equations. All the benchmarks are tabulated in Table II.

Benchmark . | Domain conditions . |
---|---|

I | $ \u2207 \xb7 u \u2260 0$ |

II | $ \u2207 \xb7 u = 0$ with singularity |

III | $ \u2207 \xb7 u = 0$ without singularity |

IV | $ u = 0$ |

Benchmark . | Domain conditions . |
---|---|

I | $ \u2207 \xb7 u \u2260 0$ |

II | $ \u2207 \xb7 u = 0$ with singularity |

III | $ \u2207 \xb7 u = 0$ without singularity |

IV | $ u = 0$ |

### B. Benchmark I

As can be seen, the one-eq scheme results in severe discrepancies between the initial interface and the final interface. These discrepancies also depend on which fluid fills the interface. In the figure, the red (blue) line represents the interface at the final time when the *red* (*blue*) fluid initially filled the inside of the dashed circle. However, for the two-eq scheme, the initial and final interface match, and changing the filling fluid from the *red* to the *blue* does not change the results.

### C. Benchmark II

*t*= 0 to $ t * = 1$, the mass term acts as a source term, while from $ t * = 1$ to $ t * = 2$, it acts as a sink. The prescribed velocity field in the polar coordinate system is given by

*r*is the distance from the reference point (which is the center of the domain) and

*θ*is the angle from the reference direction (the horizontal direction). Other properties are the same as before. Figure 7 shows the interface profile for the one-eq and two-eq schemes.

As can be seen in Fig. 7, in the case of singularity in the velocity field, the one-eq scheme leads to the wrong result with dependency on the chosen fluid. On the other hand, the two-eq scheme results in accurate results which are independent of the chosen fluid.

### D. Benchmark III

As can be seen in Fig. 8, when the velocity field is divergence-free and there is no source/sink term, the one-eq and two-eq schemes result in almost indistinguishable results.

To further investigate the differences, we set $ Pe = 500$ and $ U 0 = 0.08$. Results for the new set of variables are shown in Fig. 9. Here, we utilize contours instead of lines to show the differences.

### E. Benchmark IV

In the final benchmark, the velocity field is zero throughout the simulation. As such, Eqs. (4) reduce to diffusion equations. A square domain of the length size $ L 0 = 100$ is considered and a circle of radius $ R 0 = 20$ is placed at the center of it. Other properties are the same as Fig. 6. Figure 10 shows the interface for the two schemes. As can be seen, there is no apparent difference between the one-eq and the two-eq schemes when the filling fluid is changed.

Figures 11(a) and 11(b) plot the time history of the maximum and minimum values of $\varphi $, respectively, and Fig. 11(c) plots the summation of $ \varphi r + \varphi b$ for the two schemes when the interface is filled with the *red* fluid. It should be noted that similar results are obtained when the interface is filled with the *blue* fluid. As can be seen, both schemes produce the same results when the velocity is zero, and the convection-diffusion-like Eqs. (4) reduce to diffusion-like equations.

## V. CONCLUSION

A commonly used relation in the context of the phase-field modeling, i.e., the summation of the order parameter is equal to unity, is evaluated in numerical settings. Two different schemes are considered: in the one-equation scheme, one interface-capturing equation is solved and the order parameter of the other fluid is obtained via the summation relation. While in the two-eq scheme, two interface-capturing equations are solved for determining the order parameters of the two fluids. These two schemes were investigated in terms of accuracy and computational costs using the LBM as the numerical solvers. It was shown that using the one-eq scheme produces asymmetrical results for a symmetrical flow, and also the results depend on the fluid of choice. To further investigate the differences, we prescribed the velocity field in a way that mimics different difficulties of a real velocity field, i.e., not divergence-free velocity fields and velocity fields with singularities. Results showed that when the velocity field is not divergence-free or it has singularities, the one-eq scheme leads to inaccurate and fluid-dependent results. The problem is alleviated to some extent when the velocity field is totally divergence-free and there are not any singularities. We observe that the one-eq and two-eq schemes give the same results only if the velocity is zero. Although this paper is based on one specific interface-capturing equation, it is believed that based on the partial differential equation (PDE) of Eq. (10), we would get the same qualitative results for other methods of solving this PDE such as finite-difference or finite-volume methods and also with other phase-field equations such as the CH equation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Reza Haghani:** Conceptualization (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). **Hamidreza Erfani:** Conceptualization (supporting); Writing – review & editing (supporting). **James E. McClure:** Conceptualization (supporting); Writing – review & editing (supporting). **Carl Fredrik Berg:** Conceptualization (supporting); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: VALIDATION

^{17,33}In this study, we employ this benchmark to assess the interface-capturing equations together with the hydrodynamic equations. A 2D channel with periodic boundary conditions in the

*y*-direction and wall boundary conditions in the

*x*-direction is considered. The width of the channel is

*L*

_{0}and the

*red*fluid is filling the channel from the left wall to $ L 0 / 2$ while the

*blue*fluid is filling the channel from $ L 0 / 2$ to

*L*

_{0}. A gravity force in the

*y*-direction is imposed on the mixture which results in a velocity in the

*y*-direction. In the absence of the surface tension force, the Navier–Stokes (NS) equation simplifies to

^{17}we considered density- and viscosity-matched mixtures, i.e., $ \rho * = 1$ and $ \eta * = 1$. Figure 12 shows the results for $ L 0 = 100$,

*ξ*= 4, and $ g = 1 \xd7 10 \u2212 6$. As can be seen, there is almost no difference between the two schemes with the initial interface profile. Also in Fig. 12(a), the results of solving only the hydrodynamic equation are also included (circle symbol) for comparison in which the interface profile is prescribed as the initial condition (A2).

In Fig. 13, the results of a case with $ \rho * = 1000$ and $ \eta * = 100$ are shown with no distinguishable differences. As can be seen, the two schemes are producing almost the same results in this benchmark for different density and viscosity ratios.

## REFERENCES

*Bubbles, Drops, and Particles*

*Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science*

_{2}convection in confined aquifers: Effects of dispersion, permeability anisotropy and geochemistry

*The Lattice Boltzmann Method*