{'fraction_non_alphanumeric': 0.03184415097873935, 'fraction_numerical': 0.008491773594330492, 'mean_word_length': 4.483992467043315, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contact line dynamics in binary lattice Boltzmann simulations

We show that, when a single relaxation time lattice Boltzmann algorithm is used to solve the hydrodynamic equations of a binary fluid for which the two components have different viscosities, strong spurious velocities in the steady state lead to incorrect results for the equilibrium contact angle. We identify the origins of these spurious currents, and demonstrate how the results can be greatly improved by using a lattice Boltzmann method based on a multiple-relaxation-time algorithm. By considering capillary filling we describe the dependence of the advancing contact angle on the interface velocity.


I. INTRODUCTION
There are many different physical systems for which contact line dynamics plays an important role. These range from large scale industrial processes, such as oil recovery and coating, to micron scale phenomena, such as flow in microfluidic channels or the movement of drops across micropatterned surfaces. Despite its importance, contact line dynamics is still a challenging problem that has not been fully solved. This is partly because it is an inherently multiscale problem. The hydrodynamic flow of the fluids far away from the three phase contact line has to be matched to the microscopic details of the contact line motion [1]. In particular, a mechanism is needed by which the usual hydrodynamic, no-slip boundary condition can be violated in the vicinity of the contact line [2,3].
One solution to this problem has been to consider a class of hydrodynamic models, called diffuse interface models [4,5,6], where the interface at the contact line has a finite width. As a diffuse interface is pushed across a surface it is pulled out of equilibrium. This leads to diffusive transport of fluid across the interface and hence an effective slip at the contact line. Here we shall consider a diffuse interface model for a binary fluid, where the equilibrium properties of the fluid are described by a Landau free energy functional. The equilibrium contact angle is controlled by a surface term in the free energy.
The hydrodynamic equations of motion of diffuse interface models can be solved in many different ways, but in this paper, we shall focus on one particular mesoscopic modelling technique, the lattice Boltzmann method [7,8,9]. Lattice Boltzmann simulations have been succesfully used to study several contact line problems. Examples include the spreading of drops on chemically patterned [10] and superhydrophobic surfaces [11], modelling of contact angle hysteresis [12], and capillary filling in microfluidic channels [13,14,15].
The majority of previous work on binary fluids has concentrated on the case where the two components have equal viscosity. In this case a simple and widely used lattice Boltzmann approach, the BGK algorithm, works well, agreeing with analytical results for the equilibrium contact angle. However, if the two binary components have different viscosities, this is no longer the case: we shall show that the equilibrium contact angle is predicted incorrectly by the algorithm because there are severe spurious velocities in the steady state.
There are several important aspects of the behaviour of multicomponent fluids where it is essential or highly desirable to be able to model a two component fluid, where the components have different viscosities. These include drops moving on surfaces, where the viscosity of the surrounding fluid must be substantially smaller than that of the drop to access a rolling regime, instabilities such as viscous fingering which are driven by a viscosity difference between the two fluid components, and capillary filling, where the simple theories assume that the displaced fluid has zero viscosity.
In this paper we identify two primary reasons for the spurious currents in BGK lattice Boltzmann simulations of contact line hydrodynamics. We show how the spurious effects can be greatly suppressed by using an algorithm based on a multiple-relaxation-time lattice Boltzmann approach. We demonstrate that the new algorithm gives excellent agreement with theory, both for the equilibrium contact angle, and for the advancing contact angle, measured in capillary filling simulations.
The paper is organised in the following way: We begin, in Sec. II, by introducing the free energy model for a binary fluid system. We summarise two different lattice Boltzmann implementations that can be used to solve the binary fluid's hydrodynamic equations of motion, the BGK model and a multiple-relaxation-time method, in Secs. IIIA and B. In Sec. IV, we measure the equilibrium contact angle for a drop on a surface and find that, for the BGK approach, it deviates from the value predicted theoretically if there is a viscosity ratio between the two phases. This deviation is caused by anomalous spurious currents near to contact points. In Sec. V, we discuss the origin of the spurious currents: firstly long ranged effects and secondly non-zero velocities induced by the bounce-back boundary conditions. We propose, in Sec. VI, an algorithm based on a multiple-relaxation-time lattice Boltzmann implementation [18,19] which significantly suppresses the spurious currents at the contact point. In Sec. VII we use simulations of capillary filling in smooth channels to measure the dependence of the advancing contact angle on capillary number. Finally, in Sec. VIII, we summarise the results and conclude.

II. MODELLING THE BINARY FLUID
The equilibrium properties of a binary fluid can be described by a Landau free energy functional [6] where the bulk free energy density ψ is taken to have the form The first integral in Eq. (1) is taken over the volume of the system. ρ is the fluid density and φ is the concentration. c is a lattice velocity parameter, described below and a is a constant. This choice of ψ gives binary phase separation into two phases with φ = ±1. The κ term in Eq. (1) represents an energy contribution from gradients in φ and is related to the surface tension between the two phases through γ = 8κa/9 [6]. The second integral in Eq. (1) is over the system's surface and is used to model the fluid-solid surface tensions [16]. The parameter ξ determines the contact angle.
Taking the functional derivative of Eq. (1) with respect to φ gives the chemical potential which is constant in equilibrium. Minimisation of the free energy also shows that the gradient in φ at the boundary is ∂ ⊥ φ| b = ξ/κ [6], where the partial derivative ∂ ⊥ is taken in a direction normal to the surface. The dynamics describing how the fluid approaches equilibrium are determined by the pressure tensor The hydrodynamic equations for the binary fluid are [6] ∂ t ρ + ∂ α (ρv α ) = 0, where v is the fluid velocity, ν is the kinematic viscosity and M is a mobility. The equilibrium properties of the fluid appear in the equations of motion through the pressure tensor and the chemical potential.

III. LATTICE BOLTZMANN IMPLEMENTATIONS
The equations of motion (6)- (8) can be solved using a lattice Boltzmann algorithm. We shall consider two different lattice Boltzmann approaches in this paper and, for convenience, we summarise both here.

A. Single-relaxation-time lattice Boltzmann
To implement a lattice Boltzmann algorithm for a binary fluid in two dimensions the system is divided up into a square grid of points, with two particle distribution functions f i (r, t) and g i (r, t) on each point. The label i denotes a particular lattice velocity vector e i , defined by e 0 = (0, 0), e 1,2 = (±c, 0), e 3,4 = (0, ±c), e 5,6 = (±c, ±c), and e 7,8 = (∓c, ±c). The lattice velocity parameter c is defined by c = ∆x/∆t, where ∆x is the spacing between nearest neighbouring points and ∆t is the time step. The physical variables are obtained from the particle distribution functions using The time evolution equation for the particle distribution functions, using the standard BGK approximation [9], can be broken down into two steps. The first is a collision step This is followed by a streaming step, which moves particles along their corresponding lattice velocity vector direction f eq i and g eq i are equilibrium distribution functions, defined as a power series in the velocity. A summary of our choice of equilibria, which is motivated to help reduce spurious velocities around interfaces [17], is given in Appendix A. Note, the inter-coupling between f i and g i comes through f eq i and g eq i . In particular, the large variation in φ at an interface influences f i by ∂ α φ and ∇ 2 φ terms in f eq i . A Chapman Enskog expansion can be performed to show that the lattice Boltzmann collision (10, 11) and streaming (12,13) operations lead to the hydrodynamic equations (6)-(8) in the limit of long length and time scales [7]. The relaxation parameters τ ρ and τ φ are related to the kinematic viscosity and mobility through where Γ is a tunable parameter that appears in the equilibrium distribution.

B. Multiple-relaxation-time lattice Boltzmann
We next summarise the basic methodology behind the multiple-relaxation-time (MRT) lattice Boltzmann approach. More details are given in [18,19]. The idea behind multiple-relaxation-time lattice Boltzmann is that different relaxation parameters are used for different linear combinations of the distribution functions. In particular, the relaxation parameters responsible for generating the viscous terms in the Navier-Stokes equation (7) are set to τ ρ , those connected to conserved quantities to ∞, and all others to 1.
In multiple-relaxation-time lattice Boltzmann the collision term 1 τρ [f i − f eq i ] on the right hand side of the lattice Boltzmann equation for f i (10) is replaced by where the particle distributions f i and f eq i are written as column vectors and M is a matrix. One possible choice for M is [20 This matrix performs a change of basis. The new basis is designed to contain more physically relevant variables. For instance, when the first row is dotted with f the density ρ = i f i is obtained. Similarly, the fourth and sixth lines calculate the momentum densities ρu x and ρu y , respectively. Each of the rows in M is mutually orthogonal so the inverse follows easily as The matrix S in Eq. (16) is diagonal and has the elements where ω = 1/τ ρ now determines the fluid viscosity in Eq. (14). Note that some of the values are zero. This choice is arbitrary as these modes correspond to the conserved moments; The choice of unity for the other terms will be justified later.
It is not necessary to adopt a multi-relaxation approach for g i as there is an independent parameter Γ, which can be varied to change the mobility of particles in Eq. (15). Therefore, even when using a MRT approach, we set distribution g i to evolve according to Eq. (11) with τ φ = 1. The simulation geometry. The contact angle is defined by θ eq . The crosses denote the points which are fit to the arc of a circle to obtain a numerical estimate of the equilibrium contact angle θ eq . The grid is not to scale.

IV. MEASURING THE EQUILIBRIUM CONTACT ANGLE
In this section we check the extent to which the single relaxation time lattice Boltzmann method gives the correct equilibrium contact angle at the contact point. We find that, for a fluid with a constant viscosity modelled using τ ρ = τ α = τ β = 1, where α and β label the two coexisting bulk phases, good results are obtained. However, when we consider a difference in viscosity between the two components, τ α = 3 and τ β = 0.7, the approach does not work well. In the next section we shall explain why not and describe a method to overcome the problem. Fig. 1 shows a drop resting on a solid surface. In general, the fluid-solid (γ αs and γ βs ) and fluid-fluid (γ) surface tensions are different. At the contact point A the balance of forces is described by Young's law where θ eq defines the equilibrium contact angle. We performed simulations to verify this relation. The lattice Boltzmann system size was set to 300 × 100 lattice units and the parameters used were a = 0.04, τ φ = 1, Γ = 0.5, ∆t = ∆x = 1, and κ = 0.04, giving an interface width of W = 2 2κ/a = 2.8 lattice sites and an interfacial tension of γ = 0.038 in lattice units.
The relaxation parameter τ ρ was determined by such that it changed smoothly through the interface and had the bulk values τ α and τ β in the two bulk phases. The values τ α = 3 and τ β = 0.7 were chosen to give a viscosity ratio of R ν = (τ α − 0.5)/(τ β − 0.5) = 12.5. Non-slip boundary conditions at the walls were achieved using a standard mid-link, bounce-back method [21]. The contact angle was varied by changing the gradient of the order parameter at the solid boundary, ∂ ⊥ φ| b . Initially, a semi-circle of fluid of radius R = 35 was placed on the surface with a contact angle of 90 • . The boundary conditions were set to ∂ ⊥ φ| b = 0 and the system was evolved for 3 × 10 4 time-steps, such that it reached equilibrium. ∂ ⊥ φ| b was then quasi-statically increased over the course of 10 6 time steps and the variation in the contact angle was measured. The process was repeated, but this time decreasing ∂ ⊥ φ| b until the surface was completely wet.
Numerical measurements of contact angle were performed by matching the interface to the arc of a circle [25]. Specifically, each link between neighbouring lattice sites was examined to see if at one end it was fluid α (φ > 0) and at the other end fluid β (φ < 0). If this was the case then a linear interpolation was used to predict the point on the link where φ = 0. These crossing points are illustrated by the crosses in Fig. 1 (note that this grid is not to scale). The equilibrium contact angle was estimated by performing a least squares fit of the crossing points to a circular section and then calculating the contact angle the section made with the surface.
Results for the equilibrium contact angle for different values of τ α and τ β are shown in Fig. 2 and compared to the exact result [6], where α = arccos sin 2 θ eq . The agreement is good for τ α = τ β = 1 ( Fig. 2(a)), but there are large errors for τ α = τ β (Fig. 2(b)). In the next section we discuss the reasons behind this discrepancy.

V. THE ORIGINS OF SPURIOUS CURRENTS NEAR THE CONTACT POINT
The reason that the single relaxation time lattice Boltzmann approach gives an incorrect equilibrium contact angle for τ ρ = 1 is that near to the contact point there are strong spurious velocities which continuously push the system out of equilibrium and result in the deformation of the drop. Note that even when the surface is neutrally wetting, ∂ ⊥ φ| b = 0, the wetting angle measured numerically is ∼ 10 • larger than the expected value of θ eq = 90 • . We focus on this case and replace the drop geometry with a simpler system consisting of a stripe of component α between two neutrally wetting walls, as depicted in Fig. 3(a). This simulation was performed for a system of size 28×15 lattice units using a large relaxation parameter τ ρ = τ α = τ β = 10. Note that we do observe the correct 90 • contact angle in this case, but only because the viscosities of the two phases are the same. If we set τ α > τ β then the stripe becomes barrel shaped, consistent with the measurements reported in Fig. 2(a).
The black arrows in Fig. 3(a) show the steady-state, spurious velocity field in the system. The magnitude of the velocities is typically of order 10 −3 c. We have identified two contributions to the spurious velocities; one from long range effects, a second from the bounce back boundary conditions. We discuss each in turn.

A. Spurious velocities from long range effects
For a unbounded, system at steady-state, it is possible, by iterating the lattice Boltzmann evolution equation (10), to write the particle distribution function at any given lattice point in terms of equilibrium distributions along lines defined by the lattice velocity vector directions: Note that when τ ρ = 1 this reduces to f i (r) = f eq i (r − ∆te i ). This special case could be described as 'local', in the sense that the particle distribution function only depends on the equilibria of its neighbours.
When τ ρ > 1, contributions to the sum in Eq. (23) decay with a characteristic length This diverges as τ ρ → ∞. When τ ρ < 1, each term in the sum (23) alternates in sign and the envelope decays exponentially with length scale which diverges as τ ρ → 1 2 (This limit makes sense because it corresponds to the viscosity in Eq. (14) becoming zero.) Therefore, for τ ρ both high and low, the distribution at any given point is dependent on other nodes a long distance away.
Even without the presence of solid boundaries, it has been noted that these long range interactions can give rise to large spurious currents around curved interfaces [17]. In particular, the size of the spurious velocities scales as τ 3 ρ for large τ ρ , so choosing τ ρ large is not advisable. The problem becomes even more acute in the presence of solid boundaries. Figure 3(b) schematically shows the situation for a particular lattice node A close to the contact point. Because we focus on the case of neutrally wetting walls (contact angle 90 o ), this allows us to compare two systems we expect to have the same steady state configuration, i.e. a rectangular stripe of fluid at rest. In the first, simpler system we imagine that the walls have been replaced by periodic boundary conditions in the y-direction. In this case, contributions to Eq. (23) radiate out indefinitely in all 8 directions. This is illustrated in Fig. 3(b) by the filled circles, each of which denote a term in the sum on the right hand side of Eq. (23). (Note, the arrows indicate the direction of motion for the distribution function from which each term originates.) Because the boundaries have been removed, the symmetry of the system implies that any point above or below A at the same value of x must behave exactly the same as A. Since there can be no net flux of fluid across any given plane (as this would lead to non-conservation of mass), then in steady state, the system must be at rest. It can be numerically confirmed that no steady state spurious velocities are generated in this case.
In the second system, the periodic boundary conditions are replaced with bounce back boundary conditions, and we now discuss what implications this has. In particular, Eq. (23) needs to be modified to take into account the fact that the contributions to the sum do not extend indefinitely in all 8 directions because of the presence of the boundaries. This is illustrated at point B in Fig.  3(b), where the remaining terms in the sum, which would have originated in the direction of C, come, in fact, from a reflected branch in the direction of A. Because of the interface in the system, these two paths give significatly different contributions to the sum. Therefore, the reflection of branches at the boundary breaks the symmetry of the system in the y direction. (This can clearly be seen if we compare point A to a point on the boundary, where the incoming branches are reflected immediately.) This broken symmetry leads to the generation of the spurious, steady state velocities.
The range around the contact point over which this spurious force is active is determined by the decay lengths λ in Eqs. (24) and (25). For high or low viscosities λ is large and this is a long range effect. Correcting for it at the boundary involves extrapolating into the surface to predict, for example, the density variation down the dotted branches near to C. This is relatively simple for the 90 • case shown, as the unknown nodes can be obtained by reflecting the system across the boundary. However, in the case of arbitrary contact angle the situation is much more complicated and a general solution is far from clear.
Thus implementing solid boundaries in the multicomponent BGK lattice Boltzmann is only advisable in the case τ ρ = 1. (For the same reason it is best to choose τ φ = 1 [22].) However, to simulate phases with different viscosities it is not possible to use τ ρ = 1 in both phases. In Sec. VI we discuss how using a multiple-relaxationtime method can be used to overcome this restriction.

B. Spurious velocities from bounce-back boundary conditions
A second source of spurious velocities at the contact point is the bounce-back boundary conditions. To understand what goes wrong with the boundary conditions for the binary lattice Boltzmann model it is first important to understand why they work well when simulating a single-component fluid. Fig. 4(a) illustrates the midlink, bounce-back process for one dimension. During the streaming the population a moves down the link, is reflected, and then travels back to the node it came from but moving in the opposite direction. (Fig. 4(b) shows an alternative bounce-back scheme in which the wall lies on the boundary nodes. The arguments presented here are equally valid for this case.) Fig. 4(c) shows the distribution for a typical node in a system describing a single component fluid at equilibrium. This is a representation of the Maxwell Boltzmann distribution using a discrete number of velocity vectors. The bounce-back method acts at the boundaries of the system to switch some of the populations to travel in the opposite direction. Because this distribution is invariant under a parity transformation (that is replacing velocity vectors v by −v) then the correct distribution is preserved by the bounce-back boundary. When the system is in a state of shear the particle distribution function, in the rest frame of the fluid, is as depicted in Fig. 4(d). While this is clearly not isotropic, it still preserves invariance under the parity transformation and hence the bounce-back approach is still valid. (In fact, if non-slip boundary conditions are enforced by setting the particle distribution at the boundary to its equilibrium value at rest, an inaccurate shear velocity profile is obtained. This is because the shear-induced distortion of the distribution function is not preserved at the boundaries. ) We now return to the binary fluid case. Fig. 4(e) shows the distribution function for a typical node, lying at rest in a fluid-fluid interface. (In this particular example the interface lies perpendicular to the surface, hence the distribution has an up-down symmetry.) Note that because of the κ terms in the pressure tensor (4) the parity invariance is broken. We consider the case when the position of the node is ∆x/2 above a boundary and mid-link bounceback boundary conditions are being employed. Fig. 4(f) illustrates the situation after the bounce back step, and it clearly shows that the new distribution is not the same as in (e) (in particular, vector a is shortened and b lengthened). Therefore, bounce-back collisions result in the system continually being pushed out of equilibrium, leading to the generation of spurious velocities on or near to boundaries close to the interface between the two phases.
For τ ρ = 1, the distribution function is automatically set to its equilibrium value at each time step. While the typical distribution function in the interface is not invariant under a parity transformation, the equilibrium distribution is. Hence the spurious velocities caused by bounce back boundary conditions are suppressed in this case.

VI. SUPPRESSING THE SPURIOUS CURRENTS
We now describe how an algorithm based on the multiple-relaxation-time lattice Boltzmann approach can be used to significantly reduce the spurious currents at the contact point. The approach comprises the following four steps: Step 1: Calculate the density, concentration and velocity using the moments defined by Eq. (9).
Step 2: Set the velocity of boundary nodes to zero when calculating the equilibrium distribution function.
Step 3: Use the multiple-relaxation-time lattice Boltzmann method described in Sec. III B to perform the collision step.
Step 4: Perform the streaming step with bounce-back at the boundaries.
To justify this choice we need to consider the hydrodynamic and the non-hydrodynamic modes separately. From a hydrodynamic point of view, a non-slip boundary fixes the velocity at that boundary to zero. If the fluid is incompressible, then near to the boundary the fluid flow profile can always be approximated by a shear profile, with streamlines parallel to the surface. In section VB we argued that for a shear profile, the hydrodynamic modes, which are represented by the difference between Figs. 4(c) and (d), are invariant under a parity transformation, and so are preserved by the bounce-back method. Another interpretation of this is that, for the hydrodynamic modes, the contributions at A in Fig. 3 from the reflected branches are the same as those when periodic boundary conditions are considered (provided λ τρ is smaller than the length scale over which the shear profile approximation breaks down as we move away from the boundary). In summary, the hydrodynamic modes behave correctly for both the single and multiple-relaxation-time lattice Boltzmann models. The difference lies in the treatment of the non-hydrodynamic modes.
In the multiple-relaxation-time scheme, because we have chosen τ ρ = 1 for the non-hydrodynamic modes (this choice corresponds to the 1's in Eq. (19)), the spurious velocities generated from long ranged effects in the bulk are immediately removed. This justifies the use of the multiple-relaxation-time algorithm in Step 3. The only potential problem that remains is on the boundary nodes themselves. As discussed at the end of section VB, the distribution function in the interface between fluid phases is not invariant under a parity transformation. This generates spurious velocities at any boundary node in an interface immediately after the streaming step 4. This problem is remedied by the introduction of step 2. Note that this step is consistent with non-slip boundary conditions and does not affect the hydrodynamic correlations (i.e. the non symmetric distribution in Fig. 4(d)).
For a system with variable viscosity it would seem necessary to recalculate the collision matrix C = M −1 SM in (17) at each lattice node and at each time-step. This would be extremely slow computationally and not very practical. The approach we take is to create a lookup table with ∼ 10 4 different values of viscosity and simply pick the closest match.
We find that implementing multiple-relaxation-time lattice Boltzmann with appropriate boundary conditions leads to a significant improvement in the accuracy of the equilibrium contact angle. The results of simulations for τ α = 3, τ β = 0.7, are denoted by the crosses in Fig.  2(a), and show very good agreement with the theoretical, dashed curve. Deviation is only noticeable at small contact angles. This happens for two reasons: Firstly, the dynamics of drop wetting become very slow as the equilibrium contact angle becomes small, and so the assumption of quasi-static equilibrium, as the gradient in φ is slowly reduced, breaks down. (Tanner's law [23] states that for a completely wetting surface in two dimensions, the size of a drop spreading on the surface scales as r ∼ t 1/7 in time.) Secondly, the finite width of the interface W , which is neglected when assuming that the drop should be made up of a circular section, becomes comparable to the height of the drop.

VII. A DYNAMICAL TEST: CAPILLARY FILLING
We have shown that the equilibrium contact angle is accurately recovered in lattice Boltzmann simulations for a binary system with different viscosities only when a multiple-relaxation-time approach is employed. While this highlights a problem with the single relaxation time binary lattice Boltzmann model, contact angle measurement is a static problem. In this section we concentrate on the dynamics of a fluid penetrating a smooth microchannel to measure how the advancing contact angle changes with the velocity of the fluid interface.
Fluid is pulled into a hydrophilic capillary by the Laplace pressure across the interface. Balancing this against the viscous drag of the fluid column gives Washburn's law [24] describing the variation of the length l of fluid in the capillary with time t where h is the width of the capillary, θ a is the advancing contact angle and t 0 is an integration constant. Note that it is appropriate to use, not the static, but the advancing contact angle θ a , as this controls the curvature of the interface and hence the Laplace pressure. Moreover, Eq. (26) assumes that there is no resistance to motion from any fluid already in the capillary. Numerical results showing capillary filling of a two dimensional capillary are presented in Fig. 5   The solid lines in Fig. 5(a) and (c) are fits to the Washburn law. At later times, once inertial effects have become negligible, both the single and multiple-relaxationtime lattice Boltzmann simulations satisfy the Washburn scaling. However, the two methods provide different quantitative results, in particular, capillary filling is considerably slower for the single-relaxation-time method. The difference can be explained by considering how the advancing contact angle varies with the velocity of the interface, v I .
To lowest order in the capillary number, Ca= v I ν/γ, the advancing contact angle is related to the equilibrium angle and the capillary number by [25] cos θ a = cos θ eq − Ca log(KL/l s ), where K is a constant, L is the length scale of the system and l s is the effective slip length at the three phase contact line. Figs. 5(b) and (d) show the expected linear decrease of the measured contact angle with capillary number [25]. For the multiple-relaxation-time al-gorithm, the advancing contact angle tends to the correct value as Ca → 0. We obtain θ a |  Fig. 2, which shows that the measured equilibrium contact angle is too high in the singlerelaxation-time approach.) Since the speed of capillary filling depends on cos θ a (26), capillary filling is considerably slower using this method.
In diffuse interface models of binary fluids the contact line singularity is relieved by inter-diffusion of the two fluid components. This is governed by the mobility M ; as is apparent from Fig. 5, increasing M increases the rate of diffusion and hence the velocity of the contact line. Therefore the parameter M can be used to tune the effective slip length l s .

VIII. SUMMARY AND CONCLUSIONS
We have shown that, if the lattice Boltzmann relaxation parameter τ ρ = 1, strong spurious currents drive the contact line out of equilibrium even in a system at rest. This means that the algorithm gives incorrect values for the contact angle. In many applications it is possible to choose τ ρ = 1, thus avoiding the problem. However this parameter controls the fluid viscosity and so cannot be held at unity for both phases of a binary fluid if the two components have different viscosities.
We demonstrate that the spurious currents arise primarily from two effects. The first is a long-range contribution to the equilibrium distribution function near the contact line that effectively originates in the incorrect phase. The second is the bounce-back boundary conditions which drive the interface out of equilibrium.
Aiming to reduce the unwanted velocities we propose a revised lattice Boltzmann method, based on a multiplerelaxation-time algorithm. We show that the simulations then agree well with the analytical result for the equilibrium contact angle. Moreover the dynamic, advancing contact angle shows the expected physical behaviour, with a slip length that depends on the diffusivity of the fluid.
Using this method it will be possible to perform accurate simulations of a binary fluid where the two components have different viscosities. Examples of problems where the algorithm will prove useful include viscous fingering, rolling of viscous drops and capillary filling. Moreover, it can provide a useful approximation to a liquid-gas system in the limit that evaporationcondensation is not important.