Next Article in Journal
Investigation of Cutting Force and the Material Removal Mechanism in the Ultrasonic Vibration-Assisted Scratching of 2D-SiCf/SiC Composites
Next Article in Special Issue
Red Blood Cell Partitioning Using a Microfluidic Channel with Ladder Structure
Previous Article in Journal
Investigation of Physico-Chemical Stability and Aerodynamic Properties of Novel “Nano-in-Micro” Structured Dry Powder Inhaler System
Previous Article in Special Issue
Design and Development of a Traveling Wave Ferro-Microfluidic Device and System Rig for Potential Magnetophoretic Cell Separation and Sorting in a Water-Based Ferrofluid
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Solution of the Electric Field and Dielectrophoresis Force of Electrostatic Traveling Wave System

1
Resource Geophysics Academy, Imperial College London, London SW7 2BP, UK
2
Department of Earth Science and Engineering, Imperial College London, London SW7 2AZ, UK
3
School of Electrical Engineering and Automation, Wuhan University, Wuhan 430072, China
*
Authors to whom correspondence should be addressed.
Micromachines 2023, 14(7), 1347; https://doi.org/10.3390/mi14071347
Submission received: 19 May 2023 / Revised: 21 June 2023 / Accepted: 27 June 2023 / Published: 30 June 2023
(This article belongs to the Special Issue Microfluidic Device Fabrication and Cell Manipulation)

Abstract

:
Electrostatic traveling wave (ETW) methods have shown promising performance in dust mitigation of solar panels, particle transport and separation in in situ space resource utilization, cell manipulation, and separation in biology. The ETW field distribution is required to analyze the forces applied to particles and to evaluate ETW design parameters. This study presents the numerical results of the ETW field distribution generated by a parallel electrode array using both the charge simulation method (CSM) and the boundary element method (BEM). A low accumulated error of the CSM is achieved by properly arranging the positions and numbers of contour points and fictitious charges. The BEM can avoid the inconvenience of the charge position required in the CSM. The numerical results show extremely close agreement between the CSM and BEM. For simplification, the method of images is introduced in the implementation of the CSM and BEM. Moreover, analytical formulas are obtained for the integral of Green’s function along boundary elements. For further validation, the results are cross-checked using the finite element method (FEM). It is found that discrepancies occur at the ends of the electrode array. Finally, analyses are provided of the electric field and dielectrophoretic (DEP) components. Emphasis is given to the regions close to the electrode surfaces. These results provide guidance for the fabrication of ETW systems for various applications.

1. Introduction

An electrostatic traveling wave (ETW) field can be produced by a set of electrodes, insulated from each other and connected to AC poly-phase voltage sources. Either neutral or charged fine particles brought into such a field move due to the action of electric forces, gravitational forces, and other forces related to their different physical properties. The ETW method originates from dielectrophoresis (DEP), first defined by Pohl as the interaction between non-uniform electric fields with polarized particles and liquids [1,2]. The researchers who followed have developed a generalized theory for calculating DEP and explored its various practical applications, such as the precipitation and dispersion of liquid [3,4]. Recently, this electrokinetic phenomenon has been applied in many diverse areas, such as xerographic particle transport in electrophotography, dust mitigation of solar panels, cell manipulation and separation in biology, and lunar particle transport and separation [5,6,7,8,9,10].
A theoretical understanding of particle movement in the ETW conveyer system will allow design and operating parameter evaluation. These particles are subject to many forces, including the Coulomb force, dielectrophoretic (DEP) force, gravitational force, friction, image force, and possibly fluid drag. Of these, the Coulomb and DEP forces are predominant, and their analysis requires the accurate calculation of the electric field [11,12]. The Accurate prediction of the ETW electric field distribution is, therefore, essential.

1.1. Background of the Electric Potential Problem

The ETW field is generally produced by an array of parallel electrodes with the same width and thickness, insulated from each other, and connected to an AC, poly-phase voltage source. Figure 1 shows a diagram of a typical electrode system for particle transport and separation. A four-phase rectangular wave voltage source is often used in ETW systems [13,14,15], and is used here as an example in the calculations.
The electrode length is much larger than both the electrode width and the gap between the electrodes, and the electric field model can be simplified to a 2D problem in the xy plane, as shown in Figure 2. In phases 2 and 4, Dirichlet conditions ϕ = 0 can be designated on x = 0. In phases 1 and 3, the symmetry leads to the Neumann condition: ϕ / n = 0 . Thus, the 2D electrostatic problem can be described by Laplace’s equation Δ ϕ = 0 , associated with the boundary conditions on the electrode surfaces, which vary in time according to the system function.

1.2. Methods Development

Several approximate analytical and numerical methods have been proposed for the solution of the electric field produced by this particular arrangement of electrodes. In 1996, Wang and co-workers [16] used Green’s theorem to calculate the electric field for 2D electrode arrays. Morgan and co-workers [17,18] developed Fourier series methods and the finite element method (FEM) for calculating DEP and traveling wave forces generated by interdigitated electrode arrays. Sun and co-workers [19] completed an analytical solution using the Schwarz–Christoffel mapping method without any approximation of the boundary conditions. Gauthier and co-workers [20] developed a Fourier series method to calculate DEP force with two facing electrode arrays. However, in all these papers, either or both of the following two approximations have been applied: (1) the electrode boundaries are set as a line without considering the electrode shape, and (2) a linear potential variation between electrodes is assumed. These approximations deviate from reality in certain conditions and restrict the application scope of their methods.
In this paper, we present two alternative approaches, the charge simulation method (CSM) and the boundary element method (BEM), to calculate the electric field distribution in such an electrode arrangement. A similar approach to BEM, a boundary integral solution of a potential problem, has been obtained and has high computation efficiency using the Nyström method [21]. The CSM and BEM methods have the benefit of no requirement for any approximation of the boundary conditions or the electrode shape and can accurately predict the electric field on the electrode edges and gaps efficiently. The formulas we derived can be adapted to various designs and parameters easily.

