Lattice-Boltzmann simulations of electrowetting phenomena

TL;DR

This paper introduces a lattice-Boltzmann simulation method that models electrowetting phenomena by integrating electrostatics and hydrodynamics, validated against theoretical and experimental benchmarks.

Contribution

The authors develop a hybrid lattice-Boltzmann approach that directly incorporates electrostatic fields, enabling accurate simulation of electrowetting effects without hybrid methods.

Findings

Reproduces static droplet configurations and contact angle dependence on voltage.

Demonstrates contact angle saturation at high voltages due to electric-capillary balance.

Shows good agreement with lubrication theory for dielectric film stability.

Abstract

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.

Figures

Tables

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

Simulation parameter	Symbol	Value in §IV.2	Value in §V.2
Simulation box	$L_x\times L_y$	$512\times 288$	$[418,1256]\times 84$
Surface tension	$\gamma$	$6\times {10}^{-3}$	$8\times {10}^{-3}$
Interface width	$\ell$	$4$	$5$
Contact angle	${\theta }_0$	$[{160}^{\circ },{120}^{\circ }]$	${180}^{\circ }$
Density	$\rho$	$1$	$1$
Dynamic viscosity	${\mu }_{\text{c}}$, ${\mu }_{\text{d}}$	$1/6$, $1/6$	$1/600$, $1/3$
Mobility	$M$	$1/10$	$1/10$
Permittivity	$\epsilon$	1/6	1/6
Dielectric thickness	$d$	$2$	$2$
Initial config.		$R_0=128$	$a=1$, $H_0=20$

Equations

$$F_{\mathrm{th}}[\varphi ]:={\int }_{\Omega }\psi (\varphi ,\text{mathbfit}\nabla \varphi ){\mathrm{d}}^{3}x+{\int }_{\partial \Omega }\zeta \varphi \mathrm{d}S.$$

$$\psi (\varphi ,\text{mathbfit}\nabla \varphi ):=\frac{3\gamma }{\sqrt{8}\ell }\left[\frac{{\varphi }^4}{4}-\frac{{\varphi }^2}{2}+\frac{{\ell }^2}{2}|\text{mathbfit}\nabla \varphi {|}^2\right],$$

$${\gamma }_{\text{sd}}-{\gamma }_{\text{sc}}=\gamma \cos {\theta }_0,$$

$$\zeta =\frac{3}{2}\gamma \mathrm{sgn}({\theta }_0-\pi /2){\left\{\cos (\alpha /3)\left[1-\cos (\alpha /3)\right]\right\}}^{1/2},$$

$$\Delta p=2\gamma \kappa ,$$

$$F_{\mathrm{el}}[V]:=-\frac{1}{2}\epsilon {\int }_{\Omega }|\text{mathbfit}E{|}^2{\mathrm{d}}^{3}x,$$

$$\vartheta (\text{mathbfit}x,t):=\frac{\delta F}{\delta \varphi }=\frac{3\gamma }{\sqrt{8}\ell }\left[\varphi ({\varphi }^2-1)-{\ell }^2{\nabla }^2\varphi \right],$$

$$\text{mathbfit}{f}_{\text{cap}}=-\varphi \text{mathbfit}\nabla \vartheta ,$$

$${\varrho }_{\text{el}}(\text{mathbfit}x,t):=-\frac{\delta F}{\delta V}=-\epsilon {\nabla }^{2}V=\epsilon \text{mathbfit}\nabla \cdot \text{mathbfit}E,$$

$$\text{mathbfit}{f}_{\text{el}}={\varrho }_{\text{el}}\text{mathbfit}E,$$

$$-\text{mathbfit}\nabla \cdot \mathbf{\Pi }:=-\nabla p+\text{mathbfit}{f}_{\text{cap}}+\text{mathbfit}{f}_{\text{el}}.$$

$$\mathbf{\Pi }=\left(\varphi \vartheta -\psi \right)\mathbf{I}+\frac{3\gamma \ell }{\sqrt{8}}\text{mathbfit}\nabla \varphi \text{mathbfit}\nabla \varphi -\epsilon \left(\text{mathbfit}E\text{mathbfit}E-\frac{1}{2}|\text{mathbfit}E{|}^2\mathbf{I}\right),$$

$$\rho \left({\partial }_t+\text{mathbfit}u\cdot \text{mathbfit}\nabla \right)\text{mathbfit}u=-\nabla \cdot \left[\mathbf{\Pi }-\frac{\mu }{2}\left(\nabla \text{mathbfit}u+\nabla \text{mathbfit}{u}^{\text{T}}\right)\right],$$

$$\mu (\varphi ):=\frac{1+\varphi }{2}{\mu }_{\text{c}}+\frac{1-\varphi }{2}{\mu }_{\text{d}},$$

$${\partial }_t\varphi +\text{mathbfit}u\cdot \text{mathbfit}\nabla \varphi =M{\nabla }^2\vartheta ,$$

$$V=V_0\text{as }\varphi \to 1.$$

$${\nabla }^{2}V=0\text{as }\varphi \to -1.$$

$$\text{mathbfit}u(\text{mathbfit}{x}_{\text{b}})=0\text{forleavevmodenobreak }\text{mathbfit}{x}_{\text{b}}\in \partial \Omega .$$

$$\widehat{\text{mathbfit}n}\cdot \text{mathbfit}\nabla \varphi (\text{mathbfit}{x}_{\text{b}})=-\frac{\sqrt{8}}{3\gamma \ell }\zeta ({\theta }_0)\text{forleavevmodenobreak }\text{mathbfit}{x}_{\text{b}}\in \partial \Omega ,$$

$$V(\text{mathbfit}{x}_{\text{b}})=V_{\text{b}}\text{forleavevmodenobreak }\text{mathbfit}{x}_{\text{b}}\in \partial \Omega ,$$

$$f_q(\text{mathbfit}x+\text{mathbfit}{c}_q,t+1)=f_q(\text{mathbfit}x,t)+\mathcal{C}[f{]}_q$$

$$g_q(\text{mathbfit}x+\text{mathbfit}{c}_q,t+1)=g_q(\text{mathbfit}x,t)+\mathcal{C}[g{]}_q,$$

$$\mathcal{C}[f{]}_q:=-\sum _{r=0}^{Q-1}{\Lambda }_{qr}[f_r-f_r^{\text{eq}}](\text{mathbfit}x,t),$$

$$\mathcal{C}[g{]}_q:=-\Lambda [g_q-g_q^{\text{eq}}](\text{mathbfit}x,t),$$

$$\rho =\sum _{q=0}^{Q-1}{f}_q,$$

$$\rho \text{mathbfit}u=\sum _{q=0}^{Q-1}\text{mathbfit}{c}_q{f}_q$$

$$\varphi =\sum _{q=0}^{Q-1}{g}_q.$$

$$f_q^{\text{eq}}(\rho ,\text{mathbfit}u,\mathbf{\Pi })=w_q\left[\frac{1}{c_s}\rho \text{mathbfit}u\cdot H_q^{(1)}+\frac{1}{2c_s^2}\left(\mathbf{\Pi }+\rho \text{mathbfit}u\text{mathbfit}u\right):H_q^{(2)}\right]$$

$$f_0^{\text{eq}}(\rho ,\text{mathbfit}u,\mathbf{\Pi })=\rho -\sum _{q=1}^{Q-1}{f}_q^{\text{eq}}.$$

$$g_q^{\text{eq}}=w_q\left[\frac{1}{c_s}\varphi \text{mathbfit}u\cdot H_q^{(1)}+\frac{1}{2c_s^2}(2M\vartheta +\varphi \text{mathbfit}u\text{mathbfit}u):H_q^{(2)}\right],$$

Full text

††preprint: AIP/123-QED

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 Mugele and Baret (2005). 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 Teh et al. (2008); Ren et al. (2002), droplet actuation Baret et al. (2005) and mixing Mugele et al. (2006); Lu et al. (2017); Kudina et al. (2016), deformable optical apertures Schuhladen et al. (2016) and lenses Berge and Peseux (2000), and electronic paper displays Hayes and Feenstra (2003); You and Steckl (2010). 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 Lomax et al. (2016), and the more popular Electrowetting On Dielectric (EWOD), in which direct contact is removed by coating the electrode with a dielectric layer Nelson and Kim (2012).

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 Quilliet and Berge (2002). 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 Kuo et al. (2003); Staicu and Mugele (2006). Such a transition introduces mobile contact lines Mugele (2009), which can drastically affect the friction force acting on its overall dynamics McHale et al. (2013). 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 Buehrle et al. (2003). 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 Quilliet and Berge (2001).

Despite these important advances, the rich phenomenology of electrowetting remains to be fully understood. 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 Krüger et al. (2016). To implement electrowetting within the LBM, it has been proposed to prescribe the interaction energy of the surface Li and Fang (2009); Clime et al. (2010), 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 Aminfar and Mohammadpourfard (2009), 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 electrowetting-induced spreading of a droplet to the classical theory of Young and Lippman Lippmann (1875). 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.

II Diffuse-interface model of electrowetting phenomena

Let us consider two incompressible, immiscible fluids: a perfect conductor, corresponding to the spreading liquid, and a dielectric, corresponding to the surrounding phase. We describe the two-fluid system using a diffuse-interface model that identifies each phase using an order parameter, or phase field, $\phi({\mathbfit{x}},t)$, where $\mathbfit{x}$ denotes the position vector and $t$ denotes time. Without loss of generality, we let $\phi>0$ be the conductive phase and $\phi<0$ be the dielectric.

The Helmholtz free energy of the fluid-fluid system can be defined as Landau and Lifshitz (1980)

[TABLE]

The first term corresponds to the volumetric contribution to the free energy over the region occupied by the fluid, $\Omega$. This consists of the well-known energy density of a binary fluid Bray (1994); Cahn and Hilliard (1958),

[TABLE]

where the square-gradient term allows the coexistence of the two bulk phases, of equilibrium phase-field values $\phi=\pm 1$, separated by a diffuse interface of thickness $\ell$ and surface tension $\gamma$. The second integral in Eq. (1) corresponds to the surface interaction energy of the fluid with the solid electrode, whose boundary is denoted by $\partial\Omega$, and where the constant $\zeta$ is called the wetting potential Cahn (1977).

In equilibrium, and in the absence of an electric field, the fluid-fluid interface is expected to intersect the solid boundary at an angle $\theta_{0}$ determined by the Young-Dupré relation Cahn (1977),

[TABLE]

where $\gamma_{\text{sd}}$ and $\gamma_{\text{sc}}$ are the solid-dielectric fluid and solid-conductive 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 Briant et al. (2004):

[TABLE]

where $\alpha=\arccos(\sin^{2}\theta_{0})$. It can also be shown that, in such a limit, the pressure field, $p(\mathbfit{x})$, is uniform in each phase, but jumps across the interface satisfying the Young-Laplace relation

[TABLE]

where $\kappa$ is the interface curvature Yang et al. (1976).

To model the electrostatic behaviour of the fluid mixture we introduce the electrostatic free energy:

[TABLE]

which quantifies the potential energy density of the electric field $\mathbfit{E}({\mathbfit{x}})=-\mathbfit{\nabla}V$, where $V({\mathbfit{x}})$ is the electric potential and $\varepsilon$ is the electric permittivity Landau et al. (2013); Jackson (1999).

Out of equilibrium, local differences in the total free energy, $F=F_{\text{th}}+F_{\text{el}}$, give rise to capillary and electrostatic forces. On the one hand, changes in the phase field lead to a chemical potential field

[TABLE]

and a corresponding capillary force density

[TABLE]

which reduces to Eq. (5) in equilibrium Yang et al. (1976). On the other hand, changes in the electric potential give rise to the electric charge distribution Jackson (1999)

[TABLE]

and to the electric force density

[TABLE]

which is the Lorentz force in the absence of magnetic fields Jackson (1999).

The chemical and electrostatic force densities, Eqs. (8) and (10), together with the local pressure gradient, $-\nabla p$, change the momentum of the fluid. The resulting total force density can be written in terms of a generalised pressure tensor, $\bm{\Pi}$, i.e.,

[TABLE]

This leads to the expression

[TABLE]

where the last term in brackets is the Maxwell stress tensor Jackson (1999) and $\mathbf{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

[TABLE]

where $\mathbfit{u}(\mathbfit{x},t)$, $\rho$ and $\mu(\mathbfit{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

[TABLE]

where $\mu_{\text{c}}$ and $\mu_{\text{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 Swift et al. (1996):

[TABLE]

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.,

[TABLE]

On the other hand, for a perfect dielectric $\varrho_{\rm el}=0$, so Eq. (9) reduces to

[TABLE]

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:

[TABLE]

For the phase field, we impose the natural boundary condition

[TABLE]

where $\hat{\mathbfit{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

[TABLE]

where $V_{\text{b}}$ is the potential at the boundary.

III Lattice-Boltzmann method

In this section we formulate a lattice-Boltzmann algorithm capable of integrating Eqs. (13), (15) and (17), subject to the boundary conditions (18)-(20).

The lattice-Boltzmann method is a computational fluid dynamics solver that iterates the discretised Boltzmann equations

[TABLE]

and

[TABLE]

where $f_{q}$ and $g_{q}$ are particle distribution functions that represent the average number of fluid particles with position $\mathbfit{x}$ and velocity $\mathbfit{c}_{q}$ at time $t$. Space and time are discretised, and the velocity space is sampled by a finite set of vectors $\{\mathbfit{c}\}_{q=0}^{Q-1}$, 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_{q}^{\text{eq}}$ and $g_{q}^{\text{eq}}$. Here we use the Multi-Relaxation Time scheme (MRT) to model the collision of the $f_{q}$, i.e.,

[TABLE]

where the coefficients $\Lambda_{qr}$ determine the relaxation rate to equilibrium and are constructed using the Gram-Schmidt orthogonalisation procedure d’Humières et al. (2002). For the collision of the $g_{q}$ we use the single-relaxation time approximation,

[TABLE]

where we set $\Lambda=1$, which helps improve the stability of the numerical method without loss of generality Krüger et al. (2016).

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

[TABLE]

[TABLE]

and

[TABLE]

The equilibrium distributions, $f_{q}^{\text{eq}}$ and $g_{q}^{\text{eq}}$, 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: $\sum_{q}f_{q}^{\text{eq}}=\rho$, $\sum_{q}g_{q}^{\text{eq}}=\phi$, $\sum_{q}\mathbfit{c}_{q}f_{q}^{\text{eq}}=\rho\mathbfit{u}$, $\sum_{q}\mathbfit{c}_{q}g_{q}^{\text{eq}}=\phi\mathbfit{u}$, $\sum_{q}\mathbfit{c}_{q}\mathbfit{c}_{q}f_{q}^{\text{eq}}=\bm{\Pi}+\rho\mathbfit{u}\mathbfit{u}$ and $\sum_{q}\mathbfit{c}_{q}\mathbfit{c}_{q}g_{q}^{\text{eq}}=2M\vartheta\mathbf{I}+\phi\mathbfit{u}\mathbfit{u}$. Suitable expressions of the equilibrium distributions have been reported before Swift et al. (1996); Desplat et al. (2001). For the $f_{q}^{\rm eq}$, we use

[TABLE]

if $q\neq 0$, and

[TABLE]

For the $g_{q}^{\rm eq}$, we use

[TABLE]

if $q\neq 0$, and

[TABLE]

In these expressions, the $w_{q}$ are weighting factors determined by the geometry of the lattice, $H_{q}^{(n)}=H^{(n)}(\mathbfit{c}_{q})$ is the tensor Hermite polynomial of $n$-th degree, and $c_{s}=1/\sqrt{3}$ is a constant that represents the speed of sound Qian and Chen (2000) (see Appendix A for a list of expressions).

Using a Chapman-Enskog expansion, Eqs. (21) and (22), together with Eqs. (23)–(31), reduce to the Navier-Stokes (13) and Cahn-Hilliard (15) equations. From the expansion, the viscosity, $\mu$, is determined by the coefficients of the collision matrix, $\Lambda_{qr}$ d’Humières et al. (2002) (see Appendix A).

III.1 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 Ledesma-Aguilar et al. (2014), and introduce a third lattice-Boltzmann equation in the following form,

[TABLE]

where we use a single-relaxation-time collision operator,

[TABLE]

where $\Lambda=1$.

This new distribution function is related to the local electric potential, $V$, by the relations

[TABLE]

and

[TABLE]

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^{\text{eq}}_{q}$, using Eq. (35). This is done by writing the collision step as a differential operator acting on $h_{q}^{\text{eq}}$ (for details, see Appendix B), i.e.,

[TABLE]

Applying the summation operator, $\sum_{q}$, to Eq. (36), and using Eqs. (34) and (35), we find

[TABLE]

where we identify $\epsilon=c_{s}^{2}/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, $\mathbfit{u}$ and $\phi$.

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(\mathbfit{x})T(t)$ Arfken et al. (2013). This leads to the ordinary differential equation for the temporal part,

[TABLE]

and a partial differential equation for the spatial part,

[TABLE]

where $K=\text{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.,

[TABLE]

where the term in brackets is always negative for non-vanishing $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}\times L_{y}$. In such a case, Eq. (39) can be solved analytically Arfken et al. (2013), leading to the spectrum of eigenvalues

[TABLE]

where $l$ and $m$ are positive integers. Let us now define the transient period, $\tau_{\text{trans}}$, as the characteristic decay time associated to the smallest eigenvalue,

[TABLE]

which, for the uniform rectangular domain, is

[TABLE]

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\approx 2\pi$. Hence, from Eq. (42), the transient reduces to

[TABLE]

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.

III.2 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}\times L_{y}$. The two-phase fluid and corresponding distribution functions are defined in a simulation box of size $L_{x}\times(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

[TABLE]

everywhere in the simulation domain. The phase field, is initialised to

[TABLE]

which we specify for the specific configurations reported in §§ IV and V. The electric potential is initialised as follows.

[TABLE]

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

[TABLE]

where $\beta$ is an interpolation weight defined as

[TABLE]

where $\phi_{\text{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_{x}^{2}/8$ iterations. Nonetheless, since transient hydrodynamic flows are slow compared to the speed of sound ($|\mathbfit{u}|\ll 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 Yu et al. (2003). 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 Briant et al. (2004); Desplat et al. (2001). Finally, to implement the boundary condition on the voltage, $V_{\text{b}}$, we follow a similar approach to that of Ledesma-Aguilar, et al. Ledesma-Aguilar et al. (2014). We specify the distribution functions streaming from sites on the the solid electrode, of position vector $\mathbfit{x}_{\text{b}}$, to sites in the fluid near the solid boundary, of position vector $\mathbfit{x}_{\text{nb}}$, according to

[TABLE]

where the indices $\bar{q}$ correspond to the distribution functions that stream away from the boundary. Specifically, $\bar{q}\in\{q:\mathbfit{c}_{q}+\mathbfit{c}_{q^{\prime}}=0,\,q\in\Gamma\}$, where $\Gamma:=\{q:\mathbfit{x}_{\text{nb}}+\delta\mathbfit{c}_{q}=\mathbfit{x}_{\text{b}},\,0<\delta<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 Lippmann (1875); Mugele and Baret (2005), before comparing to our simulation results.

IV.1 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_{0}^{2}/2$, where $c$ is the capacitance per unit area of the configuration Landau et al. (2013). 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, $\mathrm{d}R$, results in a net radial force on the interface of the droplet. Hence, mechanical equilibrium is achieved when

[TABLE]

Using Eq. (3) and dividing by $\gamma$, Eq. (51) results in the Young-Lippmann relation Mugele and Baret (2005),

[TABLE]

where

[TABLE]

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 Quinn et al. (2005); Peykov et al. (2000).

IV.2 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

[TABLE]

where $\mathbfit{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 1.

We first set the potential within the conducting droplet to $V_{0}=0$ and allow the system to relax for $2\times 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, $\theta_{0}$, predicted by Eq. (4). Then, we increase the voltage by an amount $0.01\sqrt{2\gamma d/\varepsilon}$ 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\sqrt{2\gamma d/\varepsilon}$ 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 $\theta(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\rightarrow-V$). Therefore, the simulations capture the competition between electrical and capillary forces, as has been reported previously in experimental observations Buehrle et al. (2003).

Next, we carried out simulations to measure $\theta(V_{0})$ for different values of the equilibrium contact angle, $\theta_{0}$. As shown in Fig. 2b, the $\theta(V_{0})$ curves follow the same trend, with only a shift of the maximum to a value imposed by $\theta_{0}$. As shown in the inset, a plot of $\cos\theta(V_{0})-\cos\theta_{0}$ shows a linear dependence on $\eta$, which is in agreement with the theoretical prediction of Eq. (52). Fitting the simulation data to a straight line gives $c\approx 0.66$.

As the voltage in the droplet is increased, the apparent contact angle reaches a saturation value $\theta\approx 18.43^{\circ}$. The saturation effect was found to be independent of the wettability of the surface, and begins to occur when the droplet reaches $\theta\sim 50^{\circ}$. 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 Staicu and Mugele (2006).

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.

V.1 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 Oron et al. (1997),

[TABLE]

As shown by Eq. (55), the dynamics is driven by variations in the pressure within the film, $p_{\text{film}}$. This is composed of a capillary contribution, $2\gamma\kappa$, and by a contribution due to the electric stresses on the dielectric fluid, $-\frac{1}{2}\partial_{H}cV_{0}^{2}.$ For a gently curved interface, $\kappa\approx-\frac{1}{2}\partial^{2}_{x}H$. Hence,

[TABLE]

where we assume that the capacitance $c$ for a dielectric film in contact with the dielectric solid layer is given by

[TABLE]

We now study the stability the dielectric film by analysing Eq. (55) using a perturbative approach. Let us consider the initial sinusoidal interface profile

[TABLE]

where $H_{0}$ is the average height of the film, $a$ is the amplitude of the perturbation, $\lambda$ the wavelength and $\tau$ is the characteristic growth time.

Substituting Eq. (58) into Eq. (55), and assuming $a\ll H_{0}$ gives the dispersion relation

[TABLE]

where $\omega:=\mu_{\rm d}H_{0}/\gamma\tau$ is the dimensionless growth rate, and $k:=2\pi H_{0}/\lambda$ 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 $\omega=0$, corresponding to the onset of instability, gives the separatrix

[TABLE]

which gives the minimum electrowetting number for which a perturbation of given wave number leads to instability.

V.2 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

[TABLE]

with corresponds to a sinusoidal perturbation of amplitude $a=1$ and wavelength $\lambda=L_{x}$. The rest of the simulation parameters are reported in the last column of Table 1. 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 $\phi(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\pi 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 1).

Figs. 3(b) and 3(c) show the dispersion relations obtained from the simulations for $\eta=0$ and $\eta=0.03$. The data in the figures is reported in the dimensionless units of Eq. (59), where $\mu_{\rm d}$, $\gamma$, $H_{0}$ are fixed using the values reported in Table 1. For $\eta=0$, we observe the expected power-law decay, $\omega\propto-k^{4}$, predicted by the linear stability analysis. For $\eta=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\times 21$ points in the $\eta$–$k$ space. Fig. 3(d) shows the simulation results, which we present as a contour plot of $\omega$ vs $\eta$ and $k$. The separatrix, corresponding to the curve $\omega(k,\eta)=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 $\ell$. 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 ($\sim 0.08\ell$) 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}\sim 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 Jackson (1999). 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 $\theta_{0}<180^{\circ}$ 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. Andrew et al. (2016): the retracting edges collect fluid to form dewetting rims, which eventually merge to form a single circular droplet (see Fig. 4(b)). For $\theta_{0}=180^{\circ}$, 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. Staicu and Mugele (2006), 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 numerical 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.

Acknowledgements

The authors acknowledge support from EPSRC First-Grant Scheme (No. EP/P024408/1) and from EPSRC’s UK Consortium on Mesoscale Engineering Sciences (No. EP/R029598/1).

Appendix

Appendix A Parameters of the D2Q9 model

The set of velocities, $\{\mathbfit{c}_{q}\}_{q=0}^{Q-1}$, for the D2Q9 LBM is:

$\mathbfit{c}_{0}=(0,\,0),\qquad\qquad$ $\mathbfit{c}_{1}=(1,\,0),\qquad\qquad$ $\mathbfit{c}_{2}=(0,\,1),$

$\mathbfit{c}_{3}=(-1,\,0),$ $\mathbfit{c}_{4}=(0,\,-1),$ $\mathbfit{c}_{5}=(1,\,1),$

$\mathbfit{c}_{6}=(-1,\,1),$ $\mathbfit{c}_{7}=(-1,\,-1),\quad\text{ and }$ $\mathbfit{c}_{8}=(1,\,-1).$

The set of weights are: $w_{0}=4/9$, $w_{1\text{--}4}=1/9$, and $w_{5\text{--}8}=1/36$.

The tensor Hermite Polynomials are defined as Krüger et al. (2016),

[TABLE]

Explicitly, the first are,

[TABLE]

The collision matrix, $\Lambda_{kn}$ is defined as,

[TABLE]

where

[TABLE]

is the matrix of moments Krüger et al. (2016). $(M^{-1})_{kl}$, is the inverse matrix, i.e.,

[TABLE]

The matrix $S_{lm}$ contains the relaxation times, it is a diagonal matrix whose elements are given by

[TABLE]

where $\omega$ is specified by the dynamic viscosity d’Humières et al. (2002), $\mu$,

[TABLE]

Appendix B The evolution of a distribution function

Let us write the evolution of a distribution function, $f_{q}$, as

[TABLE]

where we have included $\Delta t$ as a free parameter for the present derivation that can be later set to unity. We now write the left-hand-side of Eq. (71) in a series expansion centered at $\Delta t=0$, i.e.,

[TABLE]

This implies that, the equilibrium distribution function, $f_{q}^{\text{eq}}$, can be written in terms of $f_{q}$,

[TABLE]

where we have defined the operator

[TABLE]

and $(\Lambda^{-1})_{qr}$ is the inverse of the collision matrix, which, from Eqs. (68) and (69), is given by,

[TABLE]

and, since $S$ is diagonal, the inverse is given by the reciprocal of its elements.

We can write an expression for $f_{q}$ in terms of $f_{q}^{\text{eq}}$ by inverting the relation in Eq. (73), i.e.,

[TABLE]

In this way, the collision step can be expressed in terms of the equilibrium distribution,

[TABLE]

Substituting Eq. (74) and rearranging terms, Eq. (77) becomes,

[TABLE]

Finally, we now set $\Delta t=1$, and for the case $\Lambda_{qr}=\delta_{qr}$ we obtain Eq. (36).