Cascaded lattice Boltzmann method with improved forcing scheme for large-density-ratio multiphase flow at high Reynolds and Weber numbers

A recently developed forcing scheme has allowed the pseudopotential multiphase lattice Boltzmann method to correctly reproduce coexistence curves, while expanding its range to lower surface tensions and arbitrarily high density ratios [Lycett-Brown and Luo, Phys. Rev. E 91, 023305 (2015)]. Here, a third-order Chapman-Enskog analysis is used to extend this result from the single-relaxation-time collision operator, to a multiple-relaxation-time cascaded collision operator, whose additional relaxation rates allow a significant increase in stability. Numerical results confirm that the proposed scheme enables almost independent control of density ratio, surface tension, interface width, viscosity, and the additional relaxation rates of the cascaded collision operator. This allows simulation of large density ratio flows at simultaneously high Reynolds and Weber numbers, which is demonstrated through binary collisions of water droplets in air (with density ratio up to 1000, Reynolds number 6200 and Weber number 440). This model represents a significant improvement in multiphase flow simulation by the pseudopotential lattice Boltzmann method in which real-world parameters are finally achievable.

Figure 1

Simulation results for gas densities in one dimension, for the equation of state given by Eq. (46) with ${\rho }_l^{\text{sat}}=1$. The method of Guo with $A=1/10$ (crosses) and $A=1/40$ (pluses), and the EDM with $A=1/10$ (triangles) and $A=1/40$ (diamonds) are shown compared with the expected density, ${\rho }_g^{\text{sat}}$ (black line).

Figure 2

Simulation results for gas densities in one dimension, for the equation of state given by Eq. (46) with ${\rho }_l^{\text{sat}}=1$. The present method (in the LBGK limit) with increasing interface width (decreasing $A$) using $A=1/10$ (triangles), $A=1/40$ (crosses), and $A=1/160$ (circles), is shown compared with the expected density, ${\rho }_g^{\text{sat}}$ (black line). Significant improvement over the existing methods (shown in Fig. 1) is observed.

Figure 3

Simulation results for gas densities in one dimension, for the equation of state given by Eq. (47), with ${\rho }_l^{\text{sat}}=1,{\theta }_l=1/3$ and ${\theta }_g=1/6$. The method of Guo with ${\theta }_m=-1/256$ (diamonds), the EDM with ${\theta }_m=-1/64$ (triangles) and ${\theta }_m=-1/256$ (circles), and the present method (in the LBGK limit) with ${\theta }_m=-1/64$ (crosses) and ${\theta }_m=-1/256$ (pluses), are shown compared with the expected density, ${\rho }_g^{\text{sat}}$ (black line). Decreasing ${\theta }_m$ increases interface width. The method of Guo was unstable for ${\theta }_m=-1/64$ and for ${\rho }_g^{\text{sat}}<0.01$.

Figure 4

Reduction in error in gas density with interface width, $W$, for the present forcing scheme with the equation of state given by Eq. (47). Gas density is measured far from a flat interface with $\omega =1$ in all cases. Four different density ratios are shown: ${\rho }_r=10$ (diamonds), ${\rho }_r=100$ (circles), ${\rho }_r=1000$ (crosses), and ${\rho }_r=10000$ (triangles). The solid black line shows second-order convergence.

Figure 5

Simulation results for gas densities in one dimension, for the Carnahan-Starling equation of state with $\alpha =1/16$. The results are for the method of Lycett-Brown and Luo [5], which does not include the additional correction derived here for the cascaded collision operator, in the LBGK limit with ${\omega }_b=\omega =1$ (crosses), and for $\omega =1,{\omega }_b=1.6$ (triangles), $\omega =1,{\omega }_b=0.4$ (circles), $\omega =1.8,{\omega }_b=1.2$ (diamonds), and $\omega =1.8,{\omega }_b=0.6$ (squares). They are shown compared with the Maxwell construction (black line). As in Lycett-Brown and Luo [5] the LBGK results are in very good agreement with the Maxwell construction (greatly improving on the previous methods, as shown therein), but as expected there are significant errors in the gas density for the cascaded LBM with ${\omega }_b\ne \omega$. (Liquid densities are all in close agreement with the Maxwell construction and are therefore not shown.)

Figure 6