2. Theory of the Charge Simulation Method (CSM)

2.1. Basic Principle

The charge simulation method is based on the concept of discrete charges, which has proven to be very powerful and efficient for solving many electrostatic problems [22]. Masuda and Kamimura [23] used a similar substitute charge method to calculate the electric field of parallel cylindrical electrodes.
The CSM is a numerical method for the application of the Trefftz method for the solution of the boundary value problem (BVP) in an electrostatic field where the partial differential equation is satisfied perfectly while the boundary conditions are satisfied approximately [24,25].
The basis of the CSM is the use of a group of discrete, fictitious charges as particular solutions of Poisson equations, where the distributed charges on the boundary are replaced by these discrete charges arranged outside the boundary. The fictitious charges are shown as the hollow circles in the first quadrant inside the boundary, and (xf(i), yf(i)) represents the position of the i-th fictitious charge. Using the method of images, the image fictitious charges are applied in the other three quadrants, and the positions can be easily obtained according to the symmetries. In the CSM, the solution of Laplace’s equation depends on the determination of the values of these fictitious charges. The arrangement of fictitious charges and contour points is illustrated in Figure 3. The contour points are on the boundary of electrodes with known potentials, which can be used in the equations to solve the unknown fictitious charges. The contour points are shown as the symbol × on the boundary, and (xe(i), ye(i)) represents the position of the i-th contour point. Γe and Γf represent the e-th and f-th boundaries, while Γe and Γf are the image boundaries, respectively. r1, r2, r3, and r4 are the distances between the contour point and the fictitious charge and image fictitious charges. d1, s1, s2, and δ are simulation parameters that are used to determine the positions of contour points and fictitious charges, which can be referred to in the following calculations.
The potential resulting from the superposition of discrete charges must be equal to the boundary potential ϕC on the electrode surfaces:
ϕ f = e = 1 w i = 1 n P e , f i , j Q f j
where P e , f i , j is the associated potential coefficient that can be found from the fundamental solution of Poisson’s equation, and it only depends on the related position between the calculating point and the n-th fictitious charge.
P e , f i , j = 1 2 π ε 0 ln x e i x f j 2 + y e i y f j 2
Q f j is the value of the j-th fictitious charge on the f-th electrode. n is the total number of fictitious charges in an electrode. w is the total number of electrodes. Since there is a unique solution to this boundary value problem, the potential can be found unambiguously from Equation (1), where ε0 is the permittivity of air and ( x e ( i ) , y e ( i ) ) and ( x f ( j ) , y f ( j ) ) represent the positions of the i-th fictitious charge of the e-th electrode and the j-th points on the surface of the f-th electrodes.
The application of Equations (1) and (2) leads to a system for w electrodes, each with N discrete charges:
P 11 P 12 P 1 w P 21 P 22 P 2 w   ⋮    ⋮    ⋮ P w 1 P w 2 P w w   Q 1 Q 2    ⋮ Q w = ϕ 1 ϕ 2    ⋮ ϕ w
where the elements of the submatrices are
P e f = P e , f 1 , 1 P e , f 1 , 2 P e , f 1 , N P e , f 2 , 1 P e , f 2 , 2 P e , f 2 , N   ⋮   ⋮    ⋮ P e , f N , 1 P e , f N , 2 P e , f N , N
Q f = q f 1 q f 2 . . . q f N T
ϕ f = ϕ f ϕ f . . . ϕ f T
where q f ( j )   is the value of the fictitious charge of ( x f ( j ) , y f ( j ) ). ϕ f is the potential of the f-th electrode. Equation (3) can be expressed in a simplified form:
P   Q = ϕ
where P is the potential coefficient matrix, Q is the column vector of unknown charges, and ϕ is the potential of the points on the electrode boundary.

2.2. Implementation of CSM

Considering the symmetry of the electrode array and the applied voltage, the method of images can be used to simplify the analysis. Based on the center position of o in Figure 3, the boundaries are evenly distributed in the four quadrants of the coordinate, and the fictitious charges in the first quadrant are used as the origin charges. The combination of the associated potential coefficient of the fictitious line charges and its image fictitious line charges in the other three quadrants (see Figure 3) can be expressed as follows:
P e , f i , j = 1 2 π ε 0 ln r 2 r 3 r 1 r 4        for   phase   2   and   4 1 2 π ε 0 ln 1 r 1 r 2 r 3 r 4      for   phase   1   and   3
where
r 1 = x e i x f j 2 + y e i y f j 2
r 2 = x e i + x f j 2 + y e i y f j 2
r 3 = x e i + x f j 2 + y e i + y f j 2
r 4 = x e i x f j 2 + y e i + y f j 2
The two definitions of the associated potential coefficient in (8) are automatically satisfied with Dirichlet and Neumann conditions, respectively. Formula (2) is a generalized expression for each situation.
In this way, Equation (7) can be modified as follows:
P 11 P 12 P 1 , w / 2    P 21 P 22 P 2 , w / 2   ⋮   ⋮        ⋮ P w / 2 , 1 P w / 2 , 2 P w / 2 , w / 2   Q 1 Q 2    ⋮ Q w / 2 = ϕ 1 ϕ 2    ⋮ ϕ w / 2
where the submatrices are given by
P e f = P e , f 1 , 1 P e , f 1 , 2 P e , f 1 , N / 2 P e , f 2 , 1 P e , f 2 , 2 P e , f 2 , N / 2    ⋮    ⋮       ⋮ P e , f N / 2 , 1 P e , f N / 2 , 2 P e , f N / 2 , N / 2
After obtaining the value of the fictitious charge in each position, the potential above the electrodes can be found with the following:
ϕ x , y = P x , y Q
where
P x , y = P 1 x , y P 2 x , y . . . P N / 2 x , y
Q = q 1 q 2 . . . q N / 2 T
and
P n * x , y = 1 2 π ε 0 ln r 2 x , y r 3 x , y r 1 x , y r 4 x , y          for   phase   2   and   4 1 2 π ε 0 ln 1 r 1 x , y r 2 x , y r 3 x , y r 4 x , y   for   phase   1   and   3
where
r 1 x , y = x x n 2 + y y n 2
r 2 x , y = x + x n 2 + y y n 2
r 3 x , y = x + x n 2 + y + y n 2
r 4 x , y = x x n 2 + y + y n 2
n represents the n-th fictitious charge. The electric field comes from the differentiation of the potential ϕ x , y with respect to x and y:
E x x , y = P x Q
and
E y x , y = P y Q

