3.1 Theoretical modelling

The droplet is driven by surface tension, which acts in the direction of parallel to the plane prescribed by the instantaneous contact angle and perpendicular to the contact line. Here, we consider the change in contact angle due to dynamic wetting and EWOD as stated in (2.1). The dynamic contact angle $\unicode[STIX]{x1D703}_{d}$ can be estimated by using the Hoffmann–Voinov–Tanner law (Kistler Reference Kistler and Berg1993) as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{d}=({\unicode[STIX]{x1D703}_{s}}^{3}+k\mathit{Ca})^{1/3},\end{eqnarray}$$

where $k$ is an empirical constant that depends on the experimental conditions. The term $\mathit{Ca}$ is the capillary number defined as $\mathit{Ca}=\unicode[STIX]{x1D707}_{l}U_{cl}/\unicode[STIX]{x1D70E}$ , where $\unicode[STIX]{x1D707}_{l}$ is the dynamic viscosity of the liquid and $U_{cl}$ is the speed of the contact line. We considered the generality of analysis without the assistance of numerically simulated data, so we introduced the well-used empirical relation to describe the dynamic contact angle different from the numerical model used in our simulation described later. In this modelling, the effect of unknown information such as measurement scale is put into one fitting parameter $k$ . In this study, we assume $\mathit{Ca}\sim O(10^{-3})$ , so the approximation of $k\mathit{Ca}\ll {\unicode[STIX]{x1D703}_{s}}^{3}$ can be used. We take cosine of equation (3.1) and perform Taylor expansion with $k\mathit{Ca}\ll {\unicode[STIX]{x1D703}_{s}}^{3}$ , then we have

(3.2) $$\begin{eqnarray}\cos \unicode[STIX]{x1D703}_{d}=\cos \unicode[STIX]{x1D703}_{s}-\frac{k\sin \unicode[STIX]{x1D703}_{s}}{3{\unicode[STIX]{x1D703}_{s}}^{2}}\mathit{Ca}.\end{eqnarray}$$

We substitute a front and a rear contact angle into equation (3.2) while reversing the sign of the second term on the right-hand side if required because of the sign of $U_{cl}$ :

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D703}_{dA}=\cos \unicode[STIX]{x1D703}_{sA}-\frac{k\sin \unicode[STIX]{x1D703}_{sA}}{3\unicode[STIX]{x1D703}_{sA}^{2}}\mathit{Ca}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D703}_{dR}=\cos \unicode[STIX]{x1D703}_{sR}+\frac{k\sin \unicode[STIX]{x1D703}_{sR}}{3\unicode[STIX]{x1D703}_{sR}^{2}}\mathit{Ca}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{dA}$ and $\unicode[STIX]{x1D703}_{dR}$ are the front (advancing) and rear (receding) dynamic contact angles of the drop, respectively, and $\unicode[STIX]{x1D703}_{sA}$ and $\unicode[STIX]{x1D703}_{sR}$ are the advancing and receding static contact angles, respectively.

Although real droplet geometry is varied as shown in figure 5, we approximate it by simple shape as shown in figure 2 in our theoretical model. Here, for the purpose of simplicity, we assume a racetrack shape composed of two half-circles and two straight lines. From experimental observations (Lu et al. Reference Lu, Glasner, Bertozzi and Kim2007) and numerical results shown later in § 4, the droplet shrinks quickly widthways (within 2 s) after the voltage is imposed on the electrodes and then remains as the electrode width $W$ as shown in figure 2. Note that the wetted surface area in the figure should be the same as that of the initial circular droplet ( $\unicode[STIX]{x03C0}R^{2}$ ) because of volume conservation in the channel. The length of the straight sections can be determined simply from this fact.

Figure 2.

Droplet shape assumed in the present theoretical model.

The surface tension acts on the contact line on both the upper and lower walls of the channel, as shown in figure 1. As stated in § 2, the EWOD effect changes the static contact angle appearing on the lower wall at the positive electrode. In particular, the advancing contact angle $\unicode[STIX]{x1D703}_{sA}$ in (3.3) on the lower wall should be rewritten using (2.1). Taking both EWOD and dynamic effects into consideration, the surface-tension force acting on the contact line shown in figure 2 can now be estimated. When the droplet moves while keeping the same geometrical shape as shown in figure 2, the contact-line speed $U_{cl}$ at each position on the semi-circular parts may be related simply to the droplet speed $U_{drop}$ as

