Mesoscopic Lattice Boltzmann Modeling of the Liquid-Vapor Phase Transition

Liquid-vapor phase transition is a widespread phenomenon of great importance in many natural and engineering systems. Because of its multiscale nature and macroscopic complexity Cheng et al. (2014); Onuki (2005); Boreyko and Chen (2009); Cira et al. (2015); Nikolayev and Beysens (1999); Nikolayev et al. (2006), thermodynamically consistent modeling of liquid-vapor phase transition with the underlying physics is a long-standing challenge, despite extensive studies. Physically speaking, the phase transition is a natural consequence of the molecular interaction at the microscopic level. Therefore, as a mesoscopic technique that can incorporate the underlying microscopic interaction, the lattice Boltzmann (LB) method is advocated as a promising method for modeling multiphase flows with phase transition He and Doolen (2002); He et al. (1998); Luo (1998).

The theory of the LB method for multiphase flows has been extensively studied since the early 1990s Gunstensen et al. (1991); Shan and Chen (1993); Swift et al. (1995); He et al. (1998); Luo (1998). However, most studies are inherently limited to isothermal systems, and the theory of the LB method for liquid-vapor phase transition remains largely unexplored. Recently, some liquid-vapor phase transition problems have been simulated by the LB method Hazi and Markus (2009); Gong and Cheng (2012); Li et al. (2015); Albernaz et al. (2017); Qin et al. (2019), where the popular pseudopotential LB model for isothermal systems is adopted to handle the mass and momentum conservation laws, and a supplementary macroscopic governing equation is employed to handle the energy conservation law. Because of the idea of mimicking the microscopic interaction responsible for multiphase flows, the pseudopotential LB model shows great simplicity in both concept and computation. However, it suffers from thermodynamic inconsistency He and Doolen (2002), although the coexistence densities could be numerically tuned close to the thermodynamic results. The supplementary macroscopic energy governing equation is extremely complicated Onuki (2005); Laurila et al. (2012) and it is artificially simplified with macroscopic assumptions and approximations in previous works Hazi and Markus (2009); Gong and Cheng (2012); Li et al. (2015); Albernaz et al. (2017); Qin et al. (2019). Both thermodynamic consistency and the underlying physics are sacrificed. Moreover, the simplified energy governing equation cannot be recovered from the mesoscopic level, implying that the computational simplicity is also lost.

In this Letter, we first analyze the kinetic model that combines Enskog theory for dense gases with mean-field theory for long-range molecular interaction. Guided by this kinetic model, we develop a novel mesoscopic LB model for liquid-vapor phase transition by handling the underlying microscopic molecular interaction rather than resorting to any macroscopic assumptions or approximations. The present LB model has a clear physical picture at the microscopic level and thus the conceptual and computational simplicity, and it is also kinetically and thermodynamically consistent.

The microscopic molecular interaction responsible for liquid-vapor phase transition generally consists of a short-range repulsive core and a long-range attractive tail. The short-range molecular interaction can be well modeled by Enskog theory for dense gases, and the long-range molecular interaction can be described by mean-field theory and thus modeled as a point force Rowlinson and Widom (1982). Combining Enskog theory for dense gases with mean-field theory for long-range molecular interaction, the kinetic equation for the density distribution function (DF) $f(\mathbf{x},{\bm{\xiup}},t)$ can be written as He and Doolen (2002)

where ${\bm{\xiup}}$ is the molecular velocity, $\mathbf{a}$ is the external acceleration, and $V_{\text{mean}}$ denotes the mean-field approximation of the long-range molecular potential. The Enskog collision operator $\mathit{\Omega}_{\text{Enskog}}$ is Chapman and Cowling (1970)

where $\mathit{\Omega}_{0}$ is the usual collision operator for rarefied gases, $\chi$ is the collision probability, $b=2\pi d^{3}\big{/}(3m)$ with $d$ and $m$ the molecular diameter and mass, $\mathbf{C}={\bm{\xiup}}-\mathbf{u}$, and ${\bm{\mathcal{C}}}=\mathbf{C}\big{/}\!\sqrt{2RT}$. The equilibrium DF $f^{\text{eq}}$ is

The density $\rho$ and momentum $\rho\mathbf{u}$ are calculated as

Based on the density DF, a distinct internal kinetic energy $\rho\epsilon_{k}$ and total kinetic energy $\rho e_{k}$ can be well defined:

Because of the long-range molecular interaction, the internal potential energy, defined as $\rho\epsilon_{p}=\frac{1}{2}\rho V_{\text{mean}}$, should be considered. Here, the factor $\tfrac{1}{2}$ avoids counting each interacting pair twice. Therefore, the usual internal energy and total energy are $\rho\epsilon=\rho\epsilon_{k}+\rho\epsilon_{p}$ and $\rho e=\rho e_{k}+\rho\epsilon_{p}$. Through the Chapman-Enskog (CE) analysis, the following mesomacroscopic governing equations can be derived:

where $p_{\text{\tiny BE}}=\rho RT(1+b\rho\chi)$ is the equation of state (EOS) for dense gases recovered by the Enskog collision operator, $\mathbf{F}_{\text{mean}}=-\rho\nabla V_{\text{mean}}$ is the point force for the long-range molecular interaction, $\mathbf{\Pi}$ is the viscous stress tensor, and $\mathbf{J}$ denotes the energy flux by conduction. Note that Eq. (6) should be viewed as mesomacroscopic rather than macroscopic governing equations because the involved $\mathbf{F}_{\text{mean}}$ and $\rho e_{k}$ cannot be well defined from the macroscopic viewpoint.

Equation (6c) in terms of $\rho e_{k}$ is uncommon in previous works. To derive the usual macroscopic energy governing equation, the transport equation for $\rho\epsilon_{p}$ should be first established. The mean-field approximation of the long-range molecular potential is given as Rowlinson and Widom (1982)

where $\mathbf{x}$ and $\mathbf{x}_{2}$ are the positions of two interacting molecules, $V(|\mathbf{x}_{2}-\mathbf{x}|)$ is the distance-dependent potential. Performing Taylor series expansion of $\rho(\mathbf{x}_{2})$ centered at $\mathbf{x}$, Eq. (7) can be formulated as

where $a=-\tfrac{1}{2}\int_{|\mathbf{r}|>d}V(|\mathbf{r}|)d\mathbf{r}$ and $\kappa=-\tfrac{1}{6}\int_{|\mathbf{r}|>d}|\mathbf{r}|^{2}V(|\mathbf{r}|)d\mathbf{r}$. Then, the following relation can be derived:

where $\mathbf{I}$ is the unit tensor, and $\mathbf{P}$ is the pressure tensor defined as $\nabla\cdot\mathbf{P}=\nabla p_{\text{\tiny BE}}-\mathbf{F}_{\text{mean}}$ based on Eq. (6b). Adding $\rho\epsilon_{p}$ to $\rho e_{k}$, Eq. (6c) can be rewritten in terms of $\rho e$:

The last term in Eq. (10) refers to the work done by surface tension He and Doolen (2002). Equations (6c) and (10) are physically equivalent to each other, but Eq. (6c) is much simpler than Eq. (10). This is because $\rho e_{k}$, as a moment of the density DF, is more easily calculated than $\rho e$ at the mesoscopic level, although $\rho e$ is extensively involved at the macroscopic level. Inspired by the above analysis, we will develop a mesoscopic LB model to recover Eq. (6) rather than Eq. (10), and the key points are recovering a nonideal-gas EOS like $p_{\text{\tiny BE}}$ that corresponds to the short-range molecular interaction and mimicking the long-range molecular interaction. Note that both the short- and long-range molecular interactions should be included in physically modeling liquid-vapor phase transition. Otherwise, the liquid-vapor system will suffer from density collapse or be homogenized. Before proceeding further, some discussion on the kinetic model is useful. With Eq. (8), the pressure tensor can be calculated as

where $p_{\text{\tiny EOS}}=p_{\text{\tiny BE}}-a\rho^{2}$ is the full EOS. Obviously, the above $\mathbf{P}$ is consistent with thermodynamic theory. The internal kinetic energy is $\rho\epsilon_{k}=\rho c_{v}T$ according to kinetic theory, and the total kinetic energy satisfies $\rho e_{k}=\rho\epsilon_{k}+\tfrac{1}{2}\rho|\mathbf{u}|^{2}$, where $c_{v}$ is the constant-volume specific heat. The latent heat of vaporization is $h_{lv}=h_{v}-h_{l}=a(\rho_{l}-\rho_{v})+p_{s}\big{(}\rho_{v}^{-1}-\rho_{l}^{-1}\big{)}$, where $h_{v}$ and $h_{l}$ are the specific enthalpies ($h=\epsilon+p_{\text{\tiny EOS}}/\rho$) of the saturated vapor and liquid, respectively, $\rho_{v}$ and $\rho_{l}$ are the saturated vapor and liquid densities, respectively, and $p_{s}$ is the saturation pressure.

