Ja n 20 09 Topology and shape optimization of induced-charge electro-osmotic micropumps

For a dielectric solid surrounded by an electrolyte and positioned inside an externally biased parallel-plate capacitor, we study numerically how the resulting induced-charge electro-osmotic (ICEO) flow depends on the topology and shape of the dielectric solid. In particular, we extend existing conventional electrokinetic models with an artificial design field to describe the transition from the liquid electrolyte to the solid dielectric. Using this design field, we have succeeded in applying the method of topology optimization to find system geometries with non-trivial topologies that maximize the net induced electroosmotic flow rate through the electrolytic capacitor in the direction parallel to the capacitor plates. Once found, the performance of the topology optimized geometries has been validated by transferring them to conventional electrokinetic models not relying on the artificial design field. Our results show the importance of the topology and shape of the dielectric solid in ICEO systems and point to new designs of ICEO micropumps with significantly improved performance. Submitted to: New J. Phys.


Introduction
Induced-charge electro-osmotic (ICEO) flow is generated when an external electric field polarizes a solid object placed in an electrolytic solution [1,2]. Initially, the object acquires a position-dependent potential difference ζ relative to the bulk electrolyte. However, this potential is screened out by the counterions in the electrolyte by the formation of an electric double layer of width λ D at the surface of the object. The ions in the diffusive part of the double layer then electromigrate in the resulting external electric field tangential to the surface of the object, and by viscous forces they drag the fluid along. At the outer surface of the double layer, a resulting effective slip velocity is thus established. For a review of ICEO see Squires and Bazant [3].
The ICEO effect may be utilized in microfluidic devices for fluid manipulation, as proposed in 2004 by Bazant and Squires [1]. Theoretically, various simple dielectric shapes have been analyzed for their ability to pump and mix liquids [3,4]. Experimentally ICEO was observed and the basic model validated against particle image velocimetry in 2005 [2], and later it was used in a microfluidic mixer, where a number of triangular shapes act as passive mixers [5]. However, no studies have been carried out concerning the impact of topology changes of the dielectric shapes on mixing or pumping efficiency. In this work, we focus on the application of topology optimization to ICEO systems. With this method it is possible to optimize the dielectric shapes for many purposes, such as mixing and pumping efficiency.
Our model system consists of two externally biased, parallel capacitor plates confining an electrolyte. A dielectric solid is shaped and positioned in the electrolyte, and the external bias 3 induces ICEO flow at the dielectric surfaces. In this work, we focus on optimizing the topology and shape of the dielectric solid to generate the maximal flow perpendicular to the external applied electric field. This example of establishing an optimized ICEO micropump serves as a demonstration of the implemented topology optimization method.
Following the method of Borrvall and Petersson [6] and the implementation by Olesen et al [7] of topology optimization in microfluidic systems, we introduce an artificial design field γ (r) in the governing equations. The design field varies continuously from zero to unity, and it defines to what degree a point in the design domain is occupied by dielectric solid or electrolytic fluid. Here, γ = 0 is the limit of pure solid and γ = 1 is the limit of pure fluid, while intermediate values of γ represent a mixture of solid and fluid. In this way, the discrete problem of placing and shaping the dielectric solid in the electrolytic fluid is converted into a continuous problem, where the sharp borders between solid and electrolyte are replaced by continuous transitions throughout the design domain. In some sense one can think of the solid/fluid mixture as a sort of ion-exchange membrane in the form of a sponge with varying permeability. This continuum formulation allows for an efficient gradient-based optimization of the problem.
In one important aspect our system differs from other systems previously studied by topology optimization: induced-charge electro-osmosis is a boundary effect relying on the polarization and screening charges in a nanometer-sized region around the solid/fluid interface. Previously studied systems have all been relying on bulk properties such as the distribution of solids in mechanical stress analysis [8], photonic band gap structures in optical wave guides [9], and acoustic noise reduction [10], or on the distribution of solids and liquids in viscous channel networks [6,7,11] and chemical microreactors [12]. In our case, as for most other applications of topology optimization, no mathematical proof exists that the topology optimization routine indeed will result in an optimized structure. Moreover, since the boundary effects of our problem result in a numerical analysis that is very sensitive to the initial conditions, on the meshing, and on the specific form of the design field, we take the more pragmatic approach of finding non-trivial geometries by the use of topology optimization, and then validate the optimality by transferring the geometries to conventional electrokinetic models not relying on the artificial design field.