Simulation results for gas densities in one dimension, for the Carnahan-Starling equation of state with $\alpha =1/16$. The results are for the present method with $\omega =1,{\omega }_b=1.6$ (triangles), $\omega =1,{\omega }_b=0.4$ (circles), $\omega =1.8,{\omega }_b=1.2$ (diamonds), and $\omega =1.8,{\omega }_b=0.6$ (squares). They are shown compared with the Maxwell construction (black line). All results are in very good agreement with the Maxwell construction and with each other.

Figure 7

Measured interface widths, $W$, for varying ${\omega }_b$, for the Carnahan-Starling equation of state with $\alpha =1/32$ and $T=0.0455$ giving a density ratio of approximately 1000. Results are shown for $\kappa =0$ (circles), $\kappa =1$ (crosses), $\kappa =2$ (triangles), and $\kappa =3$ (diamonds). The solid line shows the theoretical prediction for the interface based on Eq. (52) and a measured interface width of 9.46 for the case ${\omega }_b=1$ and $\kappa =0$ (filled circle). All measured interface widths are in very good agreement with the predicted values.

Figure 8

Density profiles across a flat interface for the present method (crosses), the EDM (circles), and the method of Guo (triangles), for the equation of state given by Eq. (46) with ${\rho }_l^{\text{sat}}=1,{\rho }_g^{\text{sat}}=0.2$ and $A=1/10$, compared with the theoretical solution, Eq. (53).

Figure 9

Error in gas density (compared with the predicted gas density from the Maxwell construction) with varying surface tension parameter, ${\sigma }_c$, for the present forcing scheme. The results are for the Carnahan-Starling equation of state, Eq. (45), with $T=0.0455$ giving a density ratio of approximately 1000. ${\sigma }_0$ is ${\sigma }_c$ with $\kappa =0$. Interface width is kept approximately constant at $W=9.4$. Results are shown for the LBGK limit, ${\omega }_b=\omega =1$ (crosses), and for $\omega =1$ and ${\omega }_b=0.6$ (circles), $\omega =1$ and ${\omega }_b=1.6$ (triangles), and $\omega =1.8$ and ${\omega }_b=1.2$ (diamonds).

Figure 10

Error in gas density (compared with ${\rho }_g^{\text{sat}}$) with varying surface tension parameter, ${\sigma }_c$, for the present forcing scheme. The results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}$ giving an expected density ratio of 1000, and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. ${\sigma }_0$ is ${\sigma }_c$ with $\kappa =0$. Interface width is kept approximately constant at $W=10.7$. Results are shown for the LBGK limit, ${\omega }_b=\omega =1$ (crosses), and for $\omega =1$ and ${\omega }_b=0.6$ (circles), $\omega =1$ and ${\omega }_b=1.6$ (triangles), and $\omega =1.8$ and ${\omega }_b=1.2$ (diamonds).

Figure 11

Measured surface tensions, ${\sigma }_m$, against theoretical surface tension ${\sigma }_{th}$ (both normalized by surface tension at $\kappa =0,{\sigma }_0$). Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Interface width is kept approximately constant at $W=11$. Surface tension measurements are made using the Laplace law. Results for $\omega ={\omega }_b=1$ (crosses), $\omega =1$ and ${\omega }_b=1.6$ (circles), and $\omega =1.8$ and ${\omega }_b=1.2$ (diamonds) are almost indistinguishable, and in good agreement with the theoretical prediction.

Figure 12

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with measured surface tension, ${\sigma }_m$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Interface width is kept approximately constant at $W=11$. Results are for $\omega ={\omega }_b=1$ ($\overline{u}$, black crosses; $I$, gray pluses), $\omega =1$ and ${\omega }_b=1.6$ ($\overline{u}$, black diamonds; $I$, gray triangles), and $\omega =1.8$ and ${\omega }_b=1.2$ ($\overline{u}$, black circles; $I$, gray squares).

Figure 13

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with ${\omega }_b$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension and interface width are kept approximately constant at ${\sigma }_m=0.014$ and $W=10.5$, respectively. Results are for $\omega =1$ ($\overline{u}$, black crosses; $I$, gray triangles) and $\omega =1.8$ ($\overline{u}$, black diamonds; $I$, gray circles).

Figure 14

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with ${\omega }_b$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension and interface width are kept approximately constant at ${\sigma }_m=0.00073$ and $W=11.3$, respectively. Results are for $\omega =1$ ($\overline{u}$, black crosses; $I$, gray triangles) and $\omega =1.6$ ($\overline{u}$, black diamonds; $I$, gray circles).

