1 Introduction

Wetting phenomena are widespread in nature and industrial applications^[1-4]. The essence of wetting phenomena is the complex balance between adhesion and cohesion forces. In general, the wettability of a solid surface is characterized by several macroscopic parameters^[5], such as the static contact angle, the advancing angle, the receding angle, and the roll-off angle. The systems concerned with wetting phenomena are typical multicomponent/multiphase complex fluid systems. Due to the sophisticated interactions involved in those systems, simulating the behavior of multicomponent/multiphase fluid systems is challenging. Meanwhile, it is also difficult to select a suitable and computationally efficient method for such systems.

Over the past two decades, the lattice Boltzmann method (LBM) has been under a rapid development and has evolved into a flexible and convenient numerical tool for simulating various flow behaviors^[6-7]. Distinguished from conventional computational fluid dynamics (CFD), LBM is based on the mesoscopic kinetic equations of particle distribution functions. Its main strength lies in the fact that the evolution progress of the particle distribution functions behaves like a solver for the macroscopic Navier-Stokes equations, yet its mesoscopic nature endows LBM with a viable way to incorporate the microscopic dynamics of the simulated systems. Due to such characteristics, LBM has emerged as a versatile method for investigating a series of complex multiphase/multi-component fluid systems, especially for systems involving the evolution of the fluid-fluid interfaces or systems concerning the interaction between fluid and solid particles.

At present, the following four types of LBM multiphase models are widely used in academic research: the R-K model^[8-10], the S-C model^[11-12], the free energy (FE) model^[13-14], and the models based on the phase-field theory^[15-18]. The R-K model is the first LBM multiphase model from a historical perspective^[19]. In this model, two sets of particle distribution functions are introduced to characterize two fluid components. The dynamics of fluid interfaces and the interactions between fluid components are determined by the color gradients. During each time evolution, a recoloring process is required, which costs tremendous computational resources. The S-C model is also known as the pseudo-potential model, as it introduces the pseudo-potential force between different components/phases at mesoscopic scale to mimic the microscopic intermolecular interactions. Phase separation occurs automatically, and no interface capturing technique is demanded. Although the S-C model is quite simple, it successfully captures the essence of non-ideal fluid behavior simulation. However, the S-C model neglects the thermodynamic properties of the system to a certain extent. Therefore, the FE model was proposed by Swift. Its major idea is to design the equilibrium particle distribution function according to the thermodynamics, aiming at taking the thermodynamic properties of the system into consideration. The main drawback of the FE model is that the viscous term in the corresponding macroscopic equation does not satisfy the Galilean invariance, and efforts were made to fix that issue^[20-21]. The models based on the phase-field theory have attracted extensive attention in recent years in multiphase/multicomponent system simulation. Such models usually use an independent variable (index function) to characterize the transition between different states or phases so as to track the evolution of interfaces. In general, the evolution of the index function follows the Cahn-Hilliard equation with convection terms, and the hydrodynamics of the system is described by the Navier-Stokes equations with additional surface tension terms. Essentially, the multiphase models based on the phase-field theory are much closer to designing specific numerical strategies for solving the Cahn-Hilliard equation within the LBM framework. Unlike other models based on the phase-field theory, the H-C-Z model proposed by He et al.^[22-23] originates from the kinetic theory for non-ideal gases and dense fluids. As indicated by Li et al.^[24], the H-C-Z model can be considered as the first phase-field theory based on the LBM multiphase model. Compared with the S-C model, the H-C-Z model is more physically sound. In addition, it is easy to incorporate various equations of state and adjust the surface tension independent of the density ratio^[25].

Both the S-C model and the H-C-Z model have played important roles in scientific research, yet they exhibit different characteristics. The simplicity of the S-C multicomponent/multiphase model lies in that simulating a multiphase/multicomponent fluid system can be implemented through directly adding the pseudo-potential force to the evolution equation of LBM as an external force term. Meanwhile, the pseudo-potential force that mimics the microscopic intermolecular interactions endows the S-C model with a clear physical meaning. Nevertheless, this model has suffered from some limitations and incurred certain questioning in some applications. These limitations can be summarized as the following four points^[26]: the existence of large spurious currents near phase interfaces, inconsistent with thermodynamics, limited density and viscous ratio, and the coupling between the equation of state (EOS) and the surface tension. Specific to the above four points, a number of techniques have been developed to improve the S-C model. For instance, Shan^[27] indicated that spurious current originates from the limited discretization in calculating the interaction force and can be greatly depressed through using schemes with higher order of isotropy. Besides, several techniques, including incorporating the realistic EOS into the S-C model^[28], introducing a multi-range interaction force^[29], and utilizing a multi-relaxation-time (MRT) collision operator^[30], all achieved significant improvements. The H-C-Z model is the discretized form of the discrete Boltzmann equation (DBE) proposed by He et al.^[31] for non-ideal gases and dense fluids. Following the idea of Chapman and Cowling, the volume exclusion effect and intermolecular attraction between molecules in dense fluids are expressed as a function of fluid density, and then these two effects are fully considered by introducing additional terms into the collision operator. In order to improve the numerical stability and accuracy, the H-C-Z model adopts a pressure evolution method and introduces the pressure distribution function instead of the density distribution function commonly used in LBM (i.e., the particle distribution function). He et al.^[22-23] investigated the complex flows of two phase interfaces induced by the Rayleigh-Taylor instability and Kelvin-Helmholtz instability, and the reported results agreed well with theoretical prediction. Lee and Lin^[32] further developed the H-C-Z model to achieve higher density ratio through a stable discretization scheme. They even indicated that a combination of the intermolecular force in the potential form and an appropriate isotropic discretization scheme is capable of eliminating spurious currents^[33-34].

Abundant research work has been made using the above two models^{[18, 35-38]}. As the original S-C model fails to adjust the surface tension while maintaining the density ratio of two phases, it turns out to restrict the selection of reasonable numerical parameters for practical simulations. In contrast, the H-C-Z model benefits from an adjustable surface tension, which in turn leads to an advantage in choosing parameters that can help to circumvent numerical instability. To compare these two multiphase models, the modified S-C model proposed by Kuzmin^[39] was applied, which is capable of separately adjusting the surface tension and the density ratio of two phases. Meanwhile, Yuan's strategy^[28] was utilized to incorporate realistic EOS and to improve the numerical stability as well. Simulations of the static contact angle on solid surfaces were then carried out using the two models, and in particular the wetting boundary condition was implemented with the geometric formulation proposed by Ding and Spelt^[40]. The maximum magnitude of the spurious currents (MMSC) and the convergence rate of the two models are compared in detail. The results of the current work can be viewed as useful references or supplementary when one simulates wetting phenomena using LBM multiphase models along with the geometric formulation.

2 Methodology 2.1 LBM

LBM is essentially a discrete form of the Boltzmann equation, which describes the macroscopic fluid motion through the evolution of a set of particle distribution functions in a phase space. The standard lattice Boltzmann (LB) equation with the external force term can be expressed as follows^[6]:

(1)

where f_i(x, t), e_i, and f_i^eq(x, t) represent the particle distribution function, the discrete velocity, and the equilibrium distribution function in the ith direction of the discrete velocity space, respectively. δt is the time step. τ is a single relaxation parameter that controls the viscosity of macroscopic fluid. The last term F_iδt in Eq. (1) relates to the external force, and its specific form depends on the chosen forcing scheme. In this work, we adopt the forcing scheme developed by He et al.^[31], which indicates

where a is the external force acceleration. In practice, Eq. (1) is split into a collision step and a streaming step as

(2)

(3)

The D2Q9 discrete velocity model^[41] is adopted as

(4)

where c is the lattice speed, defined as the ratio of the lattice space δx to the time step δt. The equilibrium distribution function f_i^eq takes the following form:

(5)

where

and , c_s is the lattice sound speed. The density and velocity of the fluid can be obtained from the moments of the particle distribution functions:

(6)

The macroscopic kinematic viscosity ν is related to the single relaxation parameter τ by

2.2 The S-C multiphase model

To mimic a single component system, the interaction force in the S-C multiphase model generally yields the following form^[11-12]:

(7)

where G is a parameter reflecting the interaction strength, and ψ(x, t) is the pseudo-potential function. According to the Taylor expansion, we have

(8)

leading to a total pressure tensor expressed as

(9)

For the bulk fluid far from the interface region, the above interaction force indicates a non-ideal EOS,

(10)

As inferred from Eq. (9), the surface tension and the density ratio simultaneously rely on G. To solve that issue, Kuzmin^[39] modified the S-C model through including the next-nearest neighborhood interaction as shown in Fig. 1.

Relatively low computational costs are demanded to accomplish separately tuning the surface tension and the density ratio. Meanwhile, Kuzmin utilized a specified set of weight coefficients to cancel the second term in Eq. (8), and a Navier-Stokes compliant surface tension term F_s=κρ▽△ρ was introduced into the interaction force. The weight coefficients can be expressed in a matrix form according to the corresponding discrete velocity as

(11)

Taylor expanding Eq.(8) and substituting the weight coefficients with Eq. (11), one has

(12)

Note that, to increase the maximum achievable density ratio and suppress the magnitude of the spurious currents, Kuzmin adopted the multi-relaxation collision operator and an extended equilibrium function, which further aggravated the computational load. To achieve the above effects, the strategy proposed by Yuan and Schaefer^[28] is applied in this work. At the same time, Yuan's strategy is also capable of incorporating various realistic EOS into the S-C model via the following pseudo-potential function:

(13)

Combining Eq. (7) and Eq. (13), one can notice that G no longer controls the strength of fluid-fluid interaction, and thus it is kept as -1 in the whole simulations to prevent the term inside the square root from being negative.

2.3 The H-C-Z multiphase model

As mentioned in the introduction, the H-C-Z multiphase model adopts two particle distribution functions to monitor the evolution of the pressure and the interface. The details of the H-C-Z multiphase model are briefly reviewed in this subsection. The evolution equations for the H-C-Z multiphase model read as follows^[22-23]:

(14)

(15)

where g_i is the distribution function for the pressure, f_i is the distribution function for the index function, Γ_i(u)=f_i^eq/ρ, ψ(ρ)=p-ρRT, R is the gas constant, and T is the temperature.

F_s represents the force associated with the surface tension controlled by the parameter κ, and F_e stands for the external force. The equilibrium function g_i^eq takes the following form:

(16)

Similarly, in practice, Eq. (14) and Eq. (15) are implemented through a collision step and a streaming step. Macroscopic fluid properties can be calculated via

(17)

The density and the viscosity linearly relate to the index function as

(18)

where ρ_l and ρ_h represent the lower density and higher density of the two phases, and ν_l and ν_h represent the lower viscosity and higher viscosity of the two phases, respectively. The gradient and Laplace operators are discretized according to Lee and Lin^[32],

(19)

(20)

where ζ represents an arbitrary physical variable.

2.4 Wetting boundary condition and numerical settings

In the framework of LBM, the effects of the solid surface usually can be taken into consideration in a strategy analogous to account for the interaction between fluid components. The wetting boundary condition can be sorted into two categories depending on how the density value of the solid lattice node is determined: assigning a fixed effective density to the solid lattice node or assigning a tunable value according to the morphology of contact line (i.e., the geometric formulation) to the solid lattice node. In this work, the latter one is adopted, and it takes the following discrete form^{[36, 40]}:

(21)

where ϕ represents the density or the index function, the first subscript i represents the coordinate along the solid surface while the second subscript represents the coordinate perpendicular to the solid surface, and θ_pre is the prescribed contact angle. Meanwhile, the half-way bounce-back scheme is applied when fluid particles encounter solid nodes. The surface tension term F_s requires an extra layer of ghost nodes when calculating the nearest fluid nodes next to solid nodes. The densities of those ghost nodes take the same values as those of the solid nodes.

The realistic EOS used in this work is the Carnahan-Starling EOS,

(22)

where a=4, b=4, and RT=1/3. According to Maxwell's construction, the theoretical values of the liquid and vapor phase are ρ_liq=0.251, ρ_gas=0.024 3. Other numerical parameters for both models are summarized in Table 1.

The relative residual errors of density and velocity between two consecutive time steps are defined as follows:

(23)

(24)

The whole system is assumed to reach a steady state after E₁ and E₂ are less than 10^-6.

3 Results and discussion 3.1 Laplace test

Surface tension is critical for multicomponent/multiphase fluid systems, as it significantly affects the behavior of the interface. It is necessary to figure out the relationship between the surface tension σ and the parameter κ beforehand. In this subsection, we turn to the Laplace law to obtain the relationships between σ and κ for both models. According to the Laplace law, for a two-dimensional case, the surface tension relates to the pressure difference inside and outside the droplet as δp=σ/R, where δp denotes the pressure difference, and R is the droplet radius at equilibrium. The size of the computational domain is 100×100, and periodic boundary conditions are utilized in both horizontal and vertical directions. A droplet with several initial radii (12, 15, 20, 25, 30 in lattice units (l.u.)) is deposited in the middle of the computational domain. In Fig. 2, the relationship between σ and κ for both models is illustrated, and the detailed results are summarized in Table 2.

Through the linear fitting and the cubic polynomial fitting, the following relationship between σ and κ for both models can be acquired:

(25)

(26)

It is necessary to point out that Eq. (25) and Eq. (26) merely reflect the numerical relationship between σ and κ rather than to indicate that σ in the two models depicts a different physical background. For the present S-C model, the higher order truncation error in Eq. (12) affects the magnitude of the surface tension. Also, it is worth mentioning that the shape of the liquid droplet simulated with the S-C model no longer keeps circular when the value of κ is below 0.2, and the Laplace law is not capable of estimating the surface tension under this circumstance. This unphysical phenomenon might originate from the spurious currents.

3.2 Numerical simulation of static contact angles

For liquid droplets spreading on the solid wall, interfacial tensions are in balance at the contact point of three phases when the equilibrium state is reached. To simulate such a problem, a computational domain with the size of 200×100 is used, and a semi-circular droplet with an initial radius of 25 l.u.~is deposited at the bottom of the computational domain. Periodic conditions are applied to the left and right boundaries. The top and bottom boundaries are solid surfaces. The values of κ in the S-C model and the H-C-Z model are set as 0.3 and 0.435 637, respectively, corresponding to an equal surface tension σ=4.813 31×10^-3. Snapshots of the static contact angle after reaching the equilibrium state for both models are shown in Fig. 3. It can be seen from Fig. 3 that, for both models, a wide range of static contact angles from wetting to non-wetting can be achieved via tuning the prescribed value. However, numerical instability occurs when the prescribed contact angle is too small.

As indicated by Hu et al.^[42], the reason responsible for the numerical instability can be concluded as below: the absolute value of tan(-θ) can be very large when θ gets close to 0 or π. Consequently, the model may be unstable due to the rapidly changing interaction force near the solid boundary. For the parameters used in this work, the minimum achievable contact angle obtained by the H-C-Z model is about 25º, while the S-C model can only obtain static contact angles larger than 33º. When the prescribed contact angle is larger than 160º, the droplets simulated by both models have been completely detached from the solid wall, deviated from the prescribed contact angle. The measured values of the static contact angle are summarized in Table 3. It is found that the static contact angles calculated by the two models are in good agreement with the prescribed ones. The standard deviations of the H-C-Z model and the S-C model are 1.53º and 1.71º, respectively. The maximum deviation approximately equals 2.8º.

3.3 Spurious currents

Spurious currents refer to the small-amplitude artificial structured velocity field arising near phase interfaces and pervade almost all LB models tending to simulate multiphase/multi-component systems. It is also one of the major factors that restrict the maximum achievable density ratio. As the magnitude of the spurious currents increases, the numerical stability of LBM will face severe challenges. Under certain circumstances, the spurious currents cannot even be distinguished from the real flow velocities, leading to certain questioning of the reliability of LBM. The velocity field of a single droplet (with an initial radius of 20 l.u.) immersed in its own vapor phase is shown in Fig. 4(a). It can be observed that the droplet reaching equilibrium is surrounded by eight structured vortices. A similar pattern is captured as well when simulating the static contact angle as depicted in Fig. 4(b) which exhibits the velocity field of a droplet deposited on a solid surface with a prescribed contact angle of 90º after reaching equilibrium.

During the Laplace test, we also monitor the relationship between the MMSC and a liquid droplet (with an initial radius of 25 l.u.) as depicted in Fig. 5. From Fig. 5, it can be clearly noticed that the MMSC obtained by the two models shows a contrary dependency on the value of κ. For the H-C-Z model, the MMSC increases linearly as κ increases, while for the S-C model increasing the value of κ can depress the MMSC. When implementing the static contact angle with the geometric formulation, the MMSC is also affected by the prescribed contact angle. The relationships between the MMSC of the two models and the value of the prescribed contact angle are illustrated in Fig. 6.

From the curves in Fig. 6, it can be acknowledged that for both models the value of the prescribed contact angle will lead to small fluctuations on the MMSC. Especially, when the value of θ_pre is relatively low, the MMSC tends to dramatically increase. Combining the results from Fig. 5 and Fig. 6, it can be deduced that the surface tension is the major affecting factor for determining the MMSC obtained by the two multiphase models while the prescribed contact angle has a minor impact. Also under the current value of surface tension, the MMSC of the S-C model is universally higher than that of the H-C-Z model.

3.4 Comparison of convergence rate

It should be noted that, as a natural dynamic scheme, LBM is not a method of choice for steady-state computations^[43]. Nevertheless, LBM is still a powerful and reliable tool for complex multicomponent or multiphase fluid systems due to its advantages. Thus, attention needs to be paid to the convergence rate of the chosen computational model when simulating problems in which a steady state will finally be approached. In this subsection, we compare the convergence rate of the two LBM multiphase models in detail. The steady state is assumed to be achieved after the convergence criteria, i.e., Eq. (23) and Eq. (24) are satisfied. A liquid droplet is initially deposited at the bottom of the computational domain. The density field is initialized in such a manner that the corresponding contact angle at the three phases has been set as the prescribed contact angle. Under this configuration, a rapid convergence shall be observed, and the interference of droplet spreading procedure can be alleviated. First of all, we monitor the convergence procedure of the S-C model and the H-C-Z model under three values of the prescribed contact angle (θ_pre=45º, ~90º, ~135º) using computational domains with three different sizes (200×100, ~400×200, ~600×300). The total time steps for the two models to reach the steady state are summarized in Table 4.

In general, as the size of computational domain increases, more time steps are required for the system to reach to the steady state. It can be acknowledged that the S-C model exhibits an overall faster convergence rate than the H-C-Z model, regardless of the size of the computational domain. We then move on to discuss the major factor that distinguishes the convergence rate of the two models. In Fig. 7, the relative residual error diagrams for θ_pre=90º are shown, and it can be found that the total time steps mainly depend on the relative residual error of velocity. In Fig. 8, corresponding to different prescribed contact angles, the convergence procedure of the two models is shown while the size of the computational domain is fixed as 200×100.

As the value of the prescribed contact angle decreases, the total time steps required for both models to reach the steady state gradually increase. In particular, the convergence rate of the H-C-Z model is dramatically affected by the relatively low prescribed contact angle. In order to determine the possible causes of such slow convergence, we monitor the pressure field and the density field obtained with the H-C-Z model for θ_pre=45º. It is found that the global magnitude of the pressure field is gradually increasing as the system approaching the steady state. Moreover, although the relative residual error of density continuously decreases, the size of the droplet is unexpectedly increasing. A possible explanation for the above problem is given by analyzing the collision step as well as the methods of calculating macroscopic properties in the H-C-Z model. The collision step of the H-C-Z model can be expressed as follows:

(27)

Summing Eq. (27) over all discrete velocity directions, meanwhile according to the definition of pressure given by Eq. (17), we have

(28)

Summing Eq. (28) over all fluid nodes, the following relationship is deduced:

(29)

The left-hand side of Eq. (29) is actually the summation of the global pressure field for the time step t+δt. Note that the streaming step, the periodic boundary condition and the bounce-back condition will not lead to increments in the pressure field. Thus, it can be noticed that a non-physical term will be introduced into the pressure field during each time step. This might be responsible for the slow convergence rate of the H-C-Z model when the value of θ_pre is relatively small.

4 Conclusions

In this paper, wetting phenomena on solid surfaces are investigated by two LBM multiphase models (i.e., the S-C model and the H-C-Z model). Specifically, the static contact angle is implemented with the geometric formulation, and a detailed comparison of the computational performance of the two multiphase models is made. The main conclusions can be listed as follows:

(ⅰ) Based on the geometric formulation, a wide range of static contact angles from wetting to non-wetting can be achieved for both the S-C model and the H-C-Z model. The results show that the achievable ranges of the H-C-Z model and the S-C model are 25º-160º and 33º-160º, respectively. This finding may be useful for researches concerning the effects of wettability on multiphase flow characteristics in complex porous media.

(ⅱ) The MMSC of the two models mainly depends on the value of the surface tension (i.e., the controlling parameter κ). Under the numerical parameters adopted in this paper, the MMSC of the H-C-Z model is lower than that of the S-C model. Yet, they are on the same order of magnitude, which indicates that both models can reflect the nature of wetting phenomena on solid surfaces. At the same time, it is noticed that the value of the prescribed contact angle will lead to small fluctuations on the MMSC.

(ⅲ) The convergence rate of the S-C model is superior to that of the H-C-Z model, especially for simulating surfaces with relatively low contact angles. Thus, when simulating problems in which a steady state exists, the S-C model is a suitable choice with an appropriate convergence rate.