Color-gradient-based phase-field equation for multiphase flow

In this paper, the underlying problem with the color-gradient (CG) method in handling density-contrast fluids is explored. It is shown that the CG method is not fluid invariant. Based on nondimensionalizing the CG method, a phase-field interface-capturing model is proposed which tackles the difficulty of handling density-contrast fluids. The proposed formulation is developed for incompressible, immiscible two-fluid flows without phase-change phenomena, and a solver based on the lattice Boltzmann method is proposed. Coupled with an available robust hydrodynamic solver, a binary fluid flow package that handles fluid flows with high density and viscosity contrasts is presented. The macroscopic and lattice Boltzmann equivalents of the formulation, which make the physical interpretation of it easier, are presented. In contrast to existing color-gradient models where the interface-capturing equations are coupled with the hydrodynamic ones and include the surface tension forces, the proposed formulation is in the same spirit as the other phase-field models such as the Cahn-Hilliard and the Allen-Cahn equations and is solely employed to capture the interface advected due to a flow velocity. As such, similarly to other phase-field models, a so-called mobility parameter comes into play. In contrast, the mobility is not related to the density field but a constant coefficient. This leads to a formulation that avoids individual speed of sound for the different fluids. On the lattice Boltzmann solver side, two separate distribution functions are adopted to solve the formulation, and another one is employed to solve the Navier-Stokes equations, yielding a total of three equations. Two series of numerical tests are conducted to validate the accuracy and stability of the model, where we compare simulated results with available analytical and numerical solutions, and good agreement is observed. In the first set the interfacial evolution equations are assessed, while in the second set the hydrodynamic effects are taken into account.

Figure 1

Diagonal translation of a circular interface at $\text{Pe}=62.5$ and $\text{Ch}=3/100$. (a) $t^{*}=0.25$, (b) $t^{*}=0.50$, (c) $t^{*}=0.75$, (d) $t^{*}=10.0$. The analytical profile is shown by a black dashed line at $t^{*}=10.0$; however, it is covered by the result of the introduced model (shown by the solid blue line) and the LB model [34] (shown by the red solid line).

Figure 2

Convergence rate of the relative error ${\parallel \delta \varphi \parallel }_1$ versus different grid resolutions represented by the grid cell number $L_0$ for the diagonal translation of a circular interface ($\text{Pe}=62.5$ and $\text{Ch}=3/100$) for the present model and the AC model.

Figure 3

Rotation of the Zalesak's disk at $\text{Pe}=500$ and $\text{Ch}=3/200$. (a) $t^{*}=0.5$, (b) $t^{*}=1.0$, (c) $t^{*}=1.5$, (d) $t^{*}=2.0$. The analytical profile is shown by the black dashed line at $t^{*}=2.0$ and is covered by the present result (shown by the solid blue line) and the AC model (shown by the red solid line). Note the red lines in the corners, indicating unphysical behavior for the AC model.

Figure 4

Rotation of the Zalesak's disk with the zero-gradient boundary condition. The parameters are the same as in Fig. 3. (a) $t^{*}=0.5$, (b) $t^{*}=1.0$, (c) $t^{*}=1.5$, (d) $t^{*}=2.0$. The analytical profile is shown by the black dashed line at $t^{*}=2.0$ and is covered by the present result (shown by the solid blue line) and the AC model (shown by the red solid line). Note the red lines at the boundaries, indicating unphysical behavior for the AC model, similar to what was observed in Fig. 3.

Figure 5

A circular interface in a shear flow at $\text{Pe}=500$ and $\text{Ch}=3/100$. (a) $t^{*}=0.5$, (b) $t^{*}=1.0$, (c) $t^{*}=1.5$, (d) $t^{*}=2.0$. The analytical profile is shown by the black dashed line at $t^{*}=2.0$ and is covered by the present result (shown by the solid blue line) and the AC model (shown by the red solid line).

Figure 6

Deformation of a circular interface at $\text{Pe}=1250$ and $\text{Ch}=3/500$. (a) $t^{*}=0.25$, (b) $t^{*}=0.5$, (c) $t^{*}=0.75$, (d) $t^{*}=1.0$. The analytical profile is shown by the black dashed line at $t^{*}=1$ and is covered by the present result (shown by the solid blue line) and the AC model (shown by the red solid line).

Figure 7

Circular droplet in a periodic domain for different density ratios of (a) ${\rho }^{*}=0.5$, (b) ${\rho }^{*}=1.0$, (c) ${\rho }^{*}=2.0$ for the original CG model [Eq. (1)] and (d) ${\rho }^{*}=0.5$, (e) ${\rho }^{*}=1.0$, (f) ${\rho }^{*}=2.0$ for the present model (Eq. (3)). Solid black lines show the interface profiles at the steady-state condition while the red dashed lines represent the initial interface.

Figure 8

Parasitic currents for a static droplet in a quiescent surrounding fluid (a) $\sigma =0.001$ and (b) $\sigma =0$. The velocity vectors are shown in blue color and scaled by $2\times {10}^5$ in lattice units. Solid black lines represent the interface at the steady state while the red dashed lines correspond to the initial interface.

Figure 9

Maximum kinetic energy versus dimensionless time for the case with interfacial tension $\sigma =0.001$. Note that for $\sigma =0$ the model results in the zero velocity field. The reference time is defined as $t_0={\eta }_{r}R/\sigma$.

Figure 10

Phase field profile along the horizontal centerline.

Figure 11

Laplace law verification.

Figure 12

Time evolution of the RTI for ${\rho }^{*}=3$ ($\mathrm{At}=0.5$), ${\eta }^{*}=1$, $\mathrm{Re}=3000$, $\mathrm{Pe}=0.68$, $\mathrm{Ch}=5/256$, and $\mathrm{Ca}=0.26$.

Figure 13

Time evolution of (a) the bubble front position and (b) the liquid front position.

Figure 14

Time evolution of the RTI for ${\rho }^{*}=1000$ ($\mathrm{At}=0.998$), ${\eta }^{*}=100$, $\mathrm{Re}=3000$, $\mathrm{Pe}=81.92$, $\mathrm{Ch}=5/128$, and $\mathrm{Ca}=0.87$.

Figure 15

Time evolution of the droplet splashing on a thin liquid film. The left frames are the results of the present model, the middle frames are the ones by Lee and Lin [54], and the right frames are the results of Shao and Shu [53].

Figure 16

Log-log plot of the spread factor $r/2R$ as a function of dimensionless time $t^{*}$. The straight line corresponds to the power law $r\sqrt{t}$.