Based on Eq. (6) derived from the kinetic model, we introduce double DFs: the density DF $f_{i}(\mathbf{x},t)$ for the mass and momentum conservation laws and an innovative, simple yet effective, total kinetic energy DF $g_{i}(\mathbf{x},t)$ for the energy conservation law. The standard D2Q9 lattice Qian et al. (1992) is considered here for simplicity, and the extension to three dimensions is straightforward. The LB equations for $f_{i}$ and $g_{i}$ are given as

where Eq. (12a) is the linear streaming process in velocity space with $\ell_{i}$ denoting $f_{i}$ or $g_{i}$ and the overbar denoting the post-collision state, Eqs. (12b) and (12c) are the local collision processes in moment space computed at position $\mathbf{x}$ and time $t$ with the moments $\mathbf{m}=\mathbf{M}(f_{i})^{\text{T}}$ and $\mathbf{n}=\mathbf{M}(g_{i})^{\text{T}}$, and $c=\delta_{x}/\delta_{t}$ is the lattice speed. The post-collision DFs are obtained via $(\bar{f}_{i})^{\text{T}}=\mathbf{M}^{-1}\bar{\mathbf{m}}$ and $(\bar{g}_{i})^{\text{T}}=\mathbf{M}^{-1}\bar{\mathbf{n}}$. Here, $\mathbf{M}$ is the orthogonal transformation matrix Lallemand and Luo (2000). A pairwise interaction force is introduced to mimic the long-range molecular interaction, which is given as Huang et al. (2019a)

where $G^{2}$ controls the interaction strength, $\omega(\delta_{x}^{2})=\tfrac{1}{3}$ and $\omega(2\delta_{x}^{2})=\tfrac{1}{12}$ maximize the isotropy degree of $\mathbf{F}_{\text{pair}}$. The density $\rho$, momentum $\rho\mathbf{u}$, and total kinetic energy $\rho e_{k}$ are calculated as

Here, $\mathbf{F}=\mathbf{F}_{\text{pair}}+\rho\mathbf{a}$ is the total force, and $q=\mathbf{F}_{\text{pair}}\cdot\mathbf{u}+\rho\mathbf{a}\cdot\mathbf{u}$ is the total work done by force. Note that $\tfrac{\delta_{t}}{2}\mathbf{F}$ and $\tfrac{\delta_{t}}{2}q$ in Eq. (14) are necessary to avoid the discrete lattice effect.

The technical details of the present mesoscopic LB model (including the equilibrium moments $\mathbf{m}^{\text{eq}}$ and $\mathbf{n}^{\text{eq}}$, the collision matrices $\mathbf{S}$ and $\mathbf{L}$, the discrete force $\mathbf{F}_{m}$, the discrete source $\mathbf{q}_{m}$, etc.) are given in Supplemental Material SM . Performing the second- and third-order CE analyses for the above LB model, the mesomacroscopic governing equations from the kinetic model [i.e., Eq. (6)] can be recovered once we set

where $p_{\text{\tiny LBE}}=c_{s}^{2}(\rho+\eta)$ is the recovered EOS for dense gases with $\eta$ a built-in variable in $\mathbf{m}^{\text{eq}}$, $\mathbf{R}_{\text{iso}}=\tfrac{1}{12}\delta_{x}^{2}\nabla\cdot\nabla\mathbf{F}_{\text{pair}}$ and $\mathbf{R}_{\text{add}}=-\tfrac{G^{2}\delta_{x}^{4}}{24}\nabla\cdot[2\nabla\rho\nabla\rho+(\nabla\rho\cdot\nabla\rho)\mathbf{I}]$ are the third-order terms by the third-order discrete lattice effect and by the compensation term $\mathbf{S}\mathbf{Q}_{m}$ in Eq. (12b), respectively, and $\rho h_{k}$ is the total kinetic enthalpy in $\mathbf{n}^{\text{eq}}$. The recovered viscous stress tensor and energy flux are given as $\mathbf{\Pi}=\rho\nu[\nabla\mathbf{u}+(\nabla\mathbf{u})^{\text{T}}-(\nabla\cdot\mathbf{u})\mathbf{I}]+\rho\varsigma(\nabla\cdot\mathbf{u})\mathbf{I}$ and $\mathbf{J}=\lambda\nabla T$, respectively, with the kinematic viscosity $\nu=c_{s}^{2}\delta_{t}\big{(}s_{p}^{-1}-\tfrac{1}{2}\big{)}$, the bulk viscosity $\varsigma=\varpi c_{s}^{2}\delta_{t}\big{(}s_{e}^{-1}-\tfrac{1}{2}\big{)}$, and the heat conductivity $\lambda=\tfrac{4+3\gamma_{1}+2\gamma_{2}}{6}C_{\text{ref}}c^{2}\delta_{t}\big{(}\sigma_{j}^{-1}-\tfrac{1}{2}\big{)}$. Here, $\varpi$ and $\gamma_{1,2}$ are model coefficients, $s_{e,p}$ and $\sigma_{j}$ are relaxation parameters, $C_{\text{ref}}$ is the reference volumetric heat capacity SM , and $c_{s}=c\big{/}\!\sqrt{3}$ is the lattice sound speed. Based on Eq. (15), the pressure tensor given by Eq. (11) can be derived, and there have