2.3. Accuracy Evaluation

The modified standard norm error [26] is introduced for the evaluation of accuracy, which is defined by the sum of the deviation of potential at contour checkpoints, calculated from
E r r o r = i = 1 m V ϕ i x i , y i V 2 / m
where ϕi (xi, yi) is the calculated potential at the i-th checkpoint, m is the total number of checkpoints, and V is the corresponding surface potential of the electrode (i.e., applied voltage on the electrode). This accuracy criterion is more rigorous than the maximum deviation previously used [23]. The electrode configuration parameters are as previously used by Kawamoto [27]. The width of the electrode is 300 µm, the pitch is 600 µm, and the amplitude of the voltage is 800 V. The following calculations are based on the same parameters.
The accuracy of the CSM is affected by the number and placement of fictitious charges and contour points. Normally, the collocated contour points are evenly spaced along the horizontal and vertical axes with pitches of d1 and δ1. The distance between the fictitious charge and the boundary, vertically and horizontally, are represented by s1 and s2, respectively:
s 1 = k 1 d 1
s 2 = k 2 δ 1
Here, k1 and k2 are the assignment factors deciding the position of fictitious charges. In Table 1, n1 and n2 represent D/d1 and δ/δ1, respectively, the numbers of contour points along the length and width of each electrode boundary. The total number of checkpoints is 250. Different combinations of n1, n2, k1, and k2 were tested to estimate the CSM numerical calculation accuracy. Phases 3 and 4 are symmetrical to phases 1 and 2 with different polarity, so the accumulated error are the same and only shown for phases 1 and 2 in Table 1. Evaluations and simulations were implemented on a personal computer with a 2.9 GHz processor and 16 GB RAM using Wolfram Mathematica. The calculation time for the CSM to obtain the error results is also presented.
The small accumulated errors confirm the accuracy of the CSM. By the appropriate arrangement of n1, n2, k1, and k2, the accumulated error can be decreased to as low as 0.02%. By keeping k1 and k2 constant, it is found that an increase in the simulation points can improve the calculation accuracy; however, further increases may impair this, as illustrated by the last three columns. This is because overly dense points may lead to an ill-conditioned matrix and a decrease in the solution accuracy. The test results with different combinations of k1 and k2 show that a reasonable arrangement of k1 and k2 can increase the accuracy of the CSM further.

3. Theory of the Boundary Element Method (BEM)

3.1. Formulation

