We present a highly efficient simulation method for the calculation of three-dimensional quasi-static interface shapes under the influence of electric fields. The method is especially useful for the simulation of microfluidic systems driven by electrowetting on dielectrics because it accounts automatically and inherently for the highly non-trivial interface shape in the vicinity of the triple-phase contact. In particular, the voltage independence of the local contact angle predicted based on analytical considerations is correctly reproduced in all our simulations. For the calculation of the shape of the interface, the geometry is triangulated and the mesh nodes are shifted until the system energy becomes minimal. The same mesh is also used to calculate the electric field using the boundary-element method. Therefore, only the surface of the geometry needs to be meshed, and no volume meshes are involved. The method can be used for the simulation of closed systems with a constant volume (e.g., droplet-based microfluidics) while preserving the volume very precisely as well as open systems (e.g., the liquid–air interface within micro-cavities or capillaries). Additional effects, such as the influence of gravitational forces, can easily be taken into account. In contrast to other efficient simulations, such as the volume-of-fluid, level-set, or phase-field methods, ideally, sharp interfaces are obtained. We calculate interface shapes for exemplary systems and compare with analytical as well as experimental results.

## I. INTRODUCTION

Electrowetting on dielectrics (EWOD) is a widely used technique for manipulating liquids and other soft matter at the microscale.^{1–9} The electrowetting effect, in general, describes the contact-angle change upon applying an electrical voltage.^{10–12} Displays,^{13} tunable lenses,^{14–16} or optical switches^{17,18} based on the EWOD effect have been reported in the literature. Additionally, the EWOD effect can be used to improve droplet coalescence^{19} or for micropumps^{20,21} due to the advantageous scaling of the surface forces. Very generally speaking, for smaller and smaller systems, EWOD-based concepts become more and more efficient, whereas it gets harder, and harder to realize, e.g., conventional pump designs on a microscopic level. For this reason, the EWOD effect forms a good basis for building efficient microfluidic systems. Unfortunately, the numerical simulation of EWOD-driven systems is very challenging due to the coupling of fluid and electrodynamics. This is especially the case for complicated, e.g., curved electrodes,^{22} partially conductive liquids or if the exact interface shape near the three-phase contact line is of interest. For this reason, a general purpose and at the same time efficient simulation method is desirable. However, often, it suffices to describe equilibrium or quasi-equilibrium states of the system, which can be simulated with a much reduced effort. A novel, highly efficient simulation method for such systems is presented in this work, which even allows for on-the-fly simulations using inexpensive hardware. The present paper starts with a short summary of EWOD basics. Next, the equations necessary for the description of closed systems with a constant volume and of open systems, where the volume is not preserved, are derived. After that, the solution method for the obtained equations is introduced and explained. Finally, the introduced simulation method is tested on examples, and the numerical results are compared with analytical predictions.

## II. FUNDAMENTALS

The reader may wish to think of a sessile droplet on top of a planar solid substrate (s) as a paradigmatic EWOD system. Aqueous salt solutions or ionic liquids are commonly used as electrically conducting liquids (l). Another immiscible non-conducting liquid, air, or vapor (v) forms the ambient medium. The equilibrium contact angle $\u03d1Y$ at the triple-phase contact line is given by Young’s equation,

where $\sigma lv$, $\sigma sv$, and $\sigma sl$ denote the respective interfacial energies. The change of the effective interfacial energy $\sigma sleff$ as a function of the applied voltage $U$ is described by the Lippmann equation,^{12}

where $\rho sl$ denotes the surface charge density at the interface. In the classical EWOD configuration, a dielectric layer insulates the droplet from the counterelectrode. Integration of Eq. (2) under the assumption that the surface charge density is given by the plate-capacitor formula $\rho sl=\epsilon 0\epsilon r,ddDU$ leads to

where $dD$ denotes the thickness of the dielectric layer, $\epsilon 0$ the vacuum permittivity, and $\epsilon r,d$ the relative permittivity of the dielectric layer. Apparently, this equation summarizes very intuitively the capacitor picture of EWOD. However, this apparent simplicity can be deceptive: One could argue incorrectly that Eq. (3) implies contact-angle saturation based on the argument that surface tensions cannot become negative.^{11,23–25} However, the effective interfacial energy is ultimately only the ratio of energy per area, and for an arbitrarily large voltage, the reduction in electrical energy is, of course, arbitrarily large.^{11} Then, the effective interfacial energy $\sigma sleff$ defined in Eq. (3) reaches arbitrarily large negative values as long as other effects do not come into play that are not accounted for in this simplified electrochemical description; e.g., trivially, the picture of an arbitrarily deformable, conductive liquid breaks down at the latest when the liquid spreads out to molecularly thin layer thickness.

Insertion of Eq. (3) into Eq. (1) and the definition of the dimensionless electrowetting number $\eta \u2261\epsilon 0\epsilon r,d2dD\sigma lvU2$ yields the well-known Young–Lippmann equation,

In a remarkable paper and in contrast to “common knowledge,” Buehrle *et al*. predicted that the Young–Lippmann equation is not valid near the triple-phase contact line:^{26} Contrary to the prediction of Eq. (4), applying a voltage does not change the local contact angle directly at the triple-phase contact line. They expressed the electric field force at a small distance $r$ away from the triple-phase contact line in terms of the Maxwell stress $\Pi el\u221dE2$ acting on the surface. They showed that even though $\Pi el$ is proportional to $\delta r\nu $ with $\u22121<\nu <0$ and, thus, diverges directly in the triple-phase contact line, the integral

vanishes and so does the electrical force directly at the triple-phase contact line. Thus, the local contact angle remains constant and is independent of the applied electrical voltage. While the local contact angle is difficult to observe and rather irrelevant in macroscopic systems, we shall see below that the contact region cannot be neglected on the micrometer scale and that, thus, the correct prediction of the local wetting properties is crucial for the description and design of many EWOD-driven micro- or nanofluidic systems. This effect is automatically taken into account in the presented method; see below in Secs. IV and V, in particular, Fig. 7 and Sec. V B.

## III. FROM PDEs TO ENERGY FUNCTIONALS

The static shape of interfaces between two immiscible liquids can be determined by solving the Young–Laplace equation.^{27} For EWOD-driven systems, this is

Here, $\u2207n$ denotes the mean curvature, where $n$ is the normal vector on the interface, $phydro$ is the hydrostatic pressure due to the influence of gravity, and $pel$ an additional pressure term caused by an external electric field. Accordingly, the exact form of the equation depends on the equation used for the calculation of the mean curvature. For the special case that the interface can be written as a function $z=f(x,y)$, the following holds for the curvature:

