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

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.


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 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.

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 x-y 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:   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 x-y 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.
Micromachines 2023, 14, x FOR PEER REVIEW 2 of 20 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.

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 x-y 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: 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. 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.

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.

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 (x f (i) , y f (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 (x e (i) , y e (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. r 1 , r 2 , r 3 , and r 4 are the distances between the contour point and the fictitious charge and image fictitious charges. d 1 , s 1 , s 2, 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: i j e f P is the associated potential coefficient that can be found from the fundamental The potential resulting from the superposition of discrete charges must be equal to the boundary potential φ C on the electrode surfaces: e, f 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.
Q j f 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 The application of Equations (1) and (2) leads to a system for w electrodes, each with N discrete charges: where the elements of the submatrices are where q (j) f is the value of the fictitious charge of (x f ). φ f is the potential of the f -th electrode. Equation (3) can be expressed in a simplified form: 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.

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: where 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: where the submatrices are given by After obtaining the value of the fictitious charge in each position, the potential above the electrodes can be found with the following: where and for phase 2 and 4 for phase 1 and 3 (15) where 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: and

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 where φ i (x i , y i ) 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 d 1 and δ 1 . The distance between the fictitious charge and the boundary, vertically and horizontally, are represented by s 1 and s 2 , respectively: Here, k 1 and k 2 are the assignment factors deciding the position of fictitious charges. In Table 1, n 1 and n 2 represent D/d 1 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 n 1 , n 2 , k 1 , and k 2 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.
Standard norm error The small accumulated errors confirm the accuracy of the CSM. By the appropriate arrangement of n 1 , n 2 , k 1 , and k 2 , the accumulated error can be decreased to as low as 0.02%. By keeping k 1 and k 2 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 k 1 and k 2 show that a reasonable arrangement of k 1 and k 2 can increase the accuracy of the CSM further.

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 P i are integration and observation points, respectively. Γ * l is the electrode surface in the first quadrant of the coordinate system. r 1 , r 2 , r 3 , and r 4 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]: where u(P) is the solution of ∇ 2 u = 0 for P ∈ Ω. and for Poutsideutheregion Ω and G(P, P i ) = In addition, G in Equation (27) Equation (29a-d) can be interpreted by the method of images, as shown in Figure 4. As a consequence, Equation (27)   The discretization of (30) yields  Equation (23) is the fundamental solution, where r 1 is the distance between points P and P i . Denoting the boundary of the electrodes as Γ 1 , Γ 2 , . . . , Γ w , it can be deduced from Equation (22) that On account of the relation in [29], where δ kl is the Kronecker delta, and noting that u(P) = u l (constant) on the boundary Γ l , it follows that Using Equations (24) and (26), this can be transformed into In addition, G in Equation (27) can be replaced with a modified fundamental solution G s , that is, where can be interpreted by the method of images, as shown in Figure 4. As a consequence, Equation (27) can be simplified to Equation (30) reduces the unknowns, remarkably, to 25% of those in the original Equation (27).
The discretization of (30) yields where for phase 1 and 3 (32) ∂u ∂n (j) l is assumed to be constant in each boundary element and needs to be solved.
Micromachines 2023, 14, 1347 Equation (33) can be written more concisely as follows: where U n 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 c i = 1: 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 Equation (40) can be simplified further by replacing G with G s as follows: Discretization of Equation (41) gives the following: where and n , u n , . . . , u Therefore, the x, y components of E can be calculated as follows: 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. (x 0 , y 0 ) is the coordinate of the segment midpoint, and (x i , y i ) is the coordinate of the observation point.  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 where γ is the length of the element. Moreover, by the substitution The integral of (48) is equivalent to   The position of P (x, y) can be denoted by where β = n y , δ = −n x . 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 r 1 , r 2 , r 3, and r 4 , respectively. By this notation, each of the integrals relevant to (32) can be written as where γ is the length of the element. Moreover, by the substitution The integral of (48) is equivalent to 1 2π where Equation (50) is valid for u = 0. When u = 0, the following result can be employed: When P 0 (x 0 , y 0 ) and P i (x i , y i ) coincide, we have κ = u = 0, and the corresponding singular integral can be obtained easily using Equation (52): 1 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.

Comparison of CSM, BEM, and FEM
The two above-mentioned approaches were used to calculate the electric field. Figures 6 and 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. 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): 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.

Comparison of CSM, BEM, and FEM
The two above-mentioned approaches were used to calculate the electric field.   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. 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): 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.

Comparison of CSM, BEM, and FEM
The two above-mentioned approaches were used to calculate the electric field. Figures 6 and 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. 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. n 1 and n 2 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. 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 Figures 9 and 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.    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 Figures 9 and 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 shear distribution of the electric field.

Distribution of Potential and Electric Field
The distribution of potential and electric field direction above eigh phases 1 and 2 are shown as contour lines and arrows in Figure 11a,b, from black bars on the bottom of the two figures signify each of the eight electro line directions can be used to predict the particle motion when the Coulom nates.
The electric field components Ex and Ey of eight electrodes in phase 2 and shown in Figures 12 and 13, which accurately include the edge effect electrode arrays. For the field of a large number of parallel electrodes, the metry of the field allows us to extend the solution with the same phase r simply repeating the periodic solution (in the middle range of the calculat positive and negative x-directions. In addition, the FEM has a large deviation from the CSM and BEM in the area of the shear distribution of the electric field.

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  In addition, the FEM has a large deviation from the CSM and BEM in the area of the shear distribution of the electric field.

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 Figures 12 and 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.  The electric field components E x and E y of eight electrodes in phase 2 are calculated and shown in Figures 12 and 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.   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.    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.   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.   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.

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]: 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 <F dep > as arrows.

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]: 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  

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].
where Ec is the dielectric strength of air and 6 3 10 V/m c E ≈ × ; 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

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].
where E c 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.
Micromachines 2023, 14, x FOR PEER REVIEW 18 of 20 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. Figure 16. Comparison between DEP and Coulomb force.

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.

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.