In the CSM, it may be challenging to arrange the fictitious charges properly, which affects the accuracy of the results. In contrast, the BEM does not rely upon fictitious charges, and only boundary discretization is required. Figure 4 illustrates the basic parameters for deploying the BEM. P and Pi are integration and observation points, respectively. Γ l * is the electrode surface in the first quadrant of the coordinate system. r1, r2, r3, and r4 are the distances between the integration point and the observation points and image fictitious charges. Using Green’s second identity with the appropriate fundamental solution, a boundary integral equation can be found as follows [28]:
c i u P i = Γ u P n G P , P i d Γ Γ G P , P i n u P d Γ
where u(P) is the solution of 2 u = 0 for P Ω . and
c i = 1    for   P inside the region Ω 1 / 2 for   P on the smooth boundary Γ 0    for   P outside u the region Ω
and
G P , P i = 1 2 π ln 1 r 1
Equation (23) is the fundamental solution, where r1 is the distance between points P and Pi. Denoting the boundary of the electrodes as Γ1, Γ2, …, Γw, it can be deduced from Equation (22) that
c i u P i = l = 1 w Γ l u P n G P , P i d Γ l = 1 w Γ l G P , P i n u P d Γ
On account of the relation in [29],
Γ l G P , P i n d Γ = 1 c i δ k l ,  for P i , P  on  the  boundary Γ k , Γ l , resp .
where δ k l is the Kronecker delta, and noting that u P = u l (constant) on the boundary Γl, it follows that
l = 1 w Γ l u P G P , P i n d Γ = l = 1 w u l 1 c i δ k l = 1 c i u k
Using Equations (24) and (26), this can be transformed into
u k = l = 1 w Γ l u n G P , P i d Γ
In addition, G in Equation (27) can be replaced with a modified fundamental solution Gs, that is,
G s P , P i = 1 2 π ln r 2 r 3 r 1 r 4      for  phase  2  and  4 1 2 π ln 1 r 1 r 2 r 3 r 4 for  phase  1  and  3
where
r 1 = x x i 2 + y y i 2
r 2 = x + x i 2 + y y i 2
r 3 = x + x i 2 + y + y i 2
r 4 = x x i 2 + y + y i 2
Equation (29a–d) can be interpreted by the method of images, as shown in Figure 4. As a consequence, Equation (27) can be simplified to
u k = l = 1 w / 2 Γ l u n G s P , P i d Γ
Equation (30) reduces the unknowns, remarkably, to 25% of those in the original Equation (27).
The discretization of (30) yields
u k i = l = 1 w / 2 j = 1 N Γ l , j u P n G s P , P i d Γ = l = 1 w j = 1 N u n l j G k , l i , j
where
G k , l i , j = 1 2 π Γ l , j ln r 2 r 3 r 1 r 4 d Γ      for  phase  2  and  4 1 2 π Γ l , j ln 1 r 1 r 2 r 3 r 4 d Γ    for  phase  1  and  3
u n l j is assumed to be constant in each boundary element and needs to be solved.
Accordingly, a linear equation system can be established with Equations (30) and (32):
G 11 G 21 G 12 G 22 . . . . G 1 , w / 2 G 2 , w / 2 . . . . . . . . . . . . . . . G w / 2 , 1 G w / 2 , 2 . . G w / 2 , w / 2 u n 1 u n 2   .   .   . u n w / 2 = s 1 s 2   .   .   . s w / 2
where
G k l = G k , l 1 , 1 G k , l 2 , 1 G k , l 1 , 2 G k , l 2 , 2 . . . . G k , l 1 , N G k , l 2 , N . . . . . . . . . . . . . . . G k , l N , 1 G k , l N , 2 . . G k , l N , N
u n l = u n l 1 u n l 2 . . . u n l N T
s l = u l u l . . . u l T
Equation (33) can be written more concisely as follows:
G U n = S
where Un is the solution vector representing the normal derivative of the potential on the electrode surfaces.
Using the obtained normal derivative u / n on the electrode surfaces, the potential of P in the region Ω can be found by Equation (24) with ci = 1:
u P i = l = 1 w Γ l u P n G P , P i d Γ l = 1 w Γ l G P , P i n u P d Γ
The second term on the RHS of (38) vanishes due to u P = u l on Γl and
Γ l G P , P i n d Γ = 0 , for   P i  inside  the  region  Ω  and  P  on  the  boundary  Γ l
Consequently, Equation (38) can be reduced to
u P i = l = 1 w Γ l u P n G P , P i d Γ
Equation (40) can be simplified further by replacing G with Gs as follows:
u P i = l = 1 w / 2 Γ l u P n G s P , P i d Γ
Discretization of Equation (41) gives the following:
u P i = G K U n
where
G K = G 1 , G 2 , , G w / 2 , G l = G l 1 , G l 2 , , G l N          l = 1 , 2 , , w / 2
with
G l j = 1 2 π Γ l , j ln r 2 r 3 r 1 r 4 d Γ ,      for  Phase  2  and  4 1 2 π Γ l , j ln 1 r 1 r 2 r 3 r 4 d Γ ,    for  Phase  1  and  3   j = 1 , 2 , , N
and
U n = u n 1 , u n 2 , , u n w / 2 T
Therefore, the x, y components of E can be calculated as follows:
E x P i = G K x i U n
and
E y P i = G K y i U n
The computation efficiency can be improved further by solving the integrals (32) analytically with the approach given in [30]. However, more convenient results can be obtained for the evaluation of the integrals. An illustration of the integral model is shown in Figure 5. Parameter s is the arclength of a linear segment, and the integral point P (x, y) moves along the linear segment. (x0, y0) is the coordinate of the segment midpoint, and (xi, yi) is the coordinate of the observation point.
The position of P (x, y) can be denoted by
x = x 0 + β s ,   y = y 0 + δ s
where β = ny, δ = −nx. It is worth mentioning that the normal vector always points to the left as s increases, and β 2 + δ 2 = 1 . The integral of (32) can be expanded as the sum form of four integrals related to r1, r2, r3, and r4, respectively. By this notation, each of the integrals relevant to (32) can be written as
1 2 π Γ ln 1 P i P d Γ = 1 2 π γ / 2 γ / 2 ln 1 x 0 + β s x i 2 + y 0 + δ s y i 2 d s
where γ is the length of the element. Moreover, by the substitution
s = t κ , κ = β x 0 x i + δ y 0 y i
The integral of (48) is equivalent to
1 2 π γ / 2 γ / 2 ln 1 x 0 + β s x i 2 + y 0 + δ s y i 2 d s = 1 4 π κ γ / 2 κ + γ / 2 ln t 2 + u 2 d t = 1 4 π t ln t 2 + u 2 2 t + 2 u arctan t / u t = κ γ / 2 t = κ + γ / 2 = 1 4 π κ ln κ + γ / 2 2 + u 2 κ γ / 2 2 + u 2 + γ 2 ln κ + γ / 2 2 + u 2 κ γ / 2 2 + u 2 2 γ + 2 u arctan κ + γ / 2 u arctan κ γ / 2 u
where
u = δ x 0 x i β y 0 y i
Equation (50) is valid for u 0 . When u = 0, the following result can be employed:
1 2 π γ / 2 γ / 2 ln 1 x 0 + β s x i 2 + y 0 + δ s y i 2 d s = 1 2 π t ln t t t = κ γ / 2 t = κ + γ / 2 = 1 2 π κ ln κ + γ / 2 κ γ / 2 + γ 2 ln κ 2 γ 2 / 4 γ
When P0(x0, y0) and Pi(xi, yi) coincide, we have κ = u = 0, and the corresponding singular integral can be obtained easily using Equation (52):
1 2 π γ / 2 γ / 2 ln 1 x 0 + β s x i 2 + y 0 + δ s y i 2 d s = γ 2 π 1 ln γ 2
Thus, the integrals of (32) can be found by proper linear combinations of the analytical solutions (50), (52), or (53). By this approach, no numerical quadrature is required, and the BEM algorithm can be implemented with very high efficiency.

3.2. Comparison of CSM, BEM, and FEM

