# Mesoscopic Lattice Boltzmann Modeling of the Liquid-Vapor Phase Transition

Rongzong Huang<EMAIL_ADDRESS>School of Energy Science and
Engineering, Central South University, 410083 Changsha, China  Huiying Wu
<EMAIL_ADDRESS>School of Mechanical Engineering, Shanghai Jiao Tong
University, 200240 Shanghai, China  Nikolaus A. Adams<EMAIL_ADDRESS>Institute of Aerodynamics and Fluid Mechanics, Technical University of Munich,
85748 Garching, Germany

(April 22, 2021)

###### Abstract

We develop a mesoscopic lattice Boltzmann model for liquid-vapor phase
transition by handling the microscopic molecular interaction. The short-range
molecular interaction is incorporated by recovering an equation of state for
dense gases, and the long-range molecular interaction is mimicked by
introducing a pairwise interaction force. Double distribution functions are
employed, with the density distribution function for the mass and momentum
conservation laws and an innovative total kinetic energy distribution function
for the energy conservation law. The recovered mesomacroscopic governing
equations are fully consistent with kinetic theory, and thermodynamic
consistency is naturally satisfied.

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)

$\partial_{t}f+{\bm{\xiup}}\cdot\nabla
f+\mathbf{a}\cdot\nabla_{\bm{\xiup}}f=\mathit{\Omega}_{\text{Enskog}}+\nabla
V_{\text{mean}}\cdot\nabla_{\bm{\xiup}}f,$ (1)

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)

$\displaystyle\mathit{\Omega}_{\text{Enskog}}=\chi\mathit{\Omega}_{0}-b\rho\chi
f^{\text{eq}}\Big{\\{}\tfrac{2}{5}\left[2{\bm{\mathcal{C}}}{\bm{\mathcal{C}}}:\nabla\mathbf{u}+\left(|{\bm{\mathcal{C}}}|^{2}-\tfrac{5}{2}\right)\nabla\cdot\mathbf{u}\right]$
(2) $\displaystyle+\mathbf{C}\cdot\left[\nabla\ln(\rho^{2}\chi
T)+\tfrac{3}{5}\left(|{\bm{\mathcal{C}}}|^{2}-\tfrac{5}{2}\right)\nabla\ln
T\right]$ $\displaystyle\Big{\\}},$

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

$f^{\text{eq}}=\dfrac{\rho}{(2\pi
RT)^{3/2}}\exp\left(-|{\bm{\mathcal{C}}}|^{2}\right).$ (3)

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

$\rho=\int fd{\bm{\xiup}},\quad\rho\mathbf{u}=\int
f{\bm{\xiup}}d{\bm{\xiup}}.$ (4)

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

$\rho\epsilon_{k}=\int
f\dfrac{|{\bm{\xiup}}-\mathbf{u}|^{2}}{2}d{\bm{\xiup}},\quad\rho e_{k}=\int
f\dfrac{|{\bm{\xiup}}|^{2}}{2}d{\bm{\xiup}}.$ (5)

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:

$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho\mathbf{u})=0,$ (6a)
$\displaystyle\partial_{t}(\rho\mathbf{u})+\nabla\cdot(\rho\mathbf{u}\mathbf{u})=-\nabla
p_{\text{\tiny
BE}}+\mathbf{F}_{\text{mean}}+\rho\mathbf{a}+\nabla\cdot\mathbf{\Pi},$ (6b)
$\displaystyle\partial_{t}(\rho e_{k})+\nabla\cdot(\rho
e_{k}\mathbf{u}+p_{\text{\tiny
BE}}\mathbf{u})=\mathbf{F}_{\text{mean}}\cdot\mathbf{u}+\rho\mathbf{a}\cdot\mathbf{u}+\nabla\cdot(\mathbf{J}+\mathbf{u}\cdot\mathbf{\Pi}),$
(6c)

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)

$V_{\text{mean}}=\int_{|\mathbf{x}_{2}-\mathbf{x}|>d}\rho(\mathbf{x}_{2})V(|\mathbf{x}_{2}-\mathbf{x}|)d\mathbf{x}_{2},$
(7)

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

$V_{\text{mean}}=-2a\rho-\kappa\nabla\cdot\nabla\rho,$ (8)

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:

$\displaystyle\partial_{t}(\rho\epsilon_{p})+\nabla\cdot(\rho\epsilon_{p}\mathbf{u})+\mathbf{F}_{\text{mean}}\cdot\mathbf{u}=\tfrac{1}{2}\rho(\partial_{t}V_{\text{mean}}-\mathbf{u}\cdot\nabla
V_{\text{mean}})$ (9)
$\displaystyle=-\nabla\cdot[\mathbf{u}\cdot(\mathbf{P}-p_{\text{\tiny
BE}}\mathbf{I})]+\nabla\mathbf{u}:\left[-\tfrac{\kappa}{2}\nabla\cdot(\rho\nabla\rho)\mathbf{I}+\kappa\nabla(\rho\nabla\rho)\right]$
,

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$:

$\displaystyle\partial_{t}(\rho e)+\nabla\cdot(\rho
e\mathbf{u}+\mathbf{u}\cdot\mathbf{P})=\rho\mathbf{a}\cdot\mathbf{u}+\nabla\cdot(\mathbf{J}+\mathbf{u}\cdot\mathbf{\Pi})$
(10)
$\displaystyle+\nabla\mathbf{u}:\left[-\tfrac{\kappa}{2}\nabla\cdot(\rho\nabla\rho)\mathbf{I}+\kappa\nabla(\rho\nabla\rho)\right]$
.

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

$\mathbf{P}=\left(p_{\text{\tiny
EOS}}-\kappa\rho\nabla\cdot\nabla\rho-\tfrac{\kappa}{2}\nabla\rho\cdot\nabla\rho\right)\mathbf{I}+\kappa\nabla\rho\nabla\rho,$
(11)

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

$\displaystyle\ell_{i}(\mathbf{x}+\mathbf{e}_{i}\delta_{t},t+\delta_{t})=\bar{\ell}_{i}(\mathbf{x},t),$
(12a)
$\displaystyle\bar{\mathbf{m}}=\mathbf{m}+\delta_{t}\mathbf{F}_{m}-\mathbf{S}\left(\mathbf{m}-\mathbf{m}^{\text{eq}}+\tfrac{\delta_{t}}{2}\mathbf{F}_{m}\right)+\mathbf{S}\mathbf{Q}_{m},$
(12b)
$\displaystyle\bar{\mathbf{n}}=\mathbf{n}+\delta_{t}\mathbf{q}_{m}-\mathbf{L}\left(\mathbf{n}-\mathbf{n}^{\text{eq}}+\tfrac{\delta_{t}}{2}\mathbf{q}_{m}\right)+c^{2}\mathbf{Y}\left(\tfrac{\mathbf{m}+\bar{\mathbf{m}}}{2}-\mathbf{m}^{\text{eq}}\right),$
(12c)

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)

$\mathbf{F}_{\text{pair}}=G^{2}\rho(\mathbf{x})\sum\nolimits_{i}\omega(|\mathbf{e}_{i}\delta_{t}|^{2})\rho(\mathbf{x}+\mathbf{e}_{i}\delta_{t})\mathbf{e}_{i}\delta_{t},$
(13)

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

$\rho=\sum\nolimits_{i}f_{i},\quad\rho\mathbf{u}=\sum\nolimits_{i}\mathbf{e}_{i}f_{i}+\tfrac{\delta_{t}}{2}\mathbf{F},\quad\rho
e_{k}=\sum\nolimits_{i}g_{i}+\tfrac{\delta_{t}}{2}q.$ (14)

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

$\begin{array}[]{l}p_{\text{\tiny BE}}=p_{\text{\tiny
LBE}},\quad\mathbf{F}_{\text{mean}}=\mathbf{F}_{\text{pair}}+\mathbf{R}_{\text{iso}}+\mathbf{R}_{\text{add}},\\\\[3.01385pt]
\rho h_{k}=\rho e_{k}+p_{\text{\tiny BE}},\end{array}$ (15)

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

$a=\tfrac{G^{2}\delta_{x}^{2}}{2},\quad\kappa=\tfrac{G^{2}\delta_{x}^{4}}{4}.$
(16)

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:

$p_{\text{\tiny EOS}}=K_{\text{\tiny EOS}}\left[\rho
RT\tfrac{1+\vartheta+\vartheta^{2}-\vartheta^{3}}{(1-\vartheta)^{3}}-\tilde{a}\rho^{2}\right],$
(17)

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

$G=K_{\text{\tiny INT}}\sqrt{2K_{\text{\tiny
EOS}}\tilde{a}\big{/}\delta_{x}^{2}}\,,$ (18)

and the lattice sound speed is chosen as

$c_{s}=K_{\text{\tiny INT}}\left.\sqrt{\left(\tfrac{\partial p_{\text{\tiny
EOS}}}{\partial\rho}\right)_{T}+2K_{\text{\tiny
EOS}}\tilde{a}\rho}\,\right|\,\raisebox{-7.3194pt}{$\scriptstyle\rho\,=\,\rho_{l}$}.$
(19)

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

$X_{i}(t)=2k\sqrt{\alpha_{v}(t+t_{0})},$ (20)

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

$\tfrac{\text{Ste}}{\exp(k^{2})\,\text{erf}(k)}=k\sqrt{\pi},$ (21)

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.

Figure 1: Distributions of (a) density $\rho$, (b) velocity $u_{x}$, (c)
temperature $T_{r}$, and (d) pressure $p_{\text{\tiny EOS}}$ at time
$t=0.401\times 10^{7}$, $0.803\times 10^{7}$, $1.204\times 10^{7}$, and
$1.606\times 10^{7}$, and (e) time evolution of the phase interface location
$X_{i}$ for the one-dimensional Stefan problem.

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.

Figure 2: Evolution of the square of the normalized droplet diameter
$(D/D_{0})^{2}$ with the normalized time
$t^{\ast}=\alpha_{v}t\big{/}D_{0}^{2}$ in the droplet evaporation process.
Snapshots of the local density (left) and temperature (right) fields are also
shown at time $t^{\ast}=12.716$, $38.147$, $63.578$, and $89.010$, where the
solid line in temperature field denotes the liquid-vapor phase interface.
Figure 3: Snapshots of the local temperature and velocity fields at different
normalized times for (a) $\text{Ste}=0.005$, (b) $\text{Ste}=0.05$, and (c)
$\text{Ste}=0.5$, and (d) evolutions of the normalized total volume of the
droplets $V/V_{0}$ with the normalized time
$t^{\ast}=\alpha_{v}t\big{/}D_{0}^{2}$ in the evaporation process of a large-
small droplet pair. The solid and dotted lines in temperature field denote the
liquid-vapor phase interface and its initial location, respectively.

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.

###### Acknowledgements.

R.H. acknowledges the support by the Alexander von Humboldt Foundation,
Germany. This work was supported by the National Natural Science Foundation of
China through Grants No. 52006244 and No. 51820105009.

## References

  * Cheng  _et al._ (2014) P. Cheng, X. Quan, S. Gong, X. Liu, and L. Yang, Adv. Heat Transf. 46, 187 (2014).EOS
  * Onuki (2005) A. Onuki, Phys. Rev. Lett. 94, 054501 (2005).EOS
  * Boreyko and Chen (2009) J. B. Boreyko and C.-H. Chen, Phys. Rev. Lett. 103, 184501 (2009).EOS
  * Cira  _et al._ (2015) N. J. Cira, A. Benusiglio, and M. Prakash, Nature 519, 446 (2015).EOS
  * Nikolayev and Beysens (1999) V. S. Nikolayev and D. A. Beysens, Europhys. Lett. 47, 345 (1999).EOS
  * Nikolayev  _et al._ (2006) V. S. Nikolayev, D. Chatain, Y. Garrabos, and D. Beysens, Phys. Rev. Lett. 97, 184503 (2006).EOS
  * He and Doolen (2002) X. He and G. D. Doolen, J. Stat. Phys. 107, 309 (2002).EOS
  * He  _et al._ (1998) X. He, X. Shan, and G. D. Doolen, Phys. Rev. E 57, R13 (1998).EOS
  * Luo (1998) L.-S. Luo, Phys. Rev. Lett. 81, 1618 (1998).EOS
  * Gunstensen  _et al._ (1991) A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Phys. Rev. A 43, 4320 (1991).EOS
  * Shan and Chen (1993) X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993).EOS
  * Swift  _et al._ (1995) M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett. 75, 830 (1995).EOS
  * Hazi and Markus (2009) G. Hazi and A. Markus, Int. J. Heat Mass Tran. 52, 1472 (2009).EOS
  * Gong and Cheng (2012) S. Gong and P. Cheng, Int. J. Heat Mass Tran. 55, 4923 (2012).EOS
  * Li  _et al._ (2015) Q. Li, Q. J. Kang, M. M. Francois, Y. L. He, and K. H. Luo, Int. J. Heat Mass Tran. 85, 787 (2015).EOS
  * Albernaz  _et al._ (2017) D. L. Albernaz, M. Do-Quang, J. C. Hermanson, and G. Amberg, J. Fluid Mech. 820, 61 (2017).EOS
  * Qin  _et al._ (2019) F. Qin, L. Del Carro, A. Mazloomi Moqaddam, Q. Kang, T. Brunschwiler, D. Derome, and J. Carmeliet, J. Fluid Mech. 866, 33 (2019).EOS
  * Laurila  _et al._ (2012) T. Laurila, A. Carlson, M. Do-Quang, T. Ala-Nissila, and G. Amberg, Phys. Rev. E 85, 026320 (2012).EOS
  * Rowlinson and Widom (1982) J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Oxford University Press, Oxford, 1982).EOS
  * Chapman and Cowling (1970) S. Chapman and T. Cowling, The Mathematical Theory of Non-Uniform Gases, 3rd ed. (Cambridge University Press, Cambridge, 1970).EOS
  * Qian  _et al._ (1992) Y. H. Qian, D. d’Humières, and P. Lallemand, Europhys. Lett. 17, 479 (1992).EOS
  * Lallemand and Luo (2000) P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546 (2000).EOS
  * Huang  _et al._ (2019a) R. Huang, H. Wu, and N. A. Adams, Phys. Rev. E 99, 023303 (2019a).EOS
  * (24) See Supplemental Material for technical details of the present mesoscopic LB model, which includes Refs. Huang and Wu (2014); Huang  _et al._ (2019b).
  * Huang and Wu (2014) R. Huang and H. Wu, J. Comput. Phys. 274, 50 (2014).EOS
  * Huang  _et al._ (2019b) R. Huang, H. Wu, and N. A. Adams, Phys. Rev. E 100, 043306 (2019b).EOS
  * Dellar (2014) P. J. Dellar, J. Comput. Phys. 259, 270 (2014).EOS
  * Geier and Pasquali (2018) M. Geier and A. Pasquali, Comput. Fluids 166, 139 (2018).EOS
  * Huang  _et al._ (2020) R. Huang, L. Lan, and Q. Li, Phys. Rev. E 102, 043304 (2020).EOS
  * Carnahan and Starling (1969) N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).EOS
  * Huang  _et al._ (2019c) R. Huang, H. Wu, and N. A. Adams, J. Comput. Phys. 392, 227 (2019c).EOS
  * Solomon (1966) A. Solomon, Math. Comp. 20, 347 (1966).EOS
  * Godsave (1953) G. A. E. Godsave, Symp. (Int.) Combust. 4, 818 (1953).EOS
  * Safari  _et al._ (2013) H. Safari, M. H. Rahimian, and M. Krafczyk, Phys. Rev. E 88, 013304 (2013).EOS
  * Wen  _et al._ (2019) Y. Wen, P. Y. Kim, S. Shi, D. Wang, X. Man, M. Doi, and T. P. Russell, Soft Matter 15, 2135 (2019).EOS
  * Man and Doi (2017) X. Man and M. Doi, Phys. Rev. Lett. 119, 044502 (2017).EOS
  * Fisher  _et al._ (1981) L. R. Fisher, R. A. Gamble, and J. Middlehurst, Nature 290, 575 (1981).EOS
  * Yang  _et al._ (2020) Q. Yang, P. Z. Sun, L. Fumagalli, Y. V. Stebunov, S. J. Haigh, Z. W. Zhou, I. V. Grigorieva, F. C. Wang, and A. K. Geim, Nature 588, 250 (2020).EOS

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