Implementation of contact angles in the pseudopotential lattice Boltzmann simulations with curved boundaries

The pseudopotential multiphase lattice Boltzmann (LB) model is a very popular model in the LB community for simulating multiphase flows. When the multiphase modeling involves a solid boundary, a numerical scheme is required to simulate the contact angle at the solid boundary. In this work, we aim at investigating the implementation of contact angles in the pseudopotential LB simulations with curved boundaries. In the pseudopotential LB model, the contact angle is usually realized by employing a solid-fluid interaction or specifying a constant virtual wall density. However, it is shown that the solid-fluid interaction scheme yields very large spurious currents in the simulations involving curved boundaries, while the virtual-density scheme produces an unphysical thick mass-transfer layer near the solid boundary although it gives much smaller spurious currents. We also extend the geometric-formulation scheme in the phase-field method to the pseudopotential LB model. Nevertheless, in comparison with the solid-fluid interaction scheme and the virtual-density scheme, the geometric-formulation scheme is relatively difficult to implement for curved boundaries and cannot be directly applied to three-dimensional space. By analyzing the features of these three schemes, we propose an improved virtual-density scheme to implement contact angles in the pseudopotential LB simulations with curved boundaries, which does not suffer from a thick mass-transfer layer near the solid boundary and retains the advantages of the original virtual-density scheme, i.e., simplicity, easiness for implementation, and low spurious currents.


I. Introduction
The lattice Boltzmann (LB) method has been developed into an efficient numerical methodology for simulating fluid flow and heat transfer in the past three decades [1][2][3][4][5][6][7][8]. Owing to its kinetic nature, the LB method has exhibited some distinct advantages over conventional numerical methods and has been widely used in modeling multiphase flows and interfacial phenomena. The existing multiphase LB models can be generally classified into four categories [1][2][3], i.e., the color-gradient LB model, the pseudopotential LB model, the free-energy LB model, and the phase-field LB model. Among these four categories, the pseudopotential LB model [9][10][11] is probably the simplest one. In this model, the intermolecular interactions are represented with an interaction force based on a density-dependent pseudopotential and the phase separation is naturally achieved by imposing a short-range attraction between different phases.
Historically, the first attempt of using the pseudopotential LB model to simulate wetting phenomena was made by Martys and Chen [12], who proposed a solid-fluid interaction scheme to describe the interaction between a fluid phase and a solid wall. Different contact angles were obtained by adjusting the interaction strength of the solid-fluid interaction. Another type of solid-fluid interactions was later developed by Raiskinmäki et al. [13,14]. In their scheme, the pseudopotential serves as a pre-sum factor, while in the solid-fluid interaction scheme of Martys and Chen the pre-sum factor is the density. Kang et al. [15,16] have also formulated a solid-fluid interaction scheme for the pseudopotential LB model and investigated the displacement of immiscible droplets subject to gravitational forces in a two-dimensional channel and a three-dimensional duct. Moreover, based on the work of Martys and Chen, Colosqui et al. [17] have proposed a modified solid-fluid interaction scheme composed of a repulsive core and an attractive tail.
According to the mechanical equilibrium of a multiphase system in the presence of a boundary condition, Benzi et al. [18] derived a formula for the contact angle of the pseudopotential LB model and presented an alternative treatment to implement wetting boundaries. They introduced a virtual wall density w  to fix the pseudopotential at a solid wall. By tuning w  from l  (density of liquid phase) to g  (density of gas phase), the contact angle in simulations can be varied from 0° to 180°. A similar scheme can also be found in the color-gradient multiphase LB model [19], which is called the fictitious-density scheme [20]. However, as shown in Ref. [20], the fictitious-density scheme leads to an unphysical thick mass-transfer layer near the solid boundary. Such a phenomenon can also be observed in the pseudopotential LB simulations using the virtual-density scheme [21].
Besides the aforementioned studies, Huang et al. [22] have investigated the wetting boundaries in the pseudopotential multi-component LB simulations and proposed a formula to determine the adhesion parameters of different components from the contact angle. In addition, the geometric-formulation scheme, which is proposed by Ding and Spelt [23] for the phase-field method, has also been employed to implement contact angles in the pseudopotential LB simulations involving flat surfaces [24,25].
Compared with the solid-fluid interaction scheme, the geometric-formulation scheme usually yields much smaller spurious currents. Moreover, it can give a slope of the liquid-gas interface that is consistent with the prescribed value of the contact angle. However, such a scheme is mainly applicable to flat surfaces and its implementation for curved boundaries is much more complicated [26] than that of the solid-fluid interaction scheme or the virtual-density scheme.
In the present work, we aim at investigating the implementation of contact angles in the pseudopotential LB simulations with curved boundaries. An improved virtual-density scheme is proposed, which retains the basic advantages of the original virtual-density scheme but does not suffer from a thick mass-transfer layer near the solid boundary. Meanwhile, it yields much smaller spurious currents than the solid-fluid interaction scheme and is easy to implement in both two-dimensional and three-dimensional space in comparison with the geometric-formulation scheme. The rest of the present paper is organized as follows. The pseudopotential multiphase LB model and the solid-fluid interaction where the constant  is utilized to realize the thermodynamic consistency [28]. For three-dimensional models (e.g., the D3Q15 and D3Q19 lattice models), readers are referred to Refs. [32,33,39].

B. Solid-fluid interaction scheme and virtual-density scheme
The intermolecular interaction force defined by Eq. (5) represents the cohesive force of a system.
When a solid wall is encountered, an adhesive force should also be considered [22]. In order to describe the interaction between a fluid and a solid wall, Martys and Chen [12] proposed the following solid-fluid interaction to mimic the adhesive force in the pseudopotential LB model: where w G is the adhesive parameter and   t s    x e is a switch function, which is equal to 1 or 0 for a solid or fluid phase, respectively. By adjusting the value of w G , different contact angles can be realized. Besides Eq. (8), some other types of solid-fluid interactions can be found in Ref. [40].
The treatment or scheme that uses a virtual density was developed by Benzi et al. [18], who introduced a constant virtual density w  to fix the pseudopotential of the solid phase, i.e.,   w   .
Then Eq. (5) can also be applied to the interaction between the fluid phase and the solid phase. Similarly, different contact angles can be obtained by tuning the value of w  . When w  varies from l  to g  , the contact angle is tuned from 0 to 180° [21]. The advantages of the virtual-density scheme lie in its simplicity and easiness for implementation, but some previous studies showed that such a scheme usually produces an unphysical mass-transfer layer near the solid boundary [7,21].

A. Curved geometric-formulation scheme
In 2007, Ding and Spelt [23] proposed a geometric-formulation scheme to implement wetting boundaries in the phase-field method. For a two-dimensional flat surface, the geometric-formulation scheme is given by ,0 , 2 a 1, 1 1, 1 tan 2 where C is the order parameter of the phase-field method, a  is an analytically prescribed contact angle, and ,0 i C is the order parameter at the ghost layer    [23] showed that the geometric-formulation scheme can give a slope of the liquid-gas interface that is consistent with the prescribed value of the contact angle.
However, Eq. (9) is only applicable to flat surfaces [24,25]. Recently, Liu and Ding [26] devised a geometric-formulation scheme for two-dimensional phase-field simulations with curved surfaces, which is also referred to as "the characteristic moving contact-line model". They considered a ghost contact-line region inside the solid phase, as illustrated in Fig In the phase-field method, the unit normal vector of the solid surface is calculated by [26] S s S where S C is the order parameter of the solid phase [26]. Since there is no such a quantity in the pseudopotential LB model, s n is evaluated as follows: where the switch function   is the same as that in Eq. (8). To improve the numerical accuracy, a high-order isotropic discretization scheme can be used to evaluate s n , such as the 8th-order isotropic scheme proposed by Sbragaglia et al. [38,41]:    When s n is determined, the unit vectors along the characteristic lines 1 l and 2 l can be obtained by the following vector rotation: where 2       . According to the unit vectors 1 n and 2 n , the intersection points 1 D and 2 D can be identified. Usually, different cases will be encountered when varying the contact angle. Figure 2 gives an example of the intersection point 2 D when the contact angle  in Fig. 1 is changed.
Obviously, the implementation of the geometric-formulation scheme is much more complex than that of the solid-fluid interaction scheme or the virtual-density scheme. More details about the determination of the points D 1 and D 2 can be found in Ref. [26].

FIG. 2. Illustration of the intersection point 2 D for different contact angles.
After identifying the intersection points 1 D and 2 D , the densities at these two points can be obtained by an interpolation of the densities at their neighboring lattice points. A quadric interpolation was used in the study of Liu and Ding [26], which involves three neighboring points around 1 D or 2 D .
Without loss of generality, one can also employ a linear interpolation. With the densities of the points 1 D and 2 D , the density at the point P can be determined by Eq. (11), and then the pseudopotential can be calculated by Eq. (6). Similar to the virtual-density scheme, the curved geometric-formulation scheme also applies Eq. (5) to the interaction between a fluid phase and a solid phase.

B. Improved virtual-density scheme
The advantage of the geometric-formulation scheme lies in that it is able to make the liquid-gas interface intersect a solid boundary at an angle in consistence with the prescribed contact angle. On the contrary, when employing the solid-fluid interaction scheme or the virtual-density scheme, we should adjust the value of w G or w  in simulations so as to achieve a required contact angle. However, as can be seen in the previous section, the implementation of the geometric-formulation scheme is very complicated in comparison with the solid-fluid interaction and virtual-density schemes. Moreover, the above curved geometric-formulation scheme cannot be directly applied to three-dimensional space due to the fact that in two-dimensional space there are only two possible characteristic lines making an angle  with s n (as shown in Fig. 1), but in three-dimensional space the characteristic lines that make an angle  with s n form a circular cone surface around s n [20]. Hence in this section we devise an improved contact angle scheme for the pseudopotential LB model, which is easy to implement in both two-dimensional and three-dimensional space.
Actually, in the geometric-formulation scheme the density at a solid point is also a virtual density, but the virtual density in the solid phase is a local quantity instead of a constant for the whole solid domain, which implies that the drawback of the original virtual-density scheme may be overcome when a local virtual density is employed. On the basis of such a consideration, we propose the following formula for the virtual density in the solid phase near a curved boundary:   (16) is larger than l  , and is taken as g  when it is smaller than g  .
We now explain why we choose  is initially placed on the circular cylinder with its center at (150, 230). The periodic boundary condition is applied in the x and y directions and the halfway bounce-back scheme [6,8,43] is used to treat the curved solid boundary, which is illustrated in Fig. 3. The kinematic viscosity is taken as   Fig. 6(c)) can be observed in the right-hand panel of Fig. 7. Moreover, it is clearly seen that the improved virtual-density scheme performs much better than the virtual-density scheme since the density variations in the results of the improved virtual-density scheme are significantly smaller than those of the virtual-density scheme. Figure 8 shows the static contact angles obtained by the geometric-formulation scheme. Some slight deviations are observed between the numerically obtained contact angles and the analytically prescribed contact angles given in Eq. (15), which may arise from the use of a linear interpolation in our simulations. Figure 9 compares the spurious currents produced by the solid-fluid interaction scheme at To quantify the numerical results, a comparison of the maximum spurious currents yielded by the four schemes is made in Fig. 10, from which we can find that the maximum spurious currents are in the order of 0.1 for the solid-fluid interaction scheme but are smaller than 0.006 for other schemes. As previously mentioned, in the geometric-formulation scheme the density within the solid phase is also a virtual density. Hence the results in Fig. 10 indicate that applying the intermolecular interaction force Eq.
(5) to the interaction between a fluid phase and a solid phase with a virtual density is better than using a solid-fluid interaction force in light of reducing the spurious currents. Moreover, Fig. 10 also shows that the maximum spurious currents yielded by the virtual-density scheme are larger than those given by the geometric-formulation scheme and the improved virtual-density scheme, which implies that the spurious currents may be further reduced by replacing a constant virtual density with a local virtual density.   can also be found in the pseudopotential LB simulations of contact angles on straight solid surfaces [40].

B. Effects of the thick mass-transfer layer
The influence of the spurious currents has been well studied in the literature. Hence in the present work we mainly reveal the adverse effects of the thick mass-transfer layer near the solid boundary caused by the virtual-density scheme. Firstly, we employ the test of Poiseuille flow between two parallel solid plates to analyze the effects of the thick mass-transfer layer. The distance between the two plates is   Fig. 13(a), although it is a little thinner than the mass-transfer layer of the case Fig. 4(a). Due to the unphysical mass-transfer layer, at 100 t t   the droplet in Fig. 13(a) has contacted the solid circular cylinder, which indicates that the three-phase contact line (reduces to contact points in 2D) appears earlier in the simulation using the virtual-density scheme.
Owing to the influences of the unphysical mass-transfer layer, the numerical results predicted by the virtual-density scheme gradually deviate from the results obtained by the geometric-formulation scheme, which can be found by comparing Fig. 13(a) with Fig. 13(b). For example, the three-phase contact points at 4000 t t   in Fig. 13(a) are much closer to the central vertical line ( 2 x x N  ) of the domain than those in Fig. 13(b). Moreover, significant deviations can be observed between the results of the virtual-density scheme and the geometric-formulation scheme at 10000 t t   . Contrarily, the improved virtual-density scheme is shown to be capable of producing numerical results consistent with those given by the geometric-formulation scheme. , 300 t  , 4000 t  , and 10000 t  .

C. Contact angles on a spherical surface
Finally, the capability of the improved virtual-density scheme for simulating three-dimensional contact angles is validated by the test of static contact angles on a spherical surface. The D3Q19 pseudopotential MRT-LB model proposed in Ref. [33] is adopted in our simulations and the lattice Initially, a solid sphere of radius 50 R  is located at (100, 100, 100) and a droplet of 45 r  is placed on the spherical surface with its center at (100, 100, 180). The periodic boundary condition is applied in all the directions and the halfway bounce-back scheme [6,8,43] is employed to treat the curved boundary. Other treatments such as the equation of state and the coexisting liquid and gas densities of the two-phase system are the same as those used in the above two-dimensional tests. Figure 14 presents the results of different three-dimensional contact angles obtained by the improved virtual-density scheme, in which the lower row displays the density contours of the x-z cross-section at 100 y  . The results clearly demonstrate that the improved virtual-density scheme is capable of modeling three-dimensional contact angles on a curved surface and does not suffer from a thick mass-transfer layer near the solid boundary, which exists in the simulations using the virtual-density scheme.

V. Summary
We have investigated the implementation of contact angles in the pseudopotential LB simulations involving curved boundaries. The solid-fluid interaction scheme and the virtual-density scheme, which are two popular schemes for the pseudopotential LB modeling of wetting phenomena, are shown to suffer from very large spurious currents and an unphysical thick mass-transfer layer near the solid boundary, respectively. A curved geometric-formulation scheme for the pseudopotential LB model has been extended from a recently developed contact angle scheme for two-dimensional phase-field simulations. Although the geometric-formulation scheme can give a slope of the liquid-gas interface that is basically consistent with the prescribed contact angle, it is rather difficult to implement (e.g., for moving solid particles) and cannot be directly applied to three-dimensional space.
Hence we have proposed an improved virtual-density scheme, which employs a local virtual density to replace the constant virtual density and therefore overcomes the drawback of the original virtual-density scheme. Meanwhile, the spurious currents produced by the improved virtual-density scheme are much smaller than those caused by the solid-fluid interaction scheme and it is much easier to implement in both two-dimensional and three-dimensional space as compared with the geometric-formulation scheme. The features of the improved virtual-density scheme have been well demonstrated by simulating contact angles on cylindrical and spherical surfaces. For simplicity, the halfway bounce-back scheme [6,8,43] is employed in the present work to treat the curved solid boundaries. In the LB community, there have been many curved boundary schemes for curved boundaries, such as the scheme proposed by Mei et al. [44], the interpolated bounce-back scheme [45], and the single-node curved boundary scheme [46]. However, it should be noted that these curved boundary scheme usually suffer from severe mass leakage in two-phase LB simulations [47].