Model system
We consider a parallel-plate capacitor externally biased with a harmonic oscillating voltage difference φ = 2V 0 cos(ωt) and enclosing an electrolyte and a dielectric solid. The two capacitor plates are positioned at z = ±H and periodic boundary conditions are applied at the open ends at x = ±L/2. The resulting bound domain of size L × 2H in the x z-plane is shown in figure 1. The system is assumed to be unbounded and translational invariant in the perpendicular y-direction. The topology and shape of the dielectric solid are determined by the numerical optimization routine acting only within a smaller rectangular, central design domain of size l × 2h. The remaining domain outside this area is filled with pure electrolyte. Double layers, or Debye screening layers, are formed in the electrolyte around each piece of dielectric solid to screen out the induced polarization charges. The pull from the external electric field on these screening layers in the design domain drives an ICEO flow in the entire domain.
The magnitude of the ICEO effect increases as the permittivity ε diel of the dielectric solid increases relative to that of the electrolyte, ε fluid = 78ε 0 . For maximum effect we therefore Figure 1. (a) A sketch of the rectangular L × 2H cross-section of the electrolytic capacitor in the x z-plane. The external voltages φ 1 and φ 2 are applied to the two infinite parallel-plate electrodes (thick black lines) at z = ±H . The voltage difference φ 1 − φ 2 induces an ICEO flow around the un-biased dielectric solid (dark gray) shaped by the topology optimization routine limited to the rectangular l × 2h design domain (light gray). The dielectric solid is surrounded by pure electrolyte (light gray and white). Periodic boundary conditions are applied at the vertical edges (dotted lines). (b) The dimensionless electric permittivity ε as a function of the design variable γ . (c) Zoom-in on the rapid convergence of ε(γ ) toward ε fluid = 1 for γ approaching unity after passing the value γ cut − γ 0.98. employ ε diel = 10 6 ε 0 throughout this work. For such a large value of the permittivity, the potential inside the dielectric solid is nearly constant as for a metal. Placed in an external electric field E, the largest potential difference, denoted the zeta potential ζ , between the dielectric solid and the surrounding electrolyte just outside the screening double layer can thus be estimated as ζ = Ea [3].
If the dielectric solid is symmetric around the x-axis, the antisymmetry of the applied external bias voltage ensures that the resulting electric potential is antisymmetric and the velocity and pressure fields are symmetric around the center plane z = 0. This symmetry appears in most of the cases studied in this paper, and when present it is exploited to obtain a significant decrease in memory requirements of the numerical calculations.
The specific goal of our analysis is to determine the topology and shape of the dielectric solid such that a maximal flow rate Q is obtained parallel to the x-axis, i.e. perpendicular to the direction of external potential field gradient.

Governing equations
We follow the conventional continuum approach to the electrokinetic modeling of the electrolytic capacitor [3]. For simplicity we consider a symmetric, binary electrolyte, where the positive and negative ions with concentrations c + and c − , respectively, have the same diffusivity D and valence charge number Z .

Bulk equations in the conventional ICEO model
Neglecting chemical reactions in the bulk of the electrolyte, the ionic transport is governed by particle conservation through the continuity equation, where J ± is the flux density of the two ionic species, respectively. Assuming a dilute electrolytic solution, the ion flux densities are governed by the Nernst-Planck equation, where the first term expresses ionic diffusion and the second term ionic electro-migration due to the electrostatic potential φ. Here, e is the elementary charge, T the absolute temperature and k B the Boltzmann constant. We note that due to the low fluid velocity v obtained in the ICEO systems under consideration, we can safely neglect the convective ion fluxes c ± v throughout this paper (see table 2). The electrostatic potential φ is determined by the charge density ρ el = Z e(c + − c − ) through Poisson's equation, where ε fluid is the fluid permittivity, which is assumed constant. The fluid velocity field v and pressure field p are governed by the continuity equation and the Navier-Stokes equation for incompressible fluids, where ρ m and η are the fluid mass density and viscosity, respectively, both assumed constant.

