Lattice Boltzmann simulation of mixtures with multicomponent van der Waals equation of state

We developed a general framework for simulating multicomponent and multiphase systems using the lattice Boltzmann framework. Despite the fact that there is no restriction on the number of components in principle, in this article we focus an application to two-component mixtures, but we also demonstrate that the algorithm works for larger numbers of components. To validate our algorithm we separately minimized this underlying free energy to generate theoretical phase diagrams for mixtures of fluids with a van der Waals–like free energy. Our phase diagrams and lattice Boltzmann simulation results are presented in a density-density plane, which best matches with lattice Boltzmann simulations performed at constant volume and temperature. All the theoretical phase diagrams are well recovered by our lattice Boltzmann method.

Figure 1

Simulation results from our baseline phase diagram (see Fig. 6) for the initial $({\rho }^A,{\rho }^B)$ pair $(0.8,0.8)$, which exhibits thermodynamically consistent three-phase equilibrium (see the text for details).

Figure 2

Global phase diagram for a binary van der Waals fluid mixture showing five regions (I–V) reproducible by the VdW EOS with pluses approximating the state of phase diagrams depicted in this paper. The five open symbols with a circular background indicate symmetric mixtures with equal molecular sizes for each component ($\xi =0$) and the five closed symbols indicate asymmetric mixtures for unequal molecular sizes. The dashed curve depicts neutral cross-component interactions ($\nu =1$) for the geometric mixing rule in Eq. (55). Azeotropy is not relevant to the current study and is not depicted here. This figure has been adapted from [39].

Figure 3

Set of points used to initialize the LB simulations to test the two-phase regions of our phase diagrams. The algorithm moves vertically from the $A$-component axis, horizontally from the $B$-component axis, diagonally from axis to axis, and diagonally from (1.0,1.0) to (1.4,1.4). The point (1.5,1.5) is also tested if mixture components are asymmetric enough to admit it.

Figure 4

Phase diagram for a van der Waals mixture of two identical components, generated using parameters of ${\theta }_{cr}^A={\theta }_{cr}^B=0.4, {\rho }_{cr}^B=1.0$, and a neutral interaction parameter $\nu =1.0(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.0)$. Overlaid on top of the theoretical diagram generated by free-energy minimization are the results of the LB simulations. Also depicted is a diagonal connecting the VdW equation discontinuities for both components.

Figure 5

This type II phase diagram is identical to the de facto single-component mixture except phase separation was induced by setting $\nu =0.7(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.3)$. The interaction parameter is not quite repulsive enough to connect all regions of phase separation or provoke three-phase behavior.

Figure 6

Baseline case of three-phase behavior in a type III-H phase diagram. (a) Phase diagram for a mixture identical to the component in Fig. 4 except $\nu =0.5(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.5)$. We also include a comparison of the LB results using the two kinds of chemical potential gradient forcing. (b) Close-up showing the three binodal intersections that define the full three-phase region and the connections between binodals that are used to inscribe the unconditionally unstable three-phase region. The minimization algorithm predicts three-phase behavior for every point in this region, and the behavior was reflected in all LB simulations in this region. (c) Crossing binodals in the vapor region of the three-phase behavior. The LB simulations of metastable points follow the binodals well after the crossing point. The LB simulations for three-phase points show a very small deviation from the point of intersection of the binodals, which is an error of $\sim {10}^{-3}$. (d) Crossing binodals in the $A$-rich liquid region of the three-phase behavior; those in the $B$-rich liquid region are similar. The LB simulations of metastable points follow the binodals well, but they show the same small deviation from the binodal intersection point as noted in (c).

Figure 7

Type III-H phase diagram showing a variation of Fig. 6 where the lattice temperature is above the critical temperatures of each component, but the interaction parameter is repulsive enough to still elicit three-phase behavior. It was generated using parameters of ${\theta }_{cr}^A={\theta }_{cr}^B=0.32, {\rho }_{cr}^B=1.0$, and $\nu =0.2(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.8)$. The maximum density ratio is $\sim 118$.

Figure 8

Type II phase diagram showing a lattice temperature that is in between the critical temperatures of the two components (${\theta }_{cr}^B<\theta <{\theta }_{cr}^A$) and asymmetric critical densities (${\rho }_{cr}^A<{\rho }_{cr}^B$); the interaction parameter is neutral. It was generated using parameters of ${\theta }_{cr}^A=0.5, {\theta }_{cr}^B=0.3, {\rho }_{cr}^B=1.2$, and $\nu =1.0(\xi =0.090909,\zeta =0.162791,\mathrm{\Lambda }=0.013339)$. Six test points defaulted to values for ${\kappa }^{c{c}^{\prime }}$ derived from single-component theory; the values were all between 3.6 and (slightly above) 3.7.

