Physical modeling of multiphase flow via lattice Boltzmann method: Numerical effects, equation of state and boundary conditions

Abstract

The aims of the present paper are threefold. First, we further study the fast Fourier transform thermal lattice Boltzmann (FFT–TLB) model for van der Waals (VDW) fluids proposed in Phys. Rev. E, 2011, 84(4): 046715. We analyze the merits of the FFT approach over the traditional finite difference scheme and investigate the effects of smoothing factors on accuracy and stability in detail. Second, we incorporate the VDW equation of state with flexible parameters into the FFT–TLB model. As a result, the revised model may be used to handle multiphase flows with various critical densities and temperatures. Third, we design appropriate boundary conditions for systems with solid walls. These improvements, from the views of numerics and physics, significantly extend the application scope of the model in science and engineering.

Graphical abstract

Keywords

van der Waals fluids / lattice Boltzmann method / FFT / equation of state

1 Introduction

Modeling and simulating multiphase flows is a challenging task for traditional computational fluid dynamics (CFD). Firstly, the CFD approach is generally based on the macroscopic continuum equations, e.g., Navier–Stokes equations (NSEs), and then discretized by finite difference, finite element, or finite volume schemes, etc. Secondly, there exists a wide range of space time scales in the multiphase system owing to the existence of complex interfaces between various phases/components, which may break the validity of the macroscopic thermohydrodynamic description [1]. Nevertheless, as a mesoscopic approach and a bridge between the microscopic molecular dynamics and the macroscopic CFD, the lattice Boltzmann (LB) method [2–16] is particularly promising in the multiphase area [2, 3, 17–45]. The reasons are also twofold. On the one hand, it should be noted that, the interparticle interaction (IPI) is the underlying physical reason for phase separation and interfacial tension in such systems. When the IPI is appropriately incorporated, the complex behaviors of the system will emerge naturally. On the other hand, the distinct advantage of the LB method, i.e., the intrinsic kinetic nature, makes the IPI be incorporated easily and flexibly into the LB framework.

Broadly speaking, to date, there are two kinds of LB models for the simulation of multiphase flows. The first is the isothermal models, where only mass and momentum conservations are kept, thermodynamic effects are not taken into account. Well-known models of this category includes the chromodynamic model by Gunstensen et al. [18], the pseudo-potential model by Shan and Chen (SC) [19, 20], and the free-energy model by Swift et al. [21, 22], etc. The chromodynamic model is a further development of the two-component lattice-gas automata model originally proposed by Rothman and Keller [17]. In this approach, red and blue colored particles are used to represent two different fluids. Phase separation is achieved through controlling the IPI based on the color gradient. In the SC-type model, nonlocal interactions between particles at neighboring lattice sites are incorporated. The interactions determine the form of the equation of state (EOS) and the interfacial dynamics of the system. In the free-energy model, besides the mass and momentum conservation constraints, additional ones are imposed on the equilibrium distribution function to make the pressure tensor consistent with that of the free-energy functional of inhomogeneous fluids.

Given the great importance of thermohydrodynamics of multiphase flows in many fields of physics and engineering, constructing LB models for nonideal fluids with thermodynamics has been attempted since the early days of LB research. This results in the second category, the thermal LB models. Existing thermal LB models can be roughly classified into four subcategories, i.e., the hybrid approach [34–37], the double-distribution-function (DDF) approach [38], the shifting approach [39], and the multispeed-extra-force (MEF) approach [40–43]. In the hybrid approach, the mass and momentum are described by an isothermal multiphase LB model, while the temperature dynamics are obtained by solving a scalar transport energy equation via various numerical schemes. To ensure that the velocity and pressure fields evolve in a thermodynamically consistent manner, an appropriate coupling between the isothermal and thermal modules of the numerical scheme is introduced through a suitably defined body force. This approach is conceptually simple and numerically robust, because of the absence of energy conservation constraint in the LB model. Alternatively, in the DDF approach, two different distribution functions are utilized. One is used to simulate the density and the velocity fields, and the other one is employed to monitor the temperature field. Actually, the DDF approach can be regarded as another version of the hybrid approach, since the only difference between them is that the energy equation in the DDF approach is solved by an auxiliary LB model, rather than by solving the energy equation at the macroscopic level as in the hybrid approach [46]. The third is the shifting approach, which implements energy conservation by introducing an additional suitable shift of the temperature field in the local equilibrium besides Shan–Chen’s velocity shift. The MEF approach is generally based on a multispeed TLB model for ideal gas system, e.g., the D2V33-TLB model proposed by Watari and Tsutahara (WT) [47]. To obtain the correct full set of thermohydrodynamics at the mesoscopic level, the multispeed model uses discrete-velocity model (DVM) with high spatial isotropy and includes higher-order of the flow velocity in the Taylor expansion of the equilibrium distribution. Although originally works only for ideal gas system, the multispeed model can be extended to multiphase system using the extra force method. The one developed by Gonnella, Lamura, and Sofonea (GLS) [41] is typical. In this model, an extra term I_{ki}, accounting for IPI, is added into the LB equation to describe the van der Waals (VDW) fluids. Chapman–Enskog (CE) procedure indicates that the model reproduces in the continuum the full set of thermohydrodynamic equations with the stress terms as recently developed by Onuki [48, 49]. Compared to TLB multiphase models listed above, the model considers all factors, e.g., viscous dissipation, interfacial tension, and compression work done by the pressure, in both the momentum and energy equations.

In a recent work [42], we improved the GLS model by using the fast Fourier transform (FFT) and its inverse (IFFT) to discretize the spatial derivatives in the convection term v_{ki}\cdot \partial f_{ki}/\partial r and the force term I_{ki}, so that the total energy conservation can be better held and the spurious velocities can be damped to a negligible scale in numerical simulations. For convenience of description, we refer to this model as the FFT–TLB model. However, the previous version is still subject to at least the following two constraints: (i) fixed critical density and temperature, and (ii) suitable only for periodic boundary conditions (PBCs). In the present study, we further analyze the merits of the FFT algorithm over traditional finite difference (FD) schemes and investigate the effects of smoothing factors on the accuracy and stability of the model in detail, then try to remove the above two limitations by inserting VDW EOS with flexible parameters and designing improved boundary conditions for nonperiodic system.

The rest of the paper is structured as follows. The FFT–TLB model is briefly reviewed in Section 2. In Section 3, we analyze the merits of the FFT scheme, incorporate new VDW EOS into the FFT–TLB model, and discuss how to implement solid wall boundary conditions when the FFT scheme is used. Finally, in Section 4, we summarize the results and suggest directions for future research.

2 Outline of the model

Here we present the FFT–TLB model for VDW fluids. The model consists of the following three parts: (i) WT model for ideal gas system, (ii) an appropriate interparticle force accounting for nonideal gas effects, and (iii) the FFT approach used to calculate the spatial derivatives.

2.1 WT model for ideal gas system: Based on the multispeed approach

By using the multispeed approach, Watari and Tsutahara proposed a TLB model for compressible flows [47]. The model uses a DVM that is at least up to seventh rank isotropy,

\begin{array}{}\text{(1)}& v_0=0,v_{ki}=v_k[\cos \frac{\text{i}\pi }{4},\sin \frac{\text{i}\pi }{4}]\end{array}

where subscript k=1,\mathrm{\cdots },4 indicates the k-th group of the particle velocities whose speed is v_k, and i=1,\mathrm{\cdots },8 denotes the direction of particle’s speed. The equilibrium distribution function f_{ki}^{\mathrm{eq}}, corresponding to a Maxwellian distribution, is expanded in polynomial up to the fourth order of fluid velocity u,

\begin{array}{}\text{(2)}& \begin{array}{cc}{f}_{ki}^{\mathrm{eq}}=\hfill & n{F}_k[(1-{\displaystyle \frac{u^2}{2T}+\frac{u^4}{8T^2})}+{\displaystyle \frac{v_{ki}\cdot u}{T}}(1-{\displaystyle \frac{u^2}{2T}})\hfill \\ \hfill & +{\displaystyle \frac{{(v_{ki}\cdot u)}^2}{2T^2}}\left(1-{\displaystyle \frac{u^2}{2T}}\right)+{\displaystyle \frac{{(v_{ki}\cdot u)}^3}{6T^3}}\hfill \\ \hfill & +{\displaystyle \frac{{(v_{ki}\cdot u)}^4}{24T^4}}]\hfill \end{array}\end{array}

with F_k the weighting factors, T the local temperature and n the density. Hydrodynamic quantities can be calculated from the following moments:

\begin{array}{}\text{(3)}& n={\displaystyle \sum _{ki}}{f}_{ki}\end{array}

\begin{array}{}\text{(4)}& nu={\displaystyle \sum _{ki}}{v}_{ki}{f}_{ki}\end{array}

\begin{array}{}\text{(5)}& nT={\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{(v_{ki}-u)}^2{f}_{ki}\end{array}

where f_{ki} is the the distribution function.

2.2 GLS model for nonideal gas system: Based on the extra force approach

Gonnella, Lamura, and Sofonea improved the WT model by introducing an appropriate force term, I_{ki}, to describe the nonideal gas effects. The GLS–LB equation reads

\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial f_{ki}}{\partial t}}+v_{ki}\cdot {\displaystyle \frac{\partial f_{ki}}{\partial r}}=-{\displaystyle \frac{1}{\tau }}\left[f_{ki}-f_{ki}^{eq}\right]+I_{ki}\end{array}

where I_{ki} takes the following form:

\begin{array}{}\text{(7)}& I_{ki}=-[A+B\cdot (v_{ki}-u)+(C+C_q){(v_{ki}-u)}^2]f_{ki}^{eq}\end{array}

The GLS-LB equation is designed to recover the following NSEs for VDW fluids:

\begin{array}{}\text{(8)}& {\partial }_{t}n+\nabla \cdot (nu)=0\end{array}

\begin{array}{}\text{(9)}& {\partial }_t(nu)+\nabla \cdot (nuu+\Pi -\sigma )=0\end{array}

\begin{array}{}\text{(10)}& {\partial }_t{e}_T+\nabla \cdot [e_{T}u+(\Pi -\sigma )\cdot u-{\kappa }_T\nabla T]=0\end{array}