The two above-mentioned approaches were used to calculate the electric field. Figure 6 and Figure 7 show the electric field at the height of 50 μm above an electrode surface with 8 and 16 electrodes, respectively. Clearly, the BEM and CSM produce highly consistent results, which confirms their validity. The code for the CSM can be found from the Supplementary Material.
In order to quantify the comparison for the accuracy of the CSM and BEM, the same process for calculating the error for the BEM was applied, and the results are shown in Table 2. n1 and n2 are the numbers on each long side and short side of boundary elements. It is clear that the CSM has slightly higher accuracy than the BEM. The CSM only needs 4.5 s to achieve the accuracy of 0.12% at phase, while the BEM needs 12 s. However, the CSM needs time to set up the positions of fictitious charges manually, which might be time-consuming for complex boundary conditions. The BEM does not need pre-processing, and the image method reduces unknown elements to 25% of those in the original equations, which greatly improves the computation efficiency of the BEM.
The finite element method (FEM) is frequently used to solve electrostatic problems [31,32]. The accuracy of the FEM is related closely to the quality of the calculating mesh used; at the edge of the conductors, extremely fine meshing is required, while further away from the electrodes, it can be coarser [33]. The area in the eight-electrode model was divided into two areas (blue and grey) by the rectangular boundary around the electrodes, as shown in Figure 8. A triangular mesh was applied. The element size outside the rectangular area is set as extremely fine, while the mesh size is set to be smaller than 4 × 10−6 m inside the rectangular area.
Comparisons between the electric field obtained by the CSM, BEM, and FEM for 8 and 16 electrodes are shown in Figure 9 and Figure 10. All methods are in good agreement in the central electrode region. However, the FEM results show significant deviation from the CSM and BEM on both sides, especially at the corners of the electrodes. It is clear that the CSM and BEM have higher accuracy than the FEM in solving this BVP. In addition, mesh generating is time-consuming with the FEM, and the impact of the increase in electrode number is more severe.
In addition, the FEM has a large deviation from the CSM and BEM in the area of the shear distribution of the electric field.

4. Electric Field and Dielectrophoretic Component Analysis

4.1. Distribution of Potential and Electric Field

The distribution of potential and electric field direction above eight electrodes in phases 1 and 2 are shown as contour lines and arrows in Figure 11a,b, from the BEM. The black bars on the bottom of the two figures signify each of the eight electrodes. The field line directions can be used to predict the particle motion when the Coulomb force dominates.
The electric field components Ex and Ey of eight electrodes in phase 2 are calculated and shown in Figure 12 and Figure 13, which accurately include the edge effect on the end of electrode arrays. For the field of a large number of parallel electrodes, the periodic symmetry of the field allows us to extend the solution with the same phase relationship by simply repeating the periodic solution (in the middle range of the calculated field) in the positive and negative x-directions.

4.2. Electric Fields with Different Electrode Thickness

Figure 14 compares the electric field magnitude with electrode thicknesses of 18 μm and 180 μm while maintaining the other model parameters the same. The electric field is evaluated at heights of 50 μm, 500 μm, and 1mm above the surface of the conveyor. At the height of 50 μm, the maximum value of the electric fields is similar. However, Figure 14b shows that electrodes with larger thickness have higher electric fields and that the difference is more distinct at higher altitudes. The effect of electrode thickness on the electric field is instructive for the design of ETW system electrode configurations.

4.3. Dielectrophoretic Component Analysis

The dielectrophoretic force has received increasing attention in particles or biological cell separation [34] and carbon nanotube manipulation [35]. The DEP force acts on particles with a dipole moment in a non-uniform electric field. The time-averaged force on the particle can be calculated as follows [17]:
F d e p = 1 4 v Re α E ˜ E ˜ 1 2 v Re α × E ˜ × E ˜
where E ˜ is a general complex amplitude of the electric field, * indicates a complex conjugate, v is the volume of the particle, and α is the effective polarizability related to particle and dielectric permittivity. In our case, the electric field is constant and calculated independently in each phase, so there is no phase variation, and the second term on the right side of Equation (54) is zero. Figure 15 shows the DEP potential E ˜ E ˜ as contour lines and vector direction <Fdep> as arrows.

4.4. Real Case of Using Dielectrophoresis

For the manipulation and transport of particles or biological cells, the theoretical analysis of different forces is crucial for the design and optimization of the conveyor system. For example, the charged particle can have both a Coulomb force and a dielectrophoretic force, and the two different forces may drive the particle to move in different directions. Figure 16 compares DEP and Coulomb forces acting on ballotini particles of various sizes and a charge of 10% of the saturation charge in the air. The saturation of particles is calculated by the relation of Pauthenier [36].
q m = 4 π r 2 ε 0 3 ε p ε p + 2 E c
where Ec is the dielectric strength of air and E c 3 × 10 6   V / m ; r is the particle size and εp is the relative dielectric permittivity of the particle. The real Clausius–Mossotti factor for the particle is assumed as 0.5. The position of the particle is fixed at the middle of the first electrode and half of the width of the electrode above the surface of the electrode, which has the largest DEP force in the y direction. Because the DEP force is volume related, it becomes more dominant as the particle size increases. The comparison is useful for the design of a system to decide the dominance force on the particle. Further analysis could be obtained, for example, on the trajectory simulation by the action of the two forces.

5. Conclusions

The CSM and BEM were explored for the numerical solution of the electric field for an interdigitated, rectangular electrode array with a specified thickness. These two methods are readily implemented with the method of images and can achieve high computational efficiency and accuracy. In addition, analytical solutions for the integral of Green’s function along the boundary elements are derived. These analytical formulas are beneficial to the efficient implementation of the BEM and can be utilized for the general 2D BEM of electrostatic problems.
The boundary conditions used in the CSM and BEM numerical calculation do not require simplification. Electric field and DEP component analysis with 8- and 16-electrode systems were provided, which showed accurate results near the electrode surface. These results can be used for the design of different ETW electrode configurations for particle transport and biological cell manipulations. The FEM results were compared with those of the CSM and BEM and showed differences at the electrode ends, indicating that the FEM may not be suitable for determining the electric field in systems with sharp boundaries. Moreover, the CSM and BEM are less time-consuming. Overall, the CSM and BEM are more general numerical methods for dealing with electrostatic field problems and can be adapted readily to more complex boundary conditions, such as a 3D model of the electrode array.
This accurate evaluation of the electric field could potentially benefit the analysis of particles and cells in transport and separation, as the estimation of particle trajectory is highly sensitive to the electric field.

