Kinetic modeling of multiphase flow based on simplified Enskog equation

Abstract

A new kinetic model for multiphase flow was presented under the framework of the discrete Boltzmann method (DBM). Significantly different from the previous DBM, a bottom-up approach was adopted in this model. The effects of molecular size and repulsion potential were described by the Enskog collision model; the attraction potential was obtained through the mean-field approximation method. The molecular interactions, which result in the non-ideal equation of state and surface tension, were directly introduced as an external force term. Several typical benchmark problems, including Couette flow, two-phase coexistence curve, the Laplace law, phase separation, and the collision of two droplets, were simulated to verify the model. Especially, for two types of droplet collisions, the strengths of two non-equilibrium effects, {\overline{D}}_2^{*} and {\overline{D}}_3^{*}, defined through the second and third order non-conserved kinetic moments of (f-f^eq), are comparatively investigated, where f (f^eq) is the (equilibrium) distribution function. It is interesting to find that during the collision process, {\overline{D}}_2^{*} is always significantly larger than {\overline{D}}_3^{*}, {\overline{D}}_2^{*} can be used to identify the different stages of the collision process and to distinguish different types of collisions. The modeling method can be directly extended to a higher-order model for the case where the non-equilibrium effect is strong, and the linear constitutive law of viscous stress is no longer valid.

Graphical abstract

Keywords

multiphase flow / discrete Boltzmann method / Enskog equation / non-equilibrium characteristics

1 Introduction

Multiphase flow is ubiquitous in nature and industrial processes, such as the formation of raindrops, fabrication of pills, oil and gas exploitation, micro-fluidic chip [1-4]. Understanding the hydrodynamic and thermodynamic characteristics of multi-phase flow is of great significance for fundamental research and engineering applications.

Generally, multiphase flow refers to the fluid flow consisting of two or more phases. It is more complicated than a single-phase flow process, due to phase interfaces and the interaction between the different phases [5]. The multiphase flow process is often accompanied by abundant interfacial dynamic behaviors, including interface generation, movement, deformation, fragmentation, and fusion, etc. These complex interfacial behaviors are challenging to measure experimentally and very difficult, if not impossible, to obtain analytic solutions.

As computation power increased in recent decades, numerical simulation has gradually become a critical alternative approach to study multiphase flow problems [6-8]. In general, the numerical simulation methods of multiphase flow can be divide into three levels: macroscopic, microscopic, and mesoscopic methods. The macroscopic approaches are mainly based on the numerical solutions of the hydrodynamic (Navier–Stokes, NS) equations combined with interface capture schemes [6]. Although those methods have been successfully applied to the large-scale multiphase flow problems, they fail to describe the micro-nano scale interface structure because the equations based on the continuum assumption may not be applicable in the micro-nano scale [8, 9]. In addition, the fluctuations are significant with the decreasing of system size and can even perceivably affect the flow characteristics at the micro-nano scale [10]. The surface tension and viscosities might need to be further considered with the effect of fluctuation. The microscopic model mainly refers to the molecular dynamic (MD) method, which is more fundamental and accurate [11]. However, the computational burden of MD simulation is often unaffordable; consequently, the time and space scales it can simulate are limited. As a bridge between macroscopic and microscopic, the mesoscopic model has developed rapidly over the past three decades [12-16]. It has been successfully applied in a wide variety of multiphase problems, including wetting, droplet collision and splashing, boiling and evaporation, phase separation, and hydrodynamic instability [6, 17-19], etc.

The mesoscopic model for multiphase flow, also known as the kinetic model, is mainly based on the Boltzmann equation. The vast majority, if not all, of the mesoscopic multiphase flow models in the literature, are developed from the lattice Boltzmann model (LBM), including the color-gradient model [20], the Shan–Chen model [21], the free energy model [22, 23], and the phase-field model [24-27], etc. It is undeniable that they are very convenient in modeling multiphase flow processes, and have achieved great success in the research of various multiphase flow problems [17]. However, it should be noted that the standard LBM, on which most multiphase flow models base, is just a solver of the incompressible Navier–Stokes (NS) equations. Although it may be more convenient to deal with the interparticle forces and the complicated boundaries [13], it can not provide kinetic information beyond the NS equations.

Based on the thermal finite different LBM of Watari and Tsutahara [28], Gonnella et al. first proposed the thermal multiphase model [29], in which the nonideal gases and the surface tension effects are introduced by an extra term I_{ki}. This multiphase flow model can correctly reproduce the transport equations of nonideal fluid established by Onuki [30], through a Chapman–Enskog multiscale expansion. Gan et al. [31] further developed a thermal multiphase model with negligible spurious velocities by introducing the Windowed Fast Fourier Transform scheme when calculating the spatial derivatives in the convection term and the external force term. Then morphological characteristics of phase separation in the thermal systems were systematically studied [32].

In the early studies, the non-equilibrium behaviors in various complex flow were all those described by hydrodynamic equations. For the convenience of description, we refer those behaviors to as Hydrodynamic Non-Equilibrium (HNE) behaviors. Besides the HNE, the most relevant Thermodynamic Non-Equilibrium (TNE) practices in various complex flow systems are attracting more and more attention with time [33-36].

It is understandable that when the TNE is very weak, the loss of considering TNE does not make a noticeable effect. But when the TNE is strong, the situation will be significantly different. For example, it has been shown that TNE’s existence directly affects the density, temperature, and pressure, as well as the magnitude and direction of flow velocity. Without considering TNE, the density, flow velocity, temperature, and pressure given will have a significant deviation [37]. The existence of TNE is the underlying cause of the appearance of heat flow and viscous stress. If the TNE is not insufficiently considered, i.e., only the linear response part is considered, the amplitude of heat flow and viscous stress obtained may be too large [38].

When the strength of TNE beyond some threshold value, it may change the directions of heat flow and viscous stress. If the TNE is not sufficiently considered (considered only the linear response part), the obtained viscous stress and heat flow maybe, even in the wrong directions [9]. In addition, in the phase separation system, both the mean TNE strength [39] and the entropy production rate increase with time in the spinnodal decomposition stage and decrease with time in the domain growth stage. Therefore, both the peak value points of the mean TNE strength and the entropy production rate can work as physical criteria to discriminate the two stages [35, 39, 40]. In the system with combustion, the TNE behaviors help to better understand the physical structures of the von Neumann peak and various non-equilibrium detonation [34, 37, 41-47]. In the system with Rayleigh–Taylor instability, the TNE behaviors around interfaces have been used to physically identify and distinguish various interfaces and design relevant interface-tracking schemes. With increasing the compressibility, more observable TNE kinetic modes appear under a given observation precision [36, 48]. The correlation between the mean density nonuniformity and mean TNE strength, the mean temperature nonuniformity and mean Non-Organized Energy Flux (NOEF) is almost 1. The correlation between the mean flow velocity nonuniformity and the mean Non-Organized Momentum Flux (NOMF) is also high, but generally less than 1 [36, 49]. The TNE effect helps to understand the effect of system dispersion (described by Knudsen number) on its kinetic behaviors [50]. In the system with Kelvin–Helmholtz instability, via some defined TNE quantity, for example, the heat flux intensity, the density interface and temperature interface can both be observed so that the material mixing and energy mixing in the Kelvin–Helmholtz instability evolution can be investigated simultaneously [9, 51]. The TNE behaviors can help to understand the mixing entropy in multi-component flows better [52].

