Lattice-Boltzmann simulations of electrowetting phenomena

We present a lattice-Boltzmann method that can simulate the coupled hydrodynamics and electrostatics equations of motion of a two-phase fluid as a means to model electrowetting phenomena. Our method has the advantage of modelling the electrostatic fields within the lattice-Boltzmann algorithm itself, eliminating the need for a hybrid method. We validate our method by reproducing the static equilibrium configuration of a droplet subject to an applied voltage and show that the apparent contact angle of the drop depends on the voltage following the Young-Lippmann equation up to contact angles of $\approx 50^\circ$. At higher voltages, we observe a saturation of the contact angle caused by the competition between electric and capillary stresses. We also study the stability of a dielectric film trapped between a conducting fluid and a solid electrode and find a good agreement with analytical predictions based on lubrication theory. Finally, we investigate the film dynamics at long times and report observations of film breakup and entrapment similar to previously reported experimental results.


I. INTRODUCTION
Electrowetting refers to the spreading of an electrically conducting liquid on a solid electrode when a voltage difference is applied between the two [1].Because of its ability to control the interaction of liquids with solid surfaces, electrowetting has triggered a number of applications, such as droplet-based microfluidic devices [2,3], droplet actuation [4] and mixing [5][6][7], deformable optical apertures [8] and lenses [9], and electronic paper displays [10,11].Broadly speaking, there are two types of electrowetting setups: Electrowetting On Conductor (EWOC), in which the conductive liquid is in direct contact with the solid electrode [12], and the more popular Electrowetting On Dielectric (EWOD), in which direct contact is removed by coating the electrode with a dielectric layer [13].
The simplest electrowetting situation, used widely in many EWOC and EWOD setups, is the spreading of a droplet of conductive liquid suspended in an ambient dielectric fluid that completely wets the solid surface [14].During the actuation, the ambient fluid forms a thin film underneath the droplet that can become unstable and break up into small "bubbles" that remain in contact with the solid [15,16].Such a transition introduces mobile contact lines [17], which can drastically affect the friction force acting on its overall dynamics [18].On the other hand, the spreading of a droplet at high voltages can reach a saturation regime, where the apparent contact angle that the droplet forms with the solid settles to a limiting value [19].At even higher voltages, the edge of the spreading droplet can become unstable, and trigger the breakup of small droplets that form coronal patterns around the mother drop [20].
Despite these important advances, the rich phenomenology of electrowetting remains to be fully under-stood.For this purpose, it is essential to develop computational methods that capture the multiphase fluid dynamics and that resolve the effect of electrostatic interactions, as these can help interpret experiments and inform theory.The Lattice-Boltzmann Method (LBM) has proved to be a powerful tool to study mulitphase fluid dynamics [21].To implement electrowetting within the LBM, it has been proposed to prescribe the interaction energy of the surface [22,23], which leads to an effective contact angle.Such an approach, however, does not capture the underlying coupling between the hydrodynamic and electrostatic fields.As a means to overcome this limitation, hybrid methods that solve the electrostatic field equations separately have been developed [24], but these come at the expense of running and coupling two numerical solvers concurrently.
Here we present a lattice-Boltzmann method capable of solving the coupled hydrodynamics-electrostatics equations that govern electrowetting phenomena within a single algorithm.We use the so-called free-energy approach as a starting point to model the multiphase fluid dynamics, and show that the effect of the electrostatic energy can be included explicitly in the corresponding energy functional.We introduce a set of lattice-Boltzmann equations, where the electrostatic potential field is determined by a new set of distribution functions.We validate this "all-in-one" method by comparing the electrowettinginduced spreading of a droplet to the classical theory of Young and Lippman [25].To illustrate the utility of the method, we present results of the stability of the thin film separating a conducting droplet and a solid electrode, considering both the linear and non-linear regimes.uid, and a dielectric, corresponding to the surrounding phase.We describe the two-fluid system using a diffuseinterface model that identifies each phase using an order parameter, or phase field, φ(x, t), where x denotes the position vector and t denotes time.Without loss of generality, we let φ > 0 be the conductive phase and φ < 0 be the dielectric.
The Helmholtz free energy of the fluid-fluid system can be defined as [26] F The first term corresponds to the volumetric contribution to the free energy over the region occupied by the fluid, Ω.This consists of the well-known energy density of a binary fluid [27,28], where the square-gradient term allows the coexistence of the two bulk phases, of equilibrium phase-field values φ = ±1, separated by a diffuse interface of thickness and surface tension γ.The second integral in Eq. ( 1) corresponds to the surface interaction energy of the fluid with the solid electrode, whose boundary is denoted by ∂Ω, and where the constant ζ is called the wetting potential [29].
In equilibrium, and in the absence of an electric field, the fluid-fluid interface is expected to intersect the solid boundary at an angle θ 0 determined by the Young-Dupré relation [29], where γ sd and γ sc are the solid-dielectric fluid and solidconductive fluid surface tensions.This is a standard result that can be obtained from Eqs. ( 1)-( 3), which yield a relation between the wetting potential and the contact angle [30]: where α = arccos(sin 2 θ 0 ).It can also be shown that, in such a limit, the pressure field, p(x), is uniform in each phase, but jumps across the interface satisfying the Young-Laplace relation where κ is the interface curvature [31].
To model the electrostatic behaviour of the fluid mixture we introduce the electrostatic free energy: which quantifies the potential energy density of the electric field E(x) = −∇V , where V (x) is the electric potential and ε is the electric permittivity [32,33].
Out of equilibrium, local differences in the total free energy, F = F th + F el , give rise to capillary and electrostatic forces.On the one hand, changes in the phase field lead to a chemical potential field and a corresponding capillary force density which reduces to Eq. ( 5) in equilibrium [31].On the other hand, changes in the electric potential give rise to the electric charge distribution [33] el (x, t) and to the electric force density which is the Lorentz force in the absence of magnetic fields [33].
The chemical and electrostatic force densities, Eqs. ( 8) and ( 10), together with the local pressure gradient, −∇p, change the momentum of the fluid.The resulting total force density can be written in terms of a generalised pressure tensor, Π, i.e., This leads to the expression where the last term in brackets is the Maxwell stress tensor [33] and I is the identity matrix.The equations of motion of the fluids are obtained as follows.First, imposing the conservation of momentum leads to the incompressible Navier-Stokes equations where u(x, t), ρ and µ(x) are the velocity field, density and dynamic viscosity of the fluid, respectively, and the superscript T denotes matrix transposition.To allow viscosity differences between the two phases we impose the local viscosity as where µ c and µ d are the bulk viscosities of the conductive and dielectric fluids.Imposing the conservation of the phase field leads to a convection-diffusion equation, often referred to as the Cahn-Hilliard equation [34]: where M is called the mobility.
To complete the formulation of the problem, we need to specify the electrostatic force density, which is a function of the potential field, V .In the following, we assume that both phases are ideal, i.e., the conductor has a vanishing electrical resistivity, while the dielectric has a vanishing electrical conductivity.It then follows that, since the electric field in the conductor is zero, the potential is constant in the bulk of that phase, i.e., On the other hand, for a perfect dielectric el = 0, so Eq. ( 9) reduces to The boundary conditions for the coupled set of PDEs, equations ( 13), ( 15) and (17), are specified as follows.
For the velocity field we impose the impenetrability and no-slip boundary conditions: For the phase field, we impose the natural boundary condition where n is the unit normal to the solid boundary, and which enforces the wetting behaviour of the fluid-fluid mixture.Finally, for the potential we impose where V b is the potential at the boundary.
The lattice-Boltzmann method is a computational fluid dynamics solver that iterates the discretised Boltzmann equations and where f q and g q are particle distribution functions that represent the average number of fluid particles with position x and velocity c q at time t.Space and time are discretised, and the velocity space is sampled by a finite set of vectors {c} Q−1 q=0 , where Q is the number of directions in which the particle populations can move.Here, we use the D2Q9 model, which consists of a two-dimensional square lattice with Q = 9 (see Appendix A).
The time evolution of the distribution functions, given by Eqs. ( 21) and ( 22), consists of a collision step and a streaming step.The collision step, performed by the second term on the right-hand-side in each equation, relaxes the distribution functions local equilibrium values, f eq q and g eq q .
Here we use the Multi-Relaxation Time scheme (MRT) to model the collision of the f q , i.e., where the coefficients Λ qr determine the relaxation rate to equilibrium and are constructed using the Gram-Schmidt orthogonalisation procedure [35].For the collision of the g q we use the single-relaxation time approximation, where we set Λ = 1, which helps improve the stability of the numerical method without loss of generality [21].
The connection between the lattice-Boltzmann equations and the hydrodynamic equations is done by relating the moments of the distribution functions to the hydrodynamic fields.The local mass, momentum and phase fields correspond to and The equilibrium distributions, f eq q and g eq q , are constructed to convey the thermodynamic behaviour of the fluid and to ensure the local conservation of mass and momentum.This is done by requiring that their moments satisfy the conditions: q f eq q = ρ, q g eq q = φ, q c q f eq q = ρu, q c q g eq q = φu, q c q c q f eq q = Π + ρuu and q c q c q g eq q = 2M ϑI + φuu.Suitable expressions of the equilibrium distributions have been reported before [34,36].For the f eq q , we use q (28) if q = 0, and For the g eq q , we use q , (30) if q = 0, and In these expressions, the w q are weighting factors determined by the geometry of the lattice, H = H (n) (c q ) is the tensor Hermite polynomial of n-th degree, and c s = 1/ √ 3 is a constant that represents the speed of sound [37] (see Appendix A for a list of expressions).

A. The electric potential
As discussed in §II, to model the effect of the electrostatic potential field, it suffices to introduce an algorithm that solves Laplace's equation in the dielectric, whilst keeping the potential to a constant value in the conductor.
Hence, we take inspiration from the diffusive dynamics which arises from the LBM itself [38], and introduce a third lattice-Boltzmann equation in the following form, where we use a single-relaxation-time collision operator, where Λ = 1.This new distribution function is related to the local electric potential, V , by the relations and Eq. ( 35) offers the advantage of setting the electric potential to a prescribed value, by fixing the right-hand side, and thus allows the modelling of a conducting liquid (for which the potential equilibrates to a constant).We now analyse the long-time, large-lengthscale behaviour of Eqs. ( 32)- (35).First, we express Eq. ( 32) in terms of the equilibrium distribution, h eq q , using Eq. ( 35).
This is done by writing the collision step as a differential operator acting on h eq q (for details, see Appendix B), i.e., ) Applying the summation operator, q , to Eq. ( 36), and using Eqs.(34) and (35), we find where we identify = c 2 s /2.During a relaxation process the first and second terms in Eq. ( 37) will asymptotically vanish, and thus, V will satisfy Eq. ( 17) at long times.In the context of electrowetting, one requires that this relaxation is faster than the typical timescales of the hydrodynamic fields, u and φ.
To quantify the transient, let us investigate the solutions of Eq. (37).Since the equation is linear, we proceed in the standard way by proposing the Ansatz V = X(x)T (t) [39].This leads to the ordinary differential equation for the temporal part, and a partial differential equation for the spatial part, where K = const., is the eigenvalue that couples the system of equations.
For the temporal part, Eq. ( 38), we look for solutions that decay at long times, i.e., where the term in brackets is always negative for nonvanishing K.
To better understand the rate of decay of the transient, which is controlled by K, let us focus on the limiting case of a uniform dielectric phase in a rectangular domain of of size L x × L y .In such a case, Eq. ( 39) can be solved analytically [39], leading to the spectrum of eigenvalues where l and m are positive integers.Let us now define the transient period, τ trans , as the characteristic decay time associated to the smallest eigenvalue, which, for the uniform rectangular domain, is The presence of a conductive phase will effectively reduce the domain of Eq. ( 39), and thus, will shift the spectrum of K to higher values.This implies that, Eq. ( 43) is an upper bound for the transient from arbitrary initial conditions to a steady state solution.
However, if the initial conditions for the electric field are close to a stationary solution, the transient number of iterations required to relax a small perturbation will be much smaller.For instance, introducing a perturbation of the order of one lattice unit to a stationary solution will lead to K ≈ 2π.Hence, from Eq. ( 42), the transient reduces to Such a fast relaxation can be particularly useful, for instance, when the bulk electrostatic potential V 0 is varied quasi statically to explore stationary wetting configurations, were a single iteration might be enough to update the electrostatic field.

