A note on the summation relation in phase-field equations

Multiphase (multi-fluid) flows are ubiquitous both in industry and in nature, including bubbles and droplets,¹ contact lines,² and flow in porous media,³ 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),⁶ 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.⁹ The two most popular models in the phase-field method context are the Cahn–Hilliard (CH)⁴ and Allen–Cahn (AC)⁵ 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.¹⁴ Since its emergence, different methods have been proposed for multiphase systems.¹⁵ In this study, a recently developed color-gradient-based phase-field model is employed as the interface-capturing equation¹⁶ 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.

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 $ϕ = ϕ r$ (assume that we solve the order parameter of the red fluid, as such $ϕ b = 1 − ϕ 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,¹⁹ or employing higher-order discretization schemes²⁰ to reduce dispersion errors. As for the AC interface-capturing equation, Ren et al.²¹ 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.²² 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.

In Sec. IV, we will show that by considering the summation relation (3), results are asymmetrical and inaccurate. It is worth mentioning that we reach the AC interface-capturing equation if we solve one equation of (4).

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.¹⁷ A linear interpolation is used to update the relaxation time across the interface, i.e., $τ = ϕ r τ r + ϕ b τ 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,²³ 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.

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 cos ( 2 π x / L 0 )$ at the interface are placed in a domain of size $L 0 × 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 = ρ r g L 0 L 0 / η r = 3000$ and the Atwood number $At = ( ρ r − ρ b ) / ( ρ r + ρ b ) = 1 / 2$⁠, where η is the dynamic viscosity and g is the gravitational acceleration. Other parameters are set to $η r / η b = 1 , Ch = 5 / 256$⁠, and $Pe = 1000$⁠.

Figure 1 shows the interface, where the order parameter is $ϕ = 1 / 2$⁠, at different dimensionless times $t * = t / L 0 / ( g × 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,²⁸ 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 (⁠ $∇ · u$⁠), calculated based on a finite difference scheme,²⁹ 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 = ∫ Ω ρ 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 $ϕ , ϕ r , ϕ b$⁠, and $ϕ r + ϕ 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 $ϕ , ϕ r$⁠, and $ϕ 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., $ϕ r + ϕ b$⁠, is not always equal to unity in numerical modelings in contrast to the theoretical relation (3). Also, the maximum and minimum values of $ϕ r + ϕ b$ reach their steady state at almost the same time of $t * ≈ 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 − 1 ) / N$ approaches unity.

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.

The difference between Eqs. (11) and (12) results in an additional term in the recovered hydrodynamic equations as pointed out by Li et al.³² 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., $∇ · u ≠ 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 $ϕ b = 1 − ϕ 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.

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.

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.

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.

As can be seen in Fig. 9, even for a divergence-free and without singularity velocity field, one-eq scheme results in an unphysical appearance of one fluid in the bulk region of the other fluid [Fig. 9(b)].

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 $ϕ$⁠, respectively, and Fig. 11(c) plots the summation of $ϕ r + ϕ 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.

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.

The authors have no conflicts to disclose.

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).

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

In Fig. 13, the results of a case with $ρ * = 1000$ and $η * = 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.