Figure 15

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with ${\omega }_3$ and ${\omega }_4$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension and interface width are kept approximately constant at ${\sigma }_m=0.014$ and $W=10.5$, respectively. Results are for varying ${\omega }_3$ (with ${\omega }_4=1$) ($\overline{u}$, black crosses; $I$, gray triangles) and varying ${\omega }_4$ (with ${\omega }_3=1$) ($\overline{u}$, black diamonds; $I$, gray circles).

Figure 16

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with ${\omega }_3$ and ${\omega }_4$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension and interface width are kept approximately constant at ${\sigma }_m=0.00073$ and $W=11.3$, respectively. Results are for varying ${\omega }_3$ (with ${\omega }_4=1$) ($\overline{u}$, black crosses; $I$, gray triangles) and varying ${\omega }_4$ (with ${\omega }_3=1$) ($\overline{u}$, black diamonds; $I$, gray circles).

Figure 17

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with measured interface width, $W$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension is kept approximately constant at ${\sigma }_m=0.0128$. Results are for $\omega ={\omega }_b=1$ ($\overline{u}$, black crosses; $I$, gray pluses), $\omega =1.0$ and ${\omega }_b=1.6$ ($\overline{u}$, black diamonds; $I$, gray triangles), and $\omega =1.8$ and ${\omega }_b=1.2$ ($\overline{u}$, black circles; $I$, gray squares).

Figure 18

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with measured interface width, $W$. Results are for the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_l=1/3$ and ${\theta }_g=1/6$. Surface tension is kept approximately constant at ${\sigma }_m=0.00074$. Results are for $\omega ={\omega }_b=1$ ($\overline{u}$, black crosses; $I$, gray pluses), $\omega =1.0$ and ${\omega }_b=1.6$ ($\overline{u}$, black diamonds; $I$, gray triangles), and $\omega =1.8$ and ${\omega }_b=1.2$ ($\overline{u}$, black circles; $I$, gray squares).

Figure 19

Variation in measured surface tension, $\sigma$, and interface width, $W$, with sound speed $c_{\mathrm{s}}=\sqrt{\theta }$ using the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_m=-1/512$. Results are for varying ${\theta }_g$ with ${\theta }_l$ fixed at $1/3$ ($\sigma$, black crosses; $W$, gray circles), and varying ${\theta }_l$ with ${\theta }_g$ fixed at $1/768$ ($\sigma$, black diamonds; $W$, gray pluses).

Figure 20

Variation in isotropy, $I$, defined in Eq. (56), and average spurious velocity, $\overline{u}$, with sound speed $c_{\mathrm{s}}=\sqrt{\theta }$ using the linear equation of state, Eq. (47), with ${\rho }_g^{\text{sat}}=0.001$ (density ratio of 1000), and ${\theta }_m=-1/512$. Results are for varying ${\theta }_g$ with ${\theta }_l$ fixed at $1/3$ ($\sigma$, black crosses; $W$, gray circles), and varying ${\theta }_l$ with ${\theta }_g$ fixed at $1/768$ ($\sigma$, black diamonds; $W$, gray pluses).

Figure 21

Head-on coalescence of equal-sized droplets. ${\rho }_r=1000,\mathrm{Re}=1380$, and $\mathrm{We}=27$.

Figure 22

Head-on separation of equal sized droplets, creating one central satellite droplet. ${\rho }_r=1000,\mathrm{Re}=1720$, and $\mathrm{We}=58$.

Figure 23

Experimental results from Pan et al. [43] for the head-on separation of equal-sized droplets, creating one central satellite droplet, at $\mathrm{Re}=1720$ and $\mathrm{We}=58$, as simulated in Fig. 22.

Figure 24

Experimental results from Pan et al. [43] for the head-on splattering of equal-sized droplets at $\mathrm{Re}=6210$ and $\mathrm{We}=440$, as simulated in Fig. 26.

Figure 25

Head-on breakup of equal-sized droplets, creating rim satellite droplets. ${\rho }_r=1000,\mathrm{Re}=4690$, and $\mathrm{We}=280$.

Figure 26

Head-on splattering of equal sized droplets. ${\rho }_r=1000,\mathrm{Re}=6210$, and $\mathrm{We}=440$.