B. Simulation setup: initial and boundary conditions
We now describe the simulation implementation to model the dynamics in an EWOD setup.The electric potential and its corresponding distribution function are defined in a simulation box of size L x × L y .The twophase fluid and corresponding distribution functions are defined in a simulation box of size L x × (L y − 2d), where d is a gap used to accommodate for a solid dielectric layer.This has the purpose of isolating the conductive phase from the bounding electrodes on the finite domain, and thus, to avoid divergences in the electric field.The permittivity of the solid dielectric is set equal to the permittivity of the dielectric fluid.
The velocity field is set to u(x, t = 0) = 0 (45) everywhere in the simulation domain.The phase field, is initialised to which we specify for the specific configurations reported in § § IV and V.The electric potential is initialised as follows. V At subsequent simulation times, and to smooth out the transition of V from the conductive to the dielectric fluid, we impose the electric potential following the interpolation scheme where β is an interpolation weight defined as where φ thr = 0.9, is a threshold value set to identify the bulk of the conductor.In this way, the potential is fixed to the prescribed value V 0 at the bulk of the conductive phase, whereas it evolves according to Eq. ( 34) in the bulk of the dielectric phase.
Using this setup, we found that the electric potential relaxes to a steady state typically after L 2 x /8 iterations.Nonetheless, since transient hydrodynamic flows are slow compared to the speed of sound (|u| c s ), we found that the distribution function h q could be updated at the same pace as f q and g q , with only one iteration required to relax the electric potential field.
We impose periodic boundary conditions along the x and y directions, and fix the solid electrode at the top and bottom boundaries of the simulation domain.To implement the no-slip boundary condition at the solid surface we use the bounce-back algorithm [40].To implement the wettability of the surface, Eq. ( 19), we compute the gradient and Laplacian of the phase field at near-boundary nodes using finite differences to then fix the corresponding incoming distribution functions from the solid surface [30,36].Finally, to implement the boundary condition on the voltage, V b , we follow a similar approach to that of Ledesma-Aguilar, et al. [38].We specify the distribution functions streaming from sites on the the solid electrode, of position vector x b , to sites in the fluid near the solid boundary, of position vector x nb , according to where the indices q correspond to the distribution functions that stream away from the boundary.Specifically, q ∈ {q : c q + c q = 0, q ∈ Γ}, where Γ := {q : x nb + δc q = x b , 0 < δ < 1} gives the indices of lattice vectors that stream towards the electrode.

