# Capillary filling using Lattice Boltzmann Equations: the case of multi-phase flows

###### Abstract

We present a systematic study of capillary filling for multi-phase flows by using mesoscopic lattice Boltzmann models describing a diffusive interface moving at a given contact angle with respect to the walls. We compare the numerical results at changing the density ratio between liquid and gas phases, and the ratio, , between the typical size of the capillary, , and the interface width, . It is shown that numerical results yield quantitative agreement with the Washburn law when both ratios are large, i./e. as the hydrodynamic limit of a infinitely thin interface is approached. We also show that in the initial stage of the filling process, transient behaviour induced by inertial effects and “vena contracta” mechanisms, may induce significant departure from the Washburn law. Both effects are under control in our lattice Boltzmann equation and in good agreement with the phenomenology of capillary filling.

###### pacs:

83.50.Rp, and 68.03.Cd^{†}

^{†}offprints: biferale@roma2.infn.it

## 1 Introduction

The physics of capillary filling is an old problem, originating with the pioneering works of Washburn washburn () and Lucas Lucas (). Recently, with the explosion of theoretical, experimental and numerical works on microphysics and nanophysics, the problem attracted a renewed interest degennes (); dussain (); washburn_rec (); tas (). Capillary filling is a typical “contact line” problem, where the subtle non-hydrodynamic effects taking place at the contact point between liquid-gas and solid phase allows the interface to move, pulled by capillary forces and contrasted by viscous forces. Usually, only the late asymptotic stage is studied, leading to the well-known Washburn law, which predicts the following relation for the position of the interface inside the capillary:

(1) |

where is the surface tension between liquid and gas, is the static contact angle, is the liquid viscosity, is the channel height and the factor depends on the geometry of the channel (here a two dimensional geometry given by two infinite parallel plates separated by a distance – see fig. 1). The above expression can be recasted in dimensionless variables and , being the capillary time . This leads to the universal law:

(2) |

As already remarked in many works washburn_rec (), the asymptotic behaviour (2) is obtained under the assumption that (i) the inertial terms in the Navier-Stokes evolution are negligible, (ii) the instantaneous bulk profile is given by the Poiseuille flow, (iii) the microscopic slip mechanism which allow for the movement of the interface is not relevant for bulk quantities (as the overall position of the interface inside the channel), (iv) inlet and outlet phenomena can be neglected (limit of infinitely long channels); (v) the liquid is filling in a capillary, either empty or filled with gas whose total mass is negligible with respect to the liquid one. In the following, we will address all these effects and show to which extent they can described by using a mesoscopic model for multiphase flows based on the discretized version of Boltzmann Equations in a lattice. The model here used is a suitable adaptation of the Shan-Chen pseudo-potential LBE shanchen () with hydrophobic/hydrophilic boundaries conditions, as developed in prenoi1 (); prenoi2 (). Other models with different boundary conditions and/or non-ideal interactions have been also used in kopo ().

## 2 LBE for capillary filling

The geometry we are going to investigate is depicted in fig. (1). The bottom and top surface is coated only in the rigth half of the channel with a boundary condition imposing a given static contact angle prenoi1 (); in the left half we impose periodic boundary conditions at top and bottom surface in order to have a flat liquid-gas interface which should mimic an “infinite reservoir”. Periodic boundary conditions are also imposed at the two lateral sides such as to ensure total conservation of mass inside the system.

### 2.1 LBE algorithm for multi-phase flows

We start from the usual lattice Boltzmann equation with a single-time relaxation Gladrow (); Saurobook ():

(3) |

where is the kinetic probability density function associated with a mesoscopic velocity , is a mean collision time (with a time lapse), the equilibrium distribution, corresponding to the Maxwellian distribution in the continuum limit. From the kinetic distributions we can define macroscopic density and momentum fields as Gladrow (); Saurobook ():

(4) |

For technical details and numerical simulations we shall refer to the nine-speed, two-dimensional model Gladrow (). The equilibrium distribution in the lattice Boltzmann equations is obtained via a low Mach number expansion of the equilibrium Maxwellian Gladrow (); Saurobook (). In order to study non-ideal effects we need to supplement the previous description with an interparticle forcing. This is done by adding a suitable in (3). In the original model shanchen (), the bulk interparticle interaction is proportional to a free parameter (the ratio of potential to thermal energy), , entering the equation for the momentum balance:

(5) |

being the static weights for the standard case of 2DQ9 Gladrow () and the pseudo-potential function which describes the fluid-fluid interactions triggered by inhomogeneities of the density profile (see shanchen (); prenoi1 (); prenoi2 () for details).

One may show prenoi1 (); prenoi2 () that the above pseudo-potential, leads to a non-ideal pressure tensor given by (upon Taylor expanding the forcing term):

(6) | |||||

where
is the sound speed.
This approach allows the definition of a static contact angle , by introducing at the walls a
suitable value for the pseudo-potential prenoi1 (),
which can span the range .
Moreover, it also defines a specific value for the surface tension,
, via the usual integration of the offset between normal and transverse
components of the pressure tensor along the liquid-gas interface
allows forshanchen (); prenoi1 (); prenoi2 ().

As to the boundary conditions on the Boltzmann populations,
the standard bounce-back rule is imposed. One can show that the
bounce-back rule gives no-slip boundary conditions up to second order
in the Knudsen number
in the hydrodynamical limit of single phase flows noijfm ().
In presence of strong density variation, close to the walls across the
interface, the velocity parallel to the wall may develop a small slip
length (of the order of the interface thickness, ) which in turn, allows for the interface to move. It is difficult to control exactly the phenomenon,
because even imposing an exact no-slip boundary conditions at the wall yeomans (), the model will develop non
trivial dynamics at the first node away from the wall, where both condensation/evaporation phenomena and/or spurious currents may conspire, leading to
an overall non-zero slip velocity. For the scope of controlling the capillary
filling, one may reabsorb all these effects within the usual Maxwell slip
boundary conditions: .
It is easy to show that in presence of a slip velocity, the Poiseuille profile becomes:

(7) |

where the velocity of the front must be identified with the mean velocity, . Therefore, Washburn law (2) becomes:

(8) |

### 2.2 Corrections to the Wahsburn law

As already remarked many yeas ago bousanquet (), the Wahsburn law (2) can be valid only if inertial forces are negligible with respect to the viscous and capillary ones. This cannot be true at the beginning of the filling process, where strong acceleration drives the interface inside the capillary. However, putting reasonable numbers for a typical microdevices experiments with water (, , , ), one realizes that the transient time, , is usually very small, of the order of a few nanoseconds, and therefore negligible for most experimental purposes. Another important effect which must be kept in mind when trying to simulate capillary filling, is the unavoidable “resistance” of the gas occupying the capillary during the liquid invasion. This is a particular “sensitive” question, because reaching the typical density ratio between liquid and gas of experimental set up represents a challenge for most numerical methods, particularly for multiphase Lattice Boltzmann which typically operates with density ratios of the order or . In order to take in to account both effects, inertia and gas dynamics, one may write down the balance between the total momentum change inside the capillary and the force acting on the liquidgas system:

(9) |

where is the total mass of liquid and gas inside the capillary at any given time. The two forces in the right hand side correspond to the capillary force, , and to the viscous drag . Following the notation of fig.(1) and the expression for the velocity profile (7) one obtains the final expression (see also napoli () for a similar derivation, without considering the slip velocity):

(10) |

In the above equation for the front dynamics, the terms in the LHS take into account the fluid inertia. Being proportional either to the acceleration or to the squared velocity, they become negligible for long times. Washburn law plus the slip correction (8) is therefore correctly recovered asymptotically, for , and in the limit when . The above equation is exact, in the case where evaporation-condensation effects are negligible, i.e. when the gas is pushed out of the capillary without any interaction with the liquid. This is not always the case for most of the mesoscopic models available in the literature shanchen (); yeomans (), based on a diffusive interface dynamics jacqmin (). As we will see, only when either the limit of thin interface is reached or when the gas phase is negligible, , the dynamics given by (10) is correctly recovered. Otherwise, deviations are observed, which are induced by condensation/evaporation effects, which may result in significant departure from the Poiseuille profile inside the gas phase.

### 2.3 Numerical Results

In fig. (2) we show the behaviour of the front position, , as a function of time for a given contact angle (), a given density ration and a given surface tension (in LBE units), at varying the channel width , from up to . As one can see the numerical results tends to be in good agreement with the solutions of (10), only for large enough values of , i.e. only when the interface becomes thin enough. For small to moderate values of , the overall asymptotic trend is only qualitatively reproduced by the Washburn law (plus slip effects) (8) with deviations which may be of the order of in the prefactor (see inset of the same figure).

Similarly, increasing the surface tension and the factor, leads to a early convergence towards the asymptotic Washburn law and to the solution of (10) even for small channel width . This is shown in fig. (3) where, at fixed , we increase and the Washburn law is approached better and better.

From both figure (2-3) one may notice that for short filling time, , strong deviations from the Washburn law are detected, even considering the extra effects induced by the inertial terms of (10). This slowing down at the early stage of the filling is indeed mainly induced by a sort of vena contracta term venacontracta (), reflecting the non-trivial matching between the reservoir and the capillary dynamics at the inlet. This term, has been argued to be describable by an additive apparent-mass correction, , to the LHS of (10).

In Fig. (4) we show an enlargement of fig(2) for small filling time superposed with the results of a numerical integration of eq. (10) with a phenomenological value for the vena contracta factor. As one can see, the agreement between the numerics and the evolution of (10) is now excellent also at short times.

## 3 Conclusions

The present study shows that Lattice Boltzmann models with pseudo-potential energy interactions are capable of reproducing the basic features of capillary filling, as described within the Washburn approximation. Two conditions for quantitative agreement have been identified: i) a sufficiently high density contrast between the dense/light phase, and a sufficiently thin interface, . Both conditions can be met within the current LB methodology, although it would clearly be desirable to extend the LB scheme in such a way to achieve density contrasts in the order of (the current state-of-the-art is approximately ) and interface widths of the order of the lattice spacing (current values are about ).

The present results set the stage for future computational studies aimed at identifying optimal interface functionalization strategies, based on physical, chemical and geometrical coating processes. Work along these lines is currently underway.

## 4 Acknowledgments

Valuable discussions with E. Paganini, D. Palmieri and D. Pisignano are kindly acknowledged. Work performed under the EC contract NMP3-CT-2006-031980 (INFLUS).

## References

- (1) E.W. Washburn, Phys. Rev. 17 (1921) 273.
- (2) R. Lucas, Kooloid-Z 23 (1918) 15.
- (3) P.G. de Gennes, Rev. Mod. Phys. 57 (1985) 827.
- (4) E.B. Dussan, Ann. Rev. Fluid Mech. 11 (1979) 11.
- (5) L.J. Yang, T.J. Yao and Y.C. Tai, J. Micromech. Microeng. 14 (2004) 220.
- (6) N.R. Tas et al., Appl. Phys. Lett. 85 (2004) 3274.
- (7) X. Shan and H. Chen, Phys. Rev E. 47 (1993) 1815.
- (8) R. Benzi, L. Biferale, M. Sbragaglia, S. Succi and F. Toschi, Phys. Rev. E 74 (2006) 021509.
- (9) M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama and F. Toschi, Phys. Rev. E 75 (2007) 026702.
- (10) P. Raiskinmaki, A. Shakib-Manesh, A. Jasberg, A. Koponen, J. Merikoski and J. Timonen, J. Stat. Phys. 107 (2002) 143.
- (11) D. Wolf-Gladrow, Lattice-Gas Cellular Automata And Lattice Boltzmann Models (Springer, New York, 2000)
- (12) S. Succi, The lattice Boltzmann Equation, Oxford Science (2001).
- (13) R. Benzi, L. Biferale, M. Sbragaglia, S. Succi, and F. Toschi, Jour. Fluid Mech. 548, (2006) 257.
- (14) A. J. Briant, A. J. Wagner, and J. M. Yeomans, Phys. Rev. E 69, (2004) 031602.
- (15) C.H. Bosanquet, Philos. Mag. 45 (1923) 525.
- (16) G. Cavaccini, V. Pianese, A. Jannelli, S. Iacono and R. Fazio, Lecture Series Comp. Comput. Sci. 7 (2006) 66.
- (17) J. Szkeley, A.W. Neumann and Y.K. Chuang, J. Colloid Interf. Sci. 69 (1979) 486.
- (18) D. Jacqmin, Jour. Fluid Mech. 402 (2000) 57.