The artificial design field γ used in the topology optimization model of ICEO
To be able to apply the method of topology optimization, it is necessary to extend the conventional ICEO model with three additional terms, all functions of a position-dependent artificial design field γ (r). The design field varies continuously from zero to unity, where γ = 0 is the limit of a pure dielectric solid and γ = 1 is the limit of a pure electrolytic fluid. The intermediate values of γ represent a mixture of solid and fluid. The first additional term concerns the purely fluid dynamic part of our problem. Here, we follow Borrvall and Petersson [6] and model the dielectric solid as a porous medium giving rise to a Darcy friction force density −α(γ )v, where α(γ ) may be regarded as a local inverse permeability, which we denote the Darcy friction. We let α(γ ) be a linear function of γ of the form α(γ ) = α max (1 − γ ), where α max = η/ 2 pore is the Darcy friction of the porous dielectric material assuming a characteristic pore size pore . In the limit of a completely impenetrable solid the value of α max approaches infinity, which leads to a vanishing fluid velocity v. The modified Navier-Stokes equation extending to the entire domain, including the dielectric material, becomes The second additional term is specific to our problem. Since the Navier-Stokes equation is now extended to include also the porous dielectric medium, our model must prevent the unphysical penetration of the electrolytic ions into the solid. Following the approach of Kilic et al [13], where current densities are expressed as gradients of chemical potentials, J ∝ −∇µ, we model the ion expulsion by adding an extra free energy term κ(γ ) to the chemical potential µ ± = ±Z eφ + k B T ln(c ± /c 0 ) + κ(γ ) of the ions, where c 0 is the bulk ionic concentration for both ionic species. As above we let κ(γ ) be a linear function of γ of the form κ(γ ) = κ max (1 − γ ), where κ max is the extra energy cost for an ion to enter a point containing a pure dielectric solid as compared with a pure electrolytic fluid. The value of κ max is set to an appropriately high value to expel the ions efficiently from the porous material while still ensuring a smooth transition from dielectric solid to electrolytic fluid. The modified ion flux density becomes The third and final additional term is also specific to our problem. Electrostatically, the transition from the dielectric solid to the electrolytic fluid is described through the Poisson equation by a γ -dependent permittivity ε(γ ). This modified permittivity varies continuously between the value ε diel of the dielectric solid and ε fluid of the electrolytic fluid. As above, we would like to choose ε(γ ) to be a linear function of γ . However, during our analysis using the aforementioned large value ε diel = 10 6 ε 0 for the solid in an aqueous electrolyte with ε fluid = 78ε 0 , we found unphysical polarization phenomena in the electrolyte due to numerical rounding-off errors for γ near, but not equal to, unity. To overcome this problem we ensured a more rapid convergence toward the value ε fluid by introducing a cut-off value γ cut 0.98, a transition width γ 0.002, and the following expression for ε(γ ), For γ γ cut we obtain the linear relation ε(γ ) = ε diel + (ε fluid − ε diel )γ , while for γ γ cut we have ε(γ ) = ε fluid , see figures 1(b) and (c). For γ sufficiently close to unity (and not only when γ equals unity with numerical precision), this cut-off procedure ensures that the calculated topological break-up of the dielectric solid indeed leads to several correctly polarized solids separated by pure electrolyte. The modified Poisson equation becomes Finally, we introduce the γ -dependent quantity, the so-called objective function [γ ], to be optimized by the topology optimization routine: the flow rate in the x-direction perpendicular to the applied potential gradient. Due to incompressibility, the flow rate Q(x) is the same through cross-sections parallel to the yz-plane at any position x. Hence, we can use the numerically 7 more stable integrated flow rate as the objective function, where is the entire geometric domain (including the design domain) andn x the unit vector in the x-direction.

Dimensionless form
To prepare the numerical implementation, the governing equations are rewritten in dimensionless form, using the characteristic parameters of the system. In conventional ICEO systems the size a of the dielectric solid is the natural choice for the characteristic length scale 0 , since the generated slip velocity at the solid surface is proportional to a. However, when employing topology optimization we have no prior knowledge of this length scale, and thus we choose it to be the fixed geometric half-length 0 = H between the capacitor plates. Further characteristic parameters are the ionic concentration c 0 of the bulk electrolyte and the thermal voltage φ 0 = k B T /(Z e). The characteristic velocity u 0 is chosen as the Helmholtz-Smoluchowski slip velocity u 0 , which for the induced zeta potential ζ = E 0 in a local electric field E = φ 0 / 0 tangential to the surface is [3] The pressure scale is set by the characteristic microfluidic pressure scale p 0 = ηu 0 / 0 . Although strictly applicable only to parallel-plate capacitors, the characteristic time τ 0 of the general system in chosen as the RC charging time of the double layer in terms of the Debye length λ D of the electrolyte [14], Moreover, three characteristic numbers are connected to the γ -dependent terms in the governing equations: the characteristic free energy κ 0 , the characteristic permittivity chosen as the bulk permittivity ε fluid and the characteristic Darcy friction coefficient α 0 . In summary, The new dimensionless variables (denoted by a tilde) thus becomẽ In the following, all variables and parameters are made dimensionless using these characteristic numbers and for convenience the tilde is henceforth omitted.

Linearized and reformulated equations
To reduce the otherwise very time-and memory-consuming numerical simulations, we choose to linearize the equations. There are several nonlinearities to consider. By virtue of a low Reynolds number Re ≈ 10 −6 , see table 2, the nonlinear Navier-Stokes equation is replaced by the linear Stokes equation. Likewise, as mentioned in section 3.1, the low Péclet number Pé ≈ 10 −3 allows us to neglect the nonlinear ionic convection flux density c ± v. This approximation implies the additional simplification that the electrodynamic problem is independent of the hydrodynamics.
Finally, we employ the linear Debye-Hückel approximation, which is valid when the electric energy Z eζ of an ion traversing the double layer is smaller than the thermal energy k B T . For our system we find Z eζ /k B T 0.5, so the linear Debye-Hückel approximation is valid, and we can utilize that the ionic concentrations only deviate slightly from the bulk equilibrium ionic concentration. The governing equations are reformulated in terms of the average ion concentration c ≡ (c + + c − )/2 and half the charge density ρ ≡ (c + − c − )/2. Thus, by expanding the fields to first order as c = 1 + δc and ρ = 0 + δρ, the resulting differential equation for ρ is decoupled from that of c. Introducing complex field notation, the applied external bias voltage is φ(t) = 2V 0 cos(ωt) = Re[2V 0 exp(iωt)], yielding a corresponding response for the potential φ and charge density ρ, with the complex amplitudes (r) = φ R (r) + i I (r) and P(r) = P R (r) + iP I (r), respectively. The resulting governing equations for the electrodynamic problem are then where we have introduced the dimensionless thickness of the linear Debye layer λ = λ D / 0 . Given the electric potential φ and the charge density P, we solve for the time-averaged hydrodynamic fields v and p , where the time-averaged electric body force density f el is given by