Therefore, thermodynamic consistency naturally emerges from our mesoscopic LB model developed in accordance with the kinetic model. Note that there exist some additional cubic terms of velocity in recovering the viscous stress tensor Dellar (2014); Geier and Pasquali (2018), which are ignored with the low Mach number condition and can also be eliminated by trivial modifications Huang et al. (2019a, 2020). Moreover, the present LB model shows satisfactory numerical stability due to the separate incorporations of the short- and long-range molecular interactions and the introduction of an innovative, simple yet effective, total kinetic energy DF.

In this work, the following full EOS combining the Carnahan-Starling expression for hard spheres Carnahan and Starling (1969) with an attractive term is specified:

where $\vartheta=\tilde{b}\rho/4$, $\tilde{a}=0.4963880577294099R^{2}T_{\text{cr}}^{2}\big{/}p_{\text{cr}}$, and $\tilde{b}=0.1872945669467330RT_{\text{cr}}\big{/}p_{\text{cr}}$. Here, $T_{\text{cr}}$ and $p_{\text{cr}}$ are the critical temperature and pressure, respectively. The interaction strength is set to

and the lattice sound speed is chosen as

Note that the scaling factors $K_{\text{\tiny EOS}}$ and $K_{\text{\tiny INT}}$ are introduced to adjust the surface tension $\sigma\propto K_{\text{\tiny EOS}}K_{\text{\tiny INT}}$ and interface thickness $W\propto K_{\text{\tiny INT}}$.

To test the applicability of our mesoscopic LB model for liquid-vapor phase transition, we perform simulations with $\varpi=1/6$, $\gamma_{1}=-2$, $\gamma_{2}=2$, $\tilde{a}=1$, $\tilde{b}=4$, $R=1$, and $\delta_{x}=1$. The reduced temperature ($T_{r}=T/T_{\text{cr}}$) is set to $T_{r,0}=0.8$, and the surface tension $\sigma=0.01$ and interface thickness $W=10$, which indicate that $K_{\text{\tiny EOS}}=0.479820$ and $K_{\text{\tiny INT}}=2.294922$. The kinematic viscosities and heat conductivities of the liquid and vapor satisfy $\nu_{l}=\nu_{v}$ and $\lambda_{l}=10\lambda_{v}$, respectively. A higher temperature $T_{r,1}=0.85$, together with the outflow and constant-pressure condition, is applied to drive the phase transition. This boundary condition is treated by the improved nonequilibrium-extrapolation scheme Huang et al. (2019c). Meanwhile, eliminating the additional cubic terms of velocity is also plugged into the LB model Huang et al. (2020). Before simulating liquid-vapor phase transition, an equilibrium droplet in periodic domain is considered. The numerical results of $\sigma$ and $W$, measured by Laplace’s law and circular fitting, are $0.0101676$ and $9.961104$, respectively. Such good agreements with the prescribed values validate the present LB model. Subsequently, the one-dimensional Stefan problem is simulated on a $1024\delta_{x}\times 4\delta_{x}$ domain heated from the left side. Neglecting convection and taking the sharp-interface limit, the analytical location of liquid-vapor phase interface can be obtained Solomon (1966):

where $\alpha_{v}=\lambda_{v}\big{/}(\rho_{v}c_{v})$, $t_{0}$ shifts the initial location, and $k$ is the root of the transcendental equation

