[
Abstract
Droplets on hydrophobic surfaces are ubiquitous in microfluidic applications and there exists a number of commonly used multicomponent and multiphase lattice Boltzmann schemes to study such systems. In this paper we focus on a popular implementation of a multicomponent model as introduced by Shan and Chen. Here, interactions between different components are implemented as repulsive forces whose strength is determined by model parameters. In this paper we present simulations of a droplet on a hydrophobic surface. We investigate the dependence of the contact angle on the simulation parameters and quantitatively compare different approaches to determine it. Results show that the method is capable of modelling the whole range of contact angles. We find that the a priori determination of the contact angle is depending on the simulation parameters with an uncertainty of 10 to 20%.
Contact angles in multicomponent LB models]Contact angle determination in multicomponent lattice Boltzmann simulations
S. Schmieschek, J. Harting]Sebastian Schmieschek and Jens Harting^{1}^{1}1Corresponding author.

[

[

Institute for Computational Physics, University of Stuttgart,
Pfaffenwaldring 27, 70569 Stuttgart, Germany
Department of Applied Physics, Eindhoven University of Technology,
P.O. Box 513, 5600 MB Eindhoven, The Netherlands.

PACS: 47.55.D, 47.11.j

Key words: lattice Boltzmann, ShanChen model, contact angle, droplets, hydrophobic surface
1 Introduction
During the last few decades the miniaturization of technical devices down to submicrometric sizes has made considerable progress. In particular, during the 1980s, socalled microelectromechanical systems (MEMS) became available for chemical, biological and technical applications leading to the rise of the discipline called “microfluidics” in the 1990s [1]. In microfluidic devices the surface to volume ratio of a fluid can be large and thus a good understanding of the behavior of the fluid close to the surface is mandatory. However, the behavior of a fluid close to a solid interface is very complex and involves the interplay of many physical and chemical properties. These include the wettability of the solid, the shear rate or flow velocity, the bulk pressure, the surface charge, the surface roughness, as well as impurities and dissolved gas.
A common concept to quantify the wettability of a surface is the so called contact angle. The contact angle is the angle at which the interface between a liquid and a gas or vapor meets a solid surface. If the contact angle is larger than 90, the surface is called nonwettable (hydrophobic if the liquid is water) and if the angle is smaller than 90, it is said to be wettable (hydrophilic). Superhydrophobic surfaces are surfaces with contact angles larger than 150. Here, almost no contact between droplet and surfaces can be observed and the effect is often referred to as “Lotus effect”. Regardless of the amount of wetting, the shape of the drop can be approximated by a truncated sphere.
For a droplet on an idealised smooth surface, the contact angle can be computed using the surface tensions between liquid and gas , liquid and surface and surface and gas as given by Young’s equation [2] (see Fig. 1),
(1.1) 
The model of Young was extended by Wenzel [3] as well as Cassie and Baxter [4] in order to take the influence of surface roughness into account. While Wenzel describes a state where the surface is completely covered by the liquid, Cassie and Baxter describe a state where gas bubbles are enclosed between the liquid and the rough surface. Both states have been observed experimentally and in simulations [5, 6]. The transition between the Wenzel and the CassieBaxter state leads to the phenomenon of contact angle hysteresis as observed for droplets on a tilted surface where one has to distinguish between the advancing and the receding contact angle [7, 8, 9]. In particular the state proposed by Cassie and Baxter is of technological interest since it can be used to significantly increase the contact angle in order to generate superhydrophobic surfaces with [10, 11, 12]. Such surfaces can be utilized to increase the flow velocity and thus the mass flux in microchannels [13, 14].
While both molecular dynamics and lattice Boltzmann methods (LBM) have been employed to simulate systems with wetting properties, only LBM allow to reach experimentally relevant time and length scales. Therefore, the method has become very popular to simulate typical problems occurring in microfluidics. A particular advantage of the lattice Boltzmann approach is the availability of established multiphase or multicomponent methods [15, 16, 17, 18, 19] and a straight forward implementation of complex boundary conditions. This allows the simulation of multiphase or multicomponent fluid flow along interacting surfaces [20, 21, 22, 14]. While the free energy based multiphase model introduced by Swift et al. [17] allows to set the contact angle directly, this possibility does not exist for the model introduced by Shan and Chen. Here, the surface tension and thus the contact angle only appear indirectly by tuning the interaction between different fluid species and the surface [15, 16]. Therefore, a proper determination of the contact angle is of fundamental importance for reliable comparisons between simulation results and those obtained from theory and experiment. For the single component multiphase ShanChen model, Benzi et al. proposed an analytical ansatz to compute the contact angle [23]. However, in this paper we focus on the multicomponent model [15] and restrict ourselves to single phase components only. For such a model, Huang et al. [24] recently proposed an estimate determining the contact angle. However, a full analytical solution of the problem is still missing. Therefore, we compare and discuss different methods to quantify in dependence on the parameters of the simulation model, namely the geometrical measurement, the approach of Huang et al., as well as utilizing measurements of the surface tension to solve Young’s equation.
2 Simulation method
A set of equations can be used to represent a standard lattice Boltzmann system involving multiple species [25]
(2.1) 
with . The singleparticle velocity distribution function indicates the density of species , having velocity , at site on a Ddimensional lattice of coordination number , at timestep . The collision operator
(2.2) 
represents the change in the singleparticle distribution function due to the collisions. A popular form is the single relaxation time , linear ‘BGK’ form [26] for the collision operator. It can be shown for low Mach numbers that the LB equations correspond to a solution of the NavierStokes equation for isothermal, quasiincompressible fluid flow. The lattice Boltzmann method is an excellent candidate to exploit the possibilities of parallel computers, as the computations a lattice site require only information about quantities at nearest neighbour lattice sites [27, 28]. The local equilibrium distribution plays a fundamental role in the dynamics of the system as shown by Eq. (2.1). In this study, we use a purely kinetic approach, for which is derived by imposing certain restrictions on the microscopic processes, such as explicit mass and global momentum conservation [29]
(2.3) 
where is the fluid density and is the macroscopic bulk velocity of the fluid, given by . are the coefficients resulting from the velocity space discretization and is the speed of sound, both of which are determined by the choice of the lattice. We use a D3Q19 implementation, i.e., a three dimensional lattice with 19 discrete velocities. Immiscibility of species is introduced in the model following Shan and Chen [15, 16], where only nearest neighbour interactions among the species are considered. These interactions are described by a selfconsistently generated mean field body force
(2.4) 
where is the socalled effective mass, which can have a general form for modeling various types of fluids (we use )[15]), and is a force coupling constant whose magnitude controls the strength of the interaction between components , and is set positive to mimic repulsion. The symbol denotes the position of a nearest neighbour. The dynamical effect of the force is realized in the BGK collision operator by adding the increment
(2.5) 
to the velocity in the equilibrium distribution (Eq. (2.3)). This naturally opens the way to introduce similar interactions between each fluid species and the channel walls, where the strength of the interaction is determined by the fluid densities, free coupling constants, and a wall interaction parameter.
For the interaction of the fluid components with the channel walls we apply midgrid bounce back boundary conditions [30] and assign interaction properties to the wall which are similar to those of an additional fluid species, i.e. we specify constant values for the force coupling constant and the density for the rest vector (, ) at wall boundary nodes of the lattice. This results in a purely local force as given in Eq. 2.4 between the flow and the boundaries. Even though one could argue that a single parameter to tune the fluidwall interaction would be sufficient, we keep our approach as close as possible to the original idea of Shan and Chen in order to benefit from the experience obtained from other works using the original model. Furthermore, the additional parameter allows more flexibility to tune the interactions in a system more complex than considered here. The fluidwall interaction can be linked to a contact angle between fluid droplets and solid walls as it is often used to quantitatively describe hydrophobic interactions [31].
In the model used, the interface between domains of different fluid species has a finite width. In order to define a position of an interface we introduce the order parameter which is zero at the interface.
We perform simulations of a droplet at an interacting surface in order to investigate the influence of the droplet size, the pseudo wall density (wettability) , and the coupling constant on the resulting contact angle. The system is initialised with a spherical cap of component and density at a smooth surface. The drop is surrounded by a fluid of component and density .
This choice of densities is made without loss of generality. In the scope of this work only one coupling parameter is used. Introduction of a density contrast at initialisation therefore mainly results in a shift in droplet size and mean density. The equilibrium density contrast, however, is fixed by the Laplace law. To quantitatively describe a droplet of fluid in a gaseous medium, typically a contrast in dynamic viscosities of the order of needs to be modelled. This is well beyond the limit of numerical stability of the model employed. Despite this fact, as shown below, the phenomenological nature of the ShanChen force allows the qualitative modelling of the whole contact angle range.
At the surface midgrid bounce back boundary conditions as well as a repulsive force with pseudo wall density are applied.
3 Geometrical determination of the contact angle
Assuming a droplet has the shape of a spherical segment, the contact angle
(3.1) 
can be obtained by measuring the base , the height and the radius of the droplet (see Fig. 2). The geometrical measurement is used as a reference to compare to the approaches of contact angle determination further below. Base and height can be determined by measuring the position where the order parameter has a value of zero. The radius is then given by . Due to the fluidwall interaction there exists an interface layer in the vicinity of the wall. Determining the base by measurement of sign change of the order parameter immediately above the wall is therefore introducing an error. To avoid this, the droplet radius is calculated from the base and height relative to a reference point sufficiently far from the interface layer. For the simulation results discussed here a heightoffset of 5 lattice units proved to be sufficient. The correct base length above the wall can then be calculated from the sodetermined radius and the actual height.
4 Dependence of the contact angle on model parameters
The size of the simulated system has a strong influence on the precision of the results. For example, due to discretization effects, a droplet cannot be approximated by a sphere if the lattice resolution is too low. Further, calculations that take the curvature of the drop into account, also require a well resolved surface of the droplet.
We checked the dependence of the contact angle on the system size for , , , and lattices and initial droplet volumes of , , , and . An example system setup is shown in Fig. 3. The left side of Fig. 4 depicts the measured contact angle in a system of two immiscible components of equal density and kinematic viscosity for and . It can be seen that the contact angle increases with increasing absolute value of and that even for the largest system size the contact angle is not fully converged. The convergence depends on the wettability parameter : the stronger the repulsion between fluid and surface, the larger the system has to be. However, considering the doubling of initial droplet volume from to the relative change in the contact angle measured is as low as approximately for , for and for . Therefore, we find a compromise between optimal use of computing time and precision of the measurement and restrict ourselves to lattices of size , and . Figure 4 (right) shows the dependence of the contact angle on the wetting parameter for . The plot shows a linear dependence of the contact angle up to about and . In the vicinity of the complete dewetting limit, the dependence becomes nonlinear.
4.1 Surface tension measurements at planar interfaces
As given by Eq. 1.1, the contact angle can be calculated if the surface tensions between liquid and gas, liquid and surface, and gas and surface are known. Only the curvature of the interface between liquid and gas depends on the size of the droplet. By assuming an infinitely large droplet on a surface, the interface between liquid and gas can be approximated as planar and the surface tension can be calculated using its mechanical definition
(4.1) 
wherein the component of the pressure tensor normal to the interface is and the component transversal to the interface is . The pressure tensor is computed as
(4.2)  
Here, the first term is equivalent to the dynamic pressure. The second term describes the distribution of the mean field body force given by eq. 2.4.
For interfaces between liquid or gas and the surface, is being computed equivalently. As introduced in Ref. [32] a sized system with periodic boundaries is filled with two 64 lattice units long lamellae of different fluid components. The densities for both components are chosen as , is varied between and in steps of 0.02. For calculating the surface tension between a fluid component and the wall, half of the system is filled with a wall component with variable wetting parameter between and in steps of . is varied as for the fluidfluid case.
The surface tension obtained is being used to calculate the contact angle as given by Eq. 1.1. Figure 5 shows the deviation of the obtained values from the ones obtained by a geometrical determination of . For the geometrical measurements, a droplet with a volume of on a flat surface is used. It can be seen that the deviation is always positive and that the dependence of on the model parameters is stronger than for the geometric measurements. In fact, already and values for of cause the contact angle to reach , while for this value is not being reached for . For all simulations have produced contact angles of .
The significant differences between the geometrical determination of the contact angle and the measurements of the surface tension have a number of reasons: first, fluids diffuse into areas where the other component is the majority. Thus, in the droplet system, the volume covered by the droplet also includes up to 5% of the surrounding fluid component which has an influence on the measured surface tension. Further, since the pressure is tensorial at the interface only, merely seven discrete data points along one axis are used to calculate the surface tension. Enhancing the resolution of the interface would, however, increase the computational cost significantly. Nonetheless, the measurement might be improved by introducing better statistics interpolating over the whole droplet interface. The fluid components are slightly compressible leading to slightly different maximum and minimum values of the steady state densities of the droplet system and the planar setup for the surface tension determination. Further, the curvature of the interface is not being taken into account. In particular for small droplets, this effect has a significant influence. Therefore, we compare our results to measurements obtained using an equation for the surface tension that takes the droplet geometry into account:
(4.3) 
Here is the radius of the interface. We integrate from the center of the droplet (). The integral is evaluated until 5 l.u. before the border of the system in order to minimize any influences due to periodic boundary conditions.
The resulting contact angles are always smaller than the ones obtained from the geometrical determination (see Fig. 6). In particular for moderate values of we find strong deviations due to a higher curvature in the curve solving the YoungLaplace equation with the measured surface tensions. There is no linear dependence of the contact angle on the surface wettability as observed in the geometrical measurements.
In a recent publication, Huang et al. postulate an estimate for the contact angle within the multiphase multicomponent ShanChen model [24]. This estimate is valid for a fixed ratio of the component densities and their coupling constants. The approach of Huang et al. is based on the assumption that the surface tension at the wall is mainly determined by the local interaction. The force acting on component , where the boundary condition is given by an interacting surface, can be written as
For the components we have
with if there is a surface in the direction of motion and if not. In proportion to these forces, the surface tensions can be calculated in dependence on the density gradient as arithmetic average of minimum and maximum density.
From this we obtain the YoungLaplace law
(4.4) 
In our case we use the same coupling for fluidfluid and fluidsolid interactions. Therefore, this equation only depends on the density gradients. However, the dependence on the coupling parameter enters implicitly. Also for this method we compare the results to the geometrical measurement of the contact angle. As before, Fig. 7 shows the deviation of the contact angle (in percent) from the values observed from the reference measurement for an initial droplet size of . The deviation is proportional to the absolute value of the coupling and decreases for low already at a wettability of . This allows to assume a dependence on the interface thickness.
The deviations of Huang’s approach compared to the geometrical measurements are up to 15 percent. Since the validity was only postulated for a limited set of parameters, i.e. , there might be a range where deviations are lower. Further, due to the implicit dependence on the coupling, we expect that it should be possible to achieve a better agreement of theory and simulation if one tunes the parameters consistently. However, this is beyond the scope of the current contribution.
5 Discussion and conclusion
We studied the dependence of the contact angle of a droplet on a hydrophobic surface by means of the ShanChen multicomponent LB model and our fluidsurface interaction model.
First, geometrical measurements of the contact angle were used to measure parameter dependencies. Parameters taken into consideration here were system size, coupling parameter and wetting parameter . The influence of the system size on the simulations is caused by finite size effects only and vanishes when simulating larger systems. Discretization errors for curved surfaces diminish then, as well as effects of strictly local force incorporation, leading for instance to finite interfacial thickness.
The pseudo density of the wall component was introduced into the model as parameter of the wetting behaviour [20]. Resonably far from the extremal cases of complete (de)wetting ( and , respectively), a linear dependency of on was observed. This behaviour can be understood following the concept of Eq. (4.4). The coupling parameter of the intercomponent interaction is the same for all components (fluidfluid as well as fluidwall) and therefore cancels from the YoungLaplace law, leaving the contact angle proportional to the ratio of densities only. Nonetheless, since the coupling parameter is determining the density gradient at the interfacial area, there is still an indirect influence on the contact angle. Here, two effects can be differentiated. Given lower coupling, the interfacial area becomes more diffuse, introducing a higher uncertainty to the determination of the position of the interface. For high values of and thus strong repulsive forces, the pseudo potential of the wall can cause the droplet to hover, thereby leaving the definition range of the contact angle.
A method to calculate the expected contact angle as a function of parameters would be expedient. A first ansatz to deduct the contact angle of a single phase multicomponent system from a simple model is given by the determination of the surface tension between each two of the three components present in the droplet system. The main advantage of this approach lies in the small system size needed and the possibility to tabulate the obtained values for future use. Because of the periodic boundaries the precision of the calculation is relying only on the dimension normal to the interface [32]. The surface tension is then determined by its mechanical definition, Eq. (4.1). Comparison between the contact angles calculated by inserting these surface tension values into the YoungLaplace law and the ones measured geometrically in droplet systems yields however large discrepancies. While the range of definition is met for coupling parameters close to numerical instability, in general the contact angle values gained from the model system are much higher than those observered in the droplet system, reaching the complete dewetting limit comparably faster.
To quantify the effect of the simplifications made in the model system, mainly by neglecting the presence of the minority component in the interfacial area as well as the curvature of the interface, the principal of measurement was utilised directly in droplet systems as well. The range of definition of was met for the whole coupling parameter range. However, in the range of linear dependence found by geometrical measurement the contact angle calculated from the YoungLaplace equation was in general lower, diverging by up to .
A problem still persisting with surface tension measurements in the droplet system are discretization effects at the curved interface. Addionally, the interfacial range where there is actual tensorial pressure, is depending on the coupling chosen, 5 to 11 lattice units wide. This introduces a large uncertainty to the integration even along the interface orthogonal to lattice directions. Whether a tuning of the contact angle behaviour by introducing separate coupling for each two components is possible is yet to be determined.
Finally, to evaluate another approach of a priori contact angle determination, an approximation introduced by Huang et al. for a multiphase multicomponent Shan Chen model was adapted to our single phase multicomponent approach. Here results comparable to the surface tension measurements in the droplet system were gained. While the range of definition is met, the contact angle values are up to lower than the geometrically measured. However, since the approximation was postulated for a fixed relation of density and coupling, eventually a change in the parameter set can decrease this deviation. Because of the high calculation cost of parameter search this has been omitted.
To conclude, utilising a pseudo wall density as wetting parameter in a single phase multicomponent Shan Chen LBM it is possible to simulate the complete range of contact angles as determined by geometrical measurement. A priori determination of the contact angle based on simulation parameters is possible with an uncertainty between 10 and 20% depending on the schemes taken into consideration as well as the parameter range.
Acknowledgments
This work was supported by the DFG priority program “nano and microfluidics” and the collaborative research centre (SFB) 716. The computations were performed at the Jülich Supercomputing Centre and the Scientific Supercomputing Centre Karlsruhe.
References
 [1] P. Tabeling. Introduction to microfluidics. Oxford University Press, 2005.
 [2] T. Young. An essay on the cohesion of fluids. Phil. Trans.: Phys. Sci. Eng., 95(1):65, 1805.
 [3] R. N. Wenzel. Surface roughness and contact angle. J. Phys. Chem., 53(9):1466–1467, 1949.
 [4] A. B. D. Cassie. Wettability of porous surfaces. Transactions of the Faraday Society, 40:546, 1944.
 [5] J. Hyväluoma and J. Harting. Slip flow over structured surfaces with entrapped microbubbles. Phys. Rev. Lett., 100:246001, 2008.
 [6] K.Y. Yeh, J. C. Yeh, and J. Changg. Contact angle hysteresis on regular pillarlike hydrophobic surfaces. Langmuir, 24(1):245, 2008.
 [7] C. Dorrer and J. Rühe. Condensation and wetting transitions on microstructured ultrahydrophobic surfaces. Langmuir, 23(7):3820, 2007.
 [8] C. Pirat, M. Sbragaglia, A. M. Peters, B. M. Borkent, R. G. H. Lammertink, M. Wessling, and D. Lohse. Multiple time scale dynamics in the breakdown of superhydrophobicity. Europhys. Lett., 81(6):66002, 2008.
 [9] H. Kusumaatmaja, M. L. Blow, A. Dupuis, and J. M. Yeomans. The collapse transition on superhydrophobic surfaces. Europhys. Lett., 81(3):36003, 2008.
 [10] D. Vandembroucq and S. Roux. Conformal mapping on rough boundaries; applications to harmonic problems. Phys. Rev. E, 55(5):6171, 1997.
 [11] A. Marmur. From hygrophilic to superhygrophobic: Theoretical conditions for making highcontactangle surfaces from lowcontactangle materials. Langmuir, 24(14):7573, 2008.
 [12] P. Roach, N. Shirtcliffe, and M. Newton. Progess in superhydrophobic surface development. Soft Matter, 4(2):224, 2008.
 [13] C. I. Bouzigues, P. Tabeling, and L. Bocquet. Nanofluidics in the Debye layer at hydrophilic and hydrophobic surfaces. Phys. Rev. Lett., 101(11):114503, 2008.
 [14] C. Kunert and J. Harting. Simulation of fluid flow in hydrophobic rough micro channels. Int. J. Comp. Fluid Dyn., 22:475, 2008.
 [15] X. Shan and H. Chen. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E, 47(3):1815, 1993.
 [16] X. Shan and H. Chen. Simulation of nonideal gases and liquidgas phase transitions by the lattice Boltzmann equation. Phys. Rev. E, 49(4):2941, 1994.
 [17] M. R. Swift, W. R. Osborn, and J. M. Yeomans. Lattice Boltzmann simulation of nonideal fluids. Phys. Rev. E, 75(5):830, 1995.
 [18] M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans. LatticeBoltzmann simulations of liquidgas and binary fluid mixtures. Phys. Rev. E, 54(5):5041, 1996.
 [19] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A, 43(8):4320, 1991.
 [20] J. Harting, C. Kunert, and H. Herrmann. Lattice Boltzmann simulations of apparent slip in hydrophobic microchannels. Europhys. Lett., pp 328–334, 2006.
 [21] C. Kunert and J. Harting. Roughness induced apparent boundary slip in microchannel flows. Phys. Rev. Lett., 99:176001, 2007.
 [22] C. Kunert and J. Harting. On the effect of surfactant adsorption and viscosity change on apparent slip in hydrophobic microchannels. Progress in CFD, 8:197, 2008.
 [23] R. Benzi, L. Biferale, M. Sbragaglia, S. Succi, and F. Toschi. Mesoscopic modeling of a twophase flow in the presence of boundaries: the contact angle. Phys. Rev. E, 74:021509, 2006.
 [24] H. Huang, D. T. Thorne, M. G. Schaap, and M. C. Sukop. Proposed approximation for contact angles in ShanandChentype multicomponent multiphase lattice Boltzmann models. Phys. Rev. E, 76:066701, 2007.
 [25] F. J. Higuera, S. Succi, and R. Benzi. Lattice gas dynamics with enhanced collisions. Europhys. Lett., 9(4):345, 1989.
 [26] P. L. Bhatnagar, E. P. Gross, and M. Krook. Model for collision processes in gases. I. small amplitude processes in charged and neutral onecomponent systems. Phys. Rev., 94(3):511, 1954.
 [27] P. J. Love, M. Nekovee, P. V. Coveney, J. Chin, N. GonzálezSegredo, and J. M. R. Martin. Simulations of amphiphilic fluids using mesoscale latticeBoltzmann and latticegas methods. Comp. Phys. Comm., 153:340, 2003.
 [28] J. Harting, M. Harvey, J. Chin, M. Venturoli, and P. V. Coveney. Largescale lattice Boltzmann simulations of complex fluids: advances through the advent of computational grids. Phil. Trans. R. Soc. Lond. A, 363:1895–1915, 2005.
 [29] S. Chen, H. Chen, D. Martínez, and W. H. Matthaeus. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys. Rev. Lett., 67(27):3776, 1991.
 [30] S. Succi. The lattice Boltzmann equation for fluid dynamics and beyond. Oxford University Press, 2001.
 [31] P. G. de Gennes. Wetting: statics and dynamics. Rev. Mod. Phys., 57:827–863, 1985.
 [32] N. GonzálezSegredo and P. V. Coveney. Coarsening dynamics of ternary amphiphilic fluids and the selfassembly of the gyroid and sponge mesophases: latticeBoltzmann simulations. Phys. Rev. E, 69:061501, 2004.