Multiphysics simulation of fluid interface shapes in microfluidic systems driven by electrowetting on dielectrics

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,¹³ 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¹⁹ 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,²² 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.

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 $ϑY$ at the triple-phase contact line is given by Young’s equation,

where $σlv$⁠, $σsv$⁠, and $σsl$ denote the respective interfacial energies. The change of the effective interfacial energy $σsleff$ as a function of the applied voltage $U$ is described by the Lippmann equation,¹² 

where $ρ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 $ρsl=ε0εr,ddDU$ leads to

where $dD$ denotes the thickness of the dielectric layer, $ε0$ the vacuum permittivity, and $ε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.¹¹ Then, the effective interfacial energy $σ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 $η≡ε0εr,d2dDσ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:²⁶ 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 $Πel∝E2$ acting on the surface. They showed that even though $Πel$ is proportional to $δrν$ with $−1<ν<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.

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

Here, $∇n$ 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:

• $∙$

  The boundary on which the boundary condition (8) must apply is not known at first (free boundary problem³¹). Direct numerical solution methods as shown, for example, by Cooray et al.³² are, therefore, based on iterative methods.

• $∙$

  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.

• $∙$

  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.³⁵ 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.⁴⁴ 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,¹⁰ 

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,⁴⁵ 

Here, $Δϱ=ϱF−ϱV$ 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=ε0εrE$ the electrical displacement. The second term $ΔΠ$ 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,⁴⁶ 

The depth of penetration of the electric field into the electrolyte is given by the Debye length $λ=ε0εrkBT/(∑inibqi2)$⁠. For example, the Debye length for a 1:1 electrolyte, such as an NaCl-solution with a concentration of $0.1moll−1$⁠, is less than 1 nm at room temperature and is much smaller than 1 $μ$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φkBT≪1$ is valid. In this case, the potential can be determined using the linear Poisson–Boltzmann equation,

where $κ=λ−1$ 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.⁴⁷ 

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.²⁶ 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 $σ$ 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 $∇Asl=−∇Asv$ and Eq. (1), it takes the form

Here, $∇f$ is the multi-dimensional gradient $(∂f/∂X1,x,∂f/∂X1,y,∂f/∂X1,z,…,∂f/∂XN,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

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 $∇E−p(∇V)$ and step size $α$ 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,

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 $δW$⁠:

Thus, the gradient $∇δ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.

We present our solution method using the expression of $∇FO$⁠. However, the method for closed systems is identical except the equation for $∇F$⁠. 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 $∇Alv,∇Asl$⁠, $∇∫zdV$⁠, $V$⁠, and $∇V$ 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 $αni$ denotes the step size. In the following, the abbreviation $∇F(Xni)=∇Fni$ 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 $αni$ has to fulfill the well-known strong Wolfe conditions:⁴⁸ 

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−Fni$ 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 $∇F(Xni−tl∇Fni)$ 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=α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 $α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×10−4$ and $c2=0.9$⁠, which are suggested by Nocedal and Wright.⁴⁹ To speed up the calculation, the last value found for $α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.

The total volume $V$ as well as the occurring gradient terms $∇V$⁠, $∇A$⁠, and $∇∫zdV$ 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.

The terms $∇V$⁠, $∇A$⁠, and $∇∫zdV$ have the dimension $3N×1$⁠, where $N$ is the number of mesh vertices. The wetting properties of the liquid and the surface are completely determined by the contact angle $ϑY$ and the surface tension $σlv$⁠. To ensure that the properties are correctly taken into account, the term $∇A$ is split into the two parts $∇Alv$ and $∇Asl$ for each vertex within each triple-phase contact line. Here, $∇Asl$ denotes the gradient of the wetted surface area, and $∇Alv$ is the gradient of the liquid–vapor interface. For all other vertices within the liquid–vapor interface, only the gradient $∇Alv$ 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(μ)−ajμ$ 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×K$⁠, 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.

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⁵⁰ as

where $j$ runs over all neighboring vertices of the vertex $Xi$⁠. The angles $β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 $δk,k=x,y,z$ of the displacement vector $δ$⁠:

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 $∇V$ and $∇A$ 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.

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 $∑μ=1tjwμ=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.

The remaining task in calculating the gradient $∇F$ 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 $ϕ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 $ϕnk$ and $Enk$⁠. On an outer boundary, a Dirichlet or Neumann condition is specified and either $ϕnk$ or $Enk$ is given. On the interfaces between two domains $I$ and $II$⁠, however, the four unknowns $ϕnI$⁠, $ϕ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 $Λn$ denotes a potential user-defined surface charge density on the $n$th triangle. Otherwise, $Λ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 $ζm$ and $ηm$ are calculated according to the GQUTM (Gaussian quadrature formula for unit triangles method).⁵¹ The integration is carried out over the canonical triangle ${0≤η≤1,0≤ζ≤1−η}$ using the linear mapping $X(ζm,ηm)=(1−ζm−ηm)Xi,1+ζmXi,2+η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 $κ=0$ holds and the Poisson–Boltzmann equation becomes the Laplace equation, these integrals can be calculated analytically as follows:

Unfortunately, in the case of $κ>0$⁠, the integrals $ΔUn,nk$ can no longer be integrated in a closed form.⁵² However, using the same idea that was used to derive the equations (46), the integrals can be put into the form

while $Δ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⋅Ngrad0$ with $0<r≤1$ 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.

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 $φ$ equals the specified voltage-dependent contact angle $ϑY$ calculated using Eq. (4). Here, $Bl=σlv/(Δϱg)$ 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−π∫0ymaxx2dy)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 $Δϱ$⁠. 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 $Δϱ$ 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μ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μ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.²⁶ 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.⁵³ 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.

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,⁵⁴ vanishing interfacial energy $σsl$⁠,^23–25 the influence of the wire electrode,⁵⁵ injection of ions at the triple-phase contact line,⁵⁶ or a dielectric breakdown of the intermediate layer.⁵⁷ 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⁴⁷ 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.

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.²¹ or Hoffmann et al.,²⁰ 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⁵⁸ 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 $ϑ(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μm$ and a permittivity of $ε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×10−15$ and $tolrel=1×10−5$ were used.