A finite element method for Electrowetting on Dielectric

The mathematical model is an extension of the contact line model proposed by Ren et al. [33, 35, 34] to include the contribution of electric field in EWoD. It consists of the two-phase incompressible Navier-Stokes equations for the fluid dynamics and the Maxwell’s equation for the electrostatic potential.

First we consider the two-phase fluid dynamics in the domain $\Omega=\Omega_{1}(t)\cup\Omega_{2}(t)$. Let $\mathbf{u}(\mathbf{x},t):\;\Omega\times[0,~{}T]\to\mathbb{R}^{2}$ be the fluid velocity and $p(\mathbf{x},t):\;\Omega\times[0,T]\to\mathbb{R}$ be the pressure. The dynamic of the system is governed by the standard incompressible Navier-Stokes equations in $\Omega_{i}(t)(i=1,2)$,

where $\tau_{d}=2\mu_{i}D(\mathbf{u})$ is the viscous stress with $D(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})$.

On the fluid interface $\Gamma(t)$, we have the following conditions

where $[\cdot]_{1}^{2}$ denotes the jump from fluid 1 to fluid 2, $\mathbf{I}_{2}\in\mathbb{R}^{2\times 2}$ is the identity matrix, $\kappa$ is the curvature of the fluid interface $\Gamma$, $\Phi$ is the electrostatic potential, and $v_{n}$ denotes the normal velocity of the fluid interface. Equation (4a) states the fluid velocity is continuous across the interface. Equation (4b) states that the tangential stress is continuous across the interface and the normal stress jump is balanced by the curvature force $\gamma\kappa\mathbf{n}$ together with the electric force $\frac{\epsilon_{2}}{2}|\nabla\Phi|^{2}\mathbf{n}$. Equation (4c) is the kinetic condition for the interface, where the fluid interface evolves according to the local fluid velocity.

At the lower solid wall $\Gamma_{i}(t)\ (i=1,2)$, the fluid velocity satisfies the no-penetration condition and the Navier boundary condition

where $\beta_{i}\ (i=1,2)$ is the friction coefficient of fluid $i$ at the wall, and $u_{s}=\mathbf{u}\cdot\mathbf{t}_{w}$ is the slip velocity of the fluids.

At the contact points $x_{l}$ and $x_{r}$, the dynamic contact angles, denoted by $\theta_{d}^{l}$ and $\theta_{d}^{r}$ respectively, depend on the contact line velocity [33, 35],

where $\theta_{Y}$ is the equilibrium contact angle satisfying the Young’s relation (1), and $\beta^{*}$ is the friction coefficient of the fluid interface at the contact line. The contact line velocity is determined by the slip velocity of the fluids: $\dot{x}_{l,r}=\left.u_{s}\right|_{x=x_{l,r}}$. The contact angle condition (6) is a force balance at the contact point: the Young stress resulted from the deviation of the dynamic contact angle from the equilibrium angle is balanced by the friction force. Note that the electric force does not contribute to the force balance at the contact line.

Furthermore, we use the no-slip boundary condition on the upper wall $\Gamma_{4}$ and the periodic boundary conditions along $\Gamma_{3}$.

The applied voltage induces electrostatic fields in the fluid region $\Omega_{2}$ and the dielectric substrate $\Omega_{3}$, where the electric potential, denoted by $\Phi$, satisfies the Laplace equation,

On the interface $\Gamma\cup\Gamma_{1}$ and at the top/bottom electrode $\Gamma_{4}\cup\Gamma_{5}$, the electric potential satisfies the Dirichlet boundary conditions

where $\phi>0$ is the imposed voltage difference. Since there is no free charge density across the interface between the dielectric fluid and the dielectric substrate, we have

where $\left[\cdot\right]_{2}^{3}$ denotes the jump from $\Omega_{2}$ to $\Omega_{3}$, $\epsilon=\epsilon_{2}$ in $\Omega_{2}$ and $\epsilon=\epsilon_{3}$ in $\Omega_{3}$. Furthermore, the periodic boundary condition is prescribed along $\Gamma_{3}\cup\Gamma_{6}$.