IV. ELECTROWETTING OF A DROPLET
In this section we validate the lattice-Boltzmann algorithm by studying the electrowetting-driven spreading of a droplet in an EWOD setup.We start by reviewing the Young-Lippmann classical theory of electrowetting [1,25], before comparing to our simulation results.

A. Review of the Young-Lipmann Theory
Consider a droplet of a conductive liquid in an EWOD setup as shown in Fig. 1.As the potential difference applied between the droplet and the electrode is increased, the electric charges begin to gather at the interface of a conductive liquid with a higher density near the grounded electrode.This configuration corresponds to a capacitor.Therefore, and neglecting the charges that accumulate on the opposite side of the solid dielectric, the electrostatic energy, per unit surface area of the electrode, is cV 2 0 /2, where c is the capacitance per unit area of the configuration [32].Because the droplet's surface is compliant, the electrostatic force leads to a spreading of the liquid on the solid electrode.
The equilibrium configuration of the droplet will be determined by the balance between the work done by the electric field against the increase in surface energy.Mechanically, an infinitesimal radial displacement of the contact line, dR, results in a net radial force on the interface of the droplet.Hence, mechanical equilibrium is achieved when Using Eq. ( 3) and dividing by γ, Eq. ( 51) results in the Young-Lippmann relation [1], where is the electrowetting number.Therefore, the contact angle of a droplet is reduced with increasing applied voltage.Experimentally, Young-Lipmann's result has been verified over a range of voltages.However, it has also been observed that at high voltages the contact angle reaches a saturation value, beyond which the theory is no longer valid [41,42].

B. Lattice-Boltzmann simulations
The initial configuration of the system consists of a circular droplet of the conducting liquid suspended in the dielectric fluid.We impose the initial conditions in the simulations using Eqs.( 45), ( 46) and (47); the initial phase field reads where X 0 = (L x /2, R 0 ), is the initial position of the centre of the droplet, and R 0 its initial radius.The rest of the simulation parameters are summarised in Table I.We first set the potential within the conducting droplet to V 0 = 0 and allow the system to relax for 2 × 10 5 iterations.As the droplet relaxes, it spreads on the surface and acquires a circular-cap shape intersecting the surface with the expected equilibrium contact angle, θ 0 , predicted by Eq. ( 4).Then, we increase the voltage by an amount 0.01 2γd/ε and allow the system to relax for a further 10 4 iterations.Once the relaxation has elapsed, the stationary configuration is recorded.The increment in the applied voltage is repeated until a maximum voltage V 0 = 3 2γd/ε is reached.
Fig. 1 shows a typical equilibrium configuration of the droplet subject to a non-zero potential.The upper part of the droplet conserves a circular shape that, extrapolated, intersects the surface at an apparent contact angle θ(V 0 ).However, near the solid surface, the inclination of the interface is closer to the prescribed equilibrium contact angle.As shown in Fig. 2b, the apparent contact angle decreases with increasing |V 0 |.Note that reversing the polarity of the applied voltage leads to the same decrease in the apparent angle; this is expected, since Eq. ( 10) is invariant upon an inversion of the polarity of the electric potential (V → −V ).Therefore, the simulations capture the competition between electrical and capillary forces, as has been reported previously in experimental observations [19].
Next, we carried out simulations to measure θ(V 0 ) for different values of the equilibrium contact angle, θ 0 .As .At V0 = 0, the shape of a droplet is circular and intersects the solid dielectric at the equilibrium contact angle θ0.For |V0| > 0, the shape of the droplet close to the solid wall is distorted by the electric field, leading to an apparent contact angle, θ(V0).At high applied voltages, the droplet reaches a limiting configuration, where the main drop develops a lip that spreads away from its centre.The region around the lip shows strong fringe fields (inset) and the charge density (dark-red colour map).(b) Variation of the contact angle in response to the electric potential, V0, in lattice-Boltzmann units.The curves show a monotonic decrease in the contact angle with the increasing magnitude of the potential.The inset shows the expected universal collapse as a function of the electrowetting number, η, predicted by the Young-Lippmann relation (dotted-dashed line) at low electrowetting numbers and a later saturation.
shown in Fig. 2b, the θ(V 0 ) curves follow the same trend, with only a shift of the maximum to a value imposed by θ 0 .As shown in the inset, a plot of cos θ(V 0 ) − cos θ 0 shows a linear dependence on η, which is in agreement with the theoretical prediction of Eq. (52).Fitting the simulation data to a straight line gives c ≈ 0.66.
As the voltage in the droplet is increased, the apparent contact angle reaches a saturation value θ ≈ 18.43 • .The saturation effect was found to be independent of the wettability of the surface, and begins to occur when the droplet reaches θ ∼ 50 • .From the simulations, we observe that at the onset of saturation the droplet develops two distinct regions.Close to its centre, the capillary forces smooth out the shape of the interface, which remains circular.However, the region close to the edge is subject to strong fringe fields, and deforms to take the shape of a 'lip', spreading away from the main drop (see panel 4 in Fig. 2(a)).The result is that the bulk profile retains a limiting shape, characterised by the saturation contact angle, while an increase in the voltage results in a further growth of the edge lip.

V. DYNAMICS OF A THIN DIELECTRIC FILM
In this section we illustrate the applicability of the lattice-Boltzmann algorithm to resolve the dynamics of electrowetting liquids.Specifically, we study the stability of a thin dielectric film confined between a solid charged wall and a conductive liquid layer.This problem is relevant in many electrowetting setups, where the spreading conductive liquid often entraps a thin film of dielectric fluid.As the dielectric film becomes thinner, it breaks up into small droplets [16].
We start by formulating the problem analytically, which yields a prediction of the stability of the film in the linear regime.We then report simulation results which we validate against this prediction, and extend our study to report results of the dynamics of the film at long times, including the regime of film breakup and droplet formation.

A. Linear-stability theory
We consider a thin, two-dimensional dielectric film of local thickness H(x, t).The film lies on top of a conducting solid electrode, located at y = −d which is coated with a thin dielectric solid layer of thickness d.At its top, the film is covered by a layer of conducting liquid of negligible viscosity.
To model the dynamics of the thin dielectric layer in the presence of an electric field, we use the lubrication equation [43], As shown by Eq. ( 55), the dynamics is driven by variations in the pressure within the film, p film .This is composed of a capillary contribution, 2γκ, and by a contribution due to the electric stresses on the dielectric fluid, where we assume that the capacitance c for a dielectric film in contact with the dielectric solid layer is given by We now study the stability the dielectric film by analysing Eq. (55) using a perturbative approach.Let us consider the initial sinusoidal interface profile where H 0 is the average height of the film, a is the amplitude of the perturbation, λ the wavelength and τ is the characteristic growth time.Substituting Eq. (58) into Eq.( 55), and assuming a H 0 gives the dispersion relation where ω := µ d H 0 /γτ is the dimensionless growth rate, and k := 2πH 0 /λ is the dimensionless wave number.The first term in Eq. ( 59) corresponds to the destabilising effect of the electric field, which dominates for long-wavelength perturbations.This competes against the stabilising effect of surface tension, which dominates for short wavelengths.Setting ω = 0, corresponding to the onset of instability, gives the separatrix η = 1 2 which gives the minimum electrowetting number for which a perturbation of given wave number leads to instability.

B. LB simulations
We impose the initial conditions in the simulations using Eqs.(45), ( 46) and (47); we introduce an initial perturbation to the interface between the conductive and dielectric fluids by imposing the phase-field profile with corresponds to a sinusoidal perturbation of amplitude a = 1 and wavelength λ = L x .The rest of the simulation parameters are reported in the last column of Table I.To allow the thermodynamic relaxation of the phase field from the initial conditions, we let the simulations run for 10 3 iterations, which we disregard.Fig. 3(a) shows a typical instantaneous configuration of the film after the transient has elapsed.Henceforth, we track the evolution of the fluid-fluid interface, whose location we take as the level curve φ(x, y) = 0. Once the location of the interface is determined, the amplitude of the perturbation is found by fitting the instantaneous level curves to the sinusoidal function y(x) = c 0 + c 1 cos(2πx/L x ), where c 0 and c 1 are fitting parameters.We then fit the measured amplitude data, c 1 (t), to the exponential function A(t) = c 2 exp(t/c 3 ), where c 3 gives the characteristic growth time.To obtain the dependence of the dispersion relation, for a given electrowetting number, we repeat the simulation by varying the system length, L x (see Table I).
Figs. 3(b) and 3(c) show the dispersion relations obtained from the simulations for η = 0 and η = 0.03.The data in the figures is reported in the dimensionless units of Eq. ( 59), where µ d , γ, H 0 are fixed using the values reported in Table I.For η = 0, we observe the expected power-law decay, ω ∝ −k 4 , predicted by the linear stability analysis.For η = 0.03, the dispersion relation shows a range of unstable wave numbers.In both cases, we find a quantitative agreement with Eq. ( 59), which is superimposed to the simulation data as a dashed line.
We measured the growth rate of the perturbation for 21 × 21 points in the η-k space.Fig. 3(d ) shows the simulation results, which we present as a contour plot of ω vs η and k.The separatrix, corresponding to the curve ω(k, η) = 0, was estimated from the data using bilinear interpolation (solid line in the figure).Overall, there is a good agreement with Eq. (60) (shown as a dashed line).We attribute the small discrepancy between the theory and the simulation results to the charge distribution at the diffuse fluid-fluid interface, which is dispersed in a region of the order of the interface thickness .This effect would then alter the capacitance of the dielectric film.Indeed, by fitting the separatrix obtained from the simulations to Eq. (60), we obtain an effective value for H 0 , which is displaced by a small amount (∼ 0.08 ) into the bulk of the conductive phase.
We now turn our attention to the growth of the perturbation at long times, when a/H 0 ∼ 1.This regime, which is not accessible by the linear theory, is revealed in detail by the simulations.As shown in Fig. 4(a), at large perturbation amplitudes inhomogeneities in the electric field become apparent.The simulations capture the increase in charge density in regions where the interface curvature is higher [33].This effect leads to a stronger electrostatic attraction in regions of the interface which lie closer to the solid electrode.As a result, the perturbation grows faster than predicted by the linear theory, and the interface is deformed to an asymmetric shape.
At longer times, the troughs of the perturbation approach the solid surface.In this regime, we found that the wettability of the solid has a strong effect on the dynamics.For θ 0 < 180 • the fluid-fluid interface touches the solid surface, breaking the film into droplets.The subsequent dynamics of the fluid-fluid interface is similar to the dewetting dynamics observed by Edwards, et al. [44]: the retracting edges collect fluid to form dewetting rims, which eventually merge to form a single circular droplet (see Fig. 4(b)).For θ 0 = 180 • , the conducting fluid cannot wet the surface and, hence, the dielectric film does not break up.Therefore, the film takes the shape of a series of 'bumps' which remain connected by a thin film  (of a thickness set by the range of the wetting potential in the simulations).This situation is reminiscent of the oil entrapment regime reported by Saticu et al. [16], who used an EWOD setup to spread water droplets immersed in silicone oil on Teflon-coated electrodes.

VI. CONCLUSIONS
We have presented a lattice-Boltzmann algorithm capable of solving the coupled hydrodynamics and electrostatics equations of motion of a two-phase fluid.The main advantage of our model is its ability to solve the electrostatics equations within the lattice-Boltzmann algorithm itself, eliminating the need for concurrent methods, such as finite differences or finite element methods, to model the electric field.
We have validated our algorithm by presenting numer-ical simulations of the electrowetting of a droplet in an Electro Wetting On Dielectric (EWOD) setup.Our results reproduce the dependence of the apparent contact angle of the droplet on the applied voltage predicted by the Young-Lippman theory.We also observe a saturation of the contact angle at high voltages.The saturation of the contact angle has been reported in experiments, and remains an open question in the field of electrowetting.
In the simulations, the effect is linked to a saturation of the interface curvature, which triggers the formation of a 'lip' at the droplet's edge.Such a balance between the electric and capillary stresses in the simulations might explain the saturation effect observed in experiments, but further experimental evidence is needed to reach a conclusion in this regard.
We have also used our algorithm to study the stability and dynamics of a thin dielectric film in an EWOD setup.For small perturbations, our simulations results agree well with the prediction of lubrication theory.Beyond this small-perturbation regime accessible by theory, we studied the long-time dynamics of the film.Our simulations show that as the film is destabilised and the interface approaches the solid surface.On wettable surfaces, the film breaks up and forms droplets that dewet from the surface.On non-wettable surfaces, we observe the entrapment of the dielectric film and the stabilisation of mound-shaped structures.

Figure 1 .
Figure 1.(Colour online) 2D LB simulations of a droplet in an EWOD set-up.A droplet of conducting liquid sits on top of a dielectric solid of thickness d.The droplet is set to an electric potential V0 and, on the other side of the dielectric surface, the electric potential is set to zero.The dielectric fluid surrounds the droplet where the electric field, E is shown by the stream plot.The dashed line corresponds to the best fitting circle to the cap and intersects the solid surface at the apparent contact angle θ.

Figure 2 .
Figure 2. Simulations of droplet spreading using an EWOD setup.(a) Stationary droplet configurations at different applied voltage, V0.At V0 = 0, the shape of a droplet is circular and intersects the solid dielectric at the equilibrium contact angle θ0.For |V0| > 0, the shape of the droplet close to the solid wall is distorted by the electric field, leading to an apparent contact angle, θ(V0).At high applied voltages, the droplet reaches a limiting configuration, where the main drop develops a lip that spreads away from its centre.The region around the lip shows strong fringe fields (inset) and the charge density (dark-red colour map).(b) Variation of the contact angle in response to the electric potential, V0, in lattice-Boltzmann units.The curves show a monotonic decrease in the contact angle with the increasing magnitude of the potential.The inset shows the expected universal collapse as a function of the electrowetting number, η, predicted by the Young-Lippmann relation (dotted-dashed line) at low electrowetting numbers and a later saturation.

Figure 3 .
Figure 3. (Colour online) Simulation results of the linear stability of a thin dielectric fluid film in an EWOD setup.(a) Close-up of the interface configuration.The conducting liquid (light blue region) is kept at a constant voltage V0, whilst the solid electrode (grey rectangle) remains grounded.The thin dielectric fluid film (white region), of initial average thickness H0, is subject to a sinusoidal perturbation of amplitude a and wavelength λ.Direct contact between the dielectric fluid and the solid electrode is prevented by a thin dielectric film (black line).The stream lines depict the electric field.(b) and (c): Dispersion relations for η = 0 and η = 0.03, respectively.The solid symbols correspond to the simulation results.The dashed lines correspond to a fit to the analytical model.The shaded envelopes represent the error from the curve fitting analysis.The inset shows the expected |ω| ∼ k 4 scaling predicted by the linear theory.(d ) Colour map of the growth rate as function of η and k.The solid line corresponds to the separatrix calculated from the simulation results using linear interpolation.The dashed line corresponds to the theoretical prediction.

Figure 4 .
Figure 4. (Colour online) Entrapment and break-up of unstable dielectric films.(a) Instantaneous configuration of an unstable dielectric film at large perturbation amplitudes.The initial simulation parameters are λ = 60, H0 = 20 and η = 1.1.The charge distribution, shown in dark red, is highest in the regions closer to the solid electrode, and the equipotential curves, perpendicular to the electric field, increase in density.(b) Long-time evolution of the dielectric film for λ = 512, H0 = 20, and η = 0.1; and θ0 = 120 • (left) and θ0 = 180 • (right).On a wettable surface, the dielectric fluid breaks into isolated films that dewet to form droplets.On a non-wettable surface the wetting potential prevents the breakup of the film, leading to its entrapment.

Table I .
Parameters for the simulations of the spreading of a droplet and the dielectric film dynamics.