where \Pi =P^{w}I+\Lambda is the reversible part of stress, comprising the VDW EOS P^w and the nonideal gas interaction term \Lambda =M\nabla n\nabla n-M(n{\nabla }^{2}n+{\left|\nabla n\right|}^2/2)I-[nT\nabla n\cdot \nabla (M/T)]I, with M=K+HT,K the surface tension coefficient and H a constant, I the unit tensor. \sigma =\eta [\nabla u+{(\nabla u)}^T-(\nabla \cdot u)I]+\zeta (\nabla \cdot u)I is the usual viscous stress tensor with the shear and bulk viscosities \eta and \zeta, respectively. e_T=nT-9n^2/8+K{\left|\nabla n\right|}^2/2+n{u}^2/2 is the total energy density including the gradient contributions.

Chapman–Enskog (CE) expansion indicates that the coefficients A, B, C, and C_q in Eq. (7) can be expressed in the following forms so as to reproduce Eqs. (8)-(10):

\begin{array}{}\text{(11)}& A=-2(C+C_q)T\end{array}

\begin{array}{}\text{(12)}& B={\displaystyle \frac{1}{nT}}[\nabla (P^w-nT)+\nabla \cdot \Lambda -\nabla (\zeta \nabla \cdot u)]\end{array}

\begin{array}{}\text{(13)}& \begin{array}{cc}C=\hfill & {\displaystyle \frac{1}{2n{T}^2}}\{(P^w-nT)\nabla \cdot u+\Lambda :\nabla u\hfill \\ \hfill & -\zeta {(\nabla \cdot u)}^2+{\displaystyle \frac{9}{8}}{n}^2\nabla \cdot u+K[-{\displaystyle \frac{1}{2}}(\nabla n\cdot \nabla n)\nabla \cdot u\hfill \\ \hfill & -n\nabla n\cdot \nabla (\nabla \cdot u)\mathbf{-}\nabla n\mathbf{\cdot }\nabla u\cdot \nabla n]\}\hfill \end{array}\end{array}

\begin{array}{}\text{(14)}& C_q={\displaystyle \frac{1}{2n{T}^2}}\nabla \cdot [2qnT\nabla T]\end{array}

It is worth remarking that the introduction of the term C_q makes the thermal conductivity {\kappa }_T be tuned independently from the shear viscosity η. As a result, the model can handle flows with variable Prandtl numbers by adjusting the parameter q,

\begin{array}{}\text{(15)}& \mathrm{Pr}={\displaystyle \frac{\eta }{{\kappa }_T}}={\displaystyle \frac{\tau }{2(\tau -q)}}\end{array}

Compared to the multiple-relaxation-time [50] or the hybrid approach [34–37] that can also achieve this goal, the extra force approach is more easy to implement and contains more physical connotation.

2.3 The FFT approach

Here we briefly review the basics of the FFT approach, and conduct a further analysis on the merits of this way in detail in the next section.

This approach originates from a general theorem of derivative based on FFT which reads [51, 52]

\begin{array}{}\text{(16)}& \widetilde{f^{\prime }}(k_x)=\mathrm{i}{k}_x\times \tilde{f}(k_x)\end{array}

where i is an imaginary unit, k_x the wave number along the x -direction, and \widetilde{f^{\prime }}(k_x) the Fourier transform of f^{\prime }(x). The theorem provides a way to calculate the spatial derivative f^{\prime }(x) composed of the following three steps: (i) Transforming f(x) in real space into \tilde{f}(k_x) in reciprocal space; (ii) multiplying \tilde{f}(k_x) with \mathrm{i}{k}_x; (iii) taking the IFT of \widetilde{f^{\prime }}(k_x), and then f^{\prime }(x) can be obtained. This approach is naturally suitable for periodic system since the FFT and IFT algorithms intrinsically adapt the PBCs. Besides, via this way, higher-order and fractional-order derivatives can be conveniently computed only if we replace \mathrm{i}{k}_x with {(\mathrm{i}{k}_x)}^n, where n is the order of the derivative. This is one of the main merits of the FFT over FD scheme.

Within the condition that f(x) is infinitely differentiable, the FFT approach gives the exact result with infinite order accuracy [51–54], typically well beyond that of the standard FD schemes. Nevertheless, in many cases it is hard to ensure that the infinite differentiability condition is well satisfied. For example if f^{\prime }(x) has a discontinuity, then the discontinuity will induce oscillations, known as the Gibbs phenomenon. The Gibbs effects influence the accuracy of the FFT not only in the neighborhood of the point of singularity, but also over the entire computational domain. In the recent work [42], based on Taylor series expansion of wave number k_x, we presented a way to construct smoothing factors that can be used to refrain the Gibbs oscillations effectively. First, we expand k_x in Taylor series

\begin{array}{}\text{(17)}& \begin{array}{ccc}{k}_x\hfill & =\hfill & {\displaystyle \frac{\arcsin [\sin (k_x\mathrm{\Delta }x/2)]}{\mathrm{\Delta }x/2}}\hfill \\ \hfill & =\hfill & {\displaystyle \frac{1}{\mathrm{\Delta }x/2}}[\sin (k_x\mathrm{\Delta }x/2)+{\displaystyle \frac{1}{6}}{\sin }^3(k_x\mathrm{\Delta }x/2)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{3}{40}}{\sin }^5(k_x\mathrm{\Delta }x/2)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{5}{112}}{\sin }^7(k_x\mathrm{\Delta }x/2)+\mathrm{\cdots }]\hfill \end{array}\end{array}

Second, in order to damp the Gibbs oscillations, or, equivalently, to filter out more high frequency waves, k_x may take the form of an appropriately truncated Taylor series expansion of \sin (k_x\mathrm{\Delta }x/2), e.g., k_x may take the following forms:

\begin{array}{}\text{(18)}& k_1\mathbf{=}{\displaystyle \frac{\sin (k_x\mathrm{\Delta }x/2)}{\mathrm{\Delta }x/2}}\end{array}

\begin{array}{}\text{(19)}& k_2\mathbf{=}{k}_1+{\displaystyle \frac{{\sin }^3(k_x\mathrm{\Delta }x/2)/6}{\mathrm{\Delta }x/2}}\end{array}

\begin{array}{}\text{(20)}& k_3\mathbf{=}{k}_2+{\displaystyle \frac{3{\sin }^5(k_x\mathrm{\Delta }x/2)/40}{\mathrm{\Delta }x/2}}\end{array}

and

\begin{array}{}\text{(21)}& k_4\mathbf{=}{k}_3+{\displaystyle \frac{5{\sin }^7(k_x\mathrm{\Delta }x/2)/112}{\mathrm{\Delta }x/2}}\end{array}

Some simple derivations demonstrate that the FFT approach with k_1, k_2, k_3, and k_4 has a second-order, fourth-order, sixth-order, and eighth-order accuracy in space, respectively. Then the smoothing factor for k_1 can be formulated as:

\begin{array}{}\text{(22)}& {\sigma }_1={\displaystyle \frac{k_1}{k}}={\displaystyle \frac{\sin (n\pi /N)}{n\pi /N}},n=-N/2,\mathrm{\cdots },N/2\end{array}

Smoothing factors {\sigma }_k for k₂, k₃, and k₄ can be calculated in a similar way.

Here we only discuss how to approximate \partial f/\partial x for simplicity. In the y-direction we can do in a similar way, since the derivatives \partial f/\partial x and \partial f/\partial y are approximated dimension by dimension. When approximating \partial f/\partial y, one fixes x and use a one-dimensional approximation in the y-direction.

3 Numerical results and physical analyses

In this section, we analyze the merits of the FFT scheme with k₁ over the traditional FD scheme, e.g., the second-order central difference (2nd-order CD) scheme; discuss effects of the smoothing factors on the FFT approach; incorporate VDW EOS with flexible parameters; discuss how to adapt suitable boundary conditions for the FFT scheme when solid walls are present.

3.1 Merits of the FFT scheme over the 2nd-order CD scheme

The FFT scheme with k₁ is 2nd-order in precision, the 2nd-order CD scheme has the same precision. While, compared to the latter, the former has higher precision, better numerical stability, more advantages to guarantee energy conservation and to refrain spurious velocities. Next, we illustrate these advantages one by one.

3.1.1 Higher precision

The FFT scheme with k₁ corresponds to the following formula:

\begin{array}{}\text{(23)}& \widetilde{f^{\prime }}(k_x)=\mathrm{i}{k}_1\times \tilde{f}(k_x)\end{array}

Taking IFT of the RHS of Eq. (23) gives

\begin{array}{}\text{(24)}& \begin{array}{cc}\mathrm{IFT}[\mathrm{i}{k}_1\times \tilde{f}(k_x)]\hfill & ={\displaystyle \frac{1}{L}}{\displaystyle \sum _{n=-N/2}^{N/2-1}}{\mathrm{e}}^{\mathrm{i}{k}_x{x}_j}\times {\displaystyle \frac{\mathrm{i}}{\mathrm{\Delta }x/2}}\sin (k_x\mathrm{\Delta }x/2)\times \tilde{f}(k_x)\hfill \\ \hfill & ={\displaystyle \frac{1}{L}}{\displaystyle \sum _{n=-N/2}^{N/2-1}}{\mathrm{e}}^{\mathrm{i}{k}_x{x}_j}{\displaystyle \frac{{\mathrm{e}}^{\mathrm{i}{k}_x\mathrm{\Delta }x/2}-{\mathrm{e}}^{-\mathrm{i}{k}_x\mathrm{\Delta }x/2}}{\mathrm{\Delta }x}}\times \tilde{f}(k_x)\hfill \\ \hfill & ={\displaystyle \frac{f(x_j+\mathrm{\Delta }x/2)-f(x_j-\mathrm{\Delta }x/2)}{\mathrm{\Delta }x}}\hfill \\ \hfill & =f^{\prime }(x_j)+{\displaystyle \frac{1}{24}}\mathrm{\Delta }{x}^2{f}^{(3)}(x_j)+\mathrm{\cdots }\hfill \end{array}\end{array}

where f^{(3)}(x_j) represents the third-order derivative. Therefore, the FFT scheme with operator k₁ has a 2nd-order accuracy in space. Similarly, if k_1=\frac{\sin (k_x\mathrm{\Delta }x)}{\mathrm{\Delta }x}, then the FFT approach corresponds to the normal 2nd-order CD scheme,

\begin{array}{}\text{(25)}& \begin{array}{cc}\mathrm{IFT}\left[\mathrm{i}{\displaystyle \frac{\sin (k_x\mathrm{\Delta }x)}{\mathrm{\Delta }x}}\times \tilde{f}(k_x)\right]\hfill & ={\displaystyle \frac{f(x_j+\mathrm{\Delta }x)-f(x_j-\mathrm{\Delta }x)}{2\mathrm{\Delta }x}}\hfill \\ \hfill & =f^{\prime }(x_j)+{\displaystyle \frac{1}{6}}\mathrm{\Delta }{x}^2{f}^{(3)}(x_j)+\mathrm{\cdots }\hfill \end{array}\end{array}

