ECMath OT6: Optimization and Control of Electrowetting on Dielectric for Digital Microfluidics in Emerging Technologies

This research is carried out in the framework of MATHEON supported by Einstein Foundation Berlin.

 Project head Prof. Dr. Michael Hintermüller (1,2) Staff Dr. Soheil Hajian (1) Project period 1 June 2017 - 31 December 2018 Affiliations (1) Humboldt-Universität zu Berlin (2) Weierstrass Institute Berlin Cooperations (ECmath) C-SE1, OT8 External cooperations Prof. T. M. Surowiec (U Marburg) Prof. S. W. Walker (Louisiana State Univ.) Prof. H. Antil (GMU) Further cooperations: K. Ito (North Carolina, Raleigh), R. Nochetto (Maryland, College Park), C. Elliott (Warwick), and J.A. Sethian (Berkeley).

## Introduction

A number of emerging key technologies in microbiology, medical diagnostic devices, personal genomics, as well as next-generation low-energy OLED displays and liquid lenses make use of a phenomenon known as electrowetting on dielectric (EWOD); see, e.g., [HZK${}^+$09, WMK${}^+$04, CMK03, GFK04]. Electrowetting involves the manipulation of small (microscopic) droplets on a dielectric surface by the actuation of the underlying current. In fact, droplets in a typical EWOD device are situated between two hydrophobic surfaces, one of which contains an array of controllable electrodes. The air-liquid-solid contact angle can then by changed by varying the voltages on separate electrodes, which causes the droplets to move. Thus, the voltages are a natural choice for influencing (controlling) the motion of a droplet, cf. figure below.

Building on very recent mathematical models for the simulation of EWOD, the proposed project will provide a mathematical, algorithmic, and computational framework for the optimization and control of such processes.

Figure: A (sub)optimal control pattern for barycenter matching (of a droplet) using a regularized sharp interface model in finite-horizon model predictive control; see [AHN$^{+}$].

## Mathematical models for EWOD

The EWOD process has been successfully modelled by the following time discretization of the Hele--Shaw-based (sharp interface) model with contact line pinning in [WSN09]: \begin{align}\textstyle \alpha \frac{\mathbf{u}^{i+1} - \mathbf{u}^i}{\delta t_{i+1}} + \beta \mathbf{u}^{i+1} + \nabla p^{i+1}= 0\qquad\text{in $\Omega^i$},\qquad\mathrm{div}\ \mathbf{u}^{i+1}&= 0&\text{in $\Omega^i$},&&\text{(1a)}\label{ewod_strong_1}\\ p^{i+1} - \kappa^{i+1} - E^{i} - \lambda^{i+1} - D_{\mathrm{visc}}(\mathbf{u}^{i+1}\cdot n^i)&= 0&\text{on $\Gamma^i$},&&\text{(1b)}\label{ewod_strong_3}\\ \lambda^{i+1} - P_{\rm pin}\partial(\|\cdot\|_{L^1(\Gamma^i)})(\mathbf{u}^{i+1}\cdot n^i) &\ni 0&\text{on $\Gamma^i$}. &&\text{(1c)}\label{ewod_strong_4} \end{align} Here, $\partial$ denotes the (convex) sub-differential, $\Omega^i \subset \mathbb R^2$ represents the (smooth) droplet domain at time $t_{i}$ with $\Gamma^i := \partial \Omega^i$ its boundary, and $\delta t_i > 0$, for $i=1,\dots,T$. Moreover, $0<\alpha\ll\beta$ and $D_{\mathrm{visc}}$, $P_{\mathrm{pin}} > 0$ are device parameters and we assume that $\mathbf{u}^{i}$ is represented on $\Omega^i$. The pressure condition (1b) models the effect of the mean curvature $\kappa^{i+1}$, the forcing term $E^{i}$, the pinning represented by $\lambda^{i+1}$, and the viscosity $D_{\rm visc}(\mathbf{u}^{i+1}\cdot n^i)$, respectively, on $p^{i+1}$ at time $t_{i+1}$. Equation (1c) is a macroscopic approximation of contact line pinning, a phenomenon due to molecular adhesion at the solid-liquid-air interface and contact angle hysteresis. Approximating the curvature term, we can write this system as \begin{align}\label{eq:ewod_vi} \text{Find } \mathbf{u}^{i+1}\!\in\! \mathbb V^{i}_{\rm sol}\!: \int_{\Gamma^{i}}|C^i\mathbf{z}|\! -\! |C^i\mathbf{u}^{i+1}| ds \ge \langle F^{i}(E^{i}) - A^{i} \mathbf{u}^{i+1}, \mathbf{z} -\mathbf{u}^{i+1} \rangle,\,\,\forall \mathbf{z} \in \mathbb V^{i}_{\rm sol},&&\text{(2)} \end{align} a VI of the second kind. Here, $\mathbb V^{i}_{\rm sol}$ is an appropriate solenoidal state space, the bounded linear operator $C^i:\mathbb V^{i}_{\rm sol}\to L^2(\Gamma^i)$ maps $\mathbf z$ to $C^i\mathbf z=\mathbf z\cdot n^i|_{\Gamma^i}$, $F^{i} :\mathcal{E}^{i} \to (\mathbb V^i_{\rm sol})^*$ is a continuous affine linear operator, and $A^{i}:\mathbb V^i_{\rm sol} \to (\mathbb V^i_{\rm sol})^*$ is a coercive and symmetric bounded linear operator.

This time-discrete sharp interface model provides an accurate description of the droplet's location and the associated pressure boundary condition. On the other hand, the phase field model [NSW14] couples a Cahn-Hilliard-Navier-Stokes system to a linear elliptic PDE for the electrostatic behavior and a parabolic equation for the electric charges. It is thermo-dynamically consistent and allows for topological changes such as splitting of droplets, but at the expense of not knowing the precise location of the droplet.

## Optimal control of EWOD

The project strives to provide a framework which allows engineers to design and test the viability of a device for subsequent fabrication, and to automate pre-existing devices. We hereby pursue both sharp interface [WSN09] and phase field models [NSW14].

The associated optimal control problems over all time intervals would determine solutions that minimize a respective overall objective $\mathcal{J}$. However, due to the non-trivial dependencies on the moving interface variables in the sharp interface context, the proof of existence of an optimal control remains ellusive without further restrictive assumptions or constraints, e.g., on the geometry. Moreover, the complexity of the phase field model poses severe challenges for a fast (real-time) numerical solution as needed for EWOD devices. In order to retrieve a viable optimal control setting, we thus remove a certain amount of difficulty and employ a model predictive control (MPC) approach.

#### Model predictive control for EWOD

In [AHN17], the following single step MPC problem is treated: $$\label{eq:oc}\ \text{Minimize }\mathcal J^i\!(\mathbf u^{i+1}\!,E^i)\text{ over }(\mathbf u^{i+1}\!,E^i)\!\in \mathbb V^i_{\rm sol}\times\mathcal{E}^i\text{ subject to } (\mathbf u^{i+1}\!,E^i)\text{ solves (2)},$$ where $\mathcal{E}^i$ is the control space and the objective $\mathcal J^i$ is an "instantaneous" version of the functional $\mathcal J$. Here, dual stationarity conditions are derived by applying a regularization procedure, exploiting techniques from PDE-constrained optimization, and then passing to the limit in the regularization parameters. Moreover, a function-space-based numerical procedure is developed by following the theoretical limit argument.

In order to improve optimality properties of the computed controls, we aim to develop algorithms based on multiple time step MPC: In a time instance $i$ and for a prediction horizon defined by an integer $M\geq 1$, we solve the optimal control problem restricted to the time interval $[t_{i+1}, t_{i+M}]$. The predicted control at $t_{i+1}$ is then employed to steer the system for one time step. This procedure is repeated with $i$ increased by one.

#### Directional differentiability of VIs of the second kind

While the method in [AHN17] represents the first algorithm for optimal control problems subject to EWOD processes, it typically can only provide some type of limiting C-stationary point in a time-discrete setting due to the applied smoothing technique. In order to guarantee more selective stationarity conditions, in [HS16] a minimization scheme is introduced for control problems subject to VIs of the first kind, which exploits the characterization of conical derivatives associated with the control-to-state operator.

Aiming to develop a similar algorithm in th context of EWOD, differentiability for VIs of the second kind was investigated in [HS17]. In fact, under regularity assumptions, the directional derivative $\mathbf d=S'(f^i;h)$ of the solution operator $S$ to (2) at $f^i = F^i(E^i)\in(\mathbb V_{\rm sol}^{i})^*$ in the direction $h\in(\mathbb V_{\rm sol}^{i})^*$ is given as the solution to the following VI of the first kind: $$\label{eq:vari_ineq_dir_der} \text{ Find } \mathbf d \in \mathcal{K}(\mathbf u^{i+1},\lambda^{i+1}): \langle A^i\mathbf d-h,\mathbf z-\mathbf d\rangle \ge 0,\,\,\forall \mathbf z \in \mathcal{K}(\mathbf u^{i+1},\lambda^{i+1}),$$ where $\mathcal{K}(\mathbf u^{i+1},\lambda^{i+1})$ is a suitable cone that is useful for numerical implementation.

#### Mesh adaptivity and consistent time discretization

Adaptive finite element methods have been successfully applied to the optimal control of VIs via the goal-oriented dual-weighted residual technique. An a posteriori error estimator is developed in the context of two-phase fluid flow (which is similar to the one in [NSW14]) in [HHK$^+$18].

The question of consistent time discretization of dynamical systems pertinent to optimal control problems and the respective adjoint equations was addressed in [HHU17]. Here, total variation diminishing Runge--Kutta (TVD-RK) schemes, which guarantee stability, best possible order and convergence of the discrete adjoint to its continuous counterpart have been constructed. In fact, it was proven that a state scheme which is a strongly stable TVD-RK method is always stable for the discrete adjoint equation. On the other hand, requiring strong stability for both the discrete state and adjoint is too strong, confining one to a first-order method, regardless of the number of stages employed in the TVD-RK scheme.

These techniques can be applied to improve the performance of numerical computations in the context of optimal control of EWOD.

## Publications

 [AHN$^{+}$] H. Antil, M. Hintermüller, R.H. Nochetto, T.M. Surowiec, and D. Wegner. Finite horizon model predictive control of electrowetting on dielectric with pinning. Journal: Interfaces and Free Boundaries, 19 (2017), 1-30. [HHU17] S. Hajian, M. Hintermüller, and S. Ulbrich. Total variation diminishing schemes in optimal control of scalar conservation laws. IMA Journal of Numerical Analysis, page drx073, 2017. [HHK$^+$17] M. Hintermüller, M. Hinze, C. Kahle, and T. Keil. A goal-oriented dual-weighted adaptive finite element approach for the optimal control of a nonsmooth Cahn--Hilliard--Navier--Stokes system. Optimization and Engineering (2018). [HS17] M. Hintermüller and T. M. Surowiec. On the directional differentiability of the solution mapping for a class of variational inequalities of the second kind. Set-Valued and Variational Analysis (2017)

## References

 [CMK03] S. K. Cho, H. Moon, and C.-J. Kim. Creating, transporting, cutting, and merging liquid droplets by electrowetting-based actuation for digital microfluidic circuits. Microelectromechanical Systems, Journal of Microelectromechanical Systems, 12(1):70-80, Feb 2003. [HS16] M. Hintermüller and T. Surowiec. A bundle-free implicit programming approach for a class of elliptic mpecs in function space. Math. Prog., 160(1):271-305, 2016. [HZK$^{+}$09] J. Heikenfeld, K. Zhou, E. Kreit, B. Raj, S. Yang, B. Sun, A. Milarcik, L. Clapp, and R. Schwartz. Electrofluidic displays using Young-Laplace transposition of brilliant pigment dispersions. Nat Photon, 3(5):292-296, 05 2009. [GFK04] J. Gong, S.-K. Fan, and C.J. Kim. Portable digital microfluidics platform with active but disposable lab-on-chip. In Micro Electro Mechanical Systems, 2004. 17th IEEE International Conference on. (MEMS), pages 355-358, 2004. [NSW14] R.H. Nochetto, A.J. Salgado, and S.W. Walker. A diffuse interface model for electrowetting with moving contact lines. Math. Models & Meth. in Appl. Sci., 24(1):67-111, 2014. [WMK$^{+}$04] A.R. Wheeler, H. Moon, C.-J. Kim, J. A. Loo, and R. L. Garrell. Electrowetting- based microfluidics for analysis of peptides and proteins by matrix-assisted laser desorption/ionization mass spectrometry. Anal. Chem., 76(16):4833-4838, 08 2004. [WSN09] S. W. Walker, B. Shapiro, and R. H. Nochetto. Electrowetting with contact line pinning: Computational modeling and comparisons with experiments. Physics of Fluids (1994-present), 21(10), 2009.