Boundary conditions
For symmetric dielectric solids we exploit the symmetry around z = 0 and consider only the upper half (0 < z < 1) of the domain. As boundary condition on the driving electrode, we set the spatially constant amplitude V 0 of the applied potential. Neglecting any electrode reactions 9 taking place at the surface, there is no net ion flux in the normal direction to the boundary with unit vectorn. Finally, for the fluid velocity we set a no-slip condition, and thus at z = 1 we have On the symmetry axis (z = 0) the potential and the charge density must be zero due to the antisymmetry of the applied potential. Furthermore, there is no fluid flux in the normal direction and the shear stresses vanish. So at z = 0 we have where the dimensionless stress tensor is , andn andt are the normal and tangential unit vectors, respectively, where the latter in 2D, contrary to 3D, is uniquely defined. On the remaining vertical boundaries (x = ±L/2 0 ), periodic boundary conditions are applied to all the fields. Corresponding boundary conditions apply to the conventional ICEO model equations (1)-(4b), without the artificial design field but with a hard-wall dielectric solid. For the boundary between a dielectric solid and electrolytic fluid, the standard electrostatic conditions apply; moreover, there is no ion flux normal to the surface, and a no-slip condition is applied to the fluid velocity.

Implementation and parameter values
For our numerical analysis we use the commercial numerical finite-element modeling tool Comsol [16] controlled by scripts written in Matlab [15]. The mathematical method of topology optimization in microfluidic systems is based on the seminal paper by Borrvall and Petersson [6], while the implementation (containing the method of moving asymptotes by Svanberg [17] 5 ) is taken from Olesen, Okkels and Bruus [7].
In Comsol all equations are entered in the divergence form ∇ · = F. The tensor contains all the generalized fluxes in the governing equations, one per row, e.g. J ± from equation (1), ε fluid ∇φ from equation (3), 0 from equation (4a) and ∇v + (∇v) T from equation (4b). The source term vector F contains all other terms in the governing equations that cannot be written as the divergence of a vector. Note that the continuity equation (4a) is not interpreted as an equation for v but rather as a divergence of a pressure-related flux (which is always zero), with the term ∇ · v acting as the corresponding source term, thus appearing in the vector F. As in [7], we have used second-order Lagrange elements for all fields except the pressure and the design field, for which first-order elements suffice. The resulting algebraic FEM-equations are solved using the sparse direct linear solvers UMFPACK or PARDISO. More details concerning the numerical implementation in Comsol can be found in [7]. Due to the challenges discussed in section 4.2 of resolving all length scales in the electrokinetic model, we have chosen to study small systems, 2H = 500 nm, with a relatively large Debye length, λ D = 20 nm. Our main goal is to provide a proof of concept for the use of topology optimization in electro-hydrodynamic systems, so the difficulty in fabricating actual devices on this sub-micrometric scale is not a major concern for us in the present work. A list of the parameter values chosen for our simulations is given in table 1.
For a typical topology optimization, like the one shown in figure 4(a), approximately 5400 FEM elements are involved. In each iteration loop of the topology optimization routine, three problems are solved: the electric problem, the hydrodynamic problem, and the so-called adjunct problem for the sensitivity analysis ∂ /∂γ (used to update the design field γ , see [7] for details), involving 4 × 10 4 , 2 × 10 4 and 7 × 10 4 degrees of freedom, respectively. On an Intel Core 2 Duo 2 GHz processor with 2 GB RAM, the typical CPU time is several hours.