Although both Eqs. (24) and (25) are 2nd-order, the remainder errors are different from each other. In Eq. (24) the remainder error is \mathrm{\Delta }{x}^2{f}^{(3)}(x_j)/24, which is smaller than the one in Eq. (25). So, from this point of view, the FFT scheme with k_1\mathbf{=}\frac{\sin (k_x\mathrm{\Delta }x/2)}{\mathrm{\Delta }x/2} has higher precision than the 2nd-order CD scheme does. Similarly, the FFT scheme with k₂ has higher accuracy than the fourth-order CD scheme does.

3.1.2 Better numerical stability

Figure 1 shows patterns calculated from GLS model with the 2nd-order CD scheme in (a) and the FFT scheme with k_1 in (b) at t=1.5. The time derivative is calculated via the 2nd-order Runge–Kutta scheme. The initial conditions are set as \rho =1+\Delta, T=0.85, u_x=u_y=0.0, where \Delta is a random density with an amplitude of 0.01. The remaining parameters are set to be v_1=1.00, v_2=1.90, v_3=2.90 , v_4=4.30, \tau ={10}^{-2}, \mathrm{\Delta }x=\mathrm{\Delta }y=1/256, \mathrm{\Delta }t={10}^{-5}, K=5\times {10}^{-6}, H=0, \zeta =0, q=-0.004, and lattice size N_x=N_y=128\times 128. PBCs are imposed on x- and y-directions. We see that when the 2nd-order CD scheme is used, spurious oscillations will occur [see Fig. 1(a)]. Numerical results indicate that such oscillations will also appear if the 4th-order CD scheme is employed. These oscillations enlarge spatial derivatives in both the convection term and the extra force term. As a result, with the decrease of \tau, the simulation becomes unstable. For example, when \taudecreases to 5\times {10}^{-5}, simulation with the 2nd-order CD scheme is unstable, while the one by FFT scheme with k₁ can keep on working.

3.1.3 Better energy conservation

Figure 2 depicts the variations of the total energy for the phase separating system plotted in Fig. 1. To be seen is that the FFT scheme has more advantage to guarantee energy conservation.

3.1.4 Smaller spurious velocities

In order to clarify this point, we conduct two simulations about the liquid-vapor interface with the 2nd-order CD scheme and the FFT scheme. Simulations were carried out over a 256\times 2 domain with PBCs in both directions. The initial conditions are set as

where {\rho }_v=0.80 and {\rho }_l=1.20 are the theoretical values at T=0.99. Parameters are unchanged. Time evolution of the maximum velocity u_{\max }^s with these two schemes are given in Fig. 3. Compared to the 2nd-order CD scheme, spurious velocities obtained from the FFT approach decrease to a negligible scale when t>20.

3.2 Effects of the smoothing factors

Smoothing factors {\sigma }_k for k_{\mathrm{CD}}, k_1, k_2, k_3, and k_4 with N=128 and \mathrm{\Delta }x=1/256 are displayed in Fig. 4, where k_{\mathrm{CD}}=\frac{\sin (k_x\mathrm{\Delta }x)}{\mathrm{\Delta }x} and {\sigma }_{\mathrm{CD}}=\frac{k_{\mathrm{CD}}}{k_x}=\frac{\sin (2n\pi /N)}{2n\pi /N} denote wave number and smoothing factor for the 2nd-order CD scheme, respectively. One sees that the lower-order smoothing factors, e.g., {\sigma }_1 and {\sigma }_2, filter out more high-frequency waves, and may result in excessively smeared approximations, which are unfaithful representations of the true physics. On the other hand, the higher-order smoothing factors, such as {\sigma }_3 and {\sigma }_4, reserve more higher-frequency waves, but may not damp the Gibbs phenomenon when the discontinuity is strong enough, then cause numerical instability. The smoothing factors should survive the dilemma of stability versus accuracy.

Compared to {\sigma }_1, {\sigma }_2, {\sigma }_3, and {\sigma }_4, {\sigma }_{\mathrm{CD}} damps the highest waves to zero. Consequently, we may conclude that the oscillations shown in Fig. 1(a) is due to the overshoot of the filtering operator {\sigma }_{\mathrm{CD}} that results in a serious lack of high-frequency waves. In other words, high-frequency waves are necessary to describe the continuous sharper interface, the lack of them will lead to the lack of corresponding information. While for a shock interface, too many high-frequency waves will induce numerical instability. Therefore, smoothing factors should be minimal but make the evolution stable at the same time.

Figure 5 shows patterns calculated from the FFT–TLB model at t=2.0, where {\sigma }_1, {\sigma }_2, {\sigma }_3, and {\sigma }_4 are used in (a), (b), (c), and (d), respectively. The initial conditions and parameters are consistent with what we used in Fig. 1. We cannot find obvious differences between these patterns. To further analyze the differences, the Minkowski analysis technique [31, 42] is employed. We choose a density threshold {\rho }_{th} and pixelize the map into high-density regions (with \rho >{\rho }_{th}) and low-density regions (with ). Then according to a general theorem of integral geometry, three Minkowski measures, the white area fraction A, the boundary length L and Euler characteristic \chi, totally describe the Turing patterns [55]. Figure 6 presents the time evolution of these measures with two density thresholds {\rho }_{th}=0.80 and {\rho }_{th}=1.20. Apparently, for {\rho }_{th}=0.80 or {\rho }_{th}=1.20, in each figure, the four lines obtained from the 2nd-order, the 4th-order, the 6th-order, and the 8th-order schemes overlap with each other, indicating that the above four smoothing factors have negligible effects on the morphological characterization of the phase separating system.