Figure 9

Type III-H phase diagram showing a lattice temperature that is close enough to the critical temperatures of each component that there are three distinct regions of phase separation; the interaction parameter is repulsive. It was generated using parameters of ${\theta }_{cr}^A=0.35, {\theta }_{cr}^B=0.36, {\rho }_{cr}^B=1.2$, and $\nu =0.6(\xi =-0.090909,\zeta =0.104859,\mathrm{\Lambda }=0.403308)$.

Figure 10

Phase diagram showing a mixture with a lattice temperature that is above the asymmetric critical temperatures of the individual components (${\theta }_{cr}^A<{\theta }_{cr}^B<\theta$); however, the interaction parameter is attractive enough to still produce two-phase behavior. It was generated using parameters of ${\theta }_{cr}^A=0.3, {\theta }_{cr}^B=0.31, {\rho }_{cr}^B=1.3$, and $\nu =1.5(\xi =-0.130435,\zeta =0.146515,\mathrm{\Lambda }=-0.483813)$.

Figure 11

Mixture with a lattice temperature that is above the critical temperature of one component but below that of the other (${\theta }_{cr}^A<\theta <{\theta }_{cr}^B$) and with repulsive interactions between the species. It was generated using parameters of ${\theta }_{cr}^A=0.32, {\theta }_{cr}^B=0.45, {\rho }_{cr}^B=1.1$, and $\nu =0.5(\xi =-0.047619,\zeta =0.214724,\mathrm{\Lambda }=-0.511663)$.

Figure 12

Phase diagram with asymmetries that have produced an auxiliary binodal for the metastable region of the phase diagram. A zoomed-in view of the auxiliary binodals is shown in (b). This phase diagram was generated using parameters of ${\theta }_{cr}^A=0.4, {\theta }_{cr}^B=0.45, {\rho }_{cr}^B=1.1$, and $\nu =0.5$ ($\xi =-0.047619, \zeta =0.106145, \mathrm{\Lambda }=0.502825$). Unlike the rest of the phase diagrams shown, the metastable and three-phase LB simulations were done with ${\kappa }^{c{c}^{\prime }}=0.2$. One point (0.0,1.0) used the single-component value of ${\kappa }^{c{c}^{\prime }}=2$.

Figure 13

Phase diagram from the shield region showing three independent three-phase regions shortly after the accessory plaits from each region connect. The single-phase region in the middle is fully enclosed and contains binodal segments that are very nearly continuous. The intermediary two-phase regions contain points where metastable behavior would not be anticipated. This phase diagram was generated using parameters of ${\theta }_{cr}^A={\theta }_{cr}^B=0.427, {\rho }_{cr}^B=1.0$, and $\nu =0.5636(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.4364)$. The metastable and three-phase LB simulations were done with ${\kappa }^{c{c}^{\prime }}=0.2$, and we only show the metastable results associated with one three-phase region to more clearly depict the binodal segments elsewhere.

Figure 14

Phase diagram from the shield region progressing towards four-phase behavior. (a) Regions of unconditionally unstable three-phase behavior, now sharing the same three equilibrium phases. The intermediary two-phase regions are now metastable. (b) Binodals in the middle single-phase region separated into three segments. The parameters are ${\theta }_{cr}^A={\theta }_{cr}^B=0.45, {\rho }_{cr}^B=1.0$, and $\nu =0.5636(\xi =0.0,\zeta =0.0,\mathrm{\Lambda }=0.4364)$. Metastable and three-phase LB simulations used ${\kappa }^{c{c}^{\prime }}=0.2$. (c) Highlight of the lower left of the three-phase density of (a). The metastable simulations lose track of their respective binodals and curl back towards the binodal crossing. (d) Elucidation of this behavior by showing a series of LB simulations that follow the curve of metastable results. Three yellow (white) circles in (c) correspond to lattice site 100 shown in (d) (see the text for details).

Figure 15

Simulation results for three components with the initial $A,B,C$ component set $(1.0,1.0,1.0)$, which exhibits thermodynamically consistent four-phase equilibrium. We use parameters ${\theta }_{cr}^A={\theta }_{cr}^B={\theta }_{cr}^C=0.5, {\rho }_{cr}^B={\rho }_{cr}^C=1.0, \nu =0.05$, and ${\kappa }^{c{c}^{\prime }{c}^{\prime \prime }}=2.0$.