Analytical and numerical validation by the conventional ICEO model
We have validated our simulations in two steps. Firstly, the conventional ICEO model not involving the design field γ (r) is validated against the analytical result for the slip velocity at a simple dielectric cylinder in an infinite ac electric field given by Squires and Bazant [3]. Secondly, the design field model is compared with the conventional model. This two-step validation procedure is necessary because of the limited computer capacity. The involved length scales in the problem make a large number of mesh elements necessary for the numerical solution by the finite element method. Four different length scales appear in the gammadependent model for the problem of a cylinder placed midway between the parallel capacitor plates: the distance H from the center of the dielectric cylinder to the biased plates, the radius a of the cylinder, the Debye length λ D and the length d over which the design field γ changes from zero to unity. This last and smallest length scale d in the problem is controlled directly by the numerical mesh size set up in the finite element method. It has to be significantly smaller than λ D to model the double layer dynamics correctly, so here a maximum value of the numerical mesh size is defined. Examples of meshing are shown in figure 2.
The analytical solution of Squires and Bazant [3] is only strictly valid in the case of an infinitely thin Debye layer in an infinite electric field. So, to compare this model with the bounded numerical model the plate distance must be increased to minimize the influence on the effective slip velocity. Furthermore, it has been shown in a numerical study by Gregersen et al [18] that the Debye length λ D should be about a factor of 10 3 smaller than the cylinder radius a to approximate the solution for the infinitely thin Debye layer model. Including the demand of d being significantly smaller than λ D , we end up with a length scale difference of at least 10 5 , which is difficult to resolve by the finite element mesh, even when mesh adaption is used. Consequently, we have checked that the slip velocity for the conventional model converges toward the analytical value when the ratio a/λ D increases, see figure 2(c). Afterwards, we have compared the solutions for the conventional and gamma-dependent models in a smaller system with a ratio of a/λ D ∼ 10 and found good agreement.

Validation of the self-consistency of the topology optimization
As an example of our validation of the self-consistency of the topology optimization, we study the dependence of the objective function Q = [ω, γ (ω, r)] on the external driving Note that structure γ a indeed yields the highest flow rate Q for ω = ω a , structure γ b maximizes Q for ω = ω b and structure γ c maximizes Q for ω = ω c .
frequency ω. As shown in figures 3(a)-(c), we have calculated the topology-optimized dielectric structures γ j = γ (ω j , r) where j = a, b, c, for three increasing frequencies ω = ω a = 1.25, ω = ω b = 12.5 and ω = ω c = 62.5. In the following, we let Q j (ω) denote the flow rate calculated at the frequency ω for a structure optimized at the frequency ω j . Firstly, we note that Q j (ω j ) decreases as the frequency increases above the characteristic frequency ω 0 = 1; Q a (ω a ) = 2.95 × 10 −3 , Q b (ω b ) = 1.82 × 10 −3 and Q c (ω c ) = 0.55 × 10 −3 . This phenomenon is a general aspect of ICEO systems, where the largest effect is expected to happen at ω ≈ 1 (in units of 1/τ 0 ).
Secondly, and most significant, we see in figure 3(d) that structure γ a is indeed the optimal structure for ω = ω a since Q a (ω a ) > Q b (ω a ), Q c (ω a ). Likewise, γ b is optimal for ω = ω b and γ c is optimal for ω = ω c .
We have gained confidence in the self-consistency of our topology optimization routine by carrying out a number of tests like the one in the example above.

Topology optimization
For each choice of parameters the topology optimization routine converges to a specific distribution of dielectric solid given by γ (r). As a starting point for the investigation of the optimization results, we used the parameters listed in table 1. As discussed above, the geometric dimensions are chosen as large as possible within the computational limitations: the Debye length is set to λ D = 20 nm and the distance between the capacitor plates to 2H = 500 nm. The external bias voltage V 0 is of the order of the thermal voltage V 0 = φ 0 = 25 mV to ensure the validity of the linear Debye-Hückel approximation. We let the bulk fluid consist of water containing small ions, such as dissolved KCl, with a concentration c 0 = 0.23 mM set by the chosen Debye length. As mentioned above, the permittivity of the dielectric solid is set to ε diel = 10 6 ε 0 for maximum ICEO effect. The artificial parameters κ and α max are chosen on a pure computational basis, where they have to mimic the real physics in the limits of fluid and solid, but also support the optimization routine when the phases are mixed.
Throughout our simulations we have primarily varied the applied frequency ω and the size l × 2h of the design domain. In figure 3 we have shown examples of large design domains with l × h = 2.0 × 0.8 covering 80% of the entire domain and frequency sweeps over three orders of magnitude. However, in the following we fix the frequency to be ω = 6.25, where the ICEO response is close to maximum. Moreover we focus on a smaller design domain l × h = 0.6 × 0.8 to obtain better spatial resolution for the given amount of computer memory and to avoid getting stuck in local optima. It should be stressed that the size of the design domain has a large effect on the specific form and size of the dielectric islands produced by the topology optimization. Also, it is important if the design domain is allowed to connect to the capacitor plates or not, see the remarks in section 6.
The converged solution found by topology optimization under these conditions is shown in figure 4(a). The shape of the porous dielectric material is shown together with a streamline plot of equidistant contours of the flow rate. We notice that many streamlines extend all the way through the domain from left to right, indicating that a horizontal flow parallel to the xaxis is indeed established. The resulting flow rate is Q = 2.99 × 10 −3 . The ICEO flow of this solution, based on the design-field model, is validated by transferring the geometrical shape of the porous dielectric medium into a conventional ICEO model with a hard-walled dielectric not depending on the design field. In the latter model the sharp interface between the dielectric solid and the electrolyte is defined by the 0.95-contour of the topology-optimized design field γ (r). The resulting geometry and streamline pattern of the conventional ICEO model is shown in figure 4(b). The flow rate is now found to be Q = Q * = 1.88 × 10 −3 . There is a close resemblance between the results of the two models both qualitatively and quantitatively. It is noticed how the number and positions of the induced flow rolls match well, and also the absolute values of the induced horizontal flow rates differ only by 37%, which is a small deviation as discussed in section 6.
Based on the simulation we can now justify the linearization of our model. The largest velocity u gap is found in the gap of width gap between the two satellite pieces and the central piece. As listed in table 2 the resulting Reynolds number is Re = 2.8 × 10 −7 , the Péclet number is Pé = 1.4 × 10 −3 and the Debye-Hückel number is Hü = 0.13.

