Green's function-based control-oriented modeling of electric field for dielectrophoresis
GGreen’s function-based control-oriented modeling of electricfield for dielectrophoresis
Martin Gurtner, a) Kristian Hengster-Movric, and Zdenˇek Hur´ak b) Faculty of Electrical Engineering, Czech Technical University in Prague,Czech Republic. (Dated: 24 September 2018)
In this paper, we propose a novel approach to obtaining a reliable and sim-ple mathematical model of a dielectrophoretic force for model-based feedbackmicromanipulation. Any such model is expected to sufficiently accuratelyrelate the voltages (electric potentials) applied to the electrodes to the result-ing forces exerted on microparticles at given locations in the workspace. Thismodel also has to be computationally simple enough to be used in real time asrequired by model-based feedback control. Most existing models involve solv-ing two- or three-dimensional mixed boundary value problems. As such, theyare usually analytically intractable and have to be solved numerically instead.A numerical solution is, however, infeasible in real time, hence such modelsare not suitable for feedback control. We present a novel approximation ofthe boundary value data for which a closed-form analytical solution is feasi-ble; we solve a mixed boundary value problem numerically off-line only once,and based on this solution we approximate the mixed boundary conditionsby Dirichlet boundary conditions. This way we get an approximated bound-ary value problem allowing the application of the analytical framework ofGreen’s functions. Thus obtained closed-form analytical solution is amenableto real-time use and closely matches the numerical solution of the originalexact problem.Keywords: dielectrophoresis, micromanipulation, Green’s function
I. INTRODUCTION
Since its invention by H. Pohl in the 1950s and 1960s , dielectrophoresis (DEP) hasproved an efficient tool for transportation, separation, and characterization of microparti-cles such as e.g. biological cells (see Refs. 3 and 4 for a recent survey and comprehensiveintroduction). More often than not, DEP is used to manipulate ensembles of large numbersof microparticles. However, recently some attempts were successful to use DEP in a feedbackcontrol scheme for a high accuracy noncontact manipulation of a single microparticle ;these developments can be viewed as a reopening of the topic first started in the 1990s .The technology has also boosted development in this area; there are reported CMOS chipsintegrating both actuation and sensing and thus enabling individual and independent ma-nipulation of thousands of cells . This technology has later been commercialized by SiliconBiosystems as a commercial product called DEPArray TM . These non-contact tweezers areusually based on a feedback control scheme, typically invoking an automatic visual tracking.The feedback control, in turn, requires a sufficiently accurate mathematical model of theunderlying physical phenomenon of DEP. The relationship between the voltages applied tothe microelectrodes and the DEP force exerted on a microparticle located at a given positionneeds to be evaluated periodically as the microparticle moves around the workspace. Sam-pling periods on a time-scale of few tens of milliseconds or even a few milliseconds are notunusual. The commonly used approaches to modeling DEP—which are based on numerical a) Electronic mail: [email protected] b) http://aa4cc.dce.fel.cvut.cz. a r X i v : . [ phy s i c s . c o m p - ph ] A p r solution of the corresponding Boundary Value Problem (BVP) , typically using
Finite Ele-ments Method (FEM) or Method of Moments (MOM) —are not feasible in real time. It ispossible to precompute and store these solutions in a computer memory (as reported in Ref.5) but this approach imposes stringent requirements on the volume of data stored. Thereare approaches described in the literature that provide analytical solutions ; however,they are only usable for simple electrode arrays and fixed harmonic voltage signals appliedto the electrodes while the feedback control requires the ability to change the voltages inreal time.In this paper, we propose a modeling methodology that provides a computationally simpleyet sufficiently accurate model of a DEP force for the purposes of feedback micromanipu-lation. We propose to combine numerical and analytical approaches to modeling of DEP.Existing models of DEP usually involve a numerical solution of an analytically intractable mixed Boundary Value Problem (mBVP) for the potential in the workspace. As the nu-merical solution is infeasible in real time and might be too large for storing in a computermemory, it is desirable to find a closed-form approximate analytical expression for the po-tential. To find such an expression, we solve numerically the original mBVP. Based on thisnumerical solution, we approximate the mBVP by a BVP for which the closed-form solutioncan be found by Green’s functions. Using the approximate closed-form expression for thepotential we obtain a model of DEP force that is computationally effective and requiresalmost no memory space. The numerical solution of the mBVP needs to be computedoff-line and only once. Thus the high computational burden associated with the numericalsolution is carried out off-line and the feedback control system uses only the approximateclosed-form analytical solution in real time.The paper is organized as follows. In Section II, we briefly present the commonly useddipole model of DEP and show what prevents its direct use in feedback micromanipulation.In Section III we propose a control-oriented model derived from the dipole model by Green’sfunctions. Experimental verification of the viability of the proposed model is provided inSection IV. The paper concludes with Section V where the main contributions of this paperare discussed.
II. FEEDBACK MANIPULATION BY DIELECTROPHORESIS
A great advantage of feedback manipulation by DEP, in contrast to the conventional useof DEP, is that it allows one to manipulate individual microparticles. Nevertheless, thiscomes with the cost of higher computational demands on the control system because thevoltages applied to the electrodes cannot be precomputed anymore and have to be adjustedin real time as required by the feedback loop.We can explain this with an aid of Fig. 1 depicting the scheme of feedback manipulation byDEP. The measured position of the microparticle is subtracted from the required position.The deviation is fed to a control system that calculates a DEP force needed to reduce thisdeviation—to steer the microparticle towards the required position. Therefore, the voltageson electrodes are required to generate such a force. Such voltages are then applied to theelectrodes and the generated DEP force acts on the microparticle moving it towards therequired position. The new position is then measured and the whole cycle is repeated. Thecrucial part of this algorithm is hidden in the control system where, in order to computethe voltages generating the desired DEP force, a model relating the voltages to the DEPforce has to be used in real time .Nevertheless, exact DEP models are rather complicated to use in real time. For instance,the widely used dipole model has the following form: time-averaged dielectrophoretic forceacting on a homogeneous spherical particle in a harmonic field is (cid:104) F DEP ( t ) (cid:105) = πε m r (cid:16) Re { K ( ω ) }∇| E | + 2Im { K ( ω ) }× (cid:0) E x ∇ ϕ x + E y ∇ ϕ y + E z ∇ ϕ z (cid:1) (cid:17) , (1)where ε m is the permittivity of the surrounding medium, r is the particle’s radius, K ( ω ) Required Pos. ControlSystem DEPSensor F DEP
Measured Pos. ParticlePhysical SystemTrue PositionDeviaton Voltages
FIG. 1. A block diagram of feedback manipulation by DEP. is a frequency dependent constant known as Clausius–Mossotti factor, E = [ E x , E y , E z ] isthe amplitude of the harmonic electric field Re { E e iωt } , the phase of the electric field isdenoted by ϕ a , a ∈ { x, y, z } and finally, Re {·} and Im {·} denote real and imaginary partsof a complex number, respectively. For brevity, the spatial dependence is omitted in thenotation.According to (1), to determine the DEP force due to applied voltages to electrodes, oneneeds to know the electric field E . The electric field E is given by E = −∇ φ , where thepotential φ is calculated from Laplace equation ∇ φ = 0 with mixed boundary conditions.Orienting the reference frame so that the electrodes lie in the x - y plane and manipulatedobjects are situated above it, the domain is defined by the half-space z >
0. The boundaryconditions are given by the voltages applied to the electrodes (Dirichlet boundary condition)and a zero-flux condition in the normal direction to the electrode plane in the interveningspace between the electrodes (Neumann boundary condition) . This BVP is analyticallyintractable and can be solved only approximately by numerical solvers. Since the controlalgorithm is supposed to run in real time and the calculation of the DEP force must nottake more then a few milliseconds, solving on-line this exact BVP numerically is infeasible.A partial remedy to this issue is to express explicitly the dependence of F DEP on thevoltages applied to the electrodes. By superposition, we express the net potential φ ( x, y, z )as a weighted sum of normalized contributions from individual electrodes. That is, φ ( x, y, z ) = n (cid:88) i =1 u i φ i ( x, y, z ) , (2)where n is the number of electrodes, u i serves as a scaling factor given by voltage on i thelectrode, φ i is the contribution to the net potential from the i th electrode when 1V isapplied to it while the remaining electrodes are grounded. Now, to determine the net po-tential φ ( x, y, z ), we have to solve n BVPs ( ∇ φ i = 0 , i = 1 , . . . , n ) that are still analyticallyintractable, but that do not change with the voltages applied to the electrodes.One can solve each of these BVPs numerically on a grid of points in advance, store thesolution and use it as a look-up table in real time. Nevertheless, this lookup table growsunacceptably large. As an example, let the microparticles be manipulated within an areaof size 1500 × × µ m. If we grid this area equidistantly with points separatedby 5 µ m, we obtain 300 × ×
60 = 5 , ,
000 points. Na¨ıve implementation of thisapproach would thus require to store [ E x , E y , E z ] and their relevant partial derivatives[ ∂E x ∂x , ∂E x ∂y , ∂E x ∂z , ∂E y ∂y , ∂E y ∂z ] for each point in order to evaluate F DEP and all that is only forone electrode.The volume of the data needed to be stored can be reduced by a method introduced byKharboutly et al. . They use a so-called Boundary Element Method and instead of storingdirectly the derivatives of the potential in points spread throughout all the domain, theystore precomputed surface charge density in a grid on electrodes. In real time, when it isrequired to compute the DEP force at a point, the surface charge density is numericallyintegrated to calculate the electric field and its derivatives at that point and that allowsthe computation of the DEP force. Nevertheless, for the previous case that still meansthat a large portion—depending on what extent of the electrode plane is occupied by theelectrodes—of 300 ×
300 = 90 ,
000 points have to be stored. Furthermore, the reduction inthe volume of data comes at the cost of higher computational complexity because all thestored data points are needed for evaluation of (1).In this paper, we propose a different approach. We approximate the previously mentioned,analytically intractable boundary value data so that a closed-form approximated solutioncan be found.
III. GREEN’S FUNCTIONS FOR MODELING OF DIELECTROPHORESIS
The solution of the Laplace equation in a half-space domain, which is the case here,with Dirichlet conditions only can be transformed into an integration by use of the Green’stheorem. The derivation can be found in Ref. 13 and the resulting formulas are φ ( x, y, z ) = z π ∞ (cid:90) (cid:90) −∞ h ( x (cid:48) , y (cid:48) )[( x − x (cid:48) ) + ( y − y (cid:48) ) + z ] / dx (cid:48) dy (cid:48) (3)for a 3D case and φ ( x, z ) = zπ (cid:90) ∞−∞ h ( x (cid:48) )( x − x (cid:48) ) + z dx (cid:48) (4)for the 2D case where one axis, in our case y , is redundant, meaning the electrode array hasinfinitely long electrodes along y axis. The functions h ( x, y ) and h ( x ) are Dirichlet boundaryconditions that represent values of the potential on the electrode plane, that is φ ( x, y, φ ( x, h ( x, y ) (or h ( x )) and that meansalso the potential on the electrode plane in the intermediate space between the electrodeswhere the mixed boundary conditions impose a zero normal flux. Furthermore, functions h ( x, y ) or h ( x ) have to be such that the evaluated integral (3) or (4) is expressible as aclosed-form expression containing only elementary functions; only then is the solution forthe potential applicable in real-time feedback control.To determine the values of the potential on the electrode plane in between the electrodes,various approximations of the decay of the potential away from the electrode can be foundin the literature . Nevertheless, they are all designed only for interdigitated electrodearrays. For more complex electrode array designs, they are either inapplicable or the in-tegrals (3) and (4) are analytically intractable for the approximation of the potential andhave to be solved numerically. Then, however, formulating the solution of the BVP as anintegration loses meaning since both the new and the original problem have to be solvednumerically.In this paper, instead, we numerically solve the original BVP with mixed boundaryconditions. In order to obtain h ( x, y ) (or h ( x )), we approximate this numerical solutionon the electrode plane (i.e. on the boundary of the domain) by an analytical model. Aswe require the integrals (3) and (4) to be expressible in closed form, we restrict the classof approximating models to piece-wise polynomial models in the 2D case and piece-wiseconstant models in the 3D case. Having the approximation of the potential on the electrodeplane, an approximate closed-form expression for the potential in the half-space domain isobtained by evaluating the integral (3) (or (4)). Thus, by (1) and (2), we also get a modelof DEP suitable for feedback control. Note, that the numerical solution is needed only toderive the control-oriented DEP model; the numerical solution is computed only once andoff-line. Therefore, all the heavy computational burden is carried off-line and the controlsystem uses the computationally more efficient model in real time.It is worth mentioning that since the potential—as the solution of the Laplace equation—is a harmonic function, it is infinitely differentiable . Furthermore, due to the Maximumprinciple the error of the approximated potential diminishes as one moves further awayfrom the electrodes. Thus, the accuracy of the model can be controlled by the accuracy ofthe approximation of the potential on the boundary of the domain.In the remainder of the paper, we apply the described methodology to two electrodearrays. (a)(b) X position [µm] P o t e n t i a l [ V ] FEMPolynomial approximation
FIG. 2. Interdigitated electrode array: (a) a side-view diagram and (b) an approximation of theboundary conditions. The black rectangles represent electrodes and the shaded rectangles on theleft and right represent the possibly added grounding plates.
A. Example 1: Interdigitated Electrode Array
Let us consider an electrode array with (2 n + 1) electrodes, the single electrode width b and center-to-center distance between the electrodes d (see Fig. 2(a)). We assume that theelectrodes are infinitely long and thus drop the dependence on the y coordinate altogether.We decompose the net potential by (2) into a weighted sum of normalized contributions φ i from individual electrodes. Furthermore, we assume that the potential contribution φ i is shift-invariant for a shift d along the x -axis. That means the potential contribution φ i isidentical for each electrode up to a shift. Mathematically, φ i ( x, z ) = φ i ± ( x ± d, z ) , i = − n + 1 , . . . , n − . (5)Clearly, this assumption does not hold for the electrodes close to the perimeter of theelectrode array. For instance, the potential contribution φ n ( x, z ) (i.e. from the electrodeon the perimeter) decays more quickly towards the ( n − φ i ( x, z ), e.g. φ ( x, z ). The remaining potential contributions are determined simply byshifting, that is φ i ( x, z ) = φ ( x + id, z ) with i ∈ {− n, . . . , n } , and the net potential is thengiven by (2). Nevertheless, to compute (4) for φ ( x, z ), we need to know h ( x ) := φ ( x, φ ( x,
0) are known only on the electrodes and unknown on the rest of thebottom boundary. To overcome this problem, we solve the original Laplace equation withmixed boundary conditions numerically by
Finite Element Method (FEM) in COMSOLMultiphysics. Then, we take the values of the potential on the bottom boundary (i.e. z = 0) between the 0-th electrode and its left adjacent electrode and fit a polynomial p ( x )to these values. The polynomial approximation of h ( x ) is˜ h ( x ) = p ( x ) x ∈ (cid:2) − d − b , − b (cid:1) , x ∈ (cid:2) − b , b (cid:3) ,p ( − x ) x ∈ (cid:0) b , d + b (cid:3) , . (6)Specifically, for an electrode array with parameters d = 2 b = 200 µ m, we fitted a third-order polynomial to the FEM solution and the fitted polynomial is p ( x ) = 1 . × − x + 4 . × − x + 5 . × − x + 2 . , (7)where x is in micrometers. X position [µm] -1000 -500 0 500 1000 F x [ N ] ×10 -10 -101 0.0 V 4.0 V 7.8 V 9.2 V 0.9 V -7 2 V -7.0 V -4.8 V 6.8 V -4.9 V 0.0 V Levitation height 120 [µm]
X position [µm] -1000 -500 0 500 1000 F z [ N ] ×10 -10 -101234 0.0 V 4.0 V 7.8 V 9.2 V 0.9 V -7 2 V -7.0 V -4.8 V 6.8 V -4.9 V 0.0 V Closed-formFEM
FIG. 3. A component-wise comparison of DEP force fields computed for interdigitated electrodearray ( d = 2 b = 200 µ m) numerically by FEM and analytically based on the approximate closed-form solution for the potential. The force fields are computed for the height of 120 µ m above theelectrode array. The applied potentials are indicated above the electrodes represented by the blackrectangles. (Multimedia view) The FEM solution together with the polynomial approximation is displayed in Fig. 2(b).Apparently, the polynomial approximation describes the FEM solution very accurately.However, one should not overlook the small humps in the gaps between the electrodes,which are completely omitted by the approximation.With ˜ h ( x ) approximating φ ( x,
0) we can compute the integral (4) and obtain an approx-imate closed-form solution for φ ( x, z ). Then, by (5) and (2) we get the net potential φ ( x, z )for any choice of electrode potentials u i . Thus, we can compute the electric field intensityand the pertaining DEP force. We do not present the evaluated integral here, because it israther lengthy and it would not serve any purpose. Nevertheless, it is crucial to mentionthat the evaluated integral is indeed expressible as a closed-form expression containing onlyelementary functions and thus it is easily applicable in real time.To validate the proposed model, we compare DEP force fields computed by (1) for thepotential obtained by numerical solution of the original BVP with mixed boundary condi-tions and for the potential obtained by solution of the approximated BVP. The comparisonis shown in Fig. 3 (Multimedia view) . The comparison is carried out for an electrodearray with nine electrodes ( n = 4), with grounded perimeter (also visible in the figure)and with the single electrode width b = 100 µ m and the distance between the electrodes d = 200 µ m. The remaining parameters are: r = 25 µ m, ε m = 7 . · − F/m and K ( ω ) = − . − . i . As usual for the standing wave DEP, the harmonic signal ap-plied on all electrodes has the same frequency and phase. Since it is rather inconvenientto compare the vector fields, the comparison is done for one particular height above theelectrodes (120 µ m) and for varying potentials on the electrodes. B. Example 2: Four-Leaf Clover Electrode Array
In the second example, we show a related approach how to approximate the values of thepotential on the electrode plane for a more complex electrode array shown in Fig. 4(a). Itconsists of four quadrants and allows manipulation of microparticles in all three directionsabove the electrode array, as it was experimentally verified . The width of the electrodesis b = 50 µ m and the center-to-center distance between the electrodes is d = 100 µ m. Weassume that the electrodes extend to infinity at one end.Along similar lines as in Example 1, based on the superposition principle (2) we expressthe net potential as a weighted sum of normalized contributions φ i from individual electrodes X position [ m] -600 -400 -200 0 200 400 600 Y po s i t i on [ m ] -600-400-2000200400600
36 25 µ µ P o t e n t i a l [ V ] (a)(b) X position [µm] Y position [µm]
X position [µm] Y po s i t i on [ µ m ] FIG. 4. Four-leaf clover electrode array: (a) a top-view diagram and (b) FEM solution for thenormalized potential contribution from one electrode.
Y pos. [µm]X pos. [µm] Y pos. [µm]X pos. [µm] Y pos. [µm]X pos. [µm] (b)(a) (c)
FIG. 5. An approximation of the desired shape of the bottom boundary condition ( z = 0) for thefour-leaf clover electrode array: (a) one block approximation, (b) staircase approximation and (c)boundary condition obtained by FEM. Black rectangles represent electrodes. and make a similar assumption that the normalized contributions are identical up to a shiftand/or rotation with respect to each other; for instance, for electrodes with indexes rangingfrom 1 to 7 (the indexes are shown in Fig. 4(a)), it holds that φ i +1 ( x, y, z ) = φ i ( x − d, y + d, z ) , (8)where d is the center-to-center distance between the electrodes.Again, thanks to this assumption, we need to compute the integral (3) only for one φ i ( x, y, z ). In this case, we choose φ ( x, y, z ) because it is close to the center of the sectorand, analogously to the previous example, the assumption (8) holds best for the electrodesin the center.It remains to find the approximation of h ( x, y ) := φ ( x, y, φ ( x, y,
0) shown in Fig. 4(b). This time, however, theintegral (3) is analytically intractable for h ( x, y ) being polynomial or even linear—we areunable to express the integral in closed form for anything other than for constant boundary Y po s i t i on X position
F - FEM x Y po s i t i on Force [N]
F - FEM y F - FEM z F - Closed-form x F - Closed-form y X position X positionX positionX position X position Y po s i t i on Y po s i t i on Y po s i t i on Y po s i t i on Force [N] Force [N]
F - Closed-form z FIG. 6. A comparison of individual components of DEP force fields calculated for the four-leafclover electrode array with single electrode width b = 50 µ m and center-to-center distance betweenthe electrodes d = 100 µ m. The force fields are calculated numerically by FEM and analyticallybased on the approximate closed-form solution for the potential. The force fields are computedfor the height of 120 µ m above the electrode array. The numerical values inside the electrodesrepresent applied potentials in volts. (Multimedia view) condition. Thus, instead of using a polynomial or linear approximation, the desired shape ofthe potential is constructed from blocks. Initially, we approximate the boundary conditionin the roughest possible way; we assume that the potential between the electrodes dropsimmediately to zero as one moves away from the electrode (see Fig. 5(a)). Then summingthe “scaled” and “shifted” versions of this boundary condition (see Fig. 5(b)) approximatelygives the desired shape (see Fig. 5(c)).Let us begin with the one-block approximation. We define the block boundary conditionfor one-side infinitely long electrode as h ( x, y ) = (cid:40) x ≤ y ∈ (cid:2) − b , b (cid:3) , , (9)where b is the width of the electrode. The staircase approximation of the desired shape isthen obtained by ˜ h ( x, y ) = N (cid:88) i =1 α i h (cid:18) x − ( β i − b , yβ i (cid:19) , (10)where N is the number of blocks, α i determines the height of the block and β i is a scalingparameter, meaning that β i = 2 scales the block so that it is twice as wide as the originalelectrode. Notice, that we assumed that the potential decays identically along the x and y axes and thus the coefficients β i determine not only the width but simultaneously also theshift of the blocks along the x axis.Given the FEM solution for mixed boundary conditions, the coefficients α i and β i in (10) FIG. 7. Photo of the apparatus used for the experimental verification of the proposed control-oriented model of DEP. can be determined by solving the following optimization problemmin α i ,β i ∈ R ,i =1 ,...,N (cid:107) ˜ h ( x , y ) − φ FEM ( x , y, (cid:107) (11)subject to: N (cid:88) i =1 α i = 1 ,β i ∈ [1 , , i = 1 , . . . , N, where φ FEM ( x, y, z ) is the FEM solution. Note that the 2-norm above measures the sizeof a function of the real y variable, but in the numerical optimization we are only ableto consider samples of y , which is not encoded in the optimization problem statement forthe sake of simplicity. With the assumption that the potential decays identically along x and y axes, both α i and β i can be determined from a y - z cross-section of φ FEM ( x, y, x to be a negative constant value x . For instance, the red curve in Fig. 5(a)represents φ FEM ( x , y,
0) for x = − µ m. The coefficients α i have to sum up to onebecause only then is the height of the piled up blocks equal to one. We assume that d = 2 b and thus restrict the coefficients β i to the interval [1 ,
3] because then the blocks cannot benarrower than the electrode and they cannot interfere with other electrodes. Even thoughthe optimization task is not convex, it still provides very good results when given a goodinitial guess. For the initial guess, we let β i grow linearly from 1 to 3 and set α i to beproportional to ∂∂y φ FEM ( x , y, N = 10.Having the approximation of the boundary condition ˜ h ( x, y ), one can calculate the inte-gral (3) and obtain a closed-form solution for φ ( x, y, z ). Instead of using the boundarycondition ˜ h ( x, y ) composed of several blocks directly, due to linearity of the integral, wecan use the one-block boundary condition h ( x, y ), calculate the integral (3) and composethe closed-form approximation for φ ( x, y, z ) in the same way as ˜ h ( x, y ) is itself composed.This is exactly how we proceed. Substitution of h ( x, y ) into the integral (3) gives φ i ( x, y, z ) = z π b (cid:90) − b (cid:90) −∞ x − x (cid:48) ) +( y − y (cid:48) ) + z ) / dx (cid:48) dy (cid:48) . (12)Again, we do not present the evaluated integral since one can easily compute it in Math-ematica, Maple, Matlab or any another computer algebra package. Nevertheless, we em-phasize that the evaluated integral is expressible in closed form usable in real time. The0 Time [s] -600-400-2000200400600 X po s i t i on [ m ] -10-50510 P o t en t i a l [ V ] ReferenceMeasurement
FIG. 8. Experimental verification of the proposed model of DEP in an experiment where amicrosphere is steered along a reference trajectory. Only the transverse coordinate (i.e. x ) isshown. The colors at electrode locations show what potentials were applied to the electrodes at agiven time. Notice that x axis here is time and thus any vertical cut shows what potentials wereapplied to the electrodes at the corresponding time. approximate closed-form solution for φ ( x, y, z ) is then obtained as φ i ( x, y, z ) = N (cid:88) i =1 α i φ i (cid:18) x − ( β i − b , yβ i , z (cid:19) . (13)Similarly, as in the 2D case, we do not compare the potentials directly because whatwe are interested in are the DEP force fields derived from the potentials. Since from thevisualization point of view it is rather inconvenient to directly compare 3D force fields,we compare their components separately. The comparisons were carried out for the sameparameters as in Example 1 and they are shown in Fig. 6 (Multimedia view). Fig. 6 showsa comparison carried out for randomly varying potentials on the electrodes. Based on thiscomparison, the force field computed by the proposed approximate closed-form model isseen to match that computed based on the FEM solution. IV. EXPERIMENTAL VERIFICATION
To verify the applicability of the proposed DEP model, we used it in an experiment wherea 50 µ m polystyrene microsphere was manipulated by a control system with a feedback loop.The goal of the control system is to steer the microsphere along a reference trajectory. Themicrosphere was suspended in water contained in a pool above an interdigitated electrodearray with six electrodes and d = 2 b = 200 µ m (see Fig. 2(a)). A detailed descriptionof the control and measuring system can be found in Ref. 5 and in Ref. 18, respectively.Figure 7 displays a photo of the hardware setup. Figure 8 shows reference and measuredtrajectories of the microsphere. In addition, the figure also displays the potentials appliedto the electrodes in order to steer the microsphere along the reference trajectory. Thesepotentials were computed in real time by the control system based on the proposed DEPmodel. It is noteworthy, that in Ref. 5 a numerical solution taking approximately 1 GB ofmemory space was used to calculate DEP force in real time, whereas here the DEP modelis represented by a closed-form analytical expression that takes almost no memory spaceand is computationally efficient. V. DISCUSSION
Although the presented approach to modeling of dielectrophoresis is demonstrated onthe dipole approximation of the microparticles, it can also be applied to more complex andaccurate multipole approximations; no modification would be needed because multipole1 (a) (b)
FIG. 9. Some other electrode array designs in the literature for which the proposed modelingmethodology is also applicable. approximations only require higher derivatives of the potential and our proposed approx-imation calculates an infinitely differentiable closed-form approximation of the potential.Even though we demonstrate the approach on two concrete electrode array designs, it canalso be used for other planar designs exhibiting similar symmetry; two such examples areshown in Fig. 9. Note that full analytical solutions for exact boundary conditions are notpossible in such cases.
VI. CONCLUSION
The major benefit of our approximate modeling methodology for a dielectrophoretic forcepresented here is that in comparison with the standard analytical or FEM-based (numerical)approaches it yields a mathematical model of DEP whose application is feasible in real time(e.g. in cycles of a few milliseconds or so) on a common laboratory PC, and yet the accuracyof the model is sufficient for the purposes of feedback micromanipulation. This methodologycombines numerical and analytical models so that the computational burden associated withthe calculation of the numerical solution is carried out off-line and based on this numericalsolution an approximated closed-form and easy-to-calculate—or briefly control-oriented—model is derived. This approach allows us to derive a control-oriented model for a broadercategory of electrode array designs than other approaches used in modeling of DEP.
ACKNOWLEDGMENTS H. A. Pohl, Journal of Applied Physics , 1182 (1958). H. A. Pohl and I. Hawk, Science , 647 (1966). R. Pethig, Biomicrofluidics , 022811 (2010). T. B. Jones, Engineering in Medicine and Biology Magazine, IEEE , 33 (2003). J. Zem´anek, T. Mich´alek, and Z. Hur´ak, Electrophoresis , 1451 (2015). M. Kharboutly and M. Gauthier, in
Robotics and Automation (ICRA), 2013 IEEE International Con-ference on (IEEE, 2013) pp. 1446–1451. M. Kharboutly, A. Melis, A. Bolopion, N. Chaillet, and M. Gauthier, in
Conference on the Manipulation,Manufacturing and Measurement on the Nanoscale, 3M NANO’12. (2012) pp. 6–pages. B. Edwards and N. Engheta, New Journal of Physics , 063012 (2012). K. V. Kaler and T. B. Jones, Biophysical journal , 173 (1990). N. Manaresi, A. Romani, G. Medoro, L. Altomare, A. Leonardi, M. Tartagni, and R. Guerrieri, IEEEJournal of Solid-State Circuits , 2297 (2003). H. Morgan, A. G. Izquierdo, D. Bakewell, N. G. Green, and A. Ramos, Journal of Physics D: AppliedPhysics , 1553 (2001). D. E. Chang, S. Loire, and I. Mezi´c, Journal of Physics D: Applied Physics , 3073 (2003). X. Wang, X.-B. Wang, F. F. Becker, and P. R. Gascoyne, Journal of physics D: applied physics , 1649(1996). M. P. Hughes,
Nanoelectromechanics in engineering and biology (CRC press, 2010). M. Kharboutly, M. Gauthier, and N. Chaillet, Journal of Applied Physics , 114312 (2009). L. C. Evans,
Partial Differential Equations , 2nd ed. (American Mathematical Society, Providence, R.I,2010). J. Zem´anek, J. Drs, and Z. Hur´ak, in (IEEE, 2014) pp. 19–25. M. Gurtner and J. Zem´anek, Measurement Science and Technology27