Since the traditional hydrodynamic model is incapable of capturing enough TNE behaviors, the works mentioned above on both HNE and TNE resorted to the recently proposed discrete Boltzmann method (DBM) [33, 36]. Compared with the standard LBM, the model system’s evolution does not resort to the simple “propagation + collision" scenario of virtual particle which is inherited from the lattice-gas-automaton model. Any suitable numerical scheme can be used to solve the discrete Boltzmann equation(DBE) according to the specific case. The non-conserved kinetic moments of (f_i-f_i^{eq}) are used to describe how the system deviates from its thermodynamic equilibrium state and to fetch information on the TNE behaviors, where f_i (f_i^{eq}) is the discrete (equilibrium) distribution function and i is the index of discrete velocity. The kinetic moment relations to be based in the construction process of DBM are dependent on the research perspective and the non-equilibrium extent which the DBM aims to describe [36]. Besides by theory, results of DBM have been confirmed and supplemented by results of molecular dynamics [53-55], direct simulation Monte Carlo [9, 56] and experiment [46].

The existing multiphase flow models based on DBM is a semi-kinetic and semi-top-down approach. The discrete Boltzmann equation is modified by an extra term to be consistent with the correct hydrodynamic equations of nonideal fluid. The specific forms of the coefficients (expressed as the gradient of macroscopic quantities) in the external force term are determined by the Chapman–Enskog multiscale expansion [29, 57]. Therefore, in the modeling of multiphase flow DBM, the hydrodynamic equations of nonideal fluid need to be known first. This partly limits the ability of DBM to study the non-equilibrium complex flow, because the hydrodynamic equations are often unknown when the non-equilibrium strength is strong. In this paper, this restriction is addressed by developing a bottom-up kinetic modeling method of multiphase flow under the framework of DBM.

The rest of this paper is organized as follows. In Section 2, the simplified Enskog equation is derived, and the multiphase model based on the discrete Enskog equation is presented. Several benchmark tests of multiphase flow are simulated to verify the new model, and the non-equilibrium characteristics during the collision of two droplets are analyzed in Section 3. Section 4 concludes the present paper.

2 Models and methods

2.1 Enskog equation and its simplification

According to molecular kinetic theory, the evolution equation of molecular velocity distribution function reads

\begin{array}{}\text{(1)}& \frac{\partial f}{\partial t}+v\cdot \nabla f+a\cdot {\nabla }_{v}f={\left(\frac{\partial f}{\partial t}\right)}_c\mathrm{,}\end{array}

where f=f(r, v, t) is the molecular velocity distribution function, r and v represent the position space and velocity space coordinates, respectively, a is the acceleration generated by the total extra force. {\left(\frac{\partial f}{\partial t}\right)}_c denotes the rate of change in distribution function due to the collisions between molecules. Given that the probability of two molecules colliding is much higher than that of three or more molecules colliding simultaneously, the hypothesis of two-body collision is reasonable.

Based on the elastic molecular collision model, when the molecular volume effect is ignored, the collision term in the Boltzmann equation can be obtained as [58]

\begin{array}{}\text{(2)}& {\left(\frac{\partial f}{\partial t}\right)}_c={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\mathrm{,}\end{array}

where f^{*} and f_1^{*} represent the post-collision molecular velocity distribution function with velocity v^* and v_1^{*}, respectively, f and f₁ represent the pre-collision molecular velocity distribution function with velocity v and v₁, respectively. v_r=\left|v-v_1\right| is the value of the relative velocity, which remains unchanged pre to and after the collision. σ and Ω denote the differential collision cross-section and solid angle, respectively.

However, this assumption is not appropriate for dense gases or liquid. With increasing the number density, compared with the mean distance between neighboring molecules, the size of a gas molecule is no longer negligible anymore, the volume effects of the molecule should be taken into account, the Boltzmann collision operator should be replaced by the Enskog collision operator which reads [59]

\begin{array}{}\text{(3)}& \begin{array}{ll}{\left(\frac{\partial f}{\partial t}\right)}_E=\hfill & {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}[\chi \left(r+\frac{d_0}{2}{\widehat{e}}_r\right)f^{*}\left(r\right)f_1^{*}\left(r+d_0{\widehat{e}}_r\right)\hfill \\ \hfill & -\chi \left(r-\frac{d_0}{2}{\widehat{e}}_r\right)f\left(r\right)f_1\left(r-d_0{\widehat{e}}_r\right)]v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\mathrm{,}\hfill \end{array}\end{array}

where d₀ is the diameter of the hard-sphere molecules and χ represents the collision probability correction considering the molecular volume effect. By using the Taylor expansion and keeping to the first derivative term, one can get that

\begin{array}{}\text{(3)}& \begin{array}{l}\chi \left(r+\frac{d_0}{2}{\widehat{e}}_r\right)=\chi (r)+\frac{d_0}{2}\nabla \chi \cdot {\widehat{e}}_r\mathrm{,}\\ \chi \left(r-\frac{d_0}{2}{\widehat{e}}_r\right)=\chi (r)-\frac{d_0}{2}\nabla \chi \cdot {\widehat{e}}_r\mathrm{,}\\ f_1^{*}\left(r+d_0{\widehat{e}}_r\right)=f_1^{*}(r)+d_0\nabla f_1^{*}\cdot {\widehat{e}}_r\mathrm{,}\\ f_1\left(r-d_0{\widehat{e}}_r\right)=f_1(r)-d_0\nabla f_1\cdot {\widehat{e}}_r\mathrm{.}\end{array}\end{array}

Consequently, the collision term in Eq. (3) becomes

\begin{array}{}\text{(4)}& \begin{array}{l}{\left(\frac{\partial f}{\partial t}\right)}_E=\chi {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\\ +d_0\chi {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}\nabla f_1^{*}+f\nabla f_1\right)\cdot {\widehat{e}}_r{v}_r\sigma \text{d}\text{Ω}\text{d}{v}_1\\ +\frac{1}{2}{d}_0{\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\nabla \chi \cdot {\widehat{e}}_r\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\end{array}\end{array}

where χ, f_1^{*}, and f_1 are all at r. If f in the latter two terms are approximated by the local equilibrium distribution function f^{eq},

\begin{array}{}\text{(5)}& f^{eq}=\rho \frac{1}{{\left(2\text{π}RT\right)}^{\mathrm{3/2}}}\exp [-\frac{{(v-u)}^2}{RT}]\mathrm{.}\end{array}

Then the Enskog collision operator becomes

\begin{array}{}\text{(6)}& \begin{array}{l}{\left(\frac{\partial f}{\partial t}\right)}_E=\chi {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\\ -f^{eq}b\rho \chi \{\left(v-u\right)\cdot \left[\frac{2}{\rho }\nabla \rho +\frac{3}{5T}\nabla T\left(\frac{{\left(v-u\right)}^2}{2RT}-\frac{5}{2}\right)\right]\\ +\frac{3}{5RT}\left(v-u\right)\left(v-u\right):\nabla u+\frac{3}{5}\left[\frac{{\left(v-u\right)}^2}{2RT}-\frac{5}{2}\right]\nabla \cdot u\}\\ -f^{eq}b\rho \left(v-u\right)\cdot \nabla \chi \mathrm{,}\end{array}\end{array}

where b\rho =\frac{2}{3}\text{π}n{d}_0{}^3.

For isothermal systems, the collision term can be reduced to

\begin{array}{}\text{(7)}& \begin{array}{l}{\left(\frac{\partial f}{\partial t}\right)}_E=\chi {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\\ -f^{eq}b\rho \chi \{\left(v-u\right)\cdot \frac{2}{\rho }\nabla \rho +\frac{3}{5RT}\left(v-u\right)\left(v-u\right):\nabla u\\ +\frac{3}{5}\left[\frac{{\left(v-u\right)}^2}{2RT}-\frac{5}{2}\right]\nabla \cdot u\}-f^{eq}b\rho \left(v-u\right)\cdot \nabla \chi \mathrm{.}\end{array}\end{array}

In addition, the last two terms in the second row of Eq. (7) satisfy the following relationships:

\begin{array}{}\text{(8)}& \begin{array}{l}{\displaystyle \int }{f}^{eq}b\rho \chi \{\frac{3}{5RT}(v-u)(v-u):\nabla u\\ +\frac{3}{5}\left[\frac{{(v-u)}^2}{2RT}-\frac{5}{2}\right]\nabla \cdot u\}\text{d}v=\mathrm{0,}\end{array}\end{array}

\begin{array}{}\text{(9)}& \begin{array}{l}{\displaystyle \int }{f}^{eq}b\rho \chi \{\frac{3}{5RT}(v-u)(v-u):\nabla u\\ +\frac{3}{5}\left[\frac{{(v-u)}^2}{2RT}-\frac{5}{2}\right]\nabla \cdot u\}v\text{d}v=0.\end{array}\end{array}

The sum of these two terms have no effect on the evolution of density and momentum [60], so they can be ignored in the modeling framework of this paper. As a result, the Enskog collision term can be further simplified as

\begin{array}{}\text{(10)}& \begin{array}{l}{\left(\frac{\partial f}{\partial t}\right)}_E=\chi {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{4\text{π}}}\left(f^{*}{f}_1^{*}-f{f}_1\right)v_r\sigma \text{d}\text{Ω}\text{d}{v}_1\\ -f^{eq}b\rho \chi \left(v-u\right)\cdot \frac{2}{\rho }\nabla \rho -f^{eq}b\rho \left(v-u\right)\cdot \nabla \chi \mathrm{.}\end{array}\end{array}

Similar to the modeling of DBM, the first term on the right-hand side can be replaced by the BGK collision operator, consequently Eq. (7) becomes

\begin{array}{}\text{(11)}& {\left(\frac{\partial f}{\partial t}\right)}_E=-\frac{1}{\tau }\left(f-f^{eq}\right)-f^{eq}b\rho \chi \left(v-u\right)\cdot \nabla \ln \left({\rho }^2\chi \right)\mathrm{.}\end{array}

Now we can get the simplified Enskog model for the approximately incompressible fluid systems when the temperature is uniform, which reads

\begin{array}{}\text{(12)}& \begin{array}{l}\frac{\partial f}{\partial t}+v\cdot \nabla f+a\cdot {\nabla }_{v}f=-\frac{1}{\tau }\left(f-f^{eq}\right)\\ -f^{eq}b\rho \chi \left(v-u\right)\cdot \nabla \ln \left({\rho }^2\chi \right)\mathrm{.}\end{array}\end{array}

Similar to previous studies, in this work we consider the case where the force term a\cdot {\nabla }_{v}f can be approximated by

\begin{array}{}\text{(13)}& a\cdot {\nabla }_{v}f=a\cdot {\nabla }_v{f}^{eq}=-\frac{F\cdot (v-u)}{\rho RT}{f}^{eq}\mathrm{,}\end{array}

As a result, the simplified Enskog equation becomes

\begin{array}{}\text{(14)}& \begin{array}{l}\frac{\partial f}{\partial t}+v\cdot \nabla f-\frac{F\cdot (v-u)}{\rho RT}{f}^{eq}=\\ -\frac{1}{\tau }\left(f-f^{eq}\right)-\frac{(v-u)\cdot \nabla (b{\rho }^{2}RT\chi )}{\rho RT}{f}^{eq}\mathrm{.}\end{array}\end{array}

The last term on the right side of the equation represents the repulsion interaction of molecules.

2.2 Multiphase flow model based on Enskog equation

Now the repulsion interaction of molecules is introduced by using the Enskog collision operator instead of the Boltzmann collision operator, then the key to modeling multiphase flow is to incorporate the intermolecular attraction. According to the average filed approximation, the mutual attraction between molecules can be regarded as an average force potential [24], which reads

\begin{array}{}\text{(15)}& \text{Φ}(r_1)={\displaystyle {\int }_{r_{12}>d_0}}{\varphi }_{attr}(r_{12})\rho (r_2)\text{d}{r}_2\mathrm{,}\end{array}

where r_{12}=\left|r_1-r_2\right| is the distance between r₁ and r₂, {\varphi }_{attr}(r_{12}) denotes the attraction potential. Expanding \rho (r_2) at r₁ and ignoring terms higher than the second order, \text{Φ}(r_1) can be approximated by

\begin{array}{}\text{(16)}& \text{Φ}(r_1)=-2a\rho -K{\nabla }^2\rho \mathrm{,}\end{array}

where the first term mainly affects the equation of state (EOS) with a=-\frac{1}{2}{\displaystyle {\int }_{r>d_0}}{\varphi }_{attr}\left(r\right)\text{d}r and the second term contributes to the surface tension. The coefficient of surface tension is K=-\frac{1}{6}{\displaystyle {\int }_{r>d_0}}{r}^2{\varphi }_{attr}\left(r\right)\text{d}r, and it is assumed to be a constant in this work.

Based on the expression of \text{Φ}(r) in Eq. (16), the total external force of molecules at r can be calculated as

\begin{array}{}\text{(17)}& F=-\rho \nabla \text{Φ}=\nabla (a{\rho }^2)+\rho \nabla (K{\nabla }^2\rho )\mathrm{.}\end{array}

Substituting the expression of F into Eq. (14) gives

\begin{array}{}\text{(18)}& \frac{\partial f}{\partial t}+v\cdot \nabla f-\frac{F^{\prime }\cdot (v-u)}{\rho RT}{f}^{eq}=-\frac{1}{\tau }\left(f-f^{eq}\right)\mathrm{,}\end{array}

where F^{\prime }=-\nabla \psi +\rho \nabla \left(K{\nabla }^2\rho \right) with \psi =b{\rho }^{2}RT\chi -a{\rho }^2. There are two terms in F′, the first term is related to the EOS and the second one corresponds to the surface tension.

According to Chapman–Enskog multi-scale expansion, the hydrodynamic equations for multiphase flow can be obtained from Eq. (18) as follows:

\begin{array}{}\text{(19)}& \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)=\mathrm{0,}\end{array}

\begin{array}{}\text{(20)}& \frac{\partial (\rho u)}{\partial t}+\nabla \cdot (\rho uu+PI+\Pi )-\rho \nabla (K{\nabla }^2\rho )=\mathrm{0,}\end{array}

where P represents the EOS of real gas in the form of

\begin{array}{}\text{(21)}& P=\rho RT\left(1+b\rho \chi \right)-a{\rho }^2\mathrm{.}\end{array}

When \chi =\frac{1-b\rho /8}{{\left(1-b\rho /4\right)}^3}, the Carnahan–Starling EOS can be obtained as

\begin{array}{}\text{(22)}& P^c=\rho RT\frac{1+\eta +{\eta }^2-{\eta }^3}{{(1-\eta )}^3}-a{\rho }^2\mathrm{,}\end{array}

with \eta =\frac{b\rho }{4}. The van der Waals (VDW) EOS

\begin{array}{}\text{(23)}& P^v=\frac{\rho RT}{1-b\rho }-a{\rho }^2\mathrm{,}\end{array}

can also be derived from Eq. (21) when \chi =\frac{1}{1-b\rho }. The \Pi in Eq. (20) represents the viscous stress which has an expression as

\begin{array}{}\text{(24)}& \Pi =\mu \left[\nabla u+{\nabla }^{T}u-\frac{2}{D}(\nabla \cdot u)I\right]\mathrm{,}\end{array}

where D denotes spatial dimension and μ is the coefficient of dynamic viscosity, \mu =\tau \rho RT. The last term, \rho \nabla (K{\nabla }^2\rho ), on the left-hand side of Eq. (20) denotes the surface tension F_s, which is equivalent to the following expression [19]

\begin{array}{}\text{(25)}& F_{\text{s}}=\nabla \cdot \left[\left(\frac{K}{2}\nabla \rho \cdot \nabla \rho +K\rho {\nabla }^2\rho \right)I-K\nabla \rho \nabla \rho \right]\mathrm{.}\end{array}

2.3 Discrete Enksog model for multiphase flow

After the previous derivation, the kinetic equation of multiphase flow is obtained as Eq. (18). Correspondingly, the evolution of discrete Enskog equation for multiphase flow reads

\begin{array}{}\text{(26)}& \frac{\partial f_{ki}}{\partial t}+v_{ki}\cdot \nabla f_{ki}-\frac{(v_{ki}-u)\cdot F^{\prime }}{\rho RT}{f}_{ki}^{eq}=-\frac{1}{\tau }\left(f_{ki}-f_{ki}^{eq}\right)\mathrm{,}\end{array}

where f_{ki}=f_{ki}(r\mathrm{,}t) is the distribution function of the discrete velocity v_ki and the subscript “ki” indicates the index of discrete velocity. f_{ki}^{eq} is the discrete local equilibrium distribution function; its general expression of any order has been given in the previous literature [9]. In the previous Chapman–Enskog multi-scale expansion process mentioned above, only the 0th to 3rd kinetic moments are needed to recover Eqs. (19) and (20). Therefore, for computational efficiency, the discrete local equilibrium distribution function is approximated by the 3rd-order Hermite polynomial which reads

\begin{array}{}\text{(27)}& \begin{array}{ll}{f}_{ki}^{eq}=\hfill & \rho F_k[\left(1-\frac{u^2}{2T}\right)+\frac{1}{1!}\left(1-\frac{u^2}{2T}\right)\frac{v_{ki}\cdot u}{T}\hfill \\ \hfill & +\frac{1}{2!}\frac{{\left(v_{ki}\cdot u\right)}^2}{T^2}+\frac{1}{3!}\frac{{\left(v_{ki}\cdot u\right)}^3}{T^3}]\hfill \end{array}\end{array}

where F_k is the weight coefficient which can be expressed as

\begin{array}{}\text{(28)}& F_1=\frac{24T^3-4\left(c_2^2+c_3^2\right)T^2+c_2^2{c}_3^{2}T}{3c_1^2\left(c_1^2-c_2^2\right)\left(c_1^2-c_3^2\right)}\mathrm{,}\end{array}

\begin{array}{}\text{(29)}& F_2=\frac{24T^3-4\left(c_1^2+c_3^2\right)T^2+c_1^2{c}_3^{2}T}{3c_2^2\left(c_2^2-c_1^2\right)\left(c_2^2-c_3^2\right)}\mathrm{,}\end{array}

\begin{array}{}\text{(30)}& F_3=\frac{24T^3-4\left(c_2^2+c_1^2\right)T^2+c_2^2{c}_1^{2}T}{3c_3^2\left(c_3^2-c_1^2\right)\left(c_3^2-c_2^2\right)}\mathrm{,}\end{array}

and

\begin{array}{}\text{(31)}& F_0=1-6\left(F_1+F_2+F_3\right)\mathrm{.}\end{array}

where c₁, c₂, and c₃ are the values of three groups of discrete velocities. The scheme of the discrete velocity model is shown in Fig. 1. The value of c_k does not affect the final results as long as the calculation can be performed stably. It should be noted that the new model can be easily extended to higher-order cases only if the discrete local equilibrium distribution function satisfies more kinetic moment relations. Higher-order terms need to be kept in the Hermite polynomials of Eq. (27) to achieve this aim. Accordingly, more discrete velocities are required to satisfy the higher-order isotropy of the discrete model. More details of the higher-order model are available in Ref. [9].

Once the weight coefficients F_k are calculated, the discrete local equilibrium distribution function f_{ki}^{eq} in Eq. (27) is also known, then the discrete distribution function in t+\text{Δ}t(f_{ki}^{t+\text{Δ}t}) can be solved from f_{ki}^t by the following evolution equation

\begin{array}{}\text{(32)}& f_{ki}^{t+\text{Δ}t}=-\frac{1}{\tau }\left(f_{ki}^t-f_{ki}^{eq}\right)-v_{ki}\cdot \nabla f_{ki}^t+\frac{(v_{ki}-u)\cdot F^{\prime }}{\rho RT}{f}_{ki}^{eq}\mathrm{.}\end{array}

It involves the calculation of spatial derivatives in both the convection term \nabla f_{ki}^t and the external force term F′. Various types of different schemes, such as finite-difference, Nine-Point Stencils (NPS), Nonoscillatory Non-free-parameter and Dissipative (NND), and Fast Fourier Transform (FFT) schemes, etc. can be used depending on the specific problems.

In addition, from Eq. (32) we see that f_{ki} and its corresponding local equilibrium distribution function f_{ki}^{eq} are both known at a certain time. Consequently, the non-equilibrium quantities {\text{Δ}}_2^{*} and {\text{Δ}}_3^{*} can be calculated as

\begin{array}{}\text{(33)}& {\text{Δ}}_2^{*}={\displaystyle \sum _{ki}}\left(f_{ki}-f_{ki}^{eq}\right)\left(v_{ki}-u\right)\left(v_{ki}-u\right)\mathrm{,}\end{array}

\begin{array}{}\text{(34)}& {\text{Δ}}_3^{*}={\displaystyle \sum _{ki}}\left(f_{ki}-f_{ki}^{eq}\right)\left(v_{ki}-u\right)\left(v_{ki}-u\right)\left(v_{ki}-u\right)\mathrm{.}\end{array}

Correspondingly, the quantities of non-equilibrium strength, D_2^{*} and D_3^{*}, can be obtained, which are defined as

\begin{array}{}\text{(35)}& D_2^{*}=\sqrt{{\left|{\text{Δ}}_{\mathrm{2,}xx}^{*}\right|}^2+{\left|{\text{Δ}}_{\mathrm{2,}xy}^{*}\right|}^2+{\left|{\text{Δ}}_{\mathrm{2,}yy}^{*}\right|}^2}\end{array}

and

\begin{array}{}\text{(36)}& D_3^{*}=\sqrt{{\left|{\text{Δ}}_{\mathrm{3,}xxx}^{*}\right|}^2+{\left|{\text{Δ}}_{\mathrm{3,}xxy}^{*}\right|}^2+{\left|{\text{Δ}}_{\mathrm{3,}xyy}^{*}\right|}^2+{\left|{\text{Δ}}_{\mathrm{3,}yyy}^{*}\right|}^2}\mathrm{,}\end{array}

respectively. In previous studies [9], it has been found that some quantities of non-equilibrium strength may perform better than individual non-equilibrium components in describing some characteristics of fluid interfaces. In this work, the role of these non-equilibrium quantities in multiphase flow will be further investigated.

3 Simulation and discussion

3.1 Couette flow

As the first test, the Couette flow is simulated to examine the viscous stress calculated by the discrete model with the discrete velocity model shown in Fig. 1. Two parallel plates, filled with fluid between them, are placed in the y direction. The left plate is stationary while the right plate moves at a constant speed U_y. The fluid is driven by the plate due to the viscous stress. This flow process is simulated by using the discrete model without considering the non-ideal gas effects and surface tension.

The speed of the moving plate is U_y=0.2 and the temperature of the whole flow field is fixed at T=1.0. The computational meshes are N_x\times N_y=100\times 1 with spatial steps \text{d}x=\text{d}y=0.002 and time step \text{d}t=0.0001. Different viscosity coefficients are obtained by changing the relaxation time τ. The initial velocity in the flow filed is zero. The non-slip boundary condition is used on both left and right boundary. The first-order forward difference is used for time discretization and the NND scheme is used for spatial discretization. The simulation results are shown in Fig. 2.

Figure 2 (a) gives the profiles of velocity U_y along the x direction at several different times with a fixed relaxation time \tau =0.002, while Fig. 2 (b) shows profiles of U_y with several different relaxation time at t=1.0. The analytical solutions denoted by the solid lines are also plotted for comparison. The simulation results are in good agreement with the analytical solutions, which verifies the accuracy of the viscous stress calculated by the new discrete model.

3.2 Two-phase coexistence curve

The two-phase coexistence problem is one of the most typical multiphase flow problems that can be described by a single-component multiphase flow model. According to the VDW EOS, below the critical temperature T_c, the system allows for two coexisting fluid with different densities under the same pressure P. The high density ({\rho }_l) corresponds to the liquid phase, while the lower density ({\rho }_v) to the vapor phase. Seen from the EOS shown in Fig. 3 (a), there are many groups of coexistence points at different pressures, as long as the value of pressure is between A and B; the liquid phase corresponds to the AB segment while the vapor phase to the CD segment. However, only one set of liquid-vapor points satisfy both the mechanical equilibrium and thermodynamic equilibrium (more specifically, the chemical potential equilibrium) at a given temperature T₀. This set of liquid-vapor points can be calculated by the Maxwell equal area rule.

The liquid-vapor coexistence points at several different temperatures are calculated by using the new multiphase flow model. The computation meshes are N_x\times N_y=250\times 1 with the spatial steps \text{d}x=\text{d}y=4\times {10}^{-3}. The time step and the relaxation time are \text{d}t=5\times {10}^{-5} and \tau =0.02, respectively. The VDW EOS in Eq. (23) is adopted with a=\mathrm{9/8} and b=\mathrm{1/3}. The coefficient of surface tension is K=2\times {10}^{-5}. The first-order forward difference is used for time discretization. The NND scheme is adopted to calculate the spatial derivative in the convection term while the NPS scheme is used to calculate the spatial derivatives in the force term F′. Periodic boundary conditions are adopted on both left and right boundaries. The initial conditions are set as follows:

\begin{array}{}\text{(37)}& \{\begin{array}{l}{\left(\rho \mathrm{,}T\mathrm{,}{u}_x\mathrm{,}{u}_y\right)}_L=\left({\rho }_v\mathrm{,0.99,0,0}\right)\mathrm{,}\hfill \\ {\left(\rho \mathrm{,}T\mathrm{,}{u}_x\mathrm{,}{u}_y\right)}_M=\left({\rho }_l\mathrm{,0.99,0,0}\right)\mathrm{,}\hfill \\ {\left(\rho \mathrm{,}T\mathrm{,}{u}_x\mathrm{,}{u}_y\right)}_R=\left({\rho }_v\mathrm{,0.99,0,0}\right)\mathrm{,}\hfill \end{array}\end{array}

where the subscript “L”, “M”, and “R” indicate the regions , , and , respectively. {\rho }_v and {\rho }_l are the theoretical vapor and liquid densities at T=0.99. According to the VDW EOS combined with the equal area rule we get that {\rho }_v=0.8045 and {\rho }_l=1.2035.

The simulation continues until the equilibrium state is achieved, then the temperature drops by 0.01 and wait for the system to achieve another equilibrium state. The above process is repeated until the temperature drops to T=0.85. The vapor and liquid densities calculated by the new model at different temperatures are plotted in Fig. 3 (b). The symbols are simulated by the new multiphase flow model while the solid line represents the analytic solution. The simulation results are in excellent agreement with the theoretical solutions, which verified the accuracy of the new multiphase model in calculating the equilibrium phase transition.

To further validate the effect of surface tension, the transition curves between the liquid and vapor phases with different coefficients of surface tension, i.e., K=2\times {10}^{-5}, K=5\times {10}^{-5}, and K=1\times {10}^{-4}, at T=0.95 are calculated and the results are shown in Fig. 4. The theoretical solutions are also plotted for comparison. The theoretical solution of phase interfaces given by the modified VDW theory reads [61]

\begin{array}{}\text{(38)}& x-x_0=-\frac{1}{\sqrt{2a/K}}{\displaystyle {\int }_{{\rho }^{*}(x_0)}^{{\rho }^{*}(x)}}\frac{1}{\sqrt{{\text{Φ}}^{*}({\rho }^{*})-{\text{Φ}}^{*}({\rho }_l^{*})}}\text{d}{\rho }^{*}\mathrm{,}\end{array}

with

\begin{array}{}\text{(39)}& {\text{Φ}}^{*}({\rho }^{*})={\rho }^{*}\xi -{\rho }^{*}{T}^{*}\left[\ln (1/{\rho }^{*}-1)+1\right]-{({\rho }^{*})}^2\end{array}

and

\begin{array}{}\text{(40)}& \xi =T^{*}\ln (1/{\rho }_s^{*}-1)-\frac{{\rho }_s^{*}{T}^{*}}{1-{\rho }_s^{*}}+2{\rho }_s^{*}\mathrm{,}\end{array}

where {\rho }^{*} and T^{*} denotes the reduced density and temperature, {\rho }^{*}=b\rho and T^{*}=bT/a, {\rho }_s^{*} represents the reduced density of the liquid or vapor corresponding to the equilibrium state at T^{*}.

The theoretical solutions are shown as solid lines and the DBM results are represented by symbols. It shows that the density profiles of the interface simulated by the new multiphase model agree well with the theoretical solutions, which verifies the accuracy of the surface tension of the new model.

3.3 Droplet suspension

To further validate the effect of surface tension, a circular droplet suspended in its vapor is simulated. According to the Laplace law, under a fixed surface tension coefficient, the pressure difference \text{Δ}P inside and outside the droplet is inversely proportional to the radius of the droplet R₀, i.e., \text{Δ}P\sim 1/R_0.

Several droplets with various radii are simulated under two different coefficients of surface tension K=2\times {10}^{-5} and K=2\times {10}^{-4}, respectively. The computational meshes are N_x\times N_y=128\times 128 with the spatial steps \text{d}x=\text{d}y=0.02 and time step dt=0.00001. Periodic boundary conditions are adopted on the boundaries and corners. The first-order forward difference is used for time discretization. The FFT scheme is used to calculate the space derivative in both the convection term and the force term. As a representative, the contour maps of density and pressure for a droplet with R_0=0.5 are shown in Figs. 5(a) and (b), respectively.

The pressure difference \text{Δ}P inside and outside the circular droplet is detected during the calculation. The profile of \text{Δ}P with time is shown in Fig. 6(a) and it eventually tends to a stable value. The calculation continues until \text{Δ}P is almost unchanged, then the pressure difference \text{Δ}P and the corresponding radius R₀ is recorded. Several circular droplets with different radii are simulated. Figure 6(b) shows the relation between the \text{Δ}P and 1/R_0, in which symbols are simulation results and the solid lines are linear fitting. Clearly, \text{Δ}P is proportional to 1/R_0, which agrees well with the Laplace law. The accuracy of the surface tension calculated by the new multiphase model is well verified from this simulation.

3.4 Phase separation

In this part, a dynamic phase separation process is simulated by using the new multiphase flow model. The initial conditions are set as

\begin{array}{}\text{(41)}& (\rho \mathrm{,}T\mathrm{,}{u}_x\mathrm{,}{u}_y)=(1.0+\delta \mathrm{,0.92,0,0})\mathrm{,}\end{array}

where δ is a random noise with amplitude 0.01, which provides a starting point for phase separation. In the phase separation process, the temperature is fixed at T=0.92, which is equivalent to that the system contacts with a large heat source with a temperature T=0.92. The computational meshes are N_x\times N_y=300\times 300 with spatial steps \text{Δ}x=\text{Δ}y=0.01. The time step and relaxation time are \text{Δ}t=1\times {10}^{-4} and \tau =0.01, respectively. Periodic boundary conditions are adopted on all boundaries and corners. The first-order forward difference is used for time discretization. The NND and NPS schemes are used to calculate the spatial derivatives in the convection term and the force term, respectively. Figure 7 shows the density contour maps at several different times. The simulation results are consistent with those in previous literature [20, 40, 57].

3.5 Droplets collision

Droplet collision is another typical multiphase flow problem and has a broad application background in industrial production. It involves complex interfacial dynamics, including interface movement, fusion, and separation, etc. In this part, we simulate the head-on collision between two droplets using the new multiphase flow model. Two droplets with the opposite speeds are suspended in the same horizontal position. The droplet on the left moves to the right at a speed U₀, while the droplet on the right moves to the left at the same speed. There are two situations after the collision depending on the initial speed U₀. If U₀ is much small, the two droplets eventually fuse together, whereas if U₀ is large enough, they separate again along the vertical direction.

Two cases with U_0=0.2 and U_0=0.5, respectively, are simulated by using the new multiphase flow model. The computational meshes are N_x\times N_y=60\times 120 with the spatial step \text{Δ}x=\text{Δ}y=0.02. The time step is \text{Δ}t=1\times {10}^{-4} and the relaxation time is \tau =2\times {10}^{-3}. The coefficient of surface tension is set as K=2\times {10}^{-4} and the temperature is fixed at T=0.92. The boundary conditions and the discrete schemes are the same as the previous simulation of phase separation. The snapshots of the collision processes for two cases are shown in Figs. 8 and 9, respectively. Fig. 8 shows the case with U_0=0.2, in which two droplets fuse together after the collision and figure 9 shows the case with U_0=0.5 where they separate after the collision. The interfacial dynamics in droplet collision, fusion, and separation are well described, and the simulation results are consistent with the previous studies [62].

The evolutionary characteristics of non-equilibrium strength quantities, D_2^{*} and D_3^{*}, are shown in Fig. 10, in which the left and right subgraphs correspond to the two kinds of droplets collision with U_0=0.2 and U_0=0.5, respectively. The spatial mean values {\overline{D}}_2^{*} and {\overline{D}}_3^{*} are used to represent the total mean non-equilibrium strength, i.e.

\begin{array}{}\text{(42)}& {\overline{D}}_m^{\text{*}}=\frac{1}{l_x{l}_y}{\displaystyle \iint }{D}_m^{\text{*}}\text{d}x\text{d}y\mathrm{,}m=\mathrm{2,3,}\end{array}

where l_x{l}_y represents the total area of the simulation region.

Figure 10(a) shows that at least five stages of the collision can be identified from the non-equilibrium strength {\overline{D}}_2^{*}. In the first stage, two droplets come close to each other, and {\overline{D}}_2^{*} gradually decreases until it reaches point “A” in the {\overline{D}}_2^{*}(t) curve. In the second stage (AB segment on the profile of {\overline{D}}_2^{*}), two droplets contact and the interface begins to fuse, and {\overline{D}}_2^{*} increases with time. After the phase interface fusion is completed, the third stage (BC segment) begins, and the merged droplet elongate vertically, accompanied by a decrease of {\overline{D}}_2^{*}. Because the kinetic energy after collision is not enough to overcome the constraint of surface tension, the merged droplet can not be separated again. As a result, in the fourth stage (CD segment), the droplet contracts vertically and the overall trend of the non-equilibrium strength is upward. In the fifth stage (after the “D" point in the {\overline{D}}_2^{*}(t) curve), the droplet lengthens and contracts alternately in the horizontal and vertical directions, and finally tends to a stable state. Correspondingly, the {\overline{D}}_2^{*} shows the trend of oscillation attenuation and finally tends to a stable value. From the profile of {\overline{D}}_3^{*}, only the first three stages can be well identified, which indicates that the roles of non-equilibrium strength of different orders are much different.

For the second type of droplet collision with U_0=0.5, there are also five stages as shown in Fig. 10(b). The profile of {\overline{D}}_2^{*} in the first four stages is similar to the case with U_0=0.2. However, in the last stage, {\overline{D}}_2^{*} continues to decrease until the separated droplets leave the simulation region at t=5. The profiles after t=5 is meaningless because the droplet has left the simulation region. So these two types of collisions can be well distinguished by the last stage of {\overline{D}}_2^{*}. Likewise, from the profile of {\overline{D}}_3^{*}, only the first three stages can be well identified. As a result, it is difficult to distinguish these two types of collision cases from the profile of {\overline{D}}_3^{*}.

The evolution of the non-equilibrium strength with time can be explained by the change of the velocity gradient in the flow field. Because {\text{Δ}}_2^{*} corresponds to the viscous stress \Pi in the hydrodynamic equations, which can be expressed as Eq. (24) in the NS level. In this approximation, D_2^{*} can be expressed as

\begin{array}{}\text{(43)}& D_2^{*}\approx {\mu }^{\mathrm{1/2}}{\left(\nabla u:\nabla u+2\frac{\partial u_x}{\partial y}\frac{\partial u_y}{\partial x}\right)}^{\mathrm{1/2}}\mathrm{.}\end{array}

We further find that the last two terms on Eq. (43) have little effect on the overall evolutionary characteristics. For the sake of convenience, we will only use {\left(\nabla u:\nabla u\right)}^{\mathrm{1/2}} in the following analysis. Figure 11 shows the profiles of {\overline{D}}_2^{*} and the spatial average of {\left(\nabla u:\nabla u\right)}^{\mathrm{1/2}}. The definition of \overline{{\left(\nabla u:\nabla u\right)}^{\mathrm{1/2}}} is the same as that of {\overline{D}}_m^{*} in Eq. (42). Figures 11(a) and (b) correspond to the first and second types of collision processes, respectively. It shows that {\overline{D}}_2^{*} and \overline{{\left(\nabla u:\nabla u\right)}^{\mathrm{1/2}}} have similar evolutionary characteristics in both Figs. 11(a) and (b), which is consistent with our previous analysis.

For the first type of droplet collision, the collision process is roughly divided into five stages as shown in Fig. 11(a). In the first stage, the two droplets approach each other independently. Part of the kinetic energy of the droplets is transferred to the flow field of vapor phase, leading to a decrease of the system’s velocity gradient whole. In the second stage, two droplets collide and fuse. During the fusion process, the interface contracts rapidly due to the work done by surface tension, which results in the increase of velocity gradient in the flow field. After the droplets are completely fused, the new droplet begins to stretch in the y direction in the third stage. In this process, the surface tension does the negative work, causing the interface to slow down and the velocity gradient to decrease. In the fourth stage, because the kinetic energy of the interface motion is not enough to overcome the constraint of surface tension, the interface begins to contract in reverse when it stretches to the longest in the y direction. The surface tension does positive work, so the interface accelerates and the velocity gradient increases. Soon the interface shrinks to its minimum in the y direction and begins to stretch again in the x direction. The velocity gradient decreases in the process of stretching and increases in the process of shrinking. Finally, in the fifth stage, the droplet is alternately stretched and contracted in the x and y direction, but the maximum position it can reach gradually shortens until it reaches a stable state.

For the second type of droplet collision, the first stages are the same as those of the first type of collision. However, in the fourth stage, the fused droplet rupture and separate again. The interface of the two separated droplets shrinks rapidly, leading to an increase of velocity gradient. After that, in the fifth stage, the single droplet decelerates independently, so the velocity gradient decreases until the droplets leave the simulation area. {\text{Δ}}_3^{*} is a higher-order non-equilibrium characteristic quantity related to heat flux but its corresponding macroscopic physical significance is not yet clear. Because the temperature is set as a constant in the simulation, the heat flow in the flow field is sharply limited, so the amplitude is relatively small and its evolution characteristics are not visible.

4 Conclusion

In this work, under the framework of the discrete Boltzmann method, a new kinetic multiphase flow model based on a discrete Enskog equation was proposed. Compared with the previous multiphase flow discrete Boltzmann model, a bottom-up modeling approach was adopted in the new model. The repulsion potential was introduced by using the Enskog collision operator instead of the Boltzmann collision operator, and the attraction potential is incorporated through average filed approximation. The effect of total intermolecular forces is ultimately taken into account via the force term of the DBM. The discrete local equilibrium distribution function is solved by the Hermite polynomials. The third order Hermite polynomial is adopted in order to improve the computational efficiency as much as possible, although it is straightforward to take the higher-order ones. Several benchmarks, including the Couette flow, two-phase coexistence curve, droplet suspension, phase separation, and droplet collisions are simulated, through which the accuracy of the new multiphase flow model is verified. In addition, for two kinds of droplet collisions, the non-equilibrium characteristics are comparatively investigated via two TNE strength quantities, {\overline{D}}_2^{*} and {\overline{D}}_3^{*}. It is found that during the collision process, {\overline{D}}_2^{*} is always significantly larger than {\overline{D}}_3^{*}, {\overline{D}}_2^{*} can be used to identify the different stages of the collision process and to distinguish different types of collisions.

References

[1] Y. Chen, Q. Xie, A. Sari, P. V. Bardy, and A. Saeedi, Oil/water/rock wettability: Influencing factors and implications for low salinity water flooding in carbonate reservoirs, Fuel 215, 171 (2018)

[2] Y. Chen and Z. Deng, Hydrodynamics of a droplet passing through a microfluidic T-junction, J. Fluid Mech. 819, 401 (2017)

[3] J. Tice, H. Song, A. Lyon, and R. Ismagilov, Formation of droplets and mixing in multiphase microfluidics at low values of the reynolds and the capillary numbers, Langmuir 19(22), 9127 (2003)

[4] A. Günther and K. Jensen, Multiphase microfluidics: From flow characteristics to chemical and materials synthesis, Lab Chip 6(12), 1487 (2006)

[5] E. Christopher, Brennen, Fundamentals of Multiphase Flow, Cambridge: Cambridge University Press, 2005

[6] R. Saurel and C. Pantano, Diffuse-interface capturing methods for compressible two-phase flows, Annu. Rev. Fluid Mech. 50(1), 105 (2018)

[7] A. Frezzotti, P. Barbante, and L. Gibelli, Direct simulation Monte Carlo applications to liquid–vapor flows, Phys. Fluids 31(6), 062103 (2019)

[8] M. Wörner, Numerical modeling of multiphase flows in microfluidics and micro process engineering: A review of methods and applications, Microfluid. Nanofluidics 12(6), 841 (2012)

[9] Y. Zhang, A. Xu, G. Zhang, Z. Chen, and P. Wang, Discrete Boltzmann method for non-equilibrium flows: Based on Shakhov model, Comput. Phys. Commun. 238, 50 (2019)

[10] M. Moseler and U. Landman, Formation, stability, and breakup of nanojets, Science 289(5482), 1165 (2000)

[11] S. Zhan, Y. Su, Z. Jin, M. Zhang, W. Wang, Y. Hao, and L. Li, Study of liquid-liquid two-phase flow in hydrophilic nanochannels by molecular simulations and theoretical modeling, Chem. Eng. J. 395, 125053 (2020)

[12] S. Wolfram, Cellular automaton fluids 1: Basic theory, J. Stat. Phys. 45(3–4), 471 (1986)

[13] S. Chen and G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30(1), 329 (1998)

[14] S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford: Oxford University Press, 2001

[15] X. He and G. D. Doolen, Thermodynamic foundations of kinetic theory and lattice Boltzmann models for multiphase flows, J. Stat. Phys. 107(1–2), 309 (2002)

[16] R. Qin, Mesoscopic interparticle potentials in the lattice Boltzmann equation for multiphase fluids, Phys. Rev. E 73(6), 066703 (2006)

[17] Q. Li, K. Luo, Q. Kang, Y. He, Q. Chen, and Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Pror. Energy Combust. Sci. 52, 62 (2016)

[18] R. Qin, Thermodynamic properties of phase separation in shear flow, Comput. Fluids 117, 11 (2015)

[19] K. Timm, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. Viggen, The Lattice Boltzmann Method - Principles and Practice, Springer, 2017

[20] D. Grunau, S. Chen, and K. Eggert, A lattice Boltzmann model for multiphase fluid flows, Phys. Fluids 5(10), 2557 (1993)

[21] X. Shan and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47(3), 1815 (1993)

[22] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Lattice Boltzmann simulation of non-ideal fluids, Phys. Rev. Lett. 75(5), 830 (1995)

[23] A. Xu, G. Gonnella, and A. Lamura, Phase-separating binary fluids under oscillatory shear, Phys. Rev. E 67(5), 056105 (2003)

[24] X. He, S. Chen, and R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability, J. Comput. Phys. 152(2), 642 (1999)

[25] H. Liang, Q. Li, B. Shi, and Z. Chai, Lattice Boltzmann simulation of three-dimensional Rayleigh-Taylor instability, Phys. Rev. E 93(3), 033113 (2016)

[26] H. Wang, X. Yuan, H. Liang, Z. Chai, and B. Shi, A brief review of the phase-field-based lattice Boltzmann method for multiphase flows, Capillarity 2(3), 33 (2019)

[27] D. Sun, A discrete kinetic scheme to model anisotropic liquid–solid phase transitions, Appl. Math. Lett. 103, 106222 (2020)

[28] M. Watari and M. Tsutahara, Two-dimensional thermal model of the finite-di ference lattice Boltzmann method with high spatial isotropy, Phys. Rev. E 67(3), 036306 (2003)

[29] G. Gonnella, A. Lamura, and V. Sofonea, Lattice Boltzmann simulation of thermal non-ideal fluids, Phys. Rev. E 76(3), 036703 (2007)

[30] A. Onuki, Dynamic van der Waals theory of two-phase fluids in heat flow, Phys. Rev. Lett. 94(5), 054501 (2005)

[31] Y. Gan, A. Xu, G. Zhang, and Y. Li, FFT-LB modeling of thermal liquid-vapor system, Commum. Theor. Phys. 57(4), 681 (2012)

[32] Y. Gan, A. Xu, G. Zhang, Y. Li, and H. Li, Phase separation in thermal systems: A lattice Boltzmann study and morphological characterization, Phys. Rev. E 84(4), 046715 (2011)

[33] A. Xu, G. Zhang, Y. Gan, F. Chen, and X. Yu, Lattice Boltzmann modeling and simulation of compressible flows, Front. Phys. 7(5), 582 (2012)

[34] A. Xu, G. Zhang, and Y. Ying, Progess of discrete Boltzmann modeling and simulation of combustion system, Acta Physica Sinica 64, 184701 (2015)

[35] A. Xu, G. Zhang, and Y. Gan, Progress in studies on discrete Boltzmann modeling of phase separation process, Mech. Eng. 38, 361 (2016)

[36] A. Xu, G. Zhang, and Y. Zhang, Discrete Boltzmann modeling of compressible flows, in: G. Z. Kyzas and A.C. Mitropoulos (Eds.), Kinetic Theory, InTech, Rijeka, 2018, Ch. 02

[37] C. Lin and K. Luo, Discrete Boltzmann modeling of unsteady reactive flows with nonequilibrium effects, Phys. Rev. E 99(1), 012142 (2019)

[38] Y. Gan, A. Xu, G. Zhang, Y. Zhang, and S. Succi, Discrete Boltzmann trans-scale modeling of high-speed compressible flows, Phys. Rev. E 97(5), 053312 (2018)

[39] F Y. Gan, A. Xu, G. Zhang, and S. Succi, Discrete Boltzmann modeling of multiphase flows: hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11(26), 5336 (2015)

[40] H Y. Zhang, A. Xu, G. Zhang, Y. Gan, Z. Chen, and S. Succi, Entropy production in thermal phase separation: A kinetic-theory approach, Soft Matter 15(10), 2245 (2019)

[41] B. Yan, A. Xu, G. Zhang, Y. Ying, and H. Li, Lattice Boltzmann model for combustion and detonation, Front. Phys. 8(1), 94 (2013)

[42] A. Xu, C. Lin, G. Zhang, and Y. Li, Multiple-relaxation-time lattice Boltzmann kinetic model for combustion, Phys. Rev. E 91(4), 043306 (2015)

[43] C. Lin, A. Xu, G. Zhang, and Y. Li, Double-distribution-function discrete Boltzmann model for combustion, Combust. Flame 164, 137 (2016)

[44] Y. Zhang, A. Xu, G. Zhang, C. Zhu, and C. Lin, Kinetic modeling of detonation and effects of negative temperature coefficient, Combust. Flame 173, 483 (2016)

[45] C. Lin and K. Luo, MRT discrete Boltzmann method for compressible exothermic reactive flows, Comput. Fluids 166, 176 (2018)

[46] C. Lin, K. Luo, L. Fei, and S. Succi, A multi-component discrete Boltzmann model for nonequilibrium reactive flows, Sci. Rep. 7(1), 14580 (2017)

[47] A. Xu, G. Zhang, Y. Zhang, P. Wang, and Y. Ying, Discrete Boltzmann model for implosion and explosion related compressible ow with spherical symmetry, Front. Phys. 13(5), 135102 (2018)

[48] H. Lai, A. Xu, G. Zhang, Y. Gan, Y. Ying, and S. Succi, Non-equilibrium thermohydrodynamic effects on the Rayleigh-Taylor instability incompressible flow, Phys. Rev. E 94(2), 023106 (2016)

[49] F. Chen, A. Xu, and G. Zhang, Viscosity, heat conductivity, and Prandtl number effects in the Rayleigh Taylor instability, Front. Phys. 11(6), 114703 (2016)

[50] H. Ye, H. Lai, D. Li, Y. Gan, C. Lin, L. Chen, and A. Xu, Knudsen number effects on two-dimensional Rayleigh-Taylor instability in compressible fluid: Based on a discrete Boltzmann method, Entropy (Basel) 22(5), 500 (2020)

[51] Y. Gan, A. Xu, G. Zhang, C. Lin, H. Lai, and Z. Liu, Nonequilibrium and morphological characterizations of Kelvin–Helmholtz instability in compressible flows, Front. Phys. 14(4), 43602 (2019)

[52] C. Lin, A. Xu, G. Zhang, K. Luo, and Y. Li, Discrete Boltzmann modeling of Rayleigh–Taylor instability in two-component compressible flows, Phys. Rev. E 96(5), 053305 (2017)

[53] H. Liu, W. Kang, Q. Zhang, Y. Zhang, H. Duan, and X. He, Molecular dynamics simulations of microscopic structure of ultra strong shock waves in dense helium, Front. Phys. 11(6), 115206 (2016)

[54] H. Liu, Y. Zhang, W. Kang, P. Zhang, H. Duan, and X. He, Molecular dynamics simulation of strong shock waves propagating in dense deuterium, taking into consideration effects of excited electrons, Phys. Rev. E 95(2), 023201 (2017)

[55] H. Liu, W. Kang, H. Duan, P. Zhang, and X. He, Recent progresses on numerical investigations of microscopic structure of strong shock waves in fluid, Sci. China Phys. Mech. Astron. 47(7), 070003 (2017)

[56] J. Meng, Y. Zhang, N. Hadjiconstantinou, G. Radtke, and X. Shan, Lattice ellipsoidal statistical BGK model for thermal non-equilibrium flows, J. Fluid Mech. 718, 347 (2013)

[57] F Y. Gan, A. Xu, G. Zhang, and S. Succi, Discrete Boltzmann modeling of multiphase flows: Hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11(26), 5336 (2015)

[58] Q. Shen, Rarefied Gas Dynamics: Fundamentals, Simulations and Micro Flows, Springer, 2005

[59] S. Chapman, T. Cowling, and D. Burnett, The Mathematical Theory of Non-Uniform Gases: An Account of The Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, Cambridge: Cambridge University Press, 1990

[60] Z. Guo and C. Zheng, Theory and Applications of Lattice Boltzmann Method, Beijing: Science Press, 2008

[61] V. Bongiorno and H. T. Davis, Modified van der Waals theory of fluid interfaces, Phys. Rev. A 12(5), 2213 (1975)

[62] H. Huang, M. Sukop, and X. Lu, Multiphase Lattice Boltzmann Methods: Theory and Application, John Wiley & Sons, Inc, 2015

Rights & permissions

Higher Education Press 2020