Comparison to simple shapes
We evaluate our result for the optimized flow rate by comparing it to those obtained for more basic, simply connected, dielectric shapes, such as triangles and perturbed circles previously studied in the literature as flow inducing objects both analytically and experimentally [3]- [5].
For comparison, such simpler shapes have been fitted into the same design domain as used in the topology optimization figure 4(a), and the conventional ICEO model without the design field was solved for the same parameter set. In figure 5(a) the resulting flow for a triangle with straight  figure 5(b) the induced flow around a perturbed cylinder with radius r (θ ) = 0.24[1 + 0.5 cos(3θ)] is depicted. Again the shape has been fitted within the allowed design domain. The resulting flow rate Q = 0.46 × 10 −3 is higher than for the triangle, but still a factor of four slower than the optimized result. It is clearly advantageous to change the topology of the dielectric solid from simply to multiply connected.
For the topology-optimized shape in figure 4(a), it is noticed that only a small amount of flow is passing between the two closely placed dielectric islands in the upper left corner of the design domain. To investigate the importance of this separation, the gap between the islands was filled out with dielectric material and the flow calculated. It turns out that this topological change only lowered the flow rate slightly (15%) to a value of Q = 1.59 × 10 −3 . Thus, the important topology of the dielectric solid in the top-half domain is the appearance of one center island crossing the antisymmetry line and one satellite island near the tip of the center island.

Shape optimization
The topology-optimized solutions are found based on the extended ICEO model involving the artificial design field γ (r). To avoid the artificial design field, it is desirable to validate and further investigate the obtained topology-optimized results by the physically more correct conventional ICEO model involving hard-walled dielectric solids. We therefore extend the reasoning behind the comparison of the two models shown in figure 4 and apply a more stringent shape optimization to the various topologies presented above. With this approach we are gaining further understanding of the specific shapes comprising the overall topology of the dielectric solid. Moreover, it is possible to point out simpler shapes, which are easier to fabricate, but still perform well.
In shape optimization, the goal is to optimize the objective function , which depends on the position and shape of the boundary between the dielectric solid and the electrolytic Figure 6. (a) The streamline pattern (thick lines) for the shape-optimized rightangled triangle fixed at the symmetry line z = 0 calculated for ω = 6.25 using the conventional ICEO model with a hard-walled dielectric solid (black). In the full domain this is a triangle symmetric around z = 0. The flow rate is Q = 0.32 × 10 −3 . (b) As in panel (a) but without constraining the triangle to be right-angled. In the full domain the shape is a four-sided polygon symmetric around z = 0. The flow rate is Q = 0.76 × 10 −3 . Note that all sharp corners of the polygons have been rounded by circular arcs of radius 0.01.
fluid. This boundary is given by a line interpolation through a small number of points on the boundary. These control points are in turn given by N design variables g = (g 1 , g 2 , . . . , g N ), so the objective function of equation (9) depending on the design field γ (r) is now written as [g] depending on the design variables g, To carry out the shape optimization we use a direct bounded Nelder-Mead simplex method [19] implemented in Matlab [20] 6 . This robust method finds the optimal point g opt in the N -dimensional design variable space by initially creating a simplex in this space, e.g. an N -dimensional polyhedron spanned by N + 1 points, one of which is the initial guess. The simplex then iteratively moves toward the optimal point by updating one of the N + 1 points at a time. During the iteration, the simplex can expand along favorable directions, shrink toward the best point, or have its worst point replaced with the point obtained by reflecting the worst point through the centroid of the remaining N points. The iteration terminates once the extension of the simplex is below a given tolerance. We note that unlike topology optimization, the simplex method relies only on values of the objective function [g] and not on the sensitivity ∂ /∂ g. 7 First, we perform shape optimization on a right-angled triangle corresponding to the one shown in figure 5(a). Due to the translation invariance in the x-direction, we can without loss of generality fix one base point of the triangle (x 1 , 0) to the right end of the simulation domain, while the other (x 2 , 0) can move freely along the baseline, in contrast to the original rectangular design. To ensure a right-angled triangle only the z-coordinate of the top point (x 2 , z 2 ) may move freely. In this case the design variable becomes the two-component vector g = (x 2 , z 2 ). The optimal right-angled triangle is shown in figure 6(a). The flow rate is Q = 0.32 × 10 −3 or 1.5 times larger than that of the original right-angled triangle confined to the design domain.
If we do not constrain the triangle to be right-angled, we instead optimize a polygon shape spanned by three corner points in the upper half of the electrolytic capacitor. So, due to the symmetry of the problem, we are in fact searching for the most optimal, symmetric four-sided polygon. The three corner points are now given as (x 1 , 0), (x 2 , 0) and (x 3 , z 3 ), and again due to translation invariance, it results in a three-component design variable g = (x 2 , x 3 , z 3 ). The resulting shape-optimized polygon is shown in figure 6(b). The flow rate is Q = 0.76 × 10 −3 , which is 3.5 times larger than that of the original right-angled triangle confined to the design domain and 2.4 times better than that of the best right-angled triangle. However, this flow rate is still a factor of 0.4 lower than the topology-optimized results.
To be able to shape-optimize the more complex shapes of figure 4 we have employed two methods to obtain a suitable set of design variables. The first method, the radial method, is illustrated in figure 7. The boundary of a given dielectric solid is defined through a cubic interpolation line through N control points (x i , z i ), i = 1, 2, . . . , N , which are parameterized in terms of two co-ordinates (x c , z c ) of a center point, two global scale factors A and B, N lengths r i and N fixed angles θ i distributed in the interval from 0 to 2π, In this case the design variable becomes g = (x c , z c , r 1 , r 2 , . . . , r N , A, B). The second parametrization method involves a decomposition into harmonic components. As before we define a central point (x c , z c ) surrounded by N control points. However, now the distances r i are decomposed into M harmonic components given by where r 0 is an overall scale parameter and ϕ n is a phase shift. In this case the design variable becomes g = (x c , z c , r 0 , A 1 , A 2 , . . . , A M , ϕ 1 , ϕ 2 , . . . , ϕ M ) .

Comparing topology optimization and shape optimization
When shape-optimizing a geometry similar to the one found by topology optimization, we let the geometry consist of two pieces: (i) an elliptic island centered on the symmetry axis and fixed to the right side of the design domain and (ii) an island with a complex shape to be placed anywhere inside the design domain, but not overlapping with the elliptic island. For the ellipse we only need to specify the major axis A and the minor axis B, so these two design parameters add to the design variable listed above for either the radial model or the harmonic decomposition model. To be able to compare with the topology-optimized solution the dielectric solid is restricted to the design domain. The result of this two-piece shape optimization is shown in figure 8. Compared with the simply connected topologies, the two-piece shape-optimized systems yields much improved flow rates. For the shape optimization involving the radial method with 16 directional angles and A = B for the complex piece, the flow rate is Q = 1.92 × 10 −3 , figure 8(a), which is 2.5 times larger than that of the shape-optimized four-sided symmetric polygon. The harmonic decomposition method, figure 8(b), yields a flow rate of Q = 1.52 × 10 −3 or 2.0 times larger than that of the polygon.
All the results for the obtained flow rates are summarized in table 3. It is seen that two-piece shape-optimized systems perform as well as the topology-optimized system, when analyzed using the conventional ICEO model without the artificial design field. We also note by comparing figures 4 and 8 that the resulting geometry found using either topology optimization or shape optimization is essentially the same. The central island of the dielectric solid is a thin structure perpendicular to the symmetry axis and covering approximately 60% of the channel width. The satellite island of complex shape is situated near the tip of the central island. It has two peaks pointing toward the central island that seem to suspend a flow roll which guides the ICEO flow through the gap between the two islands.

Concluding remarks
The main result of this work is the establishment of the topology optimization method for ICEO models extended with the design field γ (r). In contrast to the conventional ICEO model with its sharply defined, impenetrable dielectric solids, the design field ensures a continuous transition between the porous dielectric solid and the electrolytic fluid, which allows for an efficient gradient-based optimization of the problem. By concrete examples we have shown how the use of topology optimization has led to non-trivial system geometries with a flow rate increase of nearly one order of magnitude, from Q = 0.22 × 10 −3 in figure 5(a) to Q = 1.92 × 10 −3 in figure 8(a). When applied to ICEO, the design field method is qualitatively but not quantitatively correct. We have found deviations of 37% when comparing design field simulations with hardwall simulations. The magnitude of the ICEO effect is sensitive to the exact configuration of the charge density and the electric field within the only 20 nm thick double layer. Even for the relatively simple hard-wall model, we have shown in another numerical analysis [18] how the magnitude of the ICEO depends on the width of the double layer. By introducing the design field γ and the associated artificial smoothing of the transition between the dielectric solid and the liquid, the electric properties within the double layer are changed by γ in several ways: the electric field is directly affected by the permittivity ε[γ (r)], while the charge density is directly affected by the chemical potential term κ[γ (r)] pushing ions away from the solid, and (to a lesser degree) by the viscous drag from the velocity field of the electrolyte, which is affected by the Darcy term −α[γ (r)]v. Quantitative agreement between the design field model and the hard-wall model can therefore only be achieved by an extremely fine (and in practice unreachable) resolution of the transition zone, the position of which is not even known a priori. Given this insight, the discrepancy of 37% may be regarded as relatively small. Moreover, in spite of the quantitative discrepancy, the design field method nevertheless produces qualitatively correct solutions with topologically non-trivial geometries, which, as we have shown, can easily be transferred to quantitatively correct hard-wall models and further improved by shape optimization.
The topology optimization algorithm of ICEO systems leads to many local optima, and we cannot be sure that the converged solution is the global optimum. The resulting shapes and the generated flow rates depend on the initial condition for the artificial γ -field. Generally, the initial condition used throughout this paper, γ = 0.99 in the entire design domain, leads to the most optimal results compared with other initial conditions. This initial value corresponds to a very weak change from the electrolytic capacitor completely void of dielectric solid. In contrast, if we let γ = 0.01 corresponding to almost pure dielectric material in the entire design region, the resulting shapes are less optimal, i.e. the topology optimization routine is more likely to be caught in a local optimum. Furthermore, the resulting shapes turn out to be mesh dependent as well. So, we cannot conclude much about global optima. Instead, we can use the topology-optimized shapes as inspiration to improve existing designs. For this purpose, shape optimization turns out to be a powerful tool. We have shown in this work, how shape optimization can be used efficiently to refine the shape of individual pieces of the dielectric solid once its topology has been identified by topology optimization.
For all three additional γ -dependent fields α(γ ) κ(γ ) and ε(γ ) we have used (nearly) linear functions. In many previous applications of topology optimization nonlinear functions have successfully been used to find global optima by gradually changing the nonlinearity into strict linearity during the iterative procedure [6]- [8], [12]. However, we did not improve our search for a global optimum by employing such schemes, and simply applied the (nearly) linear functions during the entire iteration process.
The limited size of the design domain is in some cases restricting the free formation of optimized structures. This may be avoided by enlarging the design domain. However, starting a topology optimization in a very large domain gives a huge number of degrees of freedom, and the routine is easily caught in local minima. These local minima often yield results not as optimal as those obtained for the smaller design boxes. A solution could be to increase the design domain during the optimization iteration procedure. It should be noted that increasing the box all the way up to the capacitor plates results in solution shapes where some of the dielectric material is attached to the electrode in order to extend the electrode into the capacitor and thereby maximize the electric field locally. This may be a desirable feature for some purposes. In this work we have deliberately avoided such solutions by keeping the edges of the design domain from the capacitor plates.
Throughout the paper we have only presented results obtained for dielectric solid shapes forced to be symmetric around the center plane z = 0. However, we have performed both topology optimization and shape optimization of systems not restricted to this symmetry. In general we find that the symmetric shapes are always good candidates for the optimal design. It cannot be excluded, though, that in some cases a spontaneous symmetry breaking occurs similar to the asymmetric S-turn channel studied in [7].
By studying the optimized shapes of dielectric solids, we have noted that pointed features often occur, such as those clearly seen on the dielectric satellite island in figure 8(b). The reason for these to appear seems to be that the pointed regions of the dielectric surfaces can support large gradients in the electric potential and, associated with this, also large charge densities [21,22]. As a result large electric body forces act on the electrolyte in these regions.
At the same time the surface between the pointed features curves inward, which lowers the viscous flow resistance due to the no-slip boundary condition. This effect is similar to that obtained by creating electrode arrays of different heights in ac electro-osmosis [23,24].
Another noteworthy aspect of topology-optimized structures is that the appearance of dielectric satellite islands seems to break up flow rolls that would otherwise be present and not contribute to the flow rate. This leads to a larger net flow rate, as can be seen by comparing figures 5 and 8.
Throughout the paper we have treated the design field γ as an artificial field. However, the design-field model could perhaps achieve physical applications to systems containing ionexchange membranes, as briefly mentioned in section 1. Such membranes are indeed porous structures permeated by an electrolyte.
In conclusion, our analysis points out the great potential for improving actual ICEO-based devices by changing simply connected topologies and simple shapes of the dielectric solids into multiply connected topologies of complex shapes.