where $fx$ and $fy$ denote the partial derivatives with respect to $x$ and $y$, respectively. Accordingly, the Young–Laplace equation is a highly nonlinear partial differential equation. For this reason, the number of analytical solutions is severely limited and known only for a few two-dimensional and rotational symmetric problems and only if electric fields are neglected. However, approximation formulas have been found for some special geometries, although a considerable amount of mathematics is required.^{28–30} Suitable boundary conditions must be specified for the solution of the Young–Laplace equation. These boundary conditions are given by the contact angles at the triple-phase contact line and can be defined via

where $t$ is the unit tangent to the interface and $b$ is a unit vector that is binormal to a vector tangential to the contact line and a normal vector of the solid surface. Thus, this boundary condition is equivalent to Young’s equation [see Eq. (1)]. When Eq. (6) is solved directly, taking Eq. (8) into account, several challenges must already be overcome even when no external electric fields are present:

- $\u2219$
The boundary on which the boundary condition (8) must apply is not known at first (free boundary problem

^{31}). Direct numerical solution methods as shown, for example, by Cooray*et al.*^{32}are, therefore, based on iterative methods. - $\u2219$
Mostly, one is interested in solutions, which fulfill additional constraints (e.g., a constant volume). Such additional constraints require again iterative solution methods, as shown, e.g., in Sec. IV.

- $\u2219$
It is a highly nonlinear partial differential equation. This severely limits the choice of available solvers.

For these reasons, the Young–Laplace equation is usually not solved directly.^{29,33,34} A common method of dealing with general three-dimensional problems involving the calculation of interface shapes is, therefore, to treat the problem as a time-dependent multiphase flow problem and to simulate the system over a sufficiently long period of time until the stationary state has been reached. In this case, volume-of-fluid, phase-field, or level-set methods can be used, for example. However, with these methods, the volume conservation is often problematic. Furthermore, within some formalism, the interface is not ideally sharp. Finally, time-resolved simulations are unnecessarily computationally intensive if one is only interested in the equilibrium state of the system even if efficient methods exist for the time-resolved simulations as described, for example, by Grawitter and Stark.^{35} Another group of methods is based on the use of phase-field methods to describe the interfacial motion but model the EWOD effect solely via the Young–Lippmann equation^{36–39} and do not solve for the electrical field directly. A worth mentioning method and particularly efficient method for the computation of EWOD-actuated droplets on structured surfaces was introduced by Chamakos *et al.*^{40–43} However, this method is difficult to generalize to general electrode configurations and arbitrarily shaped fluid regions. To overcome these difficulties, a new and very efficient simulation method for modeling EWOD systems based on an energy minimization method is used here to solve Eq. (6).

Both the Young–Laplace equation and the Young equation can be obtained from purely thermodynamic considerations by searching for stationary points of a suitable thermodynamic potential.^{44} A challenge with this method is that the exact form of the energy functional depends on the system under consideration and the assumptions made about the behavior of the involved fluids. For the description of EWOD-driven systems, the free energy is particularly suitable,^{10}

To write out the terms explicitly, a system of two immiscible fluids, which are in contact with a solid surface, is considered and gravity is supposed to act in the z-direction. The $Finterfaces$ term contains all contributions caused by the interfacial energies and can be written as

where $Alv$ denotes the interfacial area between the fluids, $Asl$ describes the interfacial area of the wetted walls, and $Asv$ describes the interfacial area between the solid and the vapor. The second term $Fvol$ contains the contributions of gravitational energy and, if needed, other volume energies as well as possibly Lagrange multipliers $FL$ to constraint the volume,^{45}

Here, $\Delta \u03f1=\u03f1F\u2212\u03f1V$ is the difference in density between the fluids, $V$ is the volume of the fluid, and $g$ is the gravitational acceleration. The last term $Fel$ of Eq. (9) contains the energy contributions, which are caused by electrical fields. In order to write down this term, the correct thermodynamic potential must first be chosen. For electrical systems, there are basically two possibilities for the constraints: a constant charge or a constant voltage. The resulting thermodynamic potentials are connected via a Legendre transformation. In general, the voltage is kept constant in EWOD systems. In this case, it follows:

where $E$ denotes the electrical field and $D=\epsilon 0\epsilon rE$ the electrical displacement. The second term $\Delta \Pi $ corresponds to the osmotic pressure caused by an increased concentration of ions at the interface. Whether this pressure must be taken into account depends on the properties of the liquid under consideration. Mostly, EWOD experiments are performed with aqueous salt solutions with high ion concentrations. In these cases, the liquids can be modeled as ideal conductors and the interior of the liquid can be treated as free of electrical fields. In reality, however, the fields are screened by an accumulation of ions at the interface. In this case, the electric potential depends on the ion distribution in addition to the boundary conditions. The potential can be obtained as a solution of the nonlinear Poisson–Boltzmann equation,

where the sum runs over all species of ions where each species has a bulk concentration $nib$ and a charge $qi$. Under this assumption, it follows for the osmotic pressure,^{46}

The depth of penetration of the electric field into the electrolyte is given by the Debye length $\lambda =\epsilon 0\epsilon rkBT/(\u2211inibqi2)$. For example, the Debye length for a 1:1 electrolyte, such as an NaCl-solution with a concentration of $0.1moll\u22121$, is less than 1 nm at room temperature and is much smaller than 1 $\mu $m even for pure water at room temperature due to self-ionization. Thus, the influence of the finite conductivity can be neglected in most cases. Equation (13) complicates the problem of calculating the interface shape considerably since the equation adds an additional non-linearity. However, the problem can be simplified considerably if $q\phi kBT\u226a1$ is valid. In this case, the potential can be determined using the linear Poisson–Boltzmann equation,

where $\kappa =\lambda \u22121$ denotes the inverse Debye length. Then, the calculation of the osmotic pressure simplifies to

The presence of double layers leads to a smaller decrease in the contact angle compared to an ideally conducting liquid. For the sessile droplet in the conventional EWOD setup, the electrowetting number changes as follows for a linear monovalent salt:

A detailed derivation of this expression was given by Mugele and Buehrle.^{47}

## IV. SIMULATION METHOD