Next we present the above governing equations and the boundary/interface conditions in their dimensionless form. By choosing $L$ and $U$ as the characteristic length and velocity respectively, we rescale the physical quantities as

We define the Reynolds number $Re$, the Capillary number $Ca$, the slip length $l_{s}$, the Weber number $We$ and the parameter $\eta$ as

where $\eta$ measures the relative strength of the electric force compared to the surface tension at the fluid interface.

For ease of presentation, we will drop the hats on the dimensionless parameters and variables. In the dimensionless form, the governing equations for the fluid dynamics in $\Omega_{i}\ (i=1,2)$ read

where $\mathbf{T}=p\mathbf{I}_{2}-\frac{1}{Re}\tau_{d}$ is the stress tensor. These equations are supplemented with the following boundary/interface conditions:

• (i)

  The interface conditions on $\Gamma(t)$:

  (12a)		$\displaystyle\bigl{[}\mathbf{u}\bigr{]}^{2}_{1}=\mathbf{0},$
  (12b)		$\displaystyle We\bigl{[}\mathbf{T}\bigr{]}_{1}^{2}\cdot\mathbf{n}=\left(\kappa+\epsilon_{2}\eta\left|\nabla\Phi\right|^{2}\right)\mathbf{n},$
  (12c)		$\displaystyle\kappa=\partial_{ss}\mathbf{X}\cdot\mathbf{n},$
  (12d)		$\displaystyle v_{n}=\mathbf{u}|_{\Gamma(t)}\cdot\mathbf{n},$

  where $\mathbf{X}$ denotes the fluid interface, $s$ is the arc length parameter with $s=0$ being at the left contact point.

• (ii)

  The condition for the dynamic contact angles:

  (13a)		$\displaystyle\frac{1}{Ca}(\cos\theta_{d}^{l}-\cos\theta_{Y})=\beta^{*}\dot{x}_{l}(t),$
  (13b)		$\displaystyle\frac{1}{Ca}(\cos\theta_{d}^{r}-\cos\theta_{Y})=-\beta^{*}\dot{x}_{r}(t).$

• (iii)

  The boundary conditions on $\Gamma_{1}(t)\cup\Gamma_{2}(t)$:

  (14)

  $$\mathbf{u}\cdot\mathbf{n}_{w}=0,\quad l_{s}\mathbf{t}_{w}\cdot\tau_{d}\cdot\mathbf{n}_{w}=-\beta_{i}u_{s}.$$

• (iv)

  The no-slip condition on $\Gamma_{4}$

  (15)

  $$\mathbf{u}=\mathbf{0}.$$

• (v)

  Periodic boundary conditions on $\Gamma_{3}$:

  (16)

  $$\left.\mathbf{u}\right|_{x=-L_{x}}=\left.\mathbf{u}\right|_{x=L_{x}},\quad\left.\mathbf{T}\right|_{x=-L_{x}}=\left.\mathbf{T}\right|_{x=L_{x}}.$$

The governing equations for the electric field potential $\Phi(x,t)$ in $\Omega_{i}(i=2,3)$ read

subject to the following boundary conditions:

• (i)

  The Dirichlet Boundary conditions on $\Gamma$, $\Gamma_{1}$, $\Gamma_{4}$ and $\Gamma_{5}$:

  (18a)		$\displaystyle\Phi=1,\quad{\rm on}\ \Gamma(t)\cup\Gamma_{1}(t),$
  (18b)		$\displaystyle\Phi=0,\quad{\rm on}\ \Gamma_{4},\ \Gamma_{5}.$

• (ii)

  The interface condition on $\Gamma_{2}(t)$:

  (19)

  $$[\Phi]_{2}^{3}=0,\qquad\mathbf{n}_{w}\cdot\left[\epsilon\nabla\Phi\right]_{2}^{3}=0.$$

• (iii)

  Periodic boundary conditions on $\Gamma_{3}\cup\Gamma_{6}$:

  (20)

  $$\left.\Phi\right|_{x=-L_{x}}=\left.\Phi\right|_{x=L_{x}},\quad\left.\frac{\partial\Phi}{\partial\mathbf{n}}\right|_{x=-L_{x}}=-\left.\frac{\partial\Phi}{\partial\mathbf{n}}\right|_{x=L_{x}},$$

  where $\mathbf{n}$ is the unit outward normal to $\Omega_{2}$ on $\Gamma_{3}$.

Equations (11) and (17) together with the boundary/interface conditions (12)-(16) and (18)-(20) form the complete model for the problem of EWoD. For this dynamical system, we define the total energy as

where $|\Gamma_{1}(t)|$ and $|\Gamma(t)|$ are the total arc length of $\Gamma_{1}(t)$ and $\Gamma(t)$, respectively. The four terms represents the kinetic energy of the fluids, the interfacial energy at the solid wall, the interfacial energy of the fluid interface and the electrical energy, respectively. The dynamical system obeys the energy dissipation law

where the three terms represent the viscous dissipation in the bulk fluid with $\|\cdot\|$ being the Frobenius norm, the dissipation at the solid wall due to the friction and the dissipation at the contact points, respectively. Details for the derivation of (22) are provided in the Appendix A.

To propose the weak formulation for equations (11)-(16), we define the function space for the pressure and the velocity, respectively as

with the $L^{2}$-inner product on $\Omega=\Omega_{1}(t)\cup\Omega_{2}(t)$

We parameterize the fluid interface as $\mathbf{X}(\alpha,t)=(X(\alpha,t),Y(\alpha,t))^{T}$, where $\alpha\in{\mathbb{I}}=[0,1]$, and $\alpha=0,1$ corresponding to the left and and right contact points, respectively. We define the function spaces for the interface curvature and the interface as

equipped with the $L^{2}$-inner product on $\mathbb{I}$

Taking the inner product of (11a) with $\boldsymbol{\omega}\in\mathbb{U}$ then using the boundary/interface conditions in (12)-(16) and $\nabla\cdot\mathbf{u}=0$, we obtain [3, 43]

where $\rho=\rho_{1}\chi_{{}_{\Omega_{1}(t)}}+\rho_{2}\chi_{{}_{\Omega_{2}(t)}}$ with $\chi$ being the characteristic function. With the special treatment of the inertia term in (27), the numerical scheme enjoys the discrete stability for the fluid kinetic energy in the absence of electric field. This will be discussed later in section 3.4. For the viscous term, integrating by parts then applying the boundary/interface conditions yields

where $\mu=\mu_{1}\chi_{{}_{\Omega_{1}(t)}}+\mu_{2}\chi_{{}_{\Omega_{2}(t)}}$, $\beta=\beta_{1}\chi_{{}_{\Gamma_{1}(t)}}+\beta_{2}\chi_{{}_{\Gamma_{2}(t)}}$, and $\omega_{s}=\boldsymbol{\omega}\cdot\mathbf{t}_{w}$.

Equation (12c) for the curvature can be rewritten as $\kappa\mathbf{n}=\partial_{ss}\mathbf{X}$. Taking the inner product of it with $\boldsymbol{g}=(g_{1},~{}g_{2})^{T}\in\mathbb{X}$ on $\Gamma(t)$ then applying integration by parts yields

where we have used the fact that $\partial_{s}X|_{\alpha=0}=\cos\theta_{d}^{l}$, $\partial_{s}X|_{\alpha=1}=\cos\theta_{d}^{r}$, and the dynamic contact angle conditions in (13).

Combining these results, we obtain the weak formulation for the dynamic system (11)-(16) as follows: Given the initial fluid velocity $\mathbf{u}_{0}$ and the interface $\mathbf{X}_{0}(\alpha)$, find the fluid velocity $\mathbf{u}(\cdot,~{}t)\in\mathbb{U}$, the pressure $p(\cdot,~{}t)\in\mathbb{P}$, the fluid interface $\mathbf{X}(\cdot,~{}t)\in\mathbb{X}$, and the curvature $\kappa(\cdot,~{}t)\in\mathbb{K}$ such that

Eq. (30a) is obtained from Eq. (27) and Eq. (3.1). Eq. (30b) results from the incompressibility condition (11b). Eq. (30c) is obtained from the kinetic condition (12d). Eq. (30d) is obtained from (29).

The above system (30a)-(30d) is an extension of the weak formulation introduced in Ref. [3] for two-phase flows and Ref. [43] for moving contact lines. In the current problem for EWoD, we have the additional term for the electric force in (30a). The electric force is obtained from solving (17)-(20).

We solve equations (30a)-(30d) for $\mathbf{u}$ and $p$ on the fluid domain $\Omega$ and for $\mathbf{X}$ and $\kappa$ on the reference domain $\mathbb{I}$ using the finite element method on unfitted meshes. To that end, we uniformly partition the time domain as $[0,T]=\cup_{m=1}^{M}[t_{m-1},t_{m}]$ with $t_{m}=m\tau$, $\tau=T/M$, and the reference domain as $\mathbb{I}=\cup_{j=1}^{J_{\Gamma}}\mathbb{I}_{j}$ with $\mathbb{I}_{j}=[\alpha_{j-1},\alpha_{j}]$, $\alpha_{j}=jh$ and $h=1/J_{\Gamma}$. We approximate the function spaces $\mathbb{K}$ and $\mathbb{X}$ by the finite element spaces

where $\mathcal{P}_{1}(\mathbb{I}_{j})$ denotes the space of the polynomials of degree at most 1 over $\mathbb{I}_{j}$. Denote by $\Gamma^{m}:=\mathbf{X}^{m}(\cdot)\in\mathbb{X}^{h}$ the numerical approximation to the fluid interface $\Gamma(t)$ at $t=t_{m}$. Then $\Gamma^{m}\ (0\leq m\leq M)$ are polygonal curves consisting of ordered line segments. The unit tangential vector $\mathbf{t}^{m}$ and normal vector $\mathbf{n}^{m}$ of $\Gamma^{m}$ are constant vectors on each interval $\mathbb{I}_{j}$ with possible jumps at the nodes $\alpha_{j}$, and they can be easily computed as

where $\mathbf{h}_{j}^{m}:=\mathbf{X}^{m}(\alpha_{j})-\mathbf{X}^{m}(\alpha_{j-1})$, and $(\cdot)^{\perp}$ denotes the counter-clockwise rotation by $90$ degrees.

Let $\mathcal{T}^{h}=\left\{\bar{o}_{j}\right\}_{j=1}^{N}$ be a triangulation of the bounded domain $\Omega$, and

where $k\in\mathbb{N}^{+}$, and $\mathcal{P}_{k}(o_{j})$ denotes the space of polynomials of degree $k$ on $o_{j}$. Let $\mathbb{U}^{h}$ and $\mathbb{P}^{h}$ denote the finite element spaces for the numerical solutions for the velocity and pressure, respectively. In this work, we choose them as

which satisfies the inf-sup stability condition [1, 43]

where $\lVert\cdot\rVert_{0}$ and $\lVert\cdot\rVert_{1}$ denote the $L^{2}$ and $H^{1}$-norm on $\Omega$ respectively, and $c$ is a constant.

In the simulation, the partition of the reference domain $\mathbb{I}$ for the interface and the triangulation $\mathcal{T}^{h}$ of the fluid domain $\Omega$ are both fixed in time. As a result, the triangulation of $\Omega$ is decoupled from the discretization of the interface $\Gamma^{m}$, thus may not fit the interface, i.e. the line segments comprising of $\Gamma^{m}$ are in general not the boundaries of the elements in $\mathcal{T}^{h}$. At $t=t^{m}$, the interface $\Gamma^{m}$ divides $\Omega$ into two sub-domains, $\Omega_{1}^{m}$ and $\Omega_{2}^{m}$. Correspondingly, we split $\mathcal{T}^{h}$ into three disjoint subsets: $\mathcal{T}_{1}^{m}$ being the set of elements in $\Omega_{1}^{m}$, $\mathcal{T}_{2}^{m}$ being the set of elements in $\Omega_{2}^{m}$, and $\mathcal{T}_{f}^{m}$ being the set of elements that intersect with the interface (see Fig. 2 for an illustration). The splitting can be easily done by using a recursive algorithm as follows:

• (1)

  First form $\mathcal{T}^{m}_{f}$ by locating all elements that intersect with $\Gamma^{m}$.

• (2)

  Locate one element $o_{j*}$ in $\Omega_{1}^{m}$, for example, the one that contains the point with coordinate $(\frac{1}{2}(x_{l}^{m}+x_{r}^{m}),~{}0)$, and set $\mathcal{T}_{1}^{m}=\left\{o_{j*}\right\}$.

• (3)

  Check all neighbours of the elements in $\mathcal{T}_{1}^{m}$. If a neighbor is not in $\mathcal{T}^{m}_{f}$, then add it to $\mathcal{T}_{1}^{m}$.

The last step is repeated until no element can be added to $\mathcal{T}_{1}^{m}$. This gives the set $\mathcal{T}_{1}^{m}$. The rest of the elements not belonging to $\mathcal{T}_{1}^{m}\cup\mathcal{T}_{f}^{m}$ form the set $\mathcal{T}_{2}^{m}$.

Denote by $\rho^{m}$, $\mu^{m}$ and $\beta^{m}$ the numerical approximations of the density $\rho(\cdot,t)$, the viscosity $\mu(\cdot,t)$ and the friction coefficient $\beta(\cdot,t)$ at $t=t_{m}$, respectively. We define $\rho^{m}$ and $\mu^{m}$ as

where $1\leq j\leq N$, $0\leq m\leq M$. Similarly, we denote by $\Gamma_{1}^{m}$ and $\Gamma_{2}^{m}$ the boundary of $\Omega_{1}^{m}$ and $\Omega_{2}^{m}$ at the lower wall respectively, and define $\beta^{m}$ at the lower wall as

where $\partial o_{j}$ is the boundary of $o_{j}$ at the lower wall, and $x_{l}^{m}=X^{m}|_{\alpha=0}$ and $x_{r}^{m}=X^{m}|_{\alpha=1}$ are the two contact line points.

The finite element method is given as follows. First we partition the time domain $[0,T]$, the reference domain $\mathbb{I}$ for the interface and the fluid domain $\Omega$ as described above. Let $\Gamma^{0}:=\mathbf{X}^{0}(\cdot)\in\mathbb{X}^{h}$ and $\mathbf{u}^{0}\in\mathbb{U}^{h}$ be the initial interface and velocity field, respectively. For $m\geq 0$, we compute $\mathbf{u}^{m+1}\in\mathbb{U}^{h}$, $p^{m+1}\in\mathbb{P}^{h}$, $\mathbf{X}^{m+1}\in\mathbb{X}^{h}$, and $\kappa^{m+1}\in\mathbb{K}^{h}$ by solving the linear system

Here, $\boldsymbol{g}^{h}=(g_{1}^{h},~{}g_{2}^{h})^{T}$, $\omega_{s}^{h}=\omega^{h}\cdot\mathbf{t}_{w}$, $u_{s}^{m+1}=\mathbf{u}^{m+1}\cdot\mathbf{t}_{w}$, and $x_{l}^{m}=\mathbf{X}^{m}|_{\alpha=0}$, $x_{r}^{m}=\mathbf{X}^{m}|_{\alpha=1}$. The electric force $\epsilon_{2}\eta|\nabla\Phi|^{2}$ in (37) is computed by solving equations (17)-(20) in the domain $\Omega_{2}^{m}\cup\Omega_{3}$; its computation will be presented in section 3.3. In (37d), the derivative $\partial_{s}\mathbf{X}^{m+1}$ is taken with respect to the arc length of $\Gamma^{m}$: $\partial_{s}\mathbf{X}^{m+1}=\frac{1}{|\partial_{\alpha}\mathbf{X}^{m}|}\partial_{\alpha}\mathbf{X}^{m+1}$, and similarly for $\partial_{s}\boldsymbol{g}^{h}$. At the first time step, we set $\rho^{-1}=\rho^{0}$.