where the Stefan number is defined as $\text{Ste}=\rho_{v}c_{v}T_{\text{cr}}(T_{r,1}-T_{r,0})\big{/}(\rho_{l}h_{lv})$ and set to $\text{Ste}=0.005$ to ensure that convection can be neglected. The numerical results are shown in Fig. 1. It can be seen that liquid-vapor phase transition is successfully and accurately captured by the present LB model. The vapor slowly flows to the left with its temperature gradually rising from $T_{r,0}$ to $T_{r,1}$, while the liquid stays at rest with a uniform temperature $T_{r,0}$. Across the phase interface, the density profile can be well maintained, and the pressures in vapor and liquid balance each other (the jumps of $p_{\text{\tiny EOS}}$ within the phase interface come from the nonmonotonic EOS for liquid-vapor fluids). Moreover, the location of phase interface agrees very well with the analytical result, which suggests that the latent heat of vaporization in the mesoscopic LB model is naturally consistent with thermodynamic theory.

A liquid droplet with diameter $D_{0}=256\delta_{x}$ is then simulated on a $1024\delta_{x}\times 1024\delta_{x}$ domain heated from all the four sides. The Stefan number is set to $\text{Ste}=0.005$ and thus convection in the evaporation is quite weak. Figure 2 shows the time evolution of the square of the normalized diameter $(D/D_{0})^{2}$, together with four snapshots of the local density and temperature fields. Here, the time is normalized as $t^{\ast}=\alpha_{v}t/D_{0}^{2}$. The well-known $D^{2}$ law Godsave (1953); Safari et al. (2013) can be perfectly observed during the entire droplet lifetime, and both the interface thickness and droplet temperature can be well maintained at the prescribed values. As a further application, the evaporation of a large-small droplet pair is simulated with $\text{Ste}=0.005$, $0.05$, and $0.5$, respectively. Initially, the diameters of the two droplets are $160\delta_{x}$ and $96\delta_{x}$, respectively, and the distance between the droplet centers is $256\delta_{x}$. Figure 3 shows the snapshots of the local temperature and velocity fields, and the time evolution of the normalized volume $V/V_{0}$. Here, $V$ is the total volume of the droplets, $V_{0}=\pi D_{0}^{2}\big{/}4$, and $D_{0}=128\delta_{x}$. For $\text{Ste}=0.005$, the evaporation is quite slow, and the two droplets attract each other and coalesce into a single one. This attraction-coalescence behavior is due to the nonuniform evaporation rate along droplet surface, which is induced by the other droplet and will result in an imbalanced vapor recoil force Nikolayev et al. (2006). Such unusual behavior of evaporating droplets under slow evaporation condition is consistent with the recent experimental and theoretical results Wen et al. (2019); Man and Doi (2017). Interestingly, the local temperature slightly rises [see the middle panel in Fig. 3(a)] and the normalized volume slightly increases [see the “kink” in Fig. 3(d)] when the coalescence occurs, which can be explained as follows: At the neck formed by coalescence, the phase interface changes from convex to concave, and the local saturated vapor pressure will decrease according to the Kelvin equation in thermodynamic theory Rowlinson and Widom (1982). Therefore, the vapor nearby the neck becomes supersaturated and then condenses into liquid, resulting in the release of latent heat and also the increase of droplet volume. Here, it is noteworthy that the above condensation at the neck between two merging droplets is a kind of capillary condensation in thermodynamic theory Fisher et al. (1981); Yang et al. (2020). For $\text{Ste}=0.05$ and $0.5$, evaporation becomes much faster and convection is very strong. The two droplets repulse each other rather than attract, and the droplet lifetime is much shorter than that for $\text{Ste}=0.005$. As seen in Fig. 3(b) and 3(c), the vapor outflows originating from the droplet surfaces impact each other in the middle region between the two droplets, and thus the pressure in this region obviously increases, which then pushes the two droplets away from each other against the imbalanced vapor recoil force.

In summary, we have developed a novel mesoscopic LB model for liquid-vapor phase transition, where the short- and long-range molecular interactions are incorporated by recovering an EOS for dense gases and introducing a pairwise interaction force, respectively, and an innovative, simple yet effective, total kinetic energy DF is proposed for the energy conservation law. The same mesomacroscopic governing equations as the kinetic model can be recovered, and thus thermodynamic consistency is naturally satisfied. Because of the successful modeling of the underlying microscopic molecular interaction, the present mesoscopic LB model does not rely on any macroscopic assumptions or approximations and has the potential to provide reliable physical insights into the liquid-vapor phase transition processes.

See pages 1, of Supplemental_Material.pdf See pages 2, of Supplemental_Material.pdf