(3.5) $$\begin{eqnarray}U_{cl}=U_{drop}\cos \unicode[STIX]{x1D719},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ indicates the angle along the semi-circular parts of the contact line as shown in figure 2. The surface-tension forces acting on contact lines can be obtained by integration of infinitesimal force acting on a infinitesimal arc $(W/2)\,\text{d}\unicode[STIX]{x1D719}$ at the position $\unicode[STIX]{x1D719}$ for both top and bottom walls, and its $x$ –component is obtained by multiplying by $\cos \unicode[STIX]{x1D719}$ . Using the above (2.1), (3.3), (3.4) and (3.5), the surface-tension forces on the front and rear semi-circles can be calculated from the following integrals as

(front semi-circle)

(3.6) $$\begin{eqnarray}F_{sf}=\frac{\unicode[STIX]{x1D70E}W}{2}\int _{-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x03C0}/2}\left[2\left(\cos \unicode[STIX]{x1D703}_{sA}-\frac{k\sin \unicode[STIX]{x1D703}_{sA}}{3\unicode[STIX]{x1D703}_{sA}^{2}}\mathit{Ca}_{drop}\cos \unicode[STIX]{x1D719}\right)+\frac{cV^{2}}{2\unicode[STIX]{x1D70E}}\right]\cos \unicode[STIX]{x1D719}\,\text{d}\unicode[STIX]{x1D719}\end{eqnarray}$$

(rear semi-circle)

(3.7) $$\begin{eqnarray}F_{sr}=\frac{\unicode[STIX]{x1D70E}W}{2}\int _{-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x03C0}/2}2\left(\cos \unicode[STIX]{x1D703}_{sR}+\frac{k\sin \unicode[STIX]{x1D703}_{sR}}{3\unicode[STIX]{x1D703}_{sR}^{2}}\mathit{Ca}_{drop}\cos \unicode[STIX]{x1D719}\right)\cos \unicode[STIX]{x1D719}\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

where $\mathit{Ca}_{drop}$ is the non-dimensional drop speed ( $\unicode[STIX]{x1D707}_{l}U_{drop}/\unicode[STIX]{x1D70E}$ ). We note that the advancing (receding) dynamic contact angle should be associated with the front (rear) semi-circle because the contact line always moves so as to wet (dewet) the surface (Katoh et al. Reference Katoh, Azuma, Higashine and Miyamoto2006). In (3.6) and (3.7), the integration is taken over the semi-circular parts of the contact line. The surface tension on the lateral straight lines acts normal to the droplet velocity and does not contribute to the dynamics. The sum of (3.6) and (3.7) gives the resultant of surface-tension force as

(3.8) $$\begin{eqnarray}\displaystyle F_{s} & = & \displaystyle F_{sf}+F_{sr}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D70E}W\left[\frac{cV^{2}}{2\unicode[STIX]{x1D70E}}+2\left(\cos \unicode[STIX]{x1D703}_{sA}-\cos \unicode[STIX]{x1D703}_{sR}\right)-\frac{\unicode[STIX]{x03C0}k}{6}\left(\frac{\sin \unicode[STIX]{x1D703}_{sA}}{\unicode[STIX]{x1D703}_{sA}^{2}}+\frac{\sin \unicode[STIX]{x1D703}_{sR}}{\unicode[STIX]{x1D703}_{sR}^{2}}\right)\mathit{Ca}_{drop}\right].\end{eqnarray}$$

The terms in the bracket refer to the electrowetting effect, resistance due to contact-angle hysteresis for top/bottom walls and dynamic contact-angle effect at front/rear of top/bottom walls.

Next, we consider the forces resisting the droplet motion, namely the inertial force and the viscous drag on the wall. The inertial force can be estimated by considering the moving mass per unit time and velocity as

(3.9) $$\begin{eqnarray}F_{inert}\sim \unicode[STIX]{x1D70C}_{l}WbU_{drop}^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{l}$ is the liquid density. To estimate the viscous drag, we use the time scale for the development of a viscous boundary layer as given by Stokes’ first problem (Schlichting Reference Schlichting1975). The time scale $t_{b}$ on which the boundary layer reaches the channel centre height $b/2$ can be evaluated as

(3.10) $$\begin{eqnarray}4\sqrt{\frac{\unicode[STIX]{x1D707}_{l}t_{b}}{\unicode[STIX]{x1D70C}_{l}}}=\frac{b}{2}.\end{eqnarray}$$

If $t_{b}$ is much smaller than the time scale $b/U_{drop}$ of the droplet motion, namely

(3.11) $$\begin{eqnarray}\frac{t_{b}U_{drop}}{b}\ll 1,\end{eqnarray}$$

we can safely assume that the velocity distribution in the droplet develops as a two-dimensional laminar channel flow, namely

(3.12) $$\begin{eqnarray}u(y)=\frac{6U_{drop}}{b^{2}}\left(by-y^{2}\right).\end{eqnarray}$$

Using a Reynolds number based on the gap width, namely $Re=\unicode[STIX]{x1D70C}_{l}bU_{drop}/\unicode[STIX]{x1D707}_{l}$ , the time scale condition (3.11) is rewritten as

(3.13) $$\begin{eqnarray}\frac{t_{b}U_{drop}}{b}=\frac{\mathit{Re}}{64}\ll 1,\end{eqnarray}$$

and this condition is satisfied in general Hele-Shaw systems for microfluidics ( $\mathit{Re}\lesssim O(1)$ ). The viscous drag can then be obtained over the whole wetted surface area ( $2\unicode[STIX]{x03C0}R^{2}$ ) as

(3.14) $$\begin{eqnarray}F_{v}=\frac{12\unicode[STIX]{x03C0}R^{2}\unicode[STIX]{x1D707}_{l}U_{drop}}{b}=\frac{12\unicode[STIX]{x03C0}R^{2}\unicode[STIX]{x1D70E}}{b}\mathit{Ca}_{drop}.\end{eqnarray}$$

When the inertial force is much smaller than the viscous force, namely

(3.15) $$\begin{eqnarray}\frac{F_{inert}}{F_{v}}\sim \frac{\unicode[STIX]{x1D70C}_{l}bU_{drop}}{\unicode[STIX]{x1D707}_{l}}\frac{b}{W}\sim Re\frac{b}{W}\ll 1,\end{eqnarray}$$

(this condition is also satisfied in general Hele-Shaw systems for microfluidics ( $\mathit{Re}\lesssim O(1)$ and $b\ll W$ )), we can consider a simple force balance between surface tension and viscous drag as a good approximation for the steady state. From equations (3.8) and (3.14), the non-dimensional droplet speed in the channel can be obtained as

(3.16) $$\begin{eqnarray}\mathit{Ca}_{drop}=\frac{\displaystyle \frac{cV^{2}}{2\unicode[STIX]{x1D70E}}+2\left(\cos \unicode[STIX]{x1D703}_{sA}-\cos \unicode[STIX]{x1D703}_{sR}\right)}{\displaystyle \frac{12\unicode[STIX]{x03C0}R^{2}}{bW}+\frac{\unicode[STIX]{x03C0}k}{6}\left(\frac{\sin \unicode[STIX]{x1D703}_{sA}}{\unicode[STIX]{x1D703}_{sA}^{2}}+\frac{\sin \unicode[STIX]{x1D703}_{sR}}{\unicode[STIX]{x1D703}_{sR}^{2}}\right)}.\end{eqnarray}$$

The first term in the numerator represents the electrowetting effect and the second term represents the contact-angle hysteresis. The first term in the denominator is a non-dimensional geometrical parameter representing the effect of both driving and resistance force and the second term represents the dynamic wetting. Equation (3.16) shows that the droplet speed is determined by these four effects.

Figure 3.

Droplet speed as estimated by the present theoretical model. Experimental data are from Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007).

3.2 Model evaluation using experimental results

Here, we evaluate the present theoretical model by using the experimental results of Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007). In their paper, the measurement condition of the equilibrium contact angle ( $100^{\circ }$ ) is unclear, so instead we estimate the equilibrium contact angle $\unicode[STIX]{x1D703}_{0}$ for our model from their data on the advancing and receding contact angles. From figure 6 in Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007), the steady dynamic advancing angle of the top plate (without the electrowetting effect) is $111^{\circ }$ and the steady dynamic receding angle (averaged from the top and bottom plates) is $96^{\circ }$ . Hence, we set $\unicode[STIX]{x1D703}_{0}$ as the average of those two values, namely $103.5^{\circ }$ . We also assume that the static advancing and receding contact angles are shifted equally by $\unicode[STIX]{x1D6FC}$ from $\unicode[STIX]{x1D703}_{0}$ :

(3.17a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{sA}=\unicode[STIX]{x1D703}_{0}+\unicode[STIX]{x1D6FC},\quad \unicode[STIX]{x1D703}_{sR}=\unicode[STIX]{x1D703}_{0}-\unicode[STIX]{x1D6FC}.\end{eqnarray}$$

From figure 4 in Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007), the minimum voltage at which the droplet moves can be estimated as $V_{0}=36~\text{V}$ , at which we have $\mathit{Ca}_{drop}=0$ . Hence, we obtain the following from equation (3.16):

(3.18) $$\begin{eqnarray}0=\frac{cV_{0}^{2}}{2\unicode[STIX]{x1D70E}}+2\left[\cos (\unicode[STIX]{x1D703}_{0}+\unicode[STIX]{x1D6FC})-\cos (\unicode[STIX]{x1D703}_{0}-\unicode[STIX]{x1D6FC})\right].\end{eqnarray}$$

By solving (3.18), we obtain $\unicode[STIX]{x1D6FC}=5.1^{\circ }$ , then (3.17a,b ) yields $\unicode[STIX]{x1D703}_{sA}=108.6^{\circ }$ and $\unicode[STIX]{x1D703}_{sR}=98.4^{\circ }$ . Furthermore, from figure 11 in Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007), we can estimate the capillary number $\mathit{Ca}$ at which the dynamic contact angles achieve steady values as being $\mathit{Ca}=1.23\times 10^{-3}$ . Hence, we can determine the coefficient $k$ as $k=335$ by averaging values obtained from (3.3) and (3.4).

From the above values and the experimental conditions of Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007) ( $c=3.6\times 10^{-5}~\text{F m}^{-2}$ , $\unicode[STIX]{x1D70E}=0.067~\text{N m}^{-1}$ , $R=5.5\times 10^{-3}~\text{m}$ , $b=5.3\times 10^{-4}~\text{m}$ , $W=1\times 10^{-2}~\text{m}$ , $\unicode[STIX]{x1D70C}_{l}=1.28\times 10^{3}~\text{kg m}^{-3}$ and $\unicode[STIX]{x1D707}_{l}=0.02~\text{Pa s}$ ), we obtain the relationship shown in figure 3 between average droplet speed and applied voltage. The conditions imposed by (3.11) and (3.15) are satisfied because $t_{b}U_{drop}/b=O(10^{-2})$ and $F_{inert}/F_{v}\sim O(10^{-3})$ with $U_{drop}\sim 1\times 10^{-2}~\text{m}~\text{s}^{-1}$ . Up to 60 V, the present theoretical model (with racetrack shape assumption) can reproduce the experimentally measured relation very well. Above 60 V, Young–Lippmann equation (2.1) was not realised as shown in figure 2 in Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007). Increasing of the cosine of the contact angle with voltage saturates around 60 V, and it decreases above 70 V. Such contact angle saturation may be due to the dielectric break down, which has been reported by many researchers, e.g. Vallet, Berge & Vovelle (Reference Vallet, Berge and Vovelle1996), Blake et al. (Reference Blake, Clarke and Stattersfield2000). As a natural result of the contact angle increasing above 70 V, the droplet speed decreases in the experimental observation. The behaviour with the saturation voltage is outside of the scope of our interest. We only discuss the results below the saturation voltage 60 V.

In the theoretical droplet velocity described by (3.16), some experimental parameters, i.e. $\unicode[STIX]{x1D703}_{0}$ , $\unicode[STIX]{x1D6FC}$ and $k$ in (3.18) and (3.1) were assumed by referring the results shown in Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007). To examine the effect of parameter estimation error, $\unicode[STIX]{x1D703}_{0}$ was varied with $\pm 2^{\circ }$ . For each $\unicode[STIX]{x1D703}_{0}$ , $\unicode[STIX]{x1D6FC}$ and $k$ are estimated by (3.18) and (3.1), respectively. With a $2^{\circ }$ estimation error of $\unicode[STIX]{x1D703}_{0}$ , the fitted $k$ is varied by approximately 10 %, however, the theoretically estimated droplet speed is not so affected (below 5 % for maximum speed at 60 V). On the other hand, to examine the effect of the droplet shape assumption, as shown in figure 2, we also considered a rectangular droplet with a width of $W$ as an extreme example. For the rectangular droplet case, the effect of viscous drag is not changed from the racetrack case because of the constant wetted area, however, the integral of the contact-angle’s dynamic effect term in (3.6) and (3.7) depends on the contact-line shape. Then, as a result, $\unicode[STIX]{x03C0}$ in the second term in the denominator of (3.16) for the racetrack case is replaced by 4 for the rectangular case. In figure 3, the theoretical curve with the rectangular shape assumption is also shown. The effect of the shape assumption is not so large, the difference is approximately 8 % for a maximum speed of 60 V in the present sample cases. The present theoretical model can express the characteristic well. Later, we discuss the present model in the context of numerically simulated data with various values of the shape parameters.

4.1 Front-tracking method

We use the numerical procedure first introduced by Yamamoto et al. (Reference Yamamoto, Higashida, Tanaka, Wakimoto, Ito and Katoh2016) except that we represent the change in wettability due to the applied voltage. In the front-tracking method (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992), the gas–liquid interface is tracked by marker points and the interfacial information is calculated by triangular elements connecting the markers. Using this information, the single-fluid equations are solved using the finite-difference method. The interval of advected marker points is not uniform; therefore, the interval should be adjusted by marker addition, deletion or reconnection. To avoid this intricate process, we employed the level contour reconstruction method (LCRM) (Shin & Juric Reference Shin and Juric2002), which can treat the interface elements without logical connectivity. See Shin & Juric (Reference Shin and Juric2002) for further details.

The basic equations are the incompressible Navier–Stokes and continuity equations:

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})\right]=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\unicode[STIX]{x1D707}\left(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}}\right)\right]+\boldsymbol{F}+\unicode[STIX]{x1D70C}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ and $p$ are the fluid velocity and pressure, respectively. The density $\unicode[STIX]{x1D70C}$ and dynamic viscosity $\unicode[STIX]{x1D707}$ are given using an indicator function obtained from the unit normal of the interfacial elements. The body force $\boldsymbol{F}$ captures the effect of interfacial tension.

The flow field is determined by solving (4.1) and (4.2) on fixed grid points and with the assistance of the marker information for $\unicode[STIX]{x1D70C}$ , $\unicode[STIX]{x1D707}$ and $\boldsymbol{F}$ . Equations (4.1) and (4.2) are discretised in space by a second-order central finite-difference method on a staggered grid system. A second-order Adams–Bashforth method is applied to advance the advection term in time, and a second-order Crank–Nicolson scheme is applied for the viscous term. The pressure is linked to (4.2) by the simplified marker and cell algorithm, and the Poisson equation for pressure correction is solved using a multigrid method.

4.2 Boundary condition for moving contact line

By the generalised Navier boundary condition (GNBC; Qian et al. (Reference Qian, Wang and Sheng2003)) for a moving contact line in the front-tracking procedure (Yamamoto et al. Reference Yamamoto, Higashida, Tanaka, Wakimoto, Ito and Katoh2016), the non-dimensional slip speed of the contact line, namely the capillary number $\mathit{Ca}$ , is given as

(4.3) $$\begin{eqnarray}\mathit{Ca}=\unicode[STIX]{x1D712}(\cos \unicode[STIX]{x1D703}_{s}-\cos {\unicode[STIX]{x1D703}_{d}}^{micro}),\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{s}$ is the static contact angle, ${\unicode[STIX]{x1D703}_{d}}^{micro}$ is the microscopic dynamic contact angle and $\unicode[STIX]{x1D712}$ is the non-dimensional slip parameter that represents the dynamic property of wetting (Yamamoto et al. Reference Yamamoto, Ito, Wakimoto and Katoh2013, Reference Yamamoto, Tokieda, Wakimoto, Ito and Katoh2014).

In the front-tracking method, the obtained contact angles are measured macroscopically at the grid scale $\unicode[STIX]{x1D6E5}$ using the triangular elements. We use the relationship between the macroscopic contact angle (grid scale in this case) $\unicode[STIX]{x1D703}_{d}^{grid}$ and the microscopic angle ${\unicode[STIX]{x1D703}_{d}}^{micro}$ derived by Cox using a matched asymptotic expansion (Cox Reference Cox1986)

(4.4) $$\begin{eqnarray}G_{q}(\unicode[STIX]{x1D703}_{d}^{grid})=G_{q}({\unicode[STIX]{x1D703}_{d}}^{micro})+\mathit{Ca}\ln \left(\frac{\unicode[STIX]{x1D6E5}}{l^{micro}}\right),\end{eqnarray}$$

where $l^{micro}$ is the microscopic characteristic length scale. The function $G_{q}$ is represented as

(4.5) $$\begin{eqnarray}G_{q}(\unicode[STIX]{x1D703})=\int _{0}^{\unicode[STIX]{x1D703}}\frac{\text{d}\unicode[STIX]{x1D703}}{f_{q}(\unicode[STIX]{x1D703})},\end{eqnarray}$$

where

(4.6) $$\begin{eqnarray}f_{q}(\unicode[STIX]{x1D703})=\frac{2\sin \unicode[STIX]{x1D703}\{q^{2}(\unicode[STIX]{x1D703}^{2}-\sin ^{2}\unicode[STIX]{x1D703})+2q[\unicode[STIX]{x1D703}(\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703})+\sin ^{2}\unicode[STIX]{x1D703}]+[(\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703})^{2}-\sin ^{2}\unicode[STIX]{x1D703}]\}}{q(\unicode[STIX]{x1D703}^{2}-\sin ^{2}\unicode[STIX]{x1D703})[(\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703})+\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}]+[(\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703})^{2}-\sin ^{2}\unicode[STIX]{x1D703}](\unicode[STIX]{x1D703}-\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703})},\end{eqnarray}$$

and $q$ is the viscosity ratio $\unicode[STIX]{x1D707}_{g}/\unicode[STIX]{x1D707}_{l}$ . In Cox’s analysis, the flow field is divided into microscopic, macroscopic and intermediate scales. In the macroscopic (outer) and intermediate regions, the no-slip boundary condition is considered. In the microscopic (inner) region, the slip boundary condition is allowed; however, the no-slip condition is considered at the outer limit of the inner region. At all scales, the flow field is expanded as perturbations of the $\mathit{Ca}$ and the Stokes equations for the zeroth-order field are solved. The normal stress conditions obtained from the zeroth-order field are connected to the capillary effect of interfacial deformation. By matching the inner and outer regions via the intermediate region, the macroscopic dynamic contact angle is related to the microscopic angle on the wall. Cox (Reference Cox1986) did not consider the full microscopic behaviour, only the outer limit of the inner region. Cox viewed the microscopic contact angle as being static; however, we consider the microscopic angle as being dynamic and described by the GNBC.

To date, similar modelling efforts have been proposed for representing the unresolved microscopic behaviour (Afkhami, Zaleski & Bussmann Reference Afkhami, Zaleski and Bussmann2009; Dupont & Legendre Reference Dupont and Legendre2010; Sui & Spelt Reference Sui and Spelt2013; Sui, Ding & Spelt Reference Sui, Ding and Spelt2014). However, all of them consider a static contact angle as a microscopic contact angle. Sui & Spelt (Reference Sui and Spelt2013) compared their simulations with experimental data and showed that considering the microscopic contact angle as static leads to unrealistic microscopic length scales. In the method proposed by Yamamoto et al. (Reference Yamamoto, Tokieda, Wakimoto, Ito and Katoh2014), the unresolved microscopic behaviour is treated in macroscopic-scale simulations using Cox’s relation. However, a dynamic microscopic angle is considered and the GNBC model is applied to a microscopic dynamic contact angle and slip-velocity relationship. For theoretical analysis, Ledesma-Aguilar, Hernández-Machado & Pagonabarraga (Reference Ledesma-Aguilar, Hernández-Machado and Pagonabarraga2013) has applied the same type of modelling using a combination of (4.3) and (4.4) to study fluid entrainment.

For the boundary condition of the wall, (i) the grid-scale contact angle $\unicode[STIX]{x1D703}_{d}^{grid}$ is calculated from the triangular elements contacting the wall, (ii) the microscopic dynamic contact angle ${\unicode[STIX]{x1D703}_{d}}^{micro}$ and $\mathit{Ca}$ are calculated by solving (4.3) and (4.4) simultaneously using Newton’s iteration method with numerical integration of (4.5) and (iii) the wall velocity boundary condition is given by using $\mathit{Ca}$ and a two-dimensional delta-function distribution (Yamamoto et al. Reference Yamamoto, Higashida, Tanaka, Wakimoto, Ito and Katoh2016). Although strictly speaking the new parameter $l^{micro}$ should be predetermined in each case, a value of $10^{-9}~\text{m}$ is used generally, as discussed by Yamamoto et al. (Reference Yamamoto, Tokieda, Wakimoto, Ito and Katoh2014).

4.3 Representation of electrowetting and contact-angle hysteresis

To incorporate the fact that the droplet remains static until a threshold voltage, as shown in figure 3, contact-angle hysteresis ( $\unicode[STIX]{x1D703}_{sA}\neq \unicode[STIX]{x1D703}_{sR}$ ) should be represented. In this study, we do so by simply partitioning the grid-scale contact angle as

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{s}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D703}_{sA}\quad & \text{for }\unicode[STIX]{x1D703}_{d}^{grid}>\unicode[STIX]{x1D703}_{sA}\\ \unicode[STIX]{x1D703}_{sR}\quad & \text{for }\unicode[STIX]{x1D703}_{d}^{grid}<\unicode[STIX]{x1D703}_{sR},\end{array}\right.\end{eqnarray}$$

for (4.3). For $\unicode[STIX]{x1D703}_{sR}\leqslant \unicode[STIX]{x1D703}_{d}^{grid}\leqslant \unicode[STIX]{x1D703}_{sA}$ , the element does not affect the slip velocity. If the current contact line is on the positive electrode, variations of $\unicode[STIX]{x1D703}_{sA}$ and $\unicode[STIX]{x1D703}_{sR}$ expressed as (2.1) are taken into account to represent the electrowetting effect.

4.4 Simulation conditions

The number of grid points in the wall-normal direction for the base case (corresponding to the experiment by Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007), $b=5.3\times 10^{-4}~\text{m}$ ) was set to 8, so the grid interval $\unicode[STIX]{x1D6E5}$ was $6.625\times 10^{-4}~\text{m}$ . A uniform cubic grid system was used and the grid intervals for the other cases were also fixed to the above value except for the grid-convergence tests. The numbers of grid points in the $x$ and $z$ directions were $N_{x}=456$ and $N_{z}=160$ , respectively. The domain size in the $x$ – $z$ plane, $L_{x}=N_{x}\unicode[STIX]{x1D6E5}$ by $L_{z}=N_{z}\unicode[STIX]{x1D6E5}$ , was fixed in all cases. The boundary condition for the wall was described in § 4.2 and 4.3. The gradient-free condition was applied at $x=0,L_{x}$ and $z=L_{z}$ and the symmetric condition was applied at $z=0$ . The time interval was set to $5\times 10^{-5}~\text{s}$ . To allow the initial state to achieve equilibrium, 1000 computational steps were taken prior to the voltage being applied at $t=0$ . In the half- and double-grid-size computations for the grid-convergence tests, the time increments were changed to half and twice the basic resolution, respectively. The liquid properties were as given in § 3.2, and we used a density $\unicode[STIX]{x1D70C}_{g}=1.21~\text{kg m}^{-3}$ and dynamic viscosity $\unicode[STIX]{x1D707}_{g}=1.82\times 10^{-5}~\text{Pa s}$ for the gas. The length of the electrodes on the bottom plate $L$ was fixed at 0.01 m ( $151\unicode[STIX]{x1D6E5}$ ) for all cases. For the base case of the parametric study, the supplied voltage $V$ was set to 50.4 V and the other parameters ( $\unicode[STIX]{x1D703}_{sA}$ , $\unicode[STIX]{x1D703}_{sR}$ , $c$ and $\unicode[STIX]{x1D70E}$ ) were fixed at the values given in § 3.2. Given the symmetry in the $z$ direction, we used $W=150\unicode[STIX]{x1D6E5}$ for the base case and rounded it to the closest even number of grid points in the other cases.

Figure 4.

Position of droplet centre of mass with time for various values of slip parameter $\unicode[STIX]{x1D712}$ and resolution. Experimental data are from Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007).

4.5 Slip-parameter identification and grid-convergence tests

The non-dimensional slip parameter $\unicode[STIX]{x1D712}$ used in the GNBC front tracking can be estimated approximately using (4.3) and (4.4) with experimental values of $\mathit{Ca}_{drop}$ and $\unicode[STIX]{x1D703}_{dA}$ or $\unicode[STIX]{x1D703}_{dR}$ . However, the contact-line elements of the moving droplet move at various speeds depending on their position, so we must tune the parameter $\unicode[STIX]{x1D712}$ about its estimated value. Figure 4 shows the position of the droplet centre of mass with time for three values of $\unicode[STIX]{x1D712}$ with resolution $N_{y}=8$ , and grid resolutions of $N_{y}=4$ and 16 with $\unicode[STIX]{x1D712}=0.023$ . The position as shown has been shifted by $L_{x}/2$ and normalised by $R$ , and the time as shown has been normalised by $12\unicode[STIX]{x1D707}_{l}R^{3}/(\unicode[STIX]{x1D70E}b^{2})=2.12~\text{s}$ .

For large $\unicode[STIX]{x1D712}$ , the contact angle rapidly reaches the small static angle affected by EWOD, then the droplet moves away quickly. With $\unicode[STIX]{x1D712}=0.023$ , the simulation results coincide with the experimental data, hence we use this value hereafter.

In the grid-convergence tests, a higher resolution causes the droplet to move more slowly, which may be because the drag force is being evaluated more accurately. We can see the convergence tendency with grid resolution from figure 4. Although higher-resolution computation is desirable, we must use a cubic grid of the same size along every axis to avoid the effect of surface tension being distributed anisotropically, and this imposes huge memory demands. For higher resolution, a more suitable value for $\unicode[STIX]{x1D712}$ may be found. However, the difference of the results between $N_{y}=8$ and $N_{y}=16$ with constant $\unicode[STIX]{x1D712}$ is not so large and the variation of $\unicode[STIX]{x1D712}$ is almost proportional to variation of the droplet speed. So we consider that the suitable value of $\unicode[STIX]{x1D712}$ for large $N_{y}$ is not so greatly different from the $N_{y}=8$ case. Thus, we decided to use $N_{y}=8$ hereafter, which generally yields appropriate results.

Figure 5.

Composite images of moving droplet: simulations (right half of each panel) and experimental observation (left half of each panel), the latter from Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007).

Figure 5 shows composite images of a moving droplet simulated with $\unicode[STIX]{x1D712}=0.023$ and $N_{y}=8$ and paired with experimental observations. The left half of each panel is an experimental observation by Lu et al. (Reference Lu, Glasner, Bertozzi and Kim2007) whereas the right half is our simulation result. On the right-hand side, the grey region represents the (half) liquid droplet and the red (upper side) rectangle represents the positive electrode. Not only do the present simulations reproduce the droplet motion well (figure 4), they also capture the droplet shape very realistically (figure 5).