For the simulation, the surface of the given geometry is triangulated. An energy functional describing the particular system is introduced and minimized with respect to the position of the mesh vertices. The resulting derivative terms can be analytically obtained and efficiently evaluated for any triangulation. Energy contributions resulting from the interfacial tensions, the electric field, and gravity are included explicitly in the following. Other energy terms can easily be added. A similar approach has already been proposed by Buehrle *et al*.^{26} However, only a two-dimensional system was simulated in a small region around the triple-phase contact line. In the following, the idea is extended to three-dimensional systems considering arbitrary boundary conditions and an arbitrary number of phases. Closed systems with a constant volume as well as open systems, where both energy and matter can be exchanged, are discussed below. At the beginning, the energy terms are summarized, which are present in both cases. We assume that the fluids are incompressible. For ease of notation, it is assumed that the gravitational force acts in the $z$-direction. In the following, the derivation is performed for a system consisting of a single liquid in contact with a solid and gas phase. However, the equations can also be applied to systems containing several immiscible liquids. For simplification, it is assumed for the derivation of the equations that the wetting properties are the same on all wetted surfaces. However, the procedure can be applied analogously to systems with location-dependent wetting properties. Under these assumptions, for both closed and open systems, the following terms occur in the thermodynamic potential at a given external voltage:

It should be noted that the individual terms in this equation are not independent of each other. Thus, the energy minimization must be done simultaneously, taking all contributions into account. It should also be noted that all (physical) interfacial energies $\sigma $ in this expression are voltage independent.

The actual value of the energy is irrelevant for performing the minimization. Rather, the derivative of the energy with respect to the positions of the mesh vertices must be calculated. With $\u2207Asl=\u2212\u2207Asv$ and Eq. (1), it takes the form

Here, $\u2207f$ is the multi-dimensional gradient $(\u2202f/\u2202X1,x,\u2202f/\u2202X1,y,\u2202f/\u2202X1,z,\u2026,\u2202f/\u2202XN,z)T$, and $Xj,i$ is the $i$th component of the $j$th vertex. It describes the change of the respective quantities under displacement of each mesh vertex. The last term in Eq. (19) shows the strong coupling between the electric field and the deformation of the liquid: A deformation of the liquid leads to a change of the electric field and vice versa. However, the field can be considered constant for small deformations, and therefore, the last term can be approximated as follows:

Thus, the final equation for the energy gradient is

### A. Closed systems

For the simulation of closed systems with incompressible liquids, the volume should be preserved as well as possible during the simulation. To accomplish this, a Lagrange parameter $p$ corresponding to a pressure is introduced analogously to equilibrium thermodynamics, and the optimization is performed under the constraint $V=const.$. The free energy $FC$ and the volume $V$ change under variation of the position of the mesh vertices as follows:

If we perform an iterative algorithm for the displacement of the vertices with a search direction $\u2207E\u2212p(\u2207V)$ and step size $\alpha $ for the variation of the positions of the mesh vertices

the following volume change results:

In order for the volume change to vanish, the pressure $p$ must be chosen as follows:

Thus, the Lagrangian parameter $p$ at the equilibrium state equals the Laplacian pressure $pL$. This gives the final expression for the gradient of the free energy for closed systems,

### B. Open systems

For an open system, the volume is by definition not constant. The system considered is again a liquid in contact with a gas phase. It is assumed that the system has only a single inlet for the liquid with a given pressure $pin$. Assuming that a state of equilibrium nevertheless exists (e.g., for capillary filling or confined gas volumes), the gradient of the enthalpy fulfills,

The second term corresponds to the flow work performed on the liquid. The last term in this equation is an example of an additional energy term, which is specific for the respective system. This term takes into account that the liquid may do work on an enclosed gas (e.g., capillary rise inside a sealed capillary). The form of this term depends on the description of the behavior of the gas. If we assume for simplicity that the gases can be described as ideal gases, the following holds for $\delta W$:

Thus, the gradient $\u2207\delta W$ can be calculated as follows:

The combination of Eqs. (29), (27), and (21) finally results in the following expression for the gradient of the enthalpy:

which can be used for the optimization algorithm. The generalization of our approach to systems with multiple inlets and multiple free surfaces depends on the specific properties of the system under consideration but is straightforward.

### C. Solution method

We present our solution method using the expression of $\u2207FO$. However, the method for closed systems is identical except the equation for $\u2207F$. One can use Eq. (30) to establish an iterative solution method for the optimization problem based on the following iteration rule:

We state explicitly that the field term depends only on the index $i$ because the fields are not recalculated in each inner iteration $n$. Instead, there is an outer loop over $i$, in which the fields are recalculated for the actual geometry, and an inner loop over $n$ in which only the current geometric quantities $\u2207Alv,\u2207Asl$, $\u2207\u222bzdV$, $V$, and $\u2207V$ are recalculated and the positions of the mesh vertices are updated. Thus, a two-level solution approach is used. For the optimization, the points of the triangulation inside the surfaces of the geometry are shifted according to the following equation:

where $\alpha ni$ denotes the step size. In the following, the abbreviation $\u2207F(Xni)=\u2207Fni$ is used unless otherwise specified. In order to increase the convergence speed and make the process more stable, a line-search method is used where the parameter $\alpha ni$ has to fulfill the well-known strong Wolfe conditions:^{48}

where $c1$ and $c2$ are fixed constants, which satisfy the conditions $0<c1<c2<1$. Since we only need the gradient of enthalpy to perform the optimization procedure, it is more efficient to also perform the line-search procedure using only this quantity. For this purpose, the difference $Fn+1i\u2212Fni$ is approximated using a Newton–Cotes equation as follows:

Here, $wl$ denotes the integration weights and $tl$ are the sampling points. Since the calculation of the gradient $\u2207F(Xni\u2212tl\u2207Fni)$ is a computationally intensive operation, the integral must be approximated using as few sample points as possible. Numerical experiments showed that even a simple trapezoidal rule with $w1,2=0.5$ and $t1=0,t2=\alpha n$ provides a good enough approximation of the integral for the line-search procedure. This results in the following two inequalities, which are used to determine the step size $\alpha ni$:

Equation (36) prevents the step from being chosen too large. In contrast, Eq. (37) excludes step sizes that are too small. To use the line-search algorithm, only the constants $c1$ and $c2$ need to be specified. Our experience shows that good results are achieved for the values $c1=1\xd710\u22124$ and $c2=0.9$, which are suggested by Nocedal and Wright.^{49} To speed up the calculation, the last value found for $\alpha ni$ is used as the starting value in the next iteration. If the equilibrium form is determined in the absence of an electric field, the convergence can be significantly improved with the help of the line-search algorithm presented. However, special care must be taken if an electric field is considered. Since the electric field is determined only once in each outer iteration $i$ and the approximation made in Eq. (20) is valid only for small deformations, the acceptable step size is usually overestimated. For this reason, the number of inner iterations is reduced depending on the outer iteration number as discussed below.

### D. Mesh handling

The total volume $V$ as well as the occurring gradient terms $\u2207V$, $\u2207A$, and $\u2207\u222bzdV$ can be calculated analytically for a triangular mesh. The equations are summarized in Table I. The meaning of the used geometric parameters is illustrated in Fig. 1.

Symbols . | Mesh quantities . |
---|---|

N | Total number of triangles |

M | Total number of vertices |

t_{j} | Number of neighboring triangles for the jth vertex |

X_{j} | Vector pointing to the jth vertex of the triangulation with j ∈ [1, …, M] |

$aj\mu $ | μth neighboring vertex of the jth vertex with μ ∈ [1, …, t_{j}] |

p^{ν} | Midpoint of the νth triangle with ν ∈ [1, …, N] |

m^{ν} | Normal vectors of the νth triangle with ν ∈ [1, …, N]; |

note: identical to the corresponding $nj\mu $ vectors of the vertices | |

$nj\mu =(aj\mu \u2212Xj)\xd7(ajfj(\mu )\u2212Xj)$ | Normal vector of the μth neighboring triangle of the jth vertex |

f_{j}(μ) = mod(μ, t_{j}) + 1 | The function f_{j}(μ) ensures that the summation is done cyclically over the neighbors of the jth vertex |

Geometric quantities | Equations |

Area | $A=12\u2211\nu =1N|m\nu |$ |

Volume | $V=16\u2211\nu =1Np\nu \xb0m\nu $ |

Surface area gradient | $(\u2207A)j=12\u2211\mu =1tjnj\mu |nj\mu |\xd7(ajfj(\mu )\u2212aj\mu )$ |

Volume gradient | $(\u2207V)j=16\u2211\mu =1tjnj\mu $ |

Gradient of gravitational energy | $(\u2207\u222bxkdV)j=148\u2211l=13el\u2211\mu =1tj\u2211q=13\epsilon klq(aj,q\mu \u2212aj,qfj(\mu ))[(aj,k\mu +aj,kfj(\mu ))2+(aj,k\mu +Xj,k)2+(Xj,k+aj,kfj(\mu ))2]+148\u2211\mu =1tjnj,k\mu [2aj,k\mu +2aj,kfj(\mu )+4Xj,k]\delta klek$ |

Symbols . | Mesh quantities . |
---|---|

N | Total number of triangles |

M | Total number of vertices |

t_{j} | Number of neighboring triangles for the jth vertex |

X_{j} | Vector pointing to the jth vertex of the triangulation with j ∈ [1, …, M] |

$aj\mu $ | μth neighboring vertex of the jth vertex with μ ∈ [1, …, t_{j}] |

p^{ν} | Midpoint of the νth triangle with ν ∈ [1, …, N] |

m^{ν} | Normal vectors of the νth triangle with ν ∈ [1, …, N]; |

note: identical to the corresponding $nj\mu $ vectors of the vertices | |

$nj\mu =(aj\mu \u2212Xj)\xd7(ajfj(\mu )\u2212Xj)$ | Normal vector of the μth neighboring triangle of the jth vertex |

f_{j}(μ) = mod(μ, t_{j}) + 1 | The function f_{j}(μ) ensures that the summation is done cyclically over the neighbors of the jth vertex |

Geometric quantities | Equations |

Area | $A=12\u2211\nu =1N|m\nu |$ |

Volume | $V=16\u2211\nu =1Np\nu \xb0m\nu $ |

Surface area gradient | $(\u2207A)j=12\u2211\mu =1tjnj\mu |nj\mu |\xd7(ajfj(\mu )\u2212aj\mu )$ |

Volume gradient | $(\u2207V)j=16\u2211\mu =1tjnj\mu $ |

Gradient of gravitational energy | $(\u2207\u222bxkdV)j=148\u2211l=13el\u2211\mu =1tj\u2211q=13\epsilon klq(aj,q\mu \u2212aj,qfj(\mu ))[(aj,k\mu +aj,kfj(\mu ))2+(aj,k\mu +Xj,k)2+(Xj,k+aj,kfj(\mu ))2]+148\u2211\mu =1tjnj,k\mu [2aj,k\mu +2aj,kfj(\mu )+4Xj,k]\delta klek$ |

The terms $\u2207V$, $\u2207A$, and $\u2207\u222bzdV$ have the dimension $3N\xd71$, where $N$ is the number of mesh vertices. The wetting properties of the liquid and the surface are completely determined by the contact angle $\u03d1Y$ and the surface tension $\sigma lv$. To ensure that the properties are correctly taken into account, the term $\u2207A$ is split into the two parts $\u2207Alv$ and $\u2207Asl$ for each vertex within each triple-phase contact line. Here, $\u2207Asl$ denotes the gradient of the wetted surface area, and $\u2207Alv$ is the gradient of the liquid–vapor interface. For all other vertices within the liquid–vapor interface, only the gradient $\u2207Alv$ has to be calculated. Since the connectivity of the mesh remains unchanged during the optimization, the corresponding lists with the indices and the number of the neighbor vertices $tj$ of vertex $j$ as well as the list with the neighbor edges only need to be determined once. The vectors $ajfj(\mu )\u2212aj\mu $ correspond to all edges of the mesh. Each of these edges can be assigned a unique index. The indices of the attached edges are stored for all vertices in a matrix of size $N\xd7K$, where $K=max(tj)$ corresponds to the maximum number of neighbor vertices. Since not every vertex has the same number of adjacent edges, empty entries remain in this matrix. These empty entries point to the zero vector. Using this indexing, all gradient terms can be calculated fully parallelized for the complete triangulation in a single step.

### E. Moving mesh algorithm

In each step of the optimization, the positions of the vertices in the liquid domain are changed, but the positions of the vertices in the adjacent solid and gas domains remain unchanged. Therefore, a large deformation would lead to an inversion of the mesh elements. To circumvent this limitation, a moving mesh method is used, whereby the vertices in the adjacent solid and gas domains are also shifted. Therefore, the vertices are connected by imaginary springs. Their stiffness is chosen according to the method of semi-torsion springs proposed by Blom^{50} as

where $j$ runs over all neighboring vertices of the vertex $Xi$. The angles $\beta ijn$ denote the $n$ interior angles at the vertices, which are not $Xi$ or $Xj$ of the triangles adjacent to the edge $Xi,j$. After calculating the stiffnesses, the following system of equations is solved for each component $\delta k,k=x,y,z$ of the displacement vector $\delta $:

where $cj$ is a list with the indices of all neighbor vertices, $mj$ is a list with all indices of freely movable neighbor vertices, and $hj$ is a list with the indices of all neighbor vertices with a given displacement for the $j$th vertex, respectively. Using an advanced meshing algorithm, an almost constant number of neighbor triangles can be achieved for all vertices. Thus, the calculation effort for the calculation of the gradient terms $\u2207V$ and $\u2207A$ increases linearly with the number of mesh vertices. To calculate a moving mesh step, a linear system of equations must be solved. However, the displacement must only be determined for the freely moving vertices. This minimizes the number of degrees of freedom in the moving mesh calculations. In addition, the system matrix for calculating the moving mesh step is sparse.

### F. Remeshing algorithm

In order to keep the mesh quality high even in the case of very large deformations, an additional remeshing process was implemented. In the first step, the vertices that lie inside of the triple-phase contact lines or within geometry edges are redistributed. For this purpose, the positions of the vertices inside such edges are uniformly distributed using a linear interpolation, where the positions of mesh nodes representing sharp corners remain unchanged. Subsequently, the vertices within the feature faces are moved. For this, an improved position of the selected mesh vertices is calculated using a uniform Laplace smoothing,

where $w$ are weights that satisfy the condition $\u2211\mu =1tjw\mu =1$. To keep the volume approximately constant and the geometry as identical as possible, the motion of the vertices is restricted. To do this, the positions of the mesh vertices $cj$ are back-projected onto their tangent plane by solving the equation

where $Tj$ is the vertex normal vector corresponding to the tangent plane of the $j$th vertex. These vectors are calculated by averaging the normal vectors of the adjacent triangles of each vertex $Xj$. Since this method is only applied to vertices outside geometry edges, the neighboring triangles are always contained within the same domain, and thus, the respective vertex normals are well defined.

### G. Electric field calculation

The remaining task in calculating the gradient $\u2207F$ is the calculation of the fields $E$ and $D$. In principle, all common methods, such as finite elements or finite differences, can be used for this calculation. However, for our specific problem, the boundary-element method is particularly suitable since only the surface of the geometry needs to be meshed, whereby the computing time can be significantly shortened. In addition, the algorithms for describing the movement of the mesh are simplified since only the mesh on the surface is deformed. It is assumed that the potential can be calculated using the linear Poisson–Boltzmann equation (15). Constant shape functions are used for the boundary-element simulation. This prevents singularities of the electrical field near the triple-phase contact line. Even though the potential and the charge distribution are calculated only on the boundary, the linearity of the underlying differential equation allows the solution to be calculated anywhere in the volume if desired. Thus, the boundary-element method does not introduce any additional restriction of generality. The discretized integral equation for each domain of the geometry can be formulated as follows:

where $\varphi nk$ denotes the potential value of the $n$th triangle of the $k$th domain and $Dnk$ describes the normal component of the displacement field of the $n$th triangle of the $k$th domain. $Un,ik$ and $Tn,i$ are the surface integrals of the fundamental solution and the normal derivative of the fundamental solution of the Poisson–Boltzmann equation, respectively. They can be calculated as follows:

where $rn$ describes the midpoint of the $n$th triangle and $ri$ points to all Gauss quadrature nodes on the $i$th triangle. Thus, each triangle from each domain contributes with two unknowns $\varphi nk$ and $Enk$. On an outer boundary, a Dirichlet or Neumann condition is specified and either $\varphi nk$ or $Enk$ is given. On the interfaces between two domains $I$ and $II$, however, the four unknowns $\varphi nI$, $\varphi nII$, $EnI$, and $EnII$ must be determined for each triangle. Therefore, additional continuity conditions must be fulfilled on all triangles within an interface between two adjacent domains. These can be formulated as follows:

where $\Lambda n$ denotes a potential user-defined surface charge density on the $n$th triangle. Otherwise, $\Lambda n$. The integrals in Eq. (43) shall be discussed next. They are evaluated using an efficient Gaussian quadrature rule,

where $M$ denotes the degree of the integration rule. The weights $Wm$ and integration points $\zeta m$ and $\eta m$ are calculated according to the GQUTM (Gaussian quadrature formula for unit triangles method).^{51} The integration is carried out over the canonical triangle ${0\u2264\eta \u22641,0\u2264\zeta \u22641\u2212\eta}$ using the linear mapping $X(\zeta m,\eta m)=(1\u2212\zeta m\u2212\eta m)Xi,1+\zeta mXi,2+\eta mXi,3$, where $Xi,n$ are the coordinates of the vertices of the $i$th triangle. If $i$ is equal to $n$, the integrands in Eq. (43) become singular and cannot be evaluated precisely using the Gaussian quadrature described before. For the case where $\kappa =0$ holds and the Poisson–Boltzmann equation becomes the Laplace equation, these integrals can be calculated analytically as follows:

Unfortunately, in the case of $\kappa >0$, the integrals $\Delta Un,nk$ can no longer be integrated in a closed form.^{52} However, using the same idea that was used to derive the equations (46), the integrals can be put into the form

while $\Delta Tn,n=0$ still holds. The remaining integral is no longer singular and can be solved efficiently using a standard Gaussian quadrature formula.

The used geometric quantities are shown in Fig. 2. To evaluate Eq. (20), the field at the vertices must be known. However, the solution of Eq. (42) provides the field only at the midpoints of the triangles. Therefore, the field is averaged over the neighbors of each vertices. A comparison of the potential distribution using the linear Poisson–Boltzmann theory and the assumption of an ideal conducting fluid are shown in Fig. 3. In the case of linear Poisson–Boltzmann theory, the potential within the droplet decreases continuously with increasing distance from the dielectric surface. Thus, the surface of the drop itself is no longer an isopotential surface. In comparison, in the case of an ideally conducting liquid, the droplet is completely field-free inside. It is worth pointing out that the presented method allows combining different fluid properties in a single simulation. For example, EWOD can be studied on an ideally conducting fluid surrounded by a fluid of finite conductivity. In addition, interfacial charges can be specified, e.g., to model the influence of the absorption of surfactants. Thus, the simulation method can be used to investigate and simulate a wide range of different fluidic systems.

The calculation of the electrical field is by far the computationally most intensive part of the optimization. The effort for setting up the system matrix grows quadratically with the number of triangles. However, the calculation of the entries of the system matrix can be done separately for all domains and, thus, in parallel. Since the system matrix is densely occupied, in general, $N3$ operations are needed for the solution of the linear system of equations. For an efficient calculation, the number of field calculations must, therefore, be kept as small as possible. However, the approximation made in Eq. (20) is only valid for small deformations. For this reason, a user-specified maximum deflection is introduced. If the total deflection of any vertex exceeds this limit, the field is recalculated for the new geometry and the optimization continues. If the maximum deflection is not exceeded, a predefined number of gradient steps is calculated. The number of gradient steps is reduced after each field calculation according to the rule $Ngradi+1=r\u22c5Ngrad0$ with $0<r\u22641$ and is reset to the initial value $Ngrad0$ as soon as a new voltage value is specified. Therefore, at the beginning of each voltage increment, when the surface shape is still far away from the equilibrium position, a high number of inner iterations is performed. As the optimization proceeds, the field calculation is carried out after a decreasing number of iterations. Thus, oscillations around an energy minimum can be avoided, and at the same time, the number of field calculations can be minimized. Usually, 10–20 field calculations are sufficient for the solution to converge to an acceptable equilibrium position. To avoid unnecessary iterations, the optimization is aborted for the current voltage value and continues for the subsequent value as soon as the length of the descent vector falls below a specified tolerance threshold $tolabs$. In addition, the inner iteration is aborted, and the electric field is calculated again as soon as the length of the descent vector has decreased by a relative value $tolrel$.

The described simulation method was implemented in MATLAB in an object-oriented way. In addition to the boundary conditions and the material parameters, only a surface mesh needs to be specified as input parameters. An arbitrary number of domains as well as mixed boundary conditions are supported. It should also be mentioned that the method described here can be used for the calculation of static interface shapes in a multiphase system with an arbitrary number of phases.

## V. NUMERICAL EXAMPLES

### A. Closed systems—EWOD of a sessile droplet

To test and validate the method, the system of a sessile droplet is investigated. This system has already been extensively studied theoretically and experimentally. Moreover, the droplet shape can be determined semi-analytically even under the influence of gravity and is, therefore, a good benchmark for comparison. For the semi-analytical calculation of the droplet shape, the nonlinear system of equations,

is integrated numerically for the initial conditions $x0=y0=s0=0$ until the angle $\phi $ equals the specified voltage-dependent contact angle $\u03d1Y$ calculated using Eq. (4). Here, $Bl=\sigma lv/(\Delta \u03f1g)$ denotes the capillary length, and $Rs$ denotes the principal curvature in the apex of the droplet. Then, $RS$ is varied until the droplet volume equals the specified target volume $Vt$ by minimizing the objective function $Z(RS)=(Vt\u2212\pi \u222b0ymaxx2dy)2$. The simulation geometry consists of a cube-shaped domain describing the solid domain and an hemispherical domain modeling the droplet. During the simulation, the droplet can deform freely as required. The surface of the simulation area is triangulated. We use 6574 triangles in the surface of the liquid domain and 7208 triangles in the surface of the solid domain. As a first test, the droplet shape in a gravitational field is investigated without the influence of an external electric field for different values of the difference in density $\Delta \u03f1$. The results are shown in Fig. 4. As expected, the drop flattens for an increasing density difference. It is also clear that the simulated solution agrees very well with the semi-analytical solution shown as red curves. The calculation time for each value $\Delta \u03f1$ is less than 3 min on a standard laptop, with only one core being used.

Next, the classical EWOD system of a sessile droplet under the influence of an external electric field was investigated. In this case, the liquid is considered to be ideally conductive, and thus, the inside of the liquid is free of an electrical field. The simulated shapes of the droplets are shown in Fig. 5 for different values of the applied electrical voltage. As expected, the drop spreads with increasing voltage. To test the volume conservation, the volume of the droplet was calculated in each iteration and then the volume deviation was determined, which is shown in Fig. 6(a). The maximum volume increase is less than $12pl$ for a droplet with an initial volume of $5\mu l$. Thus, the volume is very well preserved even for strong deformations. In order to test the accuracy of the simulation method, the dielectric intermediate layer was chosen to be particularly thick, with a thickness of $dD=500\mu m$, as this makes the difference between the apparent and local contact angles near the three-phase contact line especially visible. As shown by Buehrle *et al*., the prediction of the Young–Lippmann equation deviates from the actual contact angle near the triple-phase contact line.^{26} The local contact angle is voltage-independent and equal to Young’s contact angle. It can be seen that this behavior is correctly reproduced in the simulation as shown in the insets of Fig. 5. To investigate this behavior in more detail, the shapes of the droplet contour near the triple-phase contact line were calculated. For this purpose, the triangulation was intersected with a plane through the apex of the drop, and the intersection points were determined. The shapes of the resulting curves as well as the semi-analytical calculated curves are shown in Fig. 6(b). It becomes clear that for small voltages as well as large distances from the triple-phase contact line, the curves are equal. However, strong deviations become apparent near the triple-phase contact line. As predicted, the constant local contact angle leads to a significantly increased curvature of the actual interface shape in a range of the thickness of the dielectric intermediate layer around the triple-phase contact line.^{53} Based on the calculated curves, the local contact angle can be determined. The contact angles determined in this way are shown together with the theoretical curve according to Eq. (4) in Fig. 7. It becomes clear that the simulation correctly predicts the voltage independence of the local contact angle. The remaining very minor numerical deviations are most likely discretization errors caused by the finite number and size of the triangles. The correct description of the local wetting behavior is inherent to the method and one of the great strengths of numerical EWOD modeling in microfluidics.

### B. Contact-angle saturation

An experimentally observed effect in EWOD experiments is the so-called contact-angle saturation. This refers to the fact that the contact angle does not decrease further once a certain voltage value has been reached. Many different theories have been presented to explain contact-angle saturation, for example, trapped charges within the dielectric intermediate layer,^{54} vanishing interfacial energy $\sigma sl$,^{23–25} the influence of the wire electrode,^{55} injection of ions at the triple-phase contact line,^{56} or a dielectric breakdown of the intermediate layer.^{57} The variety of different theories shows that the phenomenon of contact-angle saturation is still not well understood. However, it is clear that the field distribution and the high field strengths in the vicinity of the three-phase contact line play a decisive role. In addition, there is the question of whether contact-angle saturation is an inherent property of the EWOD effect even considering ideal dielectrics and an ideally conducting liquid. Mugele and Buehrle^{47} have shown that this is not the case in two dimensions. However, it could hitherto not be clarified whether the saturation possibly occurs for ideal dielectrics in a three-dimensional system. With the help of the method presented here, this question can be answered conclusively. To do this, droplets were simulated under the influence of high voltages. No contact-angle saturation was observed in any of these simulation. Therefore, we conclude that no contact-angle saturation is observed in three dimensions as well for ideal dielectrics if only the sum of the electrostatic, the gravitational, and the interfacial energy of the system is minimized and no further effects are taken into account. The saturation effect is, therefore, not an inherent property of the EWOD effect and must be caused by another external effect.

### C. Open system—EWOD in micro-cavities

In order to be able to test the prediction accuracy of the presented simulation method also for open systems in which the volume is not preserved, the EWOD effect in closed micro-cavities is investigated. We choose this example since the EWOD effect in closed micro-cavities can be used to build micropumps as shown by Bohm *et al*.^{21} or Hoffmann *et al*.,^{20} and the accurate prediction of the volume stroke becomes particularly important for the dimensioning of the devices. In this case, the thickness of the dielectric interlayer is not negligible compared to the size of each individual micro-cavity, resulting in non-negligible discrepancies between the analytical prediction for the volume stroke based on the Young–Lippmann equation and the experimental results. It is shown that the presented simulation method overcomes these disadvantages. The voltage difference across the dielectric interface of each micro-cavity leads to a decrease of the contact angle according to Eq. (4). The decrease of the contact angle is linked to a decrease of the interface curvature and as a result to a change of the Laplace pressure. The introduced pressure gradient initiates a flow, and the interface moves inside the micro-cavity. The vapor inside the cavity is compressed by the movement of the interface. The movement of the interface stops when the Laplace pressure $pL$ is equal to the difference of the vapor pressure $pV$ of the enclosed gas minus the applied overpressure at the inlet of the cavity $pin$. Therefore, the equilibrium condition is

This system has analytically been modeled by Matsumoto and Colgate^{58} for a simple cylindrical geometry. We briefly summarize their arguments: Under the assumption that the enclosed vapor can be described as ideal gas, the pressure is linked to the interface deflection by

where $pat$ denotes the pressure within the empty cavity, $L$ is the overall length of a micro-cavity, and $l(U)$ describes the voltage-dependent equilibrium position of the interface. Assuming a cylinder-shaped micro-cavity, the voltage-dependent Laplace pressure can be calculated by

The insertion of Eqs. (50) and (51) into the equilibrium condition given by Eq. (49) leads to an equation for the voltage-dependent equilibrium position of the interface,

where the voltage-dependent contact angle $\u03d1(U)$ is given by Eq. (4). If Eq. (52) is applied to cavities with non-cylindrical cross sections, $R$ is replaced by the hydraulic radius $Rh=2A/P$, where $A$ corresponds to the cross-sectional area and $P$ to the wetted perimeter of the cavity.

In order to compare the analytical predictions of this model with our simulation method, three-dimensional models of cavities were created and meshed. The initial simulation geometry of a rectangular cavity is shown in Fig. 8(a). The geometry consists of four domains. One outer domain forms the dielectric intermediate layer with a thickness of $dD=2.5\mu m$ and a permittivity of $\epsilon r,d=3$. The interior of the micro-cavity contains three domains, with the two domains in the lower region representing the gas and the upper domain representing the liquid. The gas domain was divided into two domains to reduce the number of freely moving mesh vertices for the moving mesh calculation. We present results for a mesh used consisting in total of $11134$ triangles. Therefore, equations with $22268$ degrees of freedom need to be solved in each field calculation. The voltage was gradually increased from 0 to 100 V in 2 V steps. A maximum of 30 field calculations were performed for each voltage value. Initially, 500 inner iterations are performed and reduced by a factor of $r=0.9$ after each field calculation. To shorten the simulation time, the last equilibrium position found for the interface is used as the starting point for the calculation at the next voltage value. The termination criteria $tolabs=1\xd710\u221215$ and $tolrel=1\xd710\u22125$ were used.

First, the static interface shape was calculated for different values of the applied overpressure $pin$, and the obtained values for the deflection $li=Vi/A$ were compared with the analytical predictions of Eq. (52) for $U=0V$. The resulting equilibrium shapes of the liquid are presented for different values of the inlet pressure in Fig. 9. The results for the calculated deflection $l$ are shown in Fig. 8(b). It becomes clear that the results agree very well.

In addition, the deflection of the fluid interface was investigated for different values of the applied voltage. Figure 10 shows the final deflection of the interface for the different voltage values. At the beginning of the simulation, a maximum of 10 000 inner iterations is carried out in order to obtain a well converged equilibrium shape of the interface as a starting position. Examples for the calculated equilibrium shapes of the interfaces are illustrated in Fig. 8(c) for different values of the applied voltage. Even for cylindrical micro-cavities, the interface shape deviates substantially from a spherical segment near the triple-phase contact line. This can again be explained by the “need” of the equilibrium surface to accommodate the Young condition for the local contact angle. This effect is also noticeable in the prediction of the volume stroke. The analytical model systematically overestimates the volume stroke. We emphasize again that in contrast to, e.g., the analytic model according to Eq. (52) and the direct numerical solution of the Young–Laplace partial differential equation, our method does not rely on any assumptions on the symmetry or convexity of the interfaces and, thus, can be used to simulate cavities with arbitrary cross sections and shapes.

## VI. MEASUREMENTS

The numerical results of Sec. V can, of course, be compared to experimental data. We do so for quadratic micro-cavities because the simulation results confirming the prediction of a voltage-independent local contact angle of sessile droplets in classical EWOD experiments presented in Sec. V A can at least qualitatively be compared to measurements done by Mugele and Buehrle.^{47} Therefore, we chose the deviation between the analytical and the simulated volume stroke within a closed micro-cavity as shown in Fig. 10 as an experimental benchmark for the simulation method. In order to do so, test chips featuring micro-cavities with quadratic cross sections with an edge length $aC$ were manufactured at the Center for Micro- and Nanotechnologies of the Technische Universität Ilmenau. As a starting point, highly doped (boron) silicon wafers are used. Afterwards, a mask was applied using conventional UV-lithography, and the cavities are generated by DRIE (deep reactive ion etching). The photoresist is then removed, and an $2.15\mu m$ thick oxide layer is created by wet thermal oxidation. Then, the oxide is removed on the back side in hydrofluoric acid using a single side etching device. Directly afterward, an approximately 150 nm thick layer of aluminum is deposited on the backside of the wafers by sputtering. The aluminum layer is used to establish the electrical contact to the measuring chip. Then, a protective layer of photoresist is applied to the front side, and the wafers are diced into $20\xd720mm2$ large pieces. Finally, the protective layer of photoresist is removed, and a hydrophobic layer is applied. For the hydrophobization, a layer of DuPont Teflon^{TM}-AF1600 with a thickness of approximately 330 nm is deposited by dip coating. The Teflon layer is cured for 2 h at $170\xb0C$. Figure 11 shows a three-dimensional model of the resulting chips. A milled part made of PMMA was utilized as a clamping device for the chips. The necessary seals as well as the fluidic and electric connections were integrated into this part. One of the fluidic connectors has been modified so that a thick-walled capillary can be connected directly. The inside of the capillary was hydrophobized with FDTS [perfluorodecyltrichlorosilane $CF3(CF2)7(CH2)2SiCl3$]. The capillary has an inner diameter of 1 mm. The upper end of the capillary was connected to a compressed air source so that the pressure $pin$ at the inlet of the capillary can be adjusted. The set pressure was measured with a GDH14AN digital manometer. The determination of the volume stroke was carried out for each voltage value for two different overpressure values: A low overpressure in the range of approximately 110 kPa and an increased pressure in the range of 130 kPa. Two measurements were taken for each pressure value. All mentioned results correspond to the respective mean values. A schematic representation of the experimental setup is shown in Fig. 12. A neMESYS low-pressure syringe pump from Cetoni was connected and used for filling the system. A precision glass syringe with a total volume of 2.5 ml was used in the syringe pump. All experiments were performed with DI water. The surface tension was previously determined using a DSA10 tensiometer from Krüss with the pendant drop method. The determined value for the surface tension was $\sigma lv=0.073Jm\u22122$. For the simulations and the calculations, the value of the initial contact angle $\u03d1Y$ is required. For this purpose, optical contact-angle measurements were first carried out on unstructured areas of each test chip. A value of approximately $110\xb0$ was determined. The voltage signals were generated using a DS345 function generator from SRS and amplified by a voltage amplifier. A sinusoidal voltage with a frequency of $0.2Hz$ was used. The amplitude was increased in 10 V steps every 20 s from 0 to 100 V. The resulting output voltage was recorded with a 4224 PicoScope card. The generation of the voltage signals was automated in order to achieve comparable measurement conditions. During the actuation, the movement of the liquid–air interface within the capillary was filmed. The videos were recorded utilizing a Keyence VH-Z20R zoom objective in combination with a Keyence VW-600 C camera. Afterward, the deflection of the interface was determined using the Keyence camera’s internal motion tracking software. The measured deflection of the liquid–air interface within the measuring capillary corresponds to the total liquid stroke of all micro-cavities. Since the number of micro-cavities is known, the fluid stroke per micro-cavity can be determined directly. All experiments were carried out under clean room conditions. The measured values can be directly compared with the predictions of the simulation as well as the analytical model, which is given by Eq. (52). The results are summarized in Fig. 13. It becomes clear that the theory according to Eq. (52) constantly overestimates the voltage-dependent deflection. This can be explained by the independence of the local contact angle from the applied voltage, as already discussed in Sec. V C. In addition, for each geometry as well as each value of overpressure, a simulation was carried out as described in Sec. V C. It becomes clear that the simulation results agree very well with the measurement results. As expected, the deflection of the interfacial motion increases with increasing etch depth. At the same time, an increase of the applied overpressure leads to a decreasing voltage dependent interface deflection. This can be understood by means of a mechanical equivalent circuit diagram. The enclosed vapor acts as a spring and generates the restoring force for the interface. For a decreasing etch depth or an increased overpressure, this spring becomes more and more stiff, leading to a decrease of the deflection of the fluid–air interface and, therefore, to a decreased volume stroke.

## VII. SUMMARY AND CONCLUSIONS

A novel simulation method was presented that efficiently predicts the voltage-dependent static three-dimensional interface shapes with high accuracy. The simulation method correctly predicts the nontrivial shape of the interface even close to the three-phase contact line very accurately. The presented method can be applied to a system with an arbitrary number of phases, position-dependent wetting properties in, e.g., surface modified systems, partially conductive liquids, multiple electrodes, or with multiple dielectrics with different permittivities. The simulation method was implemented in an easy-to-use program that expects only the meshed geometry as well as the boundary conditions and the material properties as input. Being an efficient and flexible simulation method for static three-dimensional fluid interface shapes, our method has a considerable potential for the rational design of EWOD-driven microfluidic systems. Future work could include implementation of more sophisticated minimization procedures. For example, the simulation time could be further reduced using a more efficient optimization algorithm like a quasi-Newton method with improved convergence properties compared to the steepest descent method. In this case, a limited-memory BFGS method^{59} is particularly suitable because the approximated Hessian matrix does not have to be stored explicitly, which is large due to the high number of vertices. Currently, the simulation time is mainly determined by the time needed to calculate the electric field. In the future, ways of compressing the system matrix, such as those presented by Zhao *et al*.^{60} or Liu,^{61} could be used, which would further reduce simulation time significantly.

## SUPPLEMENTARY MATERIAL

See the supplementary material to further facilitate understanding of the paper. The material provides detailed derivations to the equations listed in Table I. In addition, the underlying simulations of the provided videos are described. Further simulation examples are also discussed, including the simulation of the classical EWOD experiment of a sessile droplet on a thin dielectric layer, the deformation of the interface of a conductive liquid under the influence of multiple electrodes, and the deformation of a non-conductive droplet by two planar electrodes, i.e., dielectrowetting.

## ACKNOWLEDGMENTS

We thank the Federal Ministry of Economics and Energy (BMWi) for the funding of the research project within the framework of the Zentrales Innovationsprogramm Mittelstand (ZIM, Funding No. ZF4457306PO9). In addition, we would like to thank the staff of the HPC-cluster of the TU Ilmenau and, in particular, Henning Schwanbeck for providing the excellent computing environment. Furthermore, we would like to thank Thomas Boeck for fruitful discussions and valuable comments. Support by the Center of Micro- and Nanotechnologies (ZMN) (DFG RIsources Reference No. RI_00009), a DFG-funded core facility of the TU Ilmenau, is gratefully acknowledged.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Sebastian Bohm:** Conceptualization (equal); Project administration (equal); Software (lead); Visualization (lead); Writing – original draft (lead). **Erich Runge:** Conceptualization (equal); Project administration (equal); Validation (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*18th IEEE International Conference on Micro Electro Mechanical Systems, 2005. MEMS 2005*(IEEE, 2005), pp. 227–230.

*COMSOL Conference 2020*(COMSOL, 2020).

*Equilibrium Capillary Surfaces*, Grundlehren der Mathematischen Wissenschaften: A Series of Comprehensive Studies in Mathematics Vol. 284 (Springer, New York, 1986).

*Variational and Free Boundary Problems*, The IMA Volumes in Mathematics and Its Applications Vol. 53, edited by A. Friedman (Springer, New York, 1993).

*2014 4th International Conference on Engineering Technology and Technopreneuship (ICE2T)*(IEEE, 2014), pp. 289–294.

*2018 IEEE SENSORS*(IEEE, 2018), pp. 1–4.

*Iranian Conference on Electrical Engineering (ICEE)*(IEEE, 2018), pp. 48–52.

*Numerical Optimization*, 2nd ed., Springer Series in Operations Research and Financial Engineering (Springer Science, New York, 2006).

*IEEE Proceedings on Micro Electro Mechanical Systems, An Investigation of Micro Structures, Sensors, Actuators, Machines and Robots*(IEEE, 1990), pp. 105–110.