Figures 5 and 6 demonstrate that if we focus on the macroscopic thermohydrodynamic behaviors of multiphase flows, e.g., the evolution of density patterns, the rheological and morphological characterizations of patterns (growth law of domain growth and structure factors), the differences induced by different smoothing factors can be neglected. However, if we focus on the interfacial dynamics, e.g., spurious velocities, thermal effects of motion of interfaces and so on, we suggest to use the higher-order smoothing factors. Figure 7 shows the time evolutions of u_{\max }^s with the 2nd-order, the 4th-order, the 6th-order, and the 8th-order FFT schemes for a liquid-vapor interface. The initial conditions are consistent with what we used in Fig. 3, except for K={10}^{-5}. It is clear that, in each case, u_{\max }^s oscillates and decreases at the initial stage, then reduces to nearly a constant when t>30. More importantly, with the increase of precision, u_{\max }^s decreases further. Besides the advantage of refraining spurious velocities more effectively, the higher-order smoothing factors have more strengths in guaranteeing the total energy conservation, compared to the lower-order ones, as shown in Fig. 8. With the increase of precision, variations of the total energy e_T(t)-e_T(0) for the phase separating system displayed in Fig. 5, decrease consistently.

3.3 Inserting VDW EOS with flexible parameters

The FFT–TLB model uses VDW EOS which provides a simple description for the coexistence of vapor and liquid. It is based on the assumption of a long range mean field attraction a between molecules and a hard core excluded volume b, and gives the pressure P^w as a function of density n and temperature T,

\begin{array}{}\text{(27)}& P^w={\displaystyle \frac{nT}{1-bn}}-a{n}^2\end{array}

However, in this model, a=9/8 and b=1/3 are fixed parameters. Therefore, the critical density and temperature are fixed to n^c=1/(3b)=1.0 and T^c=8a/(27b)=1.0, respectively, which is an inconvenience in the use of the model to study the thermal multiphase flows. Actually, CE analysis demonstrates that this limitation can be removed by modifying the extra force term, i.e., replacing the term \frac{9}{8}{n}^2\nabla \cdot u in Eq. (13) with a{n}^2\nabla \cdot u. As a result, the improved model can recover VDW EOS with flexible critical density and temperature.

To check if the model can correctly reproduce the equilibrium thermodynamics of a VDW system at various temperatures, a series of simulations were performed. Here we test four cases with different critical densities and temperatures, as listed in Table 1. The initial conditions are described by

where {\rho }_v=0.955{\rho }^c and {\rho }_l=1.045{\rho }^c are the theoretical vapor and liquid densities at T=0.9995T^c, respectively. Parameters are set to be \tau ={10}^{-3}, K=5\times {10}^{-6}, leaving all the others unchanged. The initial temperature is set to be 0.9975T^c, but drops by a small value when the equilibrium state of the system has been achieved. Simulations are then run until the temperature has reduced to 0.85T^c. Here the FFT scheme with k₄ is used to discretize the spatial derivatives. Figure 9 illustrates phase diagrams recovered with LB simulations and Maxwell constructions for the four cases listed in Table 1, where the scattered points are for simulation results and solid lines are for analytical solutions. For each case, the two sets of results have a satisfying agreement, demonstrating that the improved model has high accuracy in capturing various equilibrium states under different EOSs and temperatures.

3.4 Design of appropriate boundary conditions in the presence of solid walls

Although the FFT approach is powerful for periodic systems, it can also be applied in systems which have walls in the x- and/or y-directions provided that the following simple treatments are made. For example, for a system with left and right walls, we can simply double the computational domain by making a mirror reflection of the original one around the right (or left) wall. Thus, for the extended computational domain, the left and the right walls are exactly the same. We may consider that the system periodic and the FFT scheme can be used in the horizontal direction. We show this via the following example.

Let us consider a system, as schematically shown in Fig. 10, with (i) PBCs in the y-direction and (ii) left and right walls with fixed temperatures T_L and T_R, respectively. To use the FFT algorithm, we first make a mirror reflection of the original computational domain around the right wall. Thus, the number of lattice in the x-direction is N_x^{{}^{\prime }}=2N_x. There are two outside boundaries and two inner boundaries in the x-direction, as shown in Fig. 10 by “a”, “b”, “c”, and “d”. Here, “b” and “c” are symmetric about the vertical line x=N_x/2+1/2, “a” and “d” are symmetric about the same line, too. Consequently, physical quantities on boundary “b” equal to the ones on boundary “c”, and so are the physical quantities on boundaries “a” and “d”. As a result, on boundary “a”, we set

\begin{array}{}\text{(29)}& \{\begin{array}{c}T(1,j)=T_L\hfill \\ u_x(1,j)=0.0\hfill \end{array}\end{array}

on boundary “b”,

\begin{array}{}\text{(30)}& \{\begin{array}{c}T(N_x,j)=T_R\hfill \\ u_x(N_x,j)=0.0\hfill \end{array}\end{array}

on boundary “c”,

\begin{array}{}\text{(31)}& \{\begin{array}{c}T(N_x+1,j)=T_R\hfill \\ u_x(N_x+1,j)=0.0\hfill \\ n(N_x+1,j)=n(N_x,j)\hfill \end{array}\end{array}

and on boundary “d”,

\begin{array}{}\text{(32)}& \{\begin{array}{c}T(2N_x,j)=T_L\hfill \\ u_x(2N_x,j)=0.0\hfill \\ n(2N_x,j)=n(1,j)\hfill \end{array}\end{array}

where j is the index of lattice node along the y direction.

It is noteworthy to point out that, the boundary conditions mentioned above, have the following advantages. Firstly, from the view point of numerical stability and accuracy, this approach makes physical quantities on the left boundary “a” equal to the ones on the right boundary “d”, therefore, the Gibbs oscillation can be avoided entirely; secondly, from the perspective of physics, the restraints on density in Eqs. (31) and (32) are equivalent to the adoption of the neutral wetting boundary conditions [36, 37, 56], i.e., {a\cdot \nabla n|}_{i=1,N_x}=0, where a is an inward normal unit vector to the boundaries. Moreover, on each boundary, u_x is set to be 0.0, then there is no density flux across the walls. So, in the region [1,N_x]\times [1,N_y] or [N_x+1,2N_x]\times [1,N_y], the mass conservation law is satisfied if the molecular diffusion is negligible.

Using the aforementioned boundary conditions, we conduct simulations on phase separation, where the system is kept in contact with walls at fixed temperatures. The initial conditions are

\begin{array}{}\text{(33)}& (n,T,u)|i,j,t=0=\{\begin{array}{cc}(n_L,0.9,0.0),& i⩽N_x\\ (n_R,0.9,0.0),& i>N_x\end{array}\end{array}

where n_L(i,j)=1.042+\Delta, n_R=n_L(2N_x-i+1,j), 1.042 is the mean density of liquid and vapor at T=0.9, \Delta is a random density with an amplitude of 0.001. Here we consider the following two cases: (I) T_L=T_R=0.9; (II) T_L=0.9, T_R=0.85. Figures 11 and 12 illustrate the density patterns at representative times for cases I and II, respectively. Similar to work of Gonnella et al. [36, 37, 41, 56], we mention that the phase separation process first occurs near the walls [see Fig. 11(a)], where the temperature is kept lower than the critical value. At about t=0.8 the spinodal decomposition pattern can be observed in the bulk. After that the bubbles (droplets) increase their sizes by the coalescence of relatively small low (high) density domains. Moreover, the left and right parts of the computational domain evolve in a symmetric way and we only need to focus on the behaviors of either of the two parts.

From Fig. 12, we find that the bubbles (droplets), especially the ones near the colder walls, e.g., walls “c” and “d”, grow much faster than that in Fig. 11 owing to the larger surface tension at lower temperatures. Most droplets are absorbed by the colder walls, and then in Fig. 12(d), we find a thicker liquid layer at the middle part of the computational domain, and two thinner liquid layers on the left and right boundaries, respectively. It is necessary to emphasize that, though this approach can cancel this limit, it doubles the computational load.

4 Conclusions

In this paper three open issues on the FFT–TLB model are re-addressed. Compared to the typical CD counterpart with the same order of accuracy, the FFT scheme has higher precision and better numerical stability, has more advantages to guarantee energy conservation and to refrain spurious velocities. The smoothing factors influence the accuracy and stability of this scheme by controlling the strength of high-frequency waves. Higher (Lower)-order smoothing factors owns higher (lower) accuracy, less (more) dissipation, but worse (better) stability. The smoothing factors should be minimal but make the evolution stable. The VDW EOS with flexible parameters makes possible simulating multiphase flows with various critical densities and temperatures. When solid walls are present, via making a mirror reflection of the original one around the wall, the FFT scheme can still be used. If the gravitational force is introduced further, the model can be used to study the hydrodynamic instabilities, e.g., Rayleigh–Taylor instability (RTI), competitions between RTI and Kelvin–Helmholtz instability [57–60].

References

[1] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases, London: Cambridge University Press, 1970

[2] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, New York: Oxford University Press, 2001

[3] Z. Guo and C. Zheng, Theory and Applications of Lattice Boltzmann Method, Beijing: Science Press, 2009 (in Chinese)

[4] A. Xu, Phys. Rev. E, 2005, 71(6): 066706

[5] A. Xu, Europhys. Lett., 2005, 69(2): 214

[6] Y. Gan, A. Xu, G. Zhang, X. Yu, and Y. Li, Physica A, 2008, 387(8–9): 1721

[7] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E, 2008, 77(3): 036707

[8] B. Shi and Z. Guo, Phys. Rev. E, 2009, 79(1): 016701

[9] Q. Li, Y. L. He, Y. Wang, and W. Q. Tao, Phys. Rev. E, 2007, 76(5): 056705

[10] Q. Li, Y. L. He, G. H. Tang, and W. Q. Tao, Phys. Rev. E, 2009, 80(3): 037702

[11] Q. Li, Y. L. He, G. H. Tang, and W. Q. Tao, Microfluid. Nanofluid., 2011, 10(3): 607

[12] Q. Li, K. H. Luo, Y. L. He, Y. J. Gao, and W. Q. Tao, Phys. Rev. E, 2012, 85(1): 016710

[13] H. Lai and C. Ma, J. Stat. Mech.: Theory Exp., 2010, 2010(4): P04011

[14] H. Lai and C. Ma, Phys. Rev. E, 2011, 84(4): 046708

[15] H. Li, X. Lu, H. Fang, and Y. Qian, Phys. Rev. E, 2004, 70(2): 026701

[16] B. Wen, H. Li, C. Zhang, and H. Fang, Phys. Rev. E, 2012, 85(1): 016704

[17] D. H. Rothman and J. M. Keller, J. Stat. Phys., 1988, 52(3–4): 1119

[18] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Phys. Rev. A, 1991, 43(8): 4320

[19] X. Shan and H. Chen, Phys. Rev. E, 1993, 47(3): 1815

[20] X. Shan and H. Chen, Phys. Rev. E, 1994, 49(4): 2941

[21] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett., 1995, 75(5): 830

[22] W. R. Osborn, E. Orlandini, M. R. Swift, J. M. Yeomans, and J. R. Banavar, Phys. Rev. Lett., 1995, 75(22): 4031

[23] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. E, 2003, 67(5): 056105

[24] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2004, 331(1–2): 10

[25] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2004, 344(3–4): 750

[26] A. Xu, G. Gonnella, A. Lamura, G. Amati, and F. Massaioli, Europhys. Lett., 2005, 71(4): 651

[27] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. E, 2006, 74(1): 011505

[28] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2006, 362(1): 42

[29] F. Corberi, G. Gonnella, and A. Lamura, Phys. Rev. Lett., 1998, 81(18): 3852

[30] A. Tiribocchi, N. Stella, G. Gonnella, and A. Lamura, Phys. Rev. E, 2009, 80(2): 026701

[31] V. Sofonea and K. R. Mecke, Eur. Phys. J. B, 1999, 8(1): 99

[32] V. Sofonea, A. Lamura, G. Gonnella, and A. Cristea, Phys. Rev. E, 2004, 70(4): 046702

[33] A. Cristea, G. Gonnella, A. Lamura, and V. Sofonea, Comm. Comp. Phys., 2010, 7(2): 350

[34] R. Zhang and H. Chen, Phys. Rev. E, 2003, 67(6): 066711

[35] P. Yuan and L. Schaefer, J. Fluid Eng., 2006, 128(1): 142

[36] G. Gonnella, A. Lamura, A. Piscitelli, and A. Tiribocchi, Phys. Rev. E, 2010, 82(4): 046302

[37] G. Gonnella, A. Lamura, and A. Tiribocchi, Phil. Trans. R. Soc. A, 2011, 369(1945): 2592

[38] A. Márkus and G. Házi, Phys. Rev. E, 2011, 83(4): 046705

[39] M. Sbragaglia, R. Benzi, L. Biferale, X. Shan, H. Chen, and S. Succi, J. Fluid Mech., 2009, 628: 299

[40] T. Seta, K. Kono, and S. Chen, Int. J. Mod. Phys. B, 2003, 17(1–2): 169

[41] G. Gonnella, A. Lamura, and V. Sofonea, Phys. Rev. E, 2007, 76(3): 036703

[42] Y. Gan, A. Xu, G. Zhang, and Y. Li, Phys. Rev. E, 2011, 84(4): 046715

[43] Y. Gan, A. Xu, G. Zhang, P. Zhang, and Y. Li, Europhys. Lett., 2012, 97(4): 44002

[44] H. Huang, D. T. Thorne, M. G. Schaap, and M. C. Sukop, Phys. Rev. E, 2007, 76(6): 066701

[45] H. Huang, M. Krafczyk, and X. Lu, Phys. Rev. E, 2011, 84(4): 046710

[46] L. Zheng, B. Shi, and Z. Guo, Phys. Rev. E, 2008, 78(2): 026705

[47] M. Watari and M. Tsutahara, Phys. Rev. E, 2003, 67(3): 036306

[48] A. Onuki, Phys. Rev. Lett., 2005, 94(5): 054501

[49] A. Onuki, Phys. Rev. E, 2007, 75(3): 036304

[50] F. Chen, A. Xu, G. Zhang, Y. Li, and S. Succi, Europhys. Lett., 2010, 90(5): 54003

[51] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, London: Springer-Verlag, 1987

[52] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation, Bristol: Adam Hilger, 1991

[53] J. P. Boyd, Chebyshev and Fourier Spectral Methods, New York: Dover Publications, 2000

[54] S. Orszag, Phys. Rev. Lett., 1971, 26(18): 1100

[55] H. Hadwiger, Math. Z., 1959, 71: 124

[56] G. Gonnella1, A. Lamura, and A. Piscitelli, J. Phys. A, 2008, 41(10): 105001

[57] Y. Gan, A. Xu, G. Zhang, and Y. Li, Phys. Rev. E, 2011, 83(5): 056704

[58] L. F. Wang, W. H. Ye, Z. F. Fan, and Y. J. Li, Europhys. Lett., 2010, 90(1): 15001

[59] W. H. Ye, L. F. Wang, and X. T. He, Phys. Plasmas, 2010, 17(12): 122704

[60] W. H. Ye, L. F. Wang, C. Xue, Z. F. Fan, and X. T. He, Phys. Plasmas, 2011, 18(2): 022704

Rights & permissions

© Higher Education Press and Springer-Verlag Berlin Heidelberg 2012