Supplementary Materials

The code for the boundary element method can be downloaded from https://notebookarchive.org/2022-12-3piktmq.

Author Contributions

Conceptualization, Y.Y., Y.L. and J.C.; methodology, Y.Y. and Y.L.; software, Y.Y. and Y.L.; validation, Y.Y., Y.L. and K.H.; formal analysis, Y.Y. and S.S.; writing—original draft preparation, Y.Y.; writing—review and editing, Y.Y., Y.L., J.C., K.H., S.S. and Y.W.; visualization, Y.Y.; supervision, J.C., K.H. and Y.W.; project administration, J.C. and Y.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data from this work will be shared and made available upon request to the corresponding author.

Acknowledgments

The authors are grateful to the sponsors at the Resource Geophysics Academy, Imperial College London, for supporting this research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Pohl, H.A. The motion and precipitation of suspensoids in divergent electric fields. J. Appl. Phys. 1951, 22, 869–871. [Google Scholar] [CrossRef]
  2. Pohl, H.A. Some effects of nonuniform fields on dielectrics. J. Appl. Phys. 1958, 29, 1182–1188. [Google Scholar] [CrossRef]
  3. Washizu, M.; Jones, T. Multipolar dielectrophoretic force calculation. J. Electrost. 1994, 33, 187–198. [Google Scholar] [CrossRef]
  4. Wang, X.-B.; Huang, Y.; Becker, F.; Gascoyne, P. A unified theory of dielectrophoresis and travelling wave dielectrophoresis. J. Phys. D Appl. Phys. 1994, 27, 1571. [Google Scholar] [CrossRef]
  5. Taniguchi, K.; Yamamoto, H.; Nakano, Y.; Sakai, T.; Morikuni, S.; Watanabe, S. A New Technique for Measuring the Distribution of Charge-to-Mass Ration for Toner Particles with On-Line Use. J. Imaging Sci. Technol. 2003, 47, 224–228. [Google Scholar] [CrossRef]
  6. Mazumder, M.; Sharma, R.; Biris, A.S.; Zhang, J.; Calle, C.; Zahn, M. Self-Cleaning Transparent Dust Shields for Protecting Solar Panels and Other Devices. Part. Sci. Technol. 2007, 25, 5–20. [Google Scholar] [CrossRef]
  7. Wang, X.-B.; Huang, Y.; Burt, J.P.H.; Markx, G.H.; Pethig, R. Selective dielectrophoretic confinement of bioparticles in potential energy wells. J. Phys. D Appl. Phys. 1993, 26, 1278–1285. [Google Scholar] [CrossRef]
  8. Kawamoto, H.; Kato, M.; Adachi, M. Electrostatic transport of regolith particles for sample return mission from asteroids. J. Electrost. 2016, 84, 42–47. [Google Scholar] [CrossRef]
  9. Gu, J.; Zhang, G.; Wang, Q.; Wang, C.; Liu, Y.; Yao, W.; Lyu, J. Experimental study on particles directed transport by an alternating travelling-wave electrostatic field. Powder Technol. 2022, 397, 117107. [Google Scholar] [CrossRef]
  10. Hewlin, R.L.; Edwards, M.; Schultz, C. Design and Development of a Traveling Wave Ferro-Microfluidic Device and System Rig for Potential Magnetophoretic Cell Separation and Sorting in a Water-Based Ferrofluid. Micromachines 2023, 14, 889. [Google Scholar] [CrossRef]
  11. Kawamoto, H.; Seki, K.; Kuromiya, N. Mechanism of travelling-wave transport of particles. J. Phys. D Appl. Phys. 2006, 39, 1249–1256. [Google Scholar] [CrossRef]
  12. Zouaghi, A.; Zouzou, N.; Dascalescu, L. Assessment of forces acting on fine particles on a traveling-wave electric field conveyor: Application to powder manipulation. Powder Technol. 2019, 343, 375–382. [Google Scholar] [CrossRef]
  13. Sayyah, A.; Eriksen, R.S.; Horenstein, M.N.; Mazumder, M.K. Performance Analysis of Electrodynamic Screens Based on Residual Particle Size Distribution. IEEE J. Photovolt. 2017, 7, 221–229. [Google Scholar] [CrossRef]
  14. Liu, G.Q.; Marshall, J.S. Effect of particle adhesion and interactions on motion by traveling waves on an electric curtain. J. Electrost. 2010, 68, 179–189. [Google Scholar] [CrossRef]
  15. Johnson, C.E.; Srirama, P.K.; Sharma, R.; Pruessner, K.; Zhang, J.; Mazumder, M.K. Effect of particle size distribution on the performance of electrodynamic screens. In Proceedings of the 40th IAS Annual Meeting, Hong Kong, China, 2–6 October 2005; pp. 341–345. [Google Scholar]
  16. Wang, X.; Wang, X.; Becker, F.F.; Gascoyne, P.R.C. A theoretical method of electrical field analysis for dielectrophoretic electrode arrays using Green’s theorem. J. Phys. D Appl. Phys. 1996, 29, 1649–1660. [Google Scholar] [CrossRef]
  17. Green, N.G.; Ramos, A.; Morgan, H. Numerical solution of the dielectrophoretic and travelling wave forces for interdigitated electrode arrays using the finite element method. J. Electrost. 2002, 56, 235–254. [Google Scholar] [CrossRef] [Green Version]
  18. Morgan, H.; Izquierdo, A.G.; Bakewell, D.; Green, N.G.; Ramos, A. The dielectrophoretic and travelling wave forces generated by interdigitated electrode arrays: Analytical solution using Fourier series. J. Phys. D Appl. Phys. 2001, 34, 1553–1561. [Google Scholar] [CrossRef] [Green Version]
  19. Sun, T.; Morgan, H.; Green, N.G. Analytical solutions of ac electrokinetics in interdigitated electrode arrays: Electric field, dielectrophoretic and traveling-wave dielectrophoretic forces. Phys. Rev. E Cover. Stat. Nonlinear Biol. Soft Matter Phys. 2007, 76, 046610. [Google Scholar] [CrossRef] [Green Version]
  20. Gauthier, V.; Bolopion, A.; Gauthier, M. Analytical Formulation of the Electric Field Induced by Electrode Arrays: Towards Automated Dielectrophoretic Cell Sorting. Micromachines 2017, 8, 253. [Google Scholar] [CrossRef] [Green Version]
  21. Chappell, D.J. Boundary integral solution of potential problems arising in the modelling of electrified oil films. J. Integral Equ. Appl. 2015, 27, 407–430. [Google Scholar] [CrossRef] [Green Version]
  22. Singer, H.; Steinbigler, H.; Weiss, P. A charge simulation method for the calculation of high voltage fields. IEEE Trans. Power Appar. Syst. 1974, 5, 1660–1668. [Google Scholar] [CrossRef] [Green Version]
  23. Masuda, S.; Kamimura, T. Approximate methods for calculating a non-uniform travelling. J. Electrost. 1975, 1, 351–370. [Google Scholar] [CrossRef]
  24. Trefftz, E. Ein Gegensttick zum Ritzschen Verfahren. In Proceedings of the 2nd International Congress on Applied Mechanics, Zürich, Switzerland, 12–17 September 1926; pp. 131–137. [Google Scholar]
  25. Kołodziej, J.A.; Grabski, J.K. Many names of the Trefftz method. Eng. Anal. Bound. Elem. 2018, 96, 169–178. [Google Scholar] [CrossRef]
  26. Yializis, A.; Kuffel, E.; Alexander, P. An Optimized Charge Simulation Method for the Calculation of High Voltage Fields. IEEE Trans. Power Appar. Syst. 1978, PAS-97, 2434–2440. [Google Scholar] [CrossRef]
  27. Kawamoto, H.; Shirai, K. Electrostatic Transport of Lunar Soil for In Situ Resource Utilization. J. Aerosp. Eng. 2012, 25, 132–138. [Google Scholar] [CrossRef]
  28. Riley, K.F.; Hobson, M.P.; Bence, S.J. Mathematical Methods for Physics and Engineering; Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
  29. Fichtenholz, G. Calculus Tutorial; Highe Education Press: Beijing, China, 2006; Volume 1. [Google Scholar]
  30. Luo, Y.; Theodoulidis, T.; Zhou, X.; Kyrgiazoglou, A. Calculation of AC resistance of single-layer coils using boundary-element method. IET Electr. Power Appl. 2020, 15, 1–12. [Google Scholar] [CrossRef]
  31. Sayyah, A.; Horenstein, M.N.; Mazumder, M.K.; Ahmadi, G. Electrostatic force distribution on an electrodynamic screen. J. Electrost. 2016, 81, 24–36. [Google Scholar] [CrossRef] [Green Version]
  32. Zouaghi, A.; Zouzou, N. Numerical modeling of particle motion in traveling wave solar panels cleaning device. J. Electrost. 2021, 110, 103552. [Google Scholar] [CrossRef]
  33. Sayyah, A.; Horenstein, M.N.; Mazumder, M.K. A comprehensive analysis of the electric field distribution in an electrodynamic screen. J. Electrost. 2015, 76, 115–126. [Google Scholar] [CrossRef] [Green Version]
  34. Nguyen, N.-V.; Manh, T.L.; Nguyen, T.S.; Le, V.T.; Hieu, N.V. Applied electric field analysis and numerical investigations of the continuous cell separation in a dielectrophoresis-based microfluidic channel. J. Sci. Adv. Mater. Devices 2021, 6, 11–18. [Google Scholar] [CrossRef]
  35. Krupke, R.; Hennrich, F.; Löhneysen, H.V.; Kappes, M. Separation of Metallic from Semiconducting Single-Walled Carbon Nanotubes. Science 2003, 301, 344–347. [Google Scholar] [CrossRef] [Green Version]
  36. Pauthenier, M.; Moreau-Hanot, M. La charge des particules sphériques dans un champ ionisé. J. De Phys. Radium 1932, 3, 590–613. [Google Scholar] [CrossRef]
Figure 1. Diagram showing the typical application system and applied voltage.
Figure 1. Diagram showing the typical application system and applied voltage.
Micromachines 14 01347 g001
Figure 2. Diagram showing the typical application system, consisting of the interdigitated electrode array. D is the electrode width, p is the electrode pitch, and δ is the electrode thickness. The voltage on the boundaries is the example at phase 2.
Figure 2. Diagram showing the typical application system, consisting of the interdigitated electrode array. D is the electrode width, p is the electrode pitch, and δ is the electrode thickness. The voltage on the boundaries is the example at phase 2.
Micromachines 14 01347 g002
Figure 3. Allocation of contour points and fictitious charges with the method of images.
Figure 3. Allocation of contour points and fictitious charges with the method of images.
Micromachines 14 01347 g003
Figure 4. Electrode geometry for the BEM and the method of images.
Figure 4. Electrode geometry for the BEM and the method of images.
Micromachines 14 01347 g004
Figure 5. Arclength parameters for the element integration.
Figure 5. Arclength parameters for the element integration.
Micromachines 14 01347 g005
Figure 6. Comparison of the electric field at the height of 50 μm above the 8 electrodes using the BEM and CSM.
Figure 6. Comparison of the electric field at the height of 50 μm above the 8 electrodes using the BEM and CSM.
Micromachines 14 01347 g006
Figure 7. Comparison of the electric field at the height of 50 μm above the 16 electrodes using the BEM and CSM.
Figure 7. Comparison of the electric field at the height of 50 μm above the 16 electrodes using the BEM and CSM.
Micromachines 14 01347 g007
Figure 8. Subdivision area illustration.
Figure 8. Subdivision area illustration.
Micromachines 14 01347 g008
Figure 9. Comparison of the electric field at the height of 32 μm above the 8 electrodes using the CSM, BEM, and FEM.
Figure 9. Comparison of the electric field at the height of 32 μm above the 8 electrodes using the CSM, BEM, and FEM.
Micromachines 14 01347 g009
Figure 10. Comparison of the electric field at the height of 50 μm above the 16 electrodes using the CSM, BEM, and FEM.
Figure 10. Comparison of the electric field at the height of 50 μm above the 16 electrodes using the CSM, BEM, and FEM.
Micromachines 14 01347 g010
Figure 11. Contour plot of potential and vector plot of electric field above electrodes.
Figure 11. Contour plot of potential and vector plot of electric field above electrodes.
Micromachines 14 01347 g011
Figure 12. Spatial distribution of the magnitude of field component Ex.
Figure 12. Spatial distribution of the magnitude of field component Ex.
Micromachines 14 01347 g012
Figure 13. Spatial distribution of the magnitude of field component Ey.
Figure 13. Spatial distribution of the magnitude of field component Ey.
Micromachines 14 01347 g013
Figure 14. Comparison of the electric field distribution with 18 μm and 180 μm electrode thicknesses (a) at the height of 50 μm above the surface of the conveyor and (b) at the height of 500 μm and 1 mm above the surface of the conveyor.
Figure 14. Comparison of the electric field distribution with 18 μm and 180 μm electrode thicknesses (a) at the height of 50 μm above the surface of the conveyor and (b) at the height of 500 μm and 1 mm above the surface of the conveyor.
Micromachines 14 01347 g014
Figure 15. Contour plot of DEP potential and vector plot of DEP above electrodes in phase 2.
Figure 15. Contour plot of DEP potential and vector plot of DEP above electrodes in phase 2.
Micromachines 14 01347 g015
Figure 16. Comparison between DEP and Coulomb force.
Figure 16. Comparison between DEP and Coulomb force.
Micromachines 14 01347 g016
Table 1. Standard norm error on each electrode in two different phases with a varying number of calculating points for the CSM.
Table 1. Standard norm error on each electrode in two different phases with a varying number of calculating points for the CSM.
Calculating parametersk1 = 1/6
k2 = 1/6
n1 = 200
n2 = 10
n1 = 150
n2 = 10
n1 = 100
n2 = 10
n1 = 100
n2 = 15
n1 = 100
n2 = 20
Standard norm errorPhase 10.03%0.04%0.07%0.08%0.12%
Phase 20.02%0.03%0.05%0.03%0.09%
Time 11.6 s7.6 s3.9 s4.1 s4.5 s
Calculating parametersn1 = 200
n2 = 10
k1 = 1/6
k2 = 1/6
k1 = 1/5
k2 = 1/5
k1 = 1/3
k2 = 1/3
k1 = 1/2
k2 = 1/2
k1 = 1/3
k2 = 1/6
Standard norm errorPhase 10.03%0.03%0.04%0.06%0.06%
Phase 20.02%0.02%0.03%0.05%0.04%
Time: 12.0 s12.1 s12.3 s12.1 s12.1 s
Table 2. Standard norm error on each electrode in two different phases with a varying number of calculating points for the BEM.
Table 2. Standard norm error on each electrode in two different phases with a varying number of calculating points for the BEM.
Calculating parameters n1 = 200
n2 = 10
n1 = 150
n2 = 10
n1 = 100
n2 = 10
n1 = 100
n2 = 15
n1 = 100
n2 = 20
Accumulated errorPhase 10.13%0.16%0.19%0.19%0.18%
Phase 20.10%0.12%0.15%0.14%0.13%
Time 12.0 s8.1 s4.9 s5.3 s5.4 s
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yu, Y.; Luo, Y.; Cilliers, J.; Hadler, K.; Starr, S.; Wang, Y. Numerical Solution of the Electric Field and Dielectrophoresis Force of Electrostatic Traveling Wave System. Micromachines 2023, 14, 1347. https://doi.org/10.3390/mi14071347

AMA Style

Yu Y, Luo Y, Cilliers J, Hadler K, Starr S, Wang Y. Numerical Solution of the Electric Field and Dielectrophoresis Force of Electrostatic Traveling Wave System. Micromachines. 2023; 14(7):1347. https://doi.org/10.3390/mi14071347

Chicago/Turabian Style

Yu, Yue, Yao Luo, Jan Cilliers, Kathryn Hadler, Stanley Starr, and Yanghua Wang. 2023. "Numerical Solution of the Electric Field and Dielectrophoresis Force of Electrostatic Traveling Wave System" Micromachines 14, no. 7: 1347. https://doi.org/10.3390/mi14071347

APA Style

Yu, Y., Luo, Y., Cilliers, J., Hadler, K., Starr, S., & Wang, Y. (2023). Numerical Solution of the Electric Field and Dielectrophoresis Force of Electrostatic Traveling Wave System. Micromachines, 14(7), 1347. https://doi.org/10.3390/mi14071347

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop