## Abstract

A three-dimensional Lattice Boltzmann two-phase model capable of dealing with large liquid and gas density ratios and with a partial wetting surface is introduced. This is based on a high density ratio model combined with a partial wetting boundary method. The predicted three-dimensional droplets at different partial wetting conditions at equilibrium are in good agreement with analytical solutions. Despite the large density ratio, the spurious velocity around the interface is not substantial, and is rather insensitive to the examined liquid and gas density and viscosity ratios. The influence of the gravitational force on the droplet shape is also examined through the variations of the Bond number, where the droplet shape migrates from spherical to flattened interface in tandem with the increase of the Bond number. The predicted interfaces under constant Bond number are also validated against measurements with good agreements.

## 1. Introduction

The Lattice Boltzmann method (LBM) [1] has achieved considerable success in simulating hydrodynamic problems, and its major advantages are that it is explicit, easy to implement and natural to parallelize. Thus, it is natural to extend the LBM to simulate multi-phase flows [2–4], and there are many multi-phase industrial applications, for example, bio-chip [5] and electro-wetting [6] based tunable micro lens system [7]. The correct simulation of these devices depends on the LBM’s capability to tackle high-phase density ratios and wetting boundaries.

The pseudo-potential by Shan & Chen [2], for example, has been widely employed to compute two-phase flows owing to its simplicity in implementation. However, its undesirable feature is the generation of spurious velocity along the phase interface and the limitation of the simulated density ratio. Both of these deficiencies can be alleviated by adopting different equation of states [8] or by adopting an extended pseudo-potential method [9].

Another popular alternative is the adoption of the free energy, where the evolution of the interface is captured using the Cahn–Hilliard equation. Similar to the pseudo-potential method, high density ratio simulation can cause stability problems, and this issue has been vigorously pursued, aiming at reducing the spurious velocity along the phase interface [10–13]. Jacqmin [10] indicated that the potential form of the surface tension can improve the elimination of the spurious velocity. On the other hand, Lee & Fischer [13] suggested that higher order finite differencing of the derivative operator is crucial to correct the balance of the interfacial force, and hence the elimination of the spurious velocity. In addition, adopting the potential form, Zheng *et al.* [14] proposed an LBM formulation that can recover correctly to the Cahn–Hilliard equation with the capability of computing large density ratio flows.

As indicated earlier, surface wettability control is crucial to the operation of the bio-chip and tunable liquid lens. By minimizing the surface free energy, Briant *et al.* [15] proposed boundary conditions to control the surface wetting behaviour to simulate contact line motion in a liquid–gas system. Surface partial wetting results were shown to agree with analytical theory, albeit with a low density ratio. To model a high density ratio with a partial wetting surface, Yan & Zu [16] proposed an LBM that is based on the existing models of Inamuro *et al.* [12] and Briant *et al.* [15]. However, the model of Inamuro *et al.* [12] requires a pressure correction to enforce the continuity at each collision-streaming step.

Here, a three-dimensional LBM methodology is adopted to simulate the liquid–gas system. This is based on the high-density model of Zheng *et al.* [14] combined with the partial wetting boundary method of Briant *et al.* [15]. High density ratio flows are used to scrutinize the capability of the present method through the predicted static droplet contact angle at different partial wetting conditions. The interaction of the gravity and surface tension on the droplet shape is also investigated. Analytical results and available measurements are used to examine the capability of the methodology.

## 2. Theory and method

For the present binary fluid simulations, the governing equations of the Navier–Stokes and Cahn–Hilliard equations can be written as [10,14]
2.1
2.2
2.3
where *n*=(*ρ*_{A}+*ρ*_{B})/2, *ϕ*=(*ρ*_{A}−*ρ*_{B})/2, *μ*_{ϕ} is chemical potential, *θ*_{M} is mobility, (*c*_{s} is the speed of sound) and *F*_{b} is the body force. *ρ*_{A} and *ρ*_{B} are the densities of fluid A and B, respectively, which can be high- and low-density fluids, *ρ*_{H} and *ρ*_{L}.

The chemical potential can be derived from the free-energy function [10,17]
2.4
where *κ* is a coefficient that is related to the surface tension and the thickness of the interface layer. *ψ*(*ϕ*) is the bulk free-energy density and is chosen as a double-well form [10,14]
2.5
where *A* is an amplitude parameter to control the interaction energy between the two phases and *ϕ**=(*ρ*_{H}−*ρ*_{L})/2. This form represents two equilibrium states, *ϕ** and −*ϕ** when the bulk free-energy density is zero.

Thus, the chemical potential depending on the free-energy function can be derived as
2.6
The surface tension *σ* and the interface width *W* are related to the *A* and *κ*, i.e. *A*=3*σ*/(4*Wϕ**^{4}) and *κ*=3*Wσ*/(8*ϕ**^{2}).

Solutions of the Navier–Stokes and Cahn–Hilliard equations are achieved using the Lattice Boltzmann equations [14,18]
2.7
and
2.8
and is defined as [19]
2.9
The discrete velocity *e*_{i} is defined as
2.10

The adoption of *q*=1/(*τ*_{ϕ}+0.5) [14] in equation (2.8) [18] would ensure the correct recovery of equation (2.3). The equilibrium distribution functions adopted for D3Q19 and D3Q7 models [14,20] are, respectively,
2.11
and
2.12
where the weighting factors *w*_{i} are *w*_{0}=1/3, *w*_{1,…,6}=1/18 and *w*_{7,…,18}=1/36. *Γ* is the mobility coefficient used to control the mobility *θ*_{M}=*q*(*qτ*_{ϕ}−0.5)Δ*tΓ*.

The distribution functions satisfy the conservation laws, 2.13 2.14 2.15

According to Young’s equation [21], the angle between the liquid–gas interface and the wall can be expressed as
2.16
where *σ*, *σ*_{sl} and *σ*_{sg} are the liquid–gas, solid–liquid and solid–gas surface tensions, respectively. By minimizing the surface free energy [15,22], the contact angle can be related to the gradient of the surface order parameter, i.e.
2.17
where and
2.18

The above equation can be used to impose the partial wetting conditions near the wall through the gradient of the order parameter and is incorporated in the equilibrium distribution function as suggested by Briant *et al.* [15].

Following Cheng *et al.* [20] to ensure a better numerical stability, the discretizations of ∇*ϕ* and ∇^{2}*ϕ* in equations (2.9) and (2.6) are adopted,
2.19
and
2.20

Near the wetting boundary, the discretization of ∇*ϕ* and ∇^{2}*ϕ* follow that proposed by Briant *et al.* [15]. Also, the no-slip boundary conditions in equation (2.14) are implemented according to those in Ho *et al.* [23] and Chang *et al.* [24].

## 3. Results

For all the simulations, the liquid and gas densities are, respectively, 1000 and 1, and have equal viscosity unless otherwise specified. The lattice size is 115×115×76. *τ*_{n} and *τ*_{ϕ} are 0.6 and 0.7 and surface tension *σ* is 0.1.

Attention here is first directed to validate the correct implementation of the wetting boundary conditions by comparing the predicted contact angles with the analytical solutions of equation (2.17). Simulations focus on droplets resting on surfaces with zero gravitational force, and figure 1 shows the predicted three-dimensional droplets at different contact angles with the solid surfaces at equilibrium. Further examination of the predicted angles compared with the analytical solutions shown in figure 2 indicates that the boundary conditions are correctly implemented.

Table 1 shows the maximum spurious velocity of stationary droplets at different density and viscosity ratios under neutral wetting boundary conditions. For the cases examined, the spurious velocity is of the order of 10^{−6}, which may not be substantial and is rather insensitive to the examined liquid and gas density and viscosity ratios. However, this also shows that the present discretization may generate a slight, albeit low, force imbalance of **∇**(*P*_{0}+*ϕμ*_{ϕ})=*μ*_{ϕ}**∇***ϕ* in equation (2.2). Thus, Lee & Fischer’s discretization methodology [13] needs to be further explored to eliminate this velocity completely.

The influence of the gravitational force on the droplet shape is also examined. Bond number (*Bo*) represents the interaction of the gravity and surface tension, i.e.
3.1
where *ρ*_{H}−*ρ*_{L} stands for the density difference between fluids and *g* is the downward acceleration of gravity. *L* is the characteristic length, and is denoted as the initial spherical droplet diameter.

A neutral wetting boundary condition is applied, which means that the contact angle is 90^{°}. The Bond numbers investigated are *Bo*=0,1,10 and 100. Figure 3 shows the contour view of the order parameter (*ϕ*), where the influence of gravity is clearly observed at *Bo*>1. Detailed examination of the influences of the Bond number is referred to in the interface profiles at different Bond numbers, where flattened interfaces are predicted at higher Bond number.

The variation of the contact angle can be used to manipulate the droplet shape through electro-wetting, where the contact angle on the surface is controllable and thus has a tunable focal length. Alternatively, electro-wetting can also be applied as a lens fabrication technique, where the shape of the lens is controlled electrostatically.

Figure 4 shows fabricated lens profiles under different applied voltages (0, 40 and 100 V) and hence contact angles (90^{°}, 79^{°} and 64^{°}) experimented by Hsieh & Chen [25]. The influence of the gravitational force is not negligible, as can be observed from the different height and radius values for the neutral wetting case, i.e. *θ*_{w}=90^{°}. The present predictions agree quite well with the measured data, indicating the capability of the present LBM model.

## 4. Conclusion

In the present study, a three-dimensional LBM is applied to compute binary fluid systems with large density ratios (1000) and partial wetting surfaces. This is based on the high-density model of Zheng *et al.* [14] combined with the partial wetting boundary method of Briant *et al.* [15]. The predicted three-dimensional droplets at different contact angles with the solid surfaces at equilibrium are in good agreement with analytical solutions. The influence of the gravitational force on the droplet shape is also examined through the variations of the Bond number (*Bo*) ranging from 0 to 100, where perceivable influence is observed for *Bo*>1, with the migration from spherical to flattened interfaces in tandem with the increase of the Bond number. The predicted interface under constant Bond number is also validated against measurements, with good agreement. On the other hand, despite the large density ratio, the maximum spurious velocity around the interface is not substantial (approx. 10^{−6}) and is rather insensitive to the examined liquid and gas density and viscosity ratios. However, Lee & Fischer’s discretization methodology [13] still needs to be further explored to eliminate this velocity. In addition, the applicability of the present methodology to simulate the dynamic wetting behaviour within the bio-chip and tunable micro lens systems should be further investigated.

## Acknowledgements

The project is supported by NSC 95-2221-E-007 -227 and NCHC of Taiwan.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society