Free-Energy-Based Discrete Unified Gas Kinetic Scheme for van der Waals Fluid

Abstract

The multiphase model based on free-energy theory has been experiencing long-term prosperity for its solid foundation and succinct implementation. To identify the main hindrance to developing a free-energy-based discrete unified gas-kinetic scheme (DUGKS), we introduced the classical lattice Boltzmann free-energy model into the DUGKS implemented with different flux reconstruction schemes. It is found that the force imbalance amplified by the reconstruction errors prevents the direct application of the free-energy model to the DUGKS. By coupling the well-balanced free-energy model with the DUGKS, the influences of the amplified force imbalance are entirely removed. Comparative results demonstrated a consistent performance of the well-balanced DUGKS despite the reconstruction schemes utilized. The capability of the DUGKS coupled with the well-balanced free-energy model was quantitatively validated by the coexisting density curves and Laplace’s law. In the quiescent droplet test, the magnitude of spurious currents is reduced to a machine accuracy of ${10}^{-15}$. Aside from the excellent performance of the well-balanced DUGKS in predicting steady-state multiphase flows, the spinodal decomposition test and the droplet coalescence test revealed its stability problems in dealing with transient flows. Further improvements are required on this point.

1. Introduction

Multiphase fluid flow characterized by the concurrent presence of multiple thermodynamic phases is frequently encountered in industrial processes and engineering applications. Insightful understanding of the multiphase flow behavior could facilitate improvements in manufacturing technology and production efficiency. Due to the restriction on measurement technology and the experimental platform, it is particularly challenging to reveal the flow details by experimental methods. Benefiting from the substantial improvements in computing power, numerical simulation technology has been developed into a powerful tool for the study of complicated behaviors arising in multiphase fluid flow. By numerically solving the set of interface capturing and hydrodynamic equations, a multitude of research studies [1,2,3,4] vividly detail the interface dynamics and flow structures from a macroscopic perspective. Essentially, the interfacial phenomenon represents the macroscopic manifestation of the microscopic interactions among fluid molecules [5]. Numerical methods based on realistic microscopic physics could offer in-depth findings regarding multiphase phenomena, but the heavy computational requirement of such methods for industry-scale multiphase problems is far beyond affordable. In recent years, numerical schemes constructed with the mesoscopic theory [6] have been emerging as a compelling methodology for resolving multiphase flow patterns as this bridges the gap between the macroscopic descriptions of multiphase dynamics and microscopic intermolecular interactions and, thus, generates insightful understandings at an affordable cost.

Among various previously proposed mesoscopic approaches [7,8,9], the lattice Boltzmann (LB) method [7] has received particular attention for its concise and intuitive way of representing intermolecular interactions. Generally, the lattice Boltzmann multiphase models developed in the past few decades can be categorized into four classifications: the color-gradient model [10], the phase-field model [11,12], the pseudopotential model [13], and the free-energy model [14]. The phase-field model employs independent sets of distribution functions to separately transfer mass and momentum, which could cause mass non-conservation problems near the interface region [15]. The pseudopotential model and the free-energy model employ a single set of distribution functions to ensure a coherent transport of mass and momentum, which conforms to the physical reality that mass and momentum are simultaneously transferred by the unique molecules. Compared to the pseudopotential model, where interactions are built heuristically, the free-energy model is constructed upon the stationary-action principle, which possesses a firm physical background. Over the last couple of decades, the free-energy lattice Boltzmann method has been successfully applied to numerically tackle a variety of flow issues including the contact line movement [16,17], multicomponent fluids’ flow [18,19], wetting boundaries [20,21], and large-density-ratio fluid flow [22,23]. The primitive free-energy multiphase model proposed by Swift et al. [14] reflects the interaction effects via a modified equilibrium distribution function, whose second-order moment incorporates a nonideal thermodynamic pressure tensor. However, this primitive model suffers from a lack of Galilean invariance due to the superfluous terms recovered in the momentum equation. Later, Swift et al. [24] tried to remedy this defect by introducing additional terms to the pressure tensor, but an analysis through the Chapman–Enskog expansion demonstrated that the lack of Galilean invariance cannot be entirely eliminated. Based on Swift et al.’s work, Inamuro et al. [25] proposed a Galilean-invariant free-energy model with the guidance of asymptotic theory. Kalarakis et al. [26] restored the Galilean invariance of the free-energy model to second-order accuracy by modifying the zero-order momentum flux tensor. Wagner and Li [27] replaced the contribution of the nonideal pressure tensor with a corrected force term and improved the Galilean invariance of the model in large velocity situations. Meanwhile, Lee and Fischer [28] reformulated the pressure form of the interaction force into a potential form and reduced the magnitude of the spurious velocity to a machine level, at the cost of including the information in next-nearest-neighbor cells. Subsequently, Guo et al. [5] spotted that the spurious velocity originates from the force imbalance at the discrete level. Based on this finding, Lou and Guo [29] applied the Lax–Wendroff scheme to the lattice Boltzmann free-energy model and successfully mitigated the effects of the force imbalance. Very recently, Guo [30] proposed a well-balanced lattice Boltzmann scheme with which the spurious velocity can be ultimately minimized to the machine accuracy. The previously mentioned improvements were carried out within the framework of the lattice Boltzmann method, which inherits its advantages such as great simplicity and high efficiency. However, the uniformity requirement on the grid types posed by the LB method prevents its application in industrial cases.

Developed in the framework of the finite volume method, the discrete unified gas-kinetic scheme (DUGKS) [31] suffers no restriction in terms of the grid types. With the information of the Knudsen number incorporated in the construction of the interface flux, the DUGKS exhibits the capability of properly modeling a wide range of fluid flows ranging from the continuum regime to the free-molecule regime [32]. Over the past decade, the DUGKS has proven its excellent performance in predicting microscale gas flows [33,34], multicomponent gas flows [35,36], turbulent flows [37,38,39], compressible flows [40,41,42], radiative heat transfer [43,44], and so forth [45]. A comparative study [46] has demonstrated the stability superiority of the DUGKS over that of the LB method in terms of nearly incompressible flows. However, the DUGKS studies centered on multiphase fluid flows remain limited [47,48] and the multiphase DUGKS has been primarily confined to the phase-field model [49]. Although Yang et al. [50] developed a pseudopotential-based DUGKS for binary fluid flow, a free parameter is needed to guarantee the isotropic property of the fluid interface. Inspired by the well-balanced LB scheme [30], Zeng et al. [51] proposed a well-balanced DUGKS for two-phase fluid flows using the free-energy model. Comparative results demonstrated the superior performance of the DUGKS over that of the LB method. Nevertheless, there is still a lack of an insightful comprehension as to the isotropic property of free-energy-based DUGKS. In this research, we elucidate the mechanism for the nonisotropic phenomena produced by the free-energy-based DUGKS using different reconstruction approaches. Then, we couple the well-balanced free-energy model with the DUGKS implemented with different reconstruction schemes to investigate practical van der Waals (vdW) fluid flows. The rest of this paper is organized as follows. In Section 2, the primitive and the well-balanced free-energy models are introduced, followed by the detailed explanation of the Strang-splitting DUGKS. The comparative numerical results, as well as brief discussions are presented in Section 3. Finally, a summary is given in Section 4.

2. Numerical Methodology

In this section, the first part theoretically introduces the free-energy model based on the vdW chemical potential and the second part exhaustively explains the Strang-splitting DUGKS implemented with different reconstruction schemes.

2.1. Free-Energy Model

Considering a multiphase system, the free-energy functional in terms of the fluid density $\rho$ can be expressed as [14,24]

$$\mathcal{F}=\int \varphi (\rho ,\nabla \rho )d{\mathsf{\Omega }}_V=\int \left(E_f\left(\rho \right)+\frac{\kappa }{2}{|\nabla \rho |}^2\right)d{\mathsf{\Omega }}_V,$$

(1)

where ${\mathsf{\Omega }}_V$ is the spatial region occupied by the system, $\varphi (\rho ,\nabla \rho )$ denotes the total free-energy density, in which $E_f\left(\rho \right)$ represents the bulk free-energy density, and $\frac{\kappa }{2}{|\nabla \rho |}^2$ signifies the interface free-energy density. The parameter $\kappa$ is a positive constant determined by the interface thickness and the surface tension coefficient. Minimization of the free-energy $\mathcal{F}$ that is subject to the constraint of a constant mass $\mathcal{M}$ evolves the system towards the equilibrium condition, where

$$\mathcal{M}=\int \rho d{\mathsf{\Omega }}_v.$$

(2)

To impose the mass constraint, a transformed free-energy functional $\mathcal{L}$ is constructed using the method of Lagrange multipliers:

$$\mathcal{L}=\mathcal{F}-\lambda \mathcal{M},$$

(3)

where $\lambda$ is the Lagrange multiplier. Minimization of the constrained free-energy demands the corresponding first variation to be zero:

$$\delta \mathcal{L}=0,$$

(4)

which yields the following Euler–Lagrange equation:

$$\frac{\partial \psi }{\partial \rho }-\nabla ·\left(\frac{\partial \psi }{\partial (\nabla \rho )}\right)=\frac{d{E}_f}{d\rho }-\kappa {\nabla }^2\rho -\lambda =0,$$

(5)

where

$$\psi (\rho ,\nabla \rho )=\varphi (\rho ,\nabla \rho )-\lambda \rho .$$

(6)

The chemical potential ${\mu }_c$ is defined as the variation of the free-energy $\mathcal{F}$ with respect to the density [52]:

$${\mu }_c=\frac{\delta \mathcal{F}}{\delta \rho }=\frac{d{E}_f}{d\rho }-\kappa {\nabla }^2\rho .$$

(7)

As the integrand of transformed free-energy $\mathcal{L}$ does not explicitly contain any spatial coordinates, it remains invariant regardless of the spatial translations [3]. Noether’s theorem [53] says that the invariance of the free-energy with respect to spatial translations corresponds to a conserved tensorial current $\mathit{J}$ satisfying [54]:

$$\nabla ·\mathit{J}=0,$$

(8)

where $\mathit{J}$ is a second-rank tensor given by

$$\mathit{J}=-\psi \mathit{I}+\nabla \rho \otimes \frac{\partial \psi (\rho ,\nabla \rho )}{\partial \left(\nabla \rho \right)},$$

(9)

in which $\mathit{I}$ is the identity matrix. Substituting Equations (5) and (6) into Equation (9) leads to

$$\mathit{J}=\left(\rho {\mu }_c-E_f-\frac{\kappa }{2}{|\nabla \rho |}^2\right)\mathit{I}+\kappa \nabla \rho \nabla \rho .$$

(10)

The bulk pressure $p_b$ is connected to the bulk free-energy density $E_f$ via the Legendre transform [54]:

$$p_b\left(\rho \right)=\rho \frac{d{E}_f}{d\rho }-E_f\left(\rho \right),$$

(11)

with which the conserved current tenor $\mathit{J}$ can be identified as the thermodynamic pressure tensor $\mathit{P}$ in such a way that

$$\mathit{P}\equiv \mathit{J}=\left(p_b-\kappa \rho {\nabla }^2\rho -\frac{\kappa }{2}{|\nabla \rho |}^2\right)\mathit{I}+\kappa \nabla \rho \nabla \rho .$$

(12)

With some basic algebraic manipulations, the divergence of the pressure tensor can be simplified as

$$\nabla ·\mathit{P}=\rho \nabla {\mu }_c.$$

(13)

In the traditional free-energy model [28], the total effects of excess pressure accounting for the phase interactions can be represented by the following interaction force

$$\mathit{F}=\nabla ·{\mathit{P}}_0-\nabla ·\mathit{P}=\nabla p_0-\rho \nabla {\mu }_c,$$

(14)

where ${\mathit{P}}_0=p_0\mathit{I}$ denotes the pressure tensor of an ideal gas. In the well-balanced free-energy model [30], the interaction force is defined as

$$\mathit{F}=-\rho \nabla {\mu }_c$$

(15)

in order to eliminate the force imbalance at the discrete level.

The only remaining task is to determine the bulk free-energy density $E_f$. In the work of Zeng et al. [51], $E_f$ takes a double-well form, which relates to no specific equation of state (EOS). In the current research, the bulk pressure is evaluated by the nonideal van der Waals EOS [55] expressed as

$$p_b=\frac{\rho RT}{1-b\rho }-a{\rho }^2,$$

(16)

where parameter a denotes the intermolecular interaction strength, parameter b indicates the volume correction, R stands for the gas constant, and T represents the temperature. The corresponding bulk free-energy density can be obtained by solving Equation (11):

$$E_f\left(\rho \right)=\rho RT\ln \left(\frac{\rho }{1-b\rho }\right)-a{\rho }^2.$$

(17)

The chemical potential can then be obtained according to Equation (7):

$${\mu }_c=RT\left[\ln \left(\frac{\rho }{1-b\rho }+\frac{1}{1-b\rho }\right)\right]-2a\rho -\kappa {\nabla }^2\rho ,$$

(18)

with which the interaction force $\mathit{F}$ can be evaluated. In the current research, the parameters in the vdW-EOS were set as [56] $a=9/392,b=2/21,R=1$. $\kappa$ was fixed at 0.02 if not otherwise specified. The critical density and temperature are given as ${\rho }_c=3.5$ and $T_c=1/14$.

2.2. Strang-Splitting DUGKS

In this subsection, the evolution process of the discrete unified gas-kinetic scheme is exhaustively clarified. Then, the Strang-splitting scheme for the incorporation of the interaction force is introduced.

2.2.1. Discrete Unified Gas-Kinetic Scheme

The investigation of multiphase flow problems in the current research was conducted by numerically solving the Boltzmann-BGK equation:

$$\frac{\partial f}{\partial t}+\mathbf{\xi }·{\nabla }_{\mathit{x}}f=\mathsf{\Omega }\equiv -\frac{f-f^{\mathrm{E}}}{\tau },$$

(19)

where $f=f(\mathit{x},\mathbf{\xi },t)$ denotes the distribution function (DF), referring to a cluster of particles residing at position $\mathit{x}$ with a velocity of $\mathbf{\xi }$ at time t, $\tau$ indicates the relaxation time, and $f^E$ represents the Maxwellian distribution function approached by f within each collision. The nondimensionalization of Equation (19) is presented in the Appendix A. The moments of distribution functions correspond to the conservative flow variables via

$$\rho =\int fd\mathbf{\xi }=\int f^{\mathrm{E}}d\mathbf{\xi },\rho \mathit{u}=\int \mathbf{\xi }fd\mathbf{\xi }=\int \mathbf{\xi }{f}^{\mathrm{E}}d\mathbf{\xi },$$

(20)

where $\mathit{u}$ denotes the velocity of the flow field. To numerically solve Equation (19), discretization of the physical and velocity space is a prerequisite. To determine the discrete velocity points along each single dimension, the three-point Gauss–Hermite quadrature is employed. The two-dimensional discrete velocity points can be derived from the tensor product of the single-dimensional velocities, which turns out to be the D2V9 velocity model commonly used in the LB community:

$$\begin{array}{cc}\hfill {\mathbf{\xi }}_i& =\sqrt{3c_s^2}\left[\begin{array}{ccccccccccccccccc}\hfill 0& \hfill & \hfill 1& \hfill & \hfill 0& \hfill & \hfill -1& \hfill & \hfill 0& \hfill & 1& & -1& & -1& & 1\\ \hfill 0& \hfill & \hfill 0& \hfill & \hfill 1& \hfill & \hfill 0& \hfill & \hfill -1& \hfill & 1& & 1& & -1& & -1\end{array}\right],\hfill \end{array}$$

where ${\mathbf{\xi }}_i$ is the ith discrete velocity and $c_s=1/\sqrt{3}$ is the model speed of sound. The ideal gas pressure $p_0$ shown in Equation (14) relates to the density $\rho$ through $p_0=\rho c_s^2$.

With the discretization of the velocity space, the Boltzmann-BGK equation turns into

$$\frac{\partial f_i}{\partial t}+{\mathbf{\xi }}_i·{\nabla }_{\mathit{x}}{f}_i={\mathsf{\Omega }}_i\equiv -\frac{f_i-f_i^{\mathrm{E}}}{\tau },$$

(21)

where the subscript i indicates the distribution function for particles possessing a velocity of ${\mathbf{\xi }}_i$. Subdividing the physical space into a set of grid cells and integrating Equation (21) over a certain cell lead to

$$\frac{d}{dt}{\int }_{V_c}{f}_i(\mathit{x},t)d\mathit{x}+{\int }_{\partial V_c}\left(\mathbf{\xi }·\mathit{n}\right)f_i(\mathit{x},t)dS={\int }_{V_c}{\mathsf{\Omega }}_i(\mathit{x},t)d\mathit{x},$$

(22)

where $V_c$ denotes the integral cell centered at position ${\mathit{x}}_c$, $\partial V_c$ denotes the surface boundary of the cell, $dS$ is the surface element, and $\mathit{n}$ is the unit vector normal to the surface element. Integrating Equation (22) over a time step of length $\Delta t=t_{n+1}-t_n$ yields

$$f_i^{n+1}-f_i^n+\frac{\Delta t}{|V_c|}{F}_i^{n+1/2}=\frac{\Delta t}{2}\left[{\mathsf{\Omega }}_i^{n+1}+{\mathsf{\Omega }}_i^n\right],$$

(23)

where $|V_c|$ measures the volume of cell $V_c$ and $f_i^n$ and ${\mathsf{\Omega }}_i^n$ approximate the cell averages of $V_c$ in such a way that

$$f_i^n=\frac{1}{|V_c|}{\int }_{V_c}{f}_i(\mathit{x},t_n)d\mathit{x},$$

(24a)

$${\mathsf{\Omega }}_i^n=\frac{1}{|V_c|}{\int }_{V_c}{\mathsf{\Omega }}_i(\mathit{x},t_n)d\mathit{x}.$$

(24b)

$F_i^{n+1/2}$ measures the kinetic flux at the mid-time $t_n+\Delta t/2$ by

$$F_i^{n+1/2}={\int }_{\partial V_c}\left({\mathbf{\xi }}_i·\mathit{n}\right)f_i(\mathit{x},t_n+\Delta t/2)dS.$$

(25)

Note that the midpoint rule is applied to compute the time integral of the kinetic flux and the trapezoidal rule is applied to evaluate the time integral of the collision term in Equation (23). To remove the implicit treatment of the collision term, two auxiliary distribution functions are introduced:

$${\tilde{f}}_i=f_i-\frac{\Delta t}{2}{\mathsf{\Omega }}_i,{\tilde{f}}_i^{+}=f_i+\frac{\Delta t}{2}{\mathsf{\Omega }}_i.$$

(26)

Substituting Equation (26) into Equation (23), we obtain a fully explicit evolution equation:

$${\tilde{f}}_i^{n+1}={\tilde{f}}_i^{+,n}+\frac{\Delta t}{2}{F}_i^{n+1/2}.$$

(27)

To obtain the kinetic flux $F_i^{n+1/2}$, the primitive distribution function $f_i({\mathit{x}}_f,t_{n+1/2})$ on the cell surface needs to be first evaluated. To this end, we integrate Equation (21) along the characteristic line over a time step length of $\delta t=\Delta t/2$:

$$f_i({\mathit{x}}_f,t_{n+1/2})-f_i({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n)=\frac{\delta t}{2}\left[{\mathsf{\Omega }}_i({\mathit{x}}_f,t_{n+1/2})+{\mathsf{\Omega }}_i({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n)\right].$$

(28)

Note that the trapezoidal rule is once again applied for the time integral of the collision term. Similar to the treatment of Equation (23), the implicitness of Equation (28) is eliminated with the help of the following auxiliary distribution functions:

$$\overline{f}=f-\frac{\delta t}{2}\mathsf{\Omega },{\overline{f}}^{+}=f+\frac{\delta t}{2}\mathsf{\Omega }.$$

(29)

Equation (28) can then be rearranged as

$${\overline{f}}_i({\mathit{x}}_f,t_{n+1/2})={\overline{f}}_i^{+}({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n).$$

(30)

The auxiliary distribution function ${\overline{f}}_i^{+}({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n)$ on the right-hand side of Equation (30) can be interpolated from the cell-centered ${\overline{f}}_i^{+}({\mathit{x}}_c,t_n)$, which could be directly constructed via Equation (29). Based on the expansion point of the Taylor series [57], the reconstruction schemes can be classified into the face-based reconstruction scheme (FRS) or the cell-based reconstruction scheme (CRS). The FRS takes the form of

$${\overline{f}}_i^{+}({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n)={\overline{f}}_i^{+}({\mathit{x}}_f,t_n)-{\mathbf{\xi }}_i\delta t·\nabla {\overline{f}}_i^{+}({\mathit{x}}_f,t_n),$$

(31)

in which the face-centered $f_i^{+}({\mathit{x}}_f,t_n)$ can be reconstructed from the cell-centered $f_i^{+}({\mathit{x}}_c,t_n)$ via the central difference (CD) scheme [31] or the weighted essentially non-oscillatory (WENO) scheme [58]. The upwind CRS takes the form of

$$\begin{array}{c}\hfill {\overline{f}}_i^{+}({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t)=\left\{\begin{array}{cc}{\overline{f}}_i^{+}\left({\mathit{x}}_l\right)+(\delta {\mathit{x}}_l-{\mathbf{\xi }}_i\delta t)·\nabla {\overline{f}}_i^{+}\left({\mathit{x}}_l\right)+{(\delta {\mathit{x}}_l-{\mathbf{\xi }}_i\delta t)}^2:{\nabla }^2{\overline{f}}_i^{+}\left({\mathit{x}}_l\right)/2,\hfill & {\mathbf{\xi }}_i·\mathit{n}\le 0,\hfill \\ {\overline{f}}_i^{+}\left({\mathit{x}}_r\right)+(\delta {\mathit{x}}_r-{\mathbf{\xi }}_i\delta t)·\nabla {\overline{f}}_i^{+}\left({\mathit{x}}_r\right)+{(\delta {\mathit{x}}_r-{\mathbf{\xi }}_i\delta t)}^2:{\nabla }^2{\overline{f}}_i^{+}\left({\mathit{x}}_r\right)/2,\hfill & {\mathbf{\xi }}_i·\mathit{n}>0,\hfill \end{array}\right.\end{array}$$

(32)

where $\delta {\mathit{x}}_l={\mathit{x}}_f-{\mathit{x}}_l$ measures the distance from the face center ${\mathit{x}}_f$ to the adjacent cell center ${\mathit{x}}_l$ on one side, while $\delta {\mathit{x}}_r={\mathit{x}}_f-{\mathit{x}}_r$ measures the distance from the face center ${\mathit{x}}_f$ to the adjacent cell center ${\mathit{x}}_r$ on the other side. An average value is used if ${\mathbf{\xi }}_i·\mathit{n}=0$. After finishing the reconstruction of ${\overline{f}}_i^{+}({\mathit{x}}_f-{\mathbf{\xi }}_i\delta t,t_n)$, the face-centered auxiliary distribution function ${\overline{f}}_i({\mathit{x}}_f,t_{n+1/2})$ can be directly obtained via Equation (30). With a straightforward transformation of Equation (29), the primitive distribution function $f_i({\mathit{x}}_f,t_{n+1/2})$ can be calculated by

$$f=\frac{2\tau }{2\tau +\delta t}\overline{f}+\frac{\delta t}{2\tau +\delta t}{f}^{\mathrm{E}}.$$

(33)

The kinetic flux $F_i^{n+1/2}$ can then be evaluated by its definition. After that, the auxiliary distribution function ${\tilde{f}}_i({\mathit{x}}_c,t_{n+1})$ at the next time step can be updated by Equation (27). Similarly, with a transformation of Equation (26), the primitive distribution function can be calculated by

$$f=\frac{2\tau }{2\tau +\Delta t}\tilde{f}+\frac{\Delta t}{2\tau +\Delta t}{f}^{\mathrm{E}}.$$

(34)

The equilibrium distribution function $f_i^{\mathrm{E}}$ for the primitive free-energy model is expressed as

$$f_i^{\mathrm{E}}={\omega }_i\rho \left[1+\frac{{\mathbf{\xi }}_i·\mathit{u}}{c_s^2}+\frac{\mathit{u}\mathit{u}:({\mathbf{\xi }}_i{\mathbf{\xi }}_i-c_s^2\mathit{I})}{2c_s^4}\right],$$

(35)

where ${\omega }_i=4/9$ for $i=0$, ${\omega }_i=1/9$ for $i=\left\{1,2,3,4\right\}$, and ${\omega }_i=1/36$ for $i=\left\{5,6,7,8\right\}$. The equilibrium distribution function $f_i^{\mathrm{E}}$ for the well-balanced free-energy model is defined as

$$\begin{array}{c}\hfill f_i^{\mathrm{E}}=\left\{\begin{array}{cc}\rho +{\omega }_0\rho s_0\left(\mathit{u}\right),\hfill & i=0,\hfill \\ {\omega }_i\rho s_i\left(\mathit{u}\right),\hfill & i\ne 0,\hfill \end{array}\right.\end{array}$$

(36)

where

$$s_i\left(\mathit{u}\right)=\left[\frac{{\mathbf{\xi }}_i·\mathit{u}}{c_s^2}+\frac{\mathit{u}\mathit{u}:({\mathbf{\xi }}_i{\mathbf{\xi }}_i-c_s^2\mathit{I})}{2c_s^4}\right].$$

(37)

Obviously, the information of macroscopic conservative variables should be first evaluated for the updating of the equilibrium distribution function. Considering the relationship between the auxiliary DF and the primitive DF presented in Equations (26) and (29), the cell-centered conservative variables are updated by

$$\rho =\sum _i{f}_i=\sum _i{\tilde{f}}_i,\rho \mathit{u}=\sum _i{\mathit{\xi }}_i{f}_i=\sum _i{\mathbf{\xi }}_i{\tilde{f}}_i,$$

(38)

and the face-centered conservative variables are updated by

$$\rho =\sum _i{f}_i=\sum _i{\overline{f}}_i,\rho \mathit{u}=\sum _i{\mathbf{\xi }}_i{f}_i=\sum _i{\mathbf{\xi }}_i{\overline{f}}_i.$$

(39)

The time step length $\Delta t$ is determined by the Courant–Friedrichs–Lewy (CFL) condition:

$$\Delta t=C\frac{\Delta x}{\sqrt{3c_s^2}},$$

(40)

where C denotes the CFL number and $\Delta x$ measures the grid spacing.

2.2.2. Strang-Splitting Scheme

To date, the evolution process of DUGKS without considering force terms has been exhaustively clarified. To incorporate the interaction effects between different phases, a source distribution function $f_i^{\mathrm{S}}$ accounting for the force effects is introduced. To correctly recover the macroscopic hydrodynamic equation, the expression of $f_i^{\mathrm{S}}$ for the primitive free-energy model is defined as

$$f_i^{\mathrm{S}}={\omega }_i\left[\frac{{\mathbf{\xi }}_i·\mathit{F}}{c_s^2}+\frac{\mathit{uF}:\left({\mathbf{\xi }}_i{\mathbf{\xi }}_i-c_s^2\mathit{I}\right)}{c_s^4}\right],$$

(41)

where $\mathit{F}$ is the interaction force defined in Equation (14). The expression of $f_i^{\mathrm{S}}$ for the well-balanced free-energy model is defined as

$$f_i^{\mathrm{S}}={\omega }_i\left[\frac{{\mathbf{\xi }}_i·\mathit{F}}{c_s^2}+\frac{\mathit{u}\left(\mathit{F}+{\mathbf{c}}_{\mathbf{s}}^{\mathbf{2}}\nabla \mathbf{\rho }\right):\left({\mathbf{\xi }}_i{\mathbf{\xi }}_i-c_s^2\mathit{I}\right)}{c_s^4}+\frac{1}{2}\left(\frac{{\mathbf{\xi }}_i^2}{c_s^2}-D\right)\left(\mathit{u}·\nabla \rho \right)\right],$$

(42)

where $D=2$ is the spatial dimension. To circumvent the calculation of the interaction force on the cell interface, the Strang-splitting scheme is employed [59]. With such a treatment, the force effects are considered before and after the evolution process of the DUGKS:

$$\frac{\partial f_i}{\partial t}=\frac{1}{2}{f}_i^{\mathrm{S}},$$

(43a)

$$\frac{\partial f_i}{\partial t}+{\mathbf{\xi }}_i·{\nabla }_{\mathit{x}}{f}_i={\mathsf{\Omega }}_i\equiv -\frac{f_i-f_i^{\mathrm{E}}}{\tau },$$

(43b)

$$\frac{\partial f_i}{\partial t}=\frac{1}{2}{f}_i^{\mathrm{S}},$$

(43c)

As Equation (43b) remains identical to Equation (21), it can be solved by the DUGKS procedure addressed previously. Equations (43a) and (43c) can be numerically solved by the forward Euler method:

$$f_i^{\ast }=f_i^n+\frac{\Delta t}{2}{f}_i^{\mathrm{S},n}.$$

(44)

The conservative variables should be accordingly updated via

$${\rho }^{\ast }={\rho }^n,{\mathit{u}}^{\ast }={\mathit{u}}^n+\frac{\Delta t}{2}\frac{{\mathit{F}}^n}{{\rho }^n}.$$

(45)

The gradient operator and Laplacian operator appearing in Equations (7), (14) and (15) are implemented via the isotropic difference scheme [60].

3. Numerical Results

In this section, several numerical tests are conducted by the Strang-splitting DUGKS to compare the performance of the primitive free-energy model and that of the well-balanced free-energy model. The nonisotropic property caused by the reconstruction procedure in the DUGKS is especially discussed. For steady tests, the iteration terminates once the $L_2$-norm error satisfies

$$E\left(\mathit{Q}\right)=\sqrt{\frac{{\sum }_{\mathit{x}}{\left|\mathit{Q}\left(\mathit{x},{\mathit{t}}_{\mathit{n}}\right)-\mathit{Q}\left(\mathit{x},{\mathit{t}}_{\mathit{n}-\mathbf{1000}}\right)\right|}^2}{{\sum }_{\mathit{x}}{\left|\mathit{Q}\left(\mathit{x},{\mathit{t}}_{\mathit{n}}\right)\right|}^2}}<e,$$

(46)

where $\mathit{Q}$ is either the flow density $\rho$ or the flow velocity $\mathit{u}$, $t_{n-1000}$ denotes the moment 1000 time steps ahead of $t_n$, and e is the error threshold.

3.1. Flat Interface

As a benchmark test, the flat interface has been widely applied to validate the performance of newly proposed models [30,55,56]. The computational domain is a $L_0\times 16L_0$ rectangular region with $L_0=16$. A uniform Cartesian mesh with a grid spacing of unity is employed to subdivide this domain. Initially, the region bounded by $y_L=4L_0$ and $y_H=12L_0$ is filled up with the liquid fluid, while the rest is occupied by the gas fluid. The periodic boundary condition is applied to all the sides. The relaxation time $\tau$ was fixed at 0.3. The CFL number was set as 0.5. The reduced temperature $T_r=T/T_c$ ranged from 0.55 to 0.95. The density field is initialized by

$$\rho (x,y)={\rho }_g+\frac{{\rho }_l-{\rho }_g}{2}\left[\tanh \frac{2(y-y_L)}{W}-\tanh \frac{2(y-y_H)}{W}\right],$$

(47)

where W measures the interface thickness and ${\rho }_l$ and ${\rho }_g$ represent the liquid density and the gas density, respectively. Three reconstruction schemes were utilized to explore the influences of varying reconstruction errors on the performance of the DUGKS coupled with different free-energy models. Figure 1a illustrates the coexisting curves predicted by the DUGKS coupled with the primitive free-energy model. It can be observed that varying the reconstruction schemes offers different coexisting results. The central difference face-based reconstruction scheme (CD-FRS) provides satisfactory results in conditions of a high reduced temperature $T_r$. As $T_r$ decreases, the results deviate apparently from the theoretical results generated by the Maxwell equal-area law [61]. The WENO-Z face-based reconstruction scheme (WENO-Z-FRS) and the upwind cell-based reconstruction scheme (CRS) produce inconsistent results in conditions of high $T_r$. As $T_r$ decreases, both of them suffer from the stability problem. The fact that different reconstruction schemes generate divergent outcomes results from the force imbalance in the primitive free-energy model [30]. As the standard LB method involves no reconstruction process, the influences of the force imbalance on the numerical results remain limited. When it is coupled with numerical methods containing a reconstruction process, the effect of the force imbalance becomes amplified by the reconstruction errors. Figure 1b illustrates the results produced by the DUGKS coupled with the well-balanced free-energy model, in which the force imbalance was entirely eliminated. It can be identified that the coexisting densities predicted by different reconstruction schemes coincide exactly with the theoretical results. Moreover, the DUGKS implemented with different reconstruction schemes performs equally well in conditions of a low reduced temperature $T_r$, which demonstrates the fundamental accuracy and stability of this method. Figure 2 illustrates the comparative chemical potential profiles produced by the DUGKS coupled with different free-energy models at $T_r=0.75$, $\tau =0.3$, $C=0.5$. Regardless of the reconstruction schemes utilized, the well-balanced free-energy-based DUGKS provides a nearly constant chemical potential profile, while the primitive free-energy-based DUGKS offers a varied chemical potential profile across the interfaces. Taking a closer look at the comparative profiles, we can identify that the chemical potential value produced by the DUGKS coupled with the primitive model varies along with the reconstruction schemes used, which should be attributed to the differences in the reconstruction errors. The chemical potential produced by the DUGKS coupled with the well-balanced model holds a nearly constant value of 0.006126, which demonstrates the excellent performance of the well-balanced DUGKS in predicting steady two-phase systems governed by free-energy theory.

3.2. Quiescent Droplet

The quiescent droplet test serves as one of the fundamental benchmarks for validating the basic capability of the newly proposed multiphase methods. A circular droplet is initially placed at the center of an $L_0\times L_0$ square domain, with $L_0=256$. A uniform Cartesian mesh is used to discretize the physical domain, with the grid spacing $\Delta x$ fixed at unity. The density field is initialized according to

$$\rho (x,y)=\frac{{\rho }_l+{\rho }_g}{2}-\frac{{\rho }_l-{\rho }_g}{2}\tanh \left[\frac{2\left(\sqrt{|x-x_c{|}^2+{\left|y-y_c\right|}^2}-R_d\right)}{W}\right],$$

(48)

where ${\rho }_l$ and ${\rho }_g$ correspond, respectively, to the coexisting liquid and gas densities, $(x_c,y_c)$ indicates the center location of the square domain, $R_d$ denotes the droplet radius, and W measures the interface thickness. The computing process terminates once the $L_2$-norm error of density evaluated by Equation (46) is below ${10}^{-10}$. Figure 3, Figure 4 and Figure 5 illustrate the density contours produced by the DUGKS coupled with different free-energy models and implemented by various reconstruction schemes at $T_r=0.9,\tau =0.6,C=0.5$. The interfaces produced with the primitive free-energy model suffer from the nonisotropic problem regardless of the reconstruction scheme utilized, which is caused by the force imbalance addressed previously. The second-order central-difference face-based reconstruction scheme (CD-FRS) evolves the initially circular interface into a roughly square interface, which should be attributed to the relatively large reconstruction errors. With a long time evolution, the fifth-order WENO-Z face-based reconstruction scheme (WENO-Z-FRS) shifts the quiescent droplet away from the center position. The interface profile deforms less than that produced by the CD-FRS, which might be attributed to the low level of reconstruction errors of WENO-Z. The interface profile generated by the third-order cell-based reconstruction scheme (CRS) is rather close to circular, which is due to the less nonisotropic reconstruction errors. A similar phenomenon can be observed in the results produced by the pseudopotential-based DUGKS. The interface profiles produced with the well-balanced free-energy model preserve a universal isotropic property across all reconstruction schemes, which demonstrates the elimination of the force imbalance. Figure 6 illustrates the contour of the velocity field produced by the DUGKS implemented with the CRS at $T_r=0.9,\tau =0.6,C=0.5$. When the steady-state is reached, the velocity field produced by the primitive model exhibits a typical patten of large spurious currents, while the velocity field obtained with the well-balanced model provides spurious currents of machine accuracy. The excellent performance of the well-balanced DUGKS is thus verified by the comparative results.

To quantitatively assess its capability, Laplace’s law is validated by the well-balanced DUGKS implemented with the CRS. Figure 7 illustrates the relations between the pressure jump $\Delta P$ and the reciprocal of radius $R_d$ obtained at $\tau =0.3$, $C=0.8$. The linear relation can be clearly identified from the results, which conforms to Laplace’s law: $\Delta P=\sigma /R_d$. The chemical potential varies along with the reduced temperature $T_r$, which results in the alteration of the surface tension coefficient $\sigma$. The CFL number was set as 0.8, at which the FRS fails to operate properly. The stability superiority of the CRS over that of the FRS in the condition of a large time step size makes it more appealing for multiphase flow simulations.

3.3. Spinodal Decomposition

Previous benchmark tests were limited to steady-state problems. Here, the spinodal decomposition test was adopted to assess the capability of DUGKS in dealing with transient problems. The computational domain is an $L_0\times L_0$ square region subdivided by the uniform Cartesian mesh. The grid spacing $\Delta x=1$, and the characteristic length $L_0=512$. The periodic boundary condition was applied to all the sides. The density field is initialized by

$$\rho (x,y)=({\rho }_l+{\rho }_g)/3+\mathrm{random}(0,1)/100,$$

(49)

where ${\rho }_l$ and ${\rho }_g$ represent the liquid density and the gas density and $\mathrm{random}(0,1)$ creates density fluctuations that induce the spinodal decomposition process. Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 illustrate the snapshots of the density distribution produced by the DUGKS coupled with the well-balanced free-energy model at $T_r=0.9$, $\tau =0.6$, $C=0.5$. In the early stages, the tiny fluctuations generate local inhomogeneities, which initialize the phase separation. As the system evolves, the inhomogeneities drive the material of the heavy fluid into small droplets and interfaces separating different phases begin to emerge. With the continual evolution of the whole system, some of these droplets gradually coalesce into large ones. Eventually, a complete quiescent droplet is formed. It can be identified that the results produced by the central difference face-based reconstruction scheme (CD-FRS) are nearly identical to those generated by the third-order cell-based reconstruction scheme (CRS), which demonstrates the consistent behaviors of the well-balanced DUGKS. The WENO-Z face-based reconstruction scheme (WENO-Z-FRS) fails to provide a converged solution in such a condition. Moreover, the well-balanced DUGKS fails to predict the evolution process of the spinodal decomposition when the reduced temperature is below 0.8. To investigate the multiphase flow dynamics by the well-balanced DUGKS, further improvements are required.

3.4. Droplet Coalescence

Simulations of the droplet coalescence phenomenon were used to further investigate the capacity of the well-balanced DUGKS for transient problems. The computational domain is a rectangle $2L_0\times L_0$ domain with $L_0=256$. The domain was subdivided into finite grid cells by a uniform Cartesian mesh with a grid spacing of unity. To avoid wall boundary influence, a periodic boundary condition was used in all directions. Initially, two circular droplets were arranged in accordance with [51]

$$\rho (x,y)=\frac{{\rho }_l+{\rho }_g}{2}+\frac{{\rho }_l-{\rho }_g}{2}\left[1-\tanh \left(\frac{2d_A}{W}\right)-\tanh \left(\frac{2d_B}{W}\right)\right],$$

(50)

where ${\rho }_l$ and ${\rho }_g$ correspond separately to the liquid and gas densities, W measures the interface thickness, and $d_A$ and $d_B$ are defined as

$$d_A=\sqrt{{(x-x_A)}^2+{(y-y_A)}^2}-R_0,d_B=\sqrt{{(x-x_B)}^2+{(y-y_B)}^2}-R_0,$$

(51)

in which $R_0$ denotes the droplet radius and $(x_A,y_A)=(L_0-R_0-W/2,L_0/2)$ and $(x_B,y_B)=(L_0+R_0+W/2,L_0/2)$ represent the central position of droplets A and B, respectively. Other parameters were set as $\kappa =0.02$, $R_0=0.2L_0$, $W=5$, and $\tau =0.3$. The initial profile of two droplets is illustrated in Figure 13. The coalescence process starts when the droplets come in contact with each other. As the process continues, a liquid bridge of radius $r_b$ that connects the two droplets is formed [51]. Previous research [62] identified the linear relation between the scaled radius $r^{\ast }$ and the dimensionless time $t^{\ast }$, with

$$r^{\ast }=r_b/R_0,t^{\ast }=t/\sqrt{{\rho }_l{R}_0^3\sigma },$$

(52)

where $\sigma$ is the surface tension coefficient. According to the validation of Laplace’s law illustrated in Figure 7, the surface tension coefficient is 0.1203 for $T_r=0.8$ and 0.0435 for $T_r=0.9$. Figure 14 presents the radius variation of the liquid bridge with regard to the dimensionless time $t^{\ast }$. The linear coefficient for the fitting result provided by the DUGKS using the primitive model is 1.4, while the linear coefficient for the fitting result produced with the well-balanced model is 1.03, which is in good agreement with the result predicted by Zeng et al. [51]. The evolution of the $L_2$-norm of the velocity field produced by the DUGKS using the well-balanced model at $T_r=0.8$ and $T_r=0.9$ is shown in Figure 15. It can be identified that the $L_2$-norm of the velocity field reaches a magnitude of ${10}^{-14}$, which is consistent with the results predicted at the steady-state. Figure 16 illustrates the density and velocity contours produced by the well-balanced DUGKS at $t=6\times {10}^6$, $T_r=0.8$, $\tau =0.3$, $C=0.8$. It can be observed that the interface maintains excellent isotropy and the velocity field holds a maximum magnitude of ${10}^{-16}$, which demonstrates the excellent ability of the well-balanced DUGKS. However, it is important to note that when the lowered temperature $T_r$ is less than 0.7, the DUGKS is unable to predict the coalescence process. More efforts are required to increase the stability of the well-balanced DUGKS.

4. Conclusions

A free-energy-based discrete unified gas-kinetic scheme (DUGKS) was developed by coupling the well-balanced free-energy model with the DUGKS to investigate the van der Waals fluid. The performance of this well-balanced scheme was compared against the counterpart of the DUGKS coupled with the primitive free-energy model. Comparative results produced with different reconstruction schemes demonstrated the force imbalance in the primitive free-energy model, which prevents its direct application to the DUGKS. By coupling the well-balanced free-energy model with the DUGKS, the amplified effects of the force imbalance are entirely eliminated and the influences of nonisotropic reconstruction errors on the fluid interfaces are totally removed. Numerical tests of a flat interface, quiescent droplet, spinodal decomposition, and droplet coalescence were adopted to assess the performance of the DUGKS coupled with the well-balanced free-energy model. Coexisting density curves and Laplace’s law were utilized to evaluate its capability quantitatively. It was proven that the well-balanced DUGKS could always produce consistent results despite the reconstruction schemes utilized in steady cases. When dealing with transient problems, the reconstruction scheme employing WENO-Z to evaluate face unknowns tends to be more unstable. When the reduced temperature is below 0.7, the DUGKS coupled with the well-balanced free-energy model suffers from stability problems. Further improvements are required to apply this scheme to predict transient multiphase fluid flows.