Three-Dimensional Interaction of a Large Number of Dense DEP Particles on a Plane Perpendicular to an AC Electrical Field ^†

Abstract

The interaction of dielectrophoresis (DEP) particles in an electric field has been observed in many experiments, known as the “particle chains phenomenon”. However, the study in 3D models (spherical particles) is rarely reported due to its complexity and significant computational cost. In this paper, we employed the iterative dipole moment (IDM) method to study the 3D interaction of a large number of dense DEP particles randomly distributed on a plane perpendicular to a uniform alternating current (AC) electric field in a bounded or unbounded space. The numerical results indicated that the particles cannot move out of the initial plane. The similar particles (either all positive or all negative DEP particles) always repelled each other, and did not form a chain. The dissimilar particles (a mixture of positive and negative DEP particles) always attracted each other, and formed particle chains consisting of alternately arranged positive and negative DEP particles. The particle chain patterns can be randomly multitudinous depending on the initial particle distribution, the electric properties of particles/fluid, the particle sizes and the number of particles. It is also found that the particle chain patterns can be effectively manipulated via tuning the frequency of the AC field and an almost uniform distribution of particles in a bounded plane chip can be achieved when all of the particles are similar, which may have potential applications in the particle manipulation of microfluidics.

1. Introduction

An electric force exerted on a non-conducting particle when a locally non-uniform electric field is applied is called dielectrophoresis (DEP) [1,2,3]. The particle is polarized in a non-uniform electric field, which arises from the spatial variation or time variation of the field (direct current (DC) or alternating current (AC)). The DEP force is dependent on the electric properties of the particle and fluid, the shape and size of the particle, and the frequency of the AC field. The particle is forced to move toward higher or lower electric field regions according to the difference in the electric properties of the particle and electrolyte, which are called positive and negative DEP, respectively. DEP can be easily applied to manipulate the particles and has become one of the most popular techniques in applications in biology and chemistry, such as cell separation [4,5,6], sorting [7], trapping [8,9], and assembly [10,11,12].

The electric polarization of a particle can also disturb the local electric field of the surrounding particles. When one particle is close to another, the local electric field around them becomes more non-uniform, even under a uniformly-applied field and, thus, this induces a DEP interaction force for each particle. The interaction force gets larger as the particles get closer. This electrical phenomenon is called a “particle chain”, which has been observed in other experiments [13] and is the basis of the particle assembly technique, which has been used in many applications [14,15,16,17]. Kevin et al. developed a new class of microwires by dielectrophoresis assembly from suspensions of metallic nanoparticles [14]. Shalini et al. co-assembled live cells and surface-functionalized synthetic microparticles on electrically controlled chips to yield permanent chains and one-cell-layer-thick membranes [15]. Zhou et al. used dielectrophoresis assembly to rapidly align ensembles of CdSe semiconductor nanowires near patterned microelectrodes for polarization-sensitive photodetection and biosensing applications [16]. Julia et al. developed a new microfluidic test system for liver toxicity, devising a microfluidic chip featuring cell culture chambers with integrated electrodes for the assembly of liver sinusoids by dielectrophoresis [17]. In recent years, a significant amount of effort has been devoted to the study of the DEP interactions of particles. Giner et al. [13] studied the characteristics of the particle chaining process in a mixture of polystyrene beads and yeast cells, applying a Monte Carlo simulation. Aubry et al. [18,19] developed a Lagrange multiplier–based numerical model to study the DEP interactions of particles. The DEP force is calculated by a simplified dipole method, which is sufficiently accurate when the particles are not too close. Ai et al. [20,21] and Xie et al. [22,23] performed direct numerical simulations on the interactions of two or three particles both in two-dimensional DC and AC fields by using the Maxwell stress tensor (MST) method and the arbitrary Lagrangian–Eulerian (ALE) method, and the effect of the frequency was studied. The numerical computation of the ALE method requires continuous grids moving and remeshing in the case of a larger number of particles. The MST method is the most rigorous approach to calculate the DEP forces. However, it is not practical in cases of 3D DEP interactions of a larger number of particles due to the numerical errors and the significant cost of computation. Kang et al. [24,25,26,27] performed direct numerical simulations to study the 2D interactions of five particles in both DC and AC fields by finite difference simulations. Hossan et al. [28] studied 2D interactions of a few particles in a DC field using a hybrid immersed interface-immersed boundary method. The above studies were all for a small number of 2D cylindrical particles in an electric field parallel to the particle plane. These results have deviated somewhat in comparison with the experimental results [29,30,31,32]. The DEP interactions of a large number of particles in a 3D electric field have not been well studied due to great difficulties in numerical computation.

Recently, the iterative dipole moment (IDM) method, first developed by Liu and Wu et al. [33,34,35,36], has been an effective method to simulate the DEP interactions of a large number of particles. IDM is an algebraic algorithm without complicated numerical computation for solving the Laplace equation and the Navier-Stokes equation. The accuracy, convergence, and simplicity of the IDM method have been demonstrated by comparing them with the results of the MST method in both DC and AC fields [33,34]. The motion trajectories and final chains of the particles are in good agreement with the results of the previous numerical studies and consistent with the experiments. This method has been successfully applied in the interaction of a large number of particles in 2D DC and AC fields [35,36].

In this work, the IDM model is used to simulate the DEP interactions of 3D spherical particles in a special 3D scenario where the particles are initially positioned in a plane. A uniform AC electric field is applied perpendicular to the particle plane, which is different to the prior 2D studies. No correlational numerical studies have been found so far. The phenomena of the interactions in this scenario are novel and notable, and they are distinguished from the results of previous studies. The mechanism of particle motions and the particle chains of DEP interactions of a large number of dense particles are investigated in situations of different sizes, different initial distributions, and different frequencies. Furthermore, the effect of a circular boundary to the DEP particle interactions is also investigated.

2. Mathematical Model

2.1. Physical Description

In practical microfluidic chips the depth is much smaller than the horizontal scale. The particle motion in the chip is almost confined to a plane in cases where the particle size is of the same order as the chip depth. The 2D DEP particle interactions on a plane parallel to the electric field have been well studied. In the present study, the 3D spherical particles are randomly distributed on a plane chip, and a uniform AC electric field perpendicular to the plane is applied as shown in Figure 1.

2.2. The Iterative Dipole Moment Method (IDM)

The particles are located at $(x_i,y_i,z_i)(i=1,2,\cdots N)$, respectively. When an AC electric field ${\tilde{E}}_0$ is applied, the particles can be polarized as dipole moments ${\tilde{p}}_i$, expressed as [2]

$${\tilde{p}}_i=4\mathsf{\pi }{a}_i^3{\mathsf{\epsilon }}_f{K}_i(\mathsf{\omega }){\tilde{E}}_0$$

(1)

where $a_i$ is the radius of the ith particle, $K_i(\mathsf{\omega })$ is the Clausius–Mossotti (C–M) factor, given by

$$K_i(\mathsf{\omega })=\frac{{\tilde{\mathsf{\epsilon }}}_{pi}-{\tilde{\mathsf{\epsilon }}}_f}{{\tilde{\mathsf{\epsilon }}}_{pi}+2{\tilde{\mathsf{\epsilon }}}_f}$$

(2)

where ${\tilde{\mathsf{\epsilon }}}_{pi}={\mathsf{\epsilon }}_{pi}-j{\mathsf{\sigma }}_{pi}/\mathsf{\omega }$ and ${\tilde{\mathsf{\epsilon }}}_f={\mathsf{\epsilon }}_f-j{\mathsf{\sigma }}_f/\mathsf{\omega }$ are the complex permittivities of the ith particle and fluid, respectively; $j=\sqrt{-1}$ and the superscript “~” denote the complex variables; ${\mathsf{\epsilon }}_{pi}$ and ${\mathsf{\epsilon }}_f$ are the permittivity of the ith particle and fluid, respectively; ${\mathsf{\sigma }}_{pi}$ and ${\mathsf{\sigma }}_f$ are the conductivities of the ith particle and fluid, respectively; $\mathsf{\omega }=2\mathsf{\pi }f$ is the angular frequency of the AC electric field.

The dipole moment will induce a local electric field, called the dipole-induced field in this work, written as follows [2]

$${\tilde{\mathsf{\phi }}}_i=\frac{{\tilde{p}}_i\cdot r_i}{4{\mathsf{\pi }\mathsf{\epsilon }}_f{r}_i^3}=\frac{a_i^3{K}_i(\mathsf{\omega })\cdot ({\tilde{E}}_{0x}(x-x_i)+{\tilde{E}}_{0y}(y-y_i)+{\tilde{E}}_{0z}(z-z_i))}{r_i^3}$$

(3)

where $r_i(r_i=\left|r_i\right|)$ is the position vector originating from the ith particle location $(x_i,y_i,{\mathrm{z}}_i)$ to any point $(x,y,\mathrm{z})$. When the particles are close to each other, the electric field can be strikingly influenced by the local electric fields. In our work, we correct the electric field with the dipole-induced fields of the neighboring particles, written as follows:

$${\tilde{E}}_j^{(1)}={\tilde{E}}_{0j}+{\displaystyle \sum _{i=1,i\ne j}^N{\tilde{E}}_{ij}^{(0)}}$$

(4)

where ${\tilde{E}}_{0j}$ is the original electric field, ${\tilde{E}}_{ij}^{(0)}$ is the dipole-induced field of the neighboring particles, and the components can be expressed as

$$\begin{array}{l}{({\tilde{E}}_x)}_{ij}=-\frac{\partial {\tilde{\mathsf{\phi }}}_i}{\partial x}=-a_i^3{K}_i(\mathsf{\omega })\\ \times \left[\frac{{\tilde{E}}_{0x}}{r_{ij}^3}-\frac{3({\tilde{E}}_{0x}(x_j-x_i)+{\tilde{E}}_{0y}(y_j-y_i)+{\tilde{E}}_{0z}(z_j-z_i))(x_j-x_i)}{r_{ij}^5}\right]\\ {({\tilde{E}}_y)}_{ij}=-\frac{\partial {\tilde{\mathsf{\phi }}}_i}{\partial y}=-a_i^3{K}_i(\mathsf{\omega })\\ \times \left[\frac{{\tilde{E}}_{0y}}{r_{ij}^3}-\frac{3({\tilde{E}}_{0x}(x_j-x_i)+{\tilde{E}}_{0y}(y_j-y_i)+{\tilde{E}}_{0z}(z_j-z_i))(y_j-y_i)}{r_{ij}^5}\right]\\ {({\tilde{E}}_z)}_{ij}=-\frac{\partial {\tilde{\mathsf{\phi }}}_i}{\partial z}=-a_i^3{K}_i(\mathsf{\omega })\\ \times \left[\frac{{\tilde{E}}_{0z}}{r_{ij}^3}-\frac{3({\tilde{E}}_{0x}(x_j-x_i)+{\tilde{E}}_{0y}(y_j-y_i)+{\tilde{E}}_{0z}(z_j-z_i))(z_j-z_i)}{r_{ij}^5}\right]\end{array}$$

(5)

where $r_{ij}=\left|r_{ij}\right|=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+{(z_j-z_i)}^2}$, $r_{ij}$ is pointed from ith particle to the jth particle.

The corrected field induces a new dipole ${\tilde{p}}_i$ and a new dipole-induced electric field ${\tilde{E}}_{ij}^{(1)}$, and then a second corrected electric field ${\tilde{E}}_j^{(2)}$ is obtained in the same way. The local electric field, dipole moment and the dipole-induced electric field will be iteratively corrected, until a converged electric field ${\tilde{E}}_j^{(n)}$ is obtained by

$$\left|\frac{{\tilde{E}}_{xj}^{(n+1)}-{\tilde{E}}_{xj}^{(n)}}{{\tilde{E}}_{xj}^{(n+1)}}\right|\le \mathsf{\tau },\left|\frac{{\tilde{E}}_{yj}^{(n+1)}-{\tilde{E}}_{yj}^{(n)}}{{\tilde{E}}_{yj}^{(n+1)}}\right|\le \mathsf{\tau },\left|\frac{{\tilde{E}}_{zj}^{(n+1)}-{\tilde{E}}_{zj}^{(n)}}{{\tilde{E}}_{zj}^{(n+1)}}\right|\le \mathsf{\tau }$$

(6)

where τ denotes the convergence criteria of the iteration error.

When the electric field is corrected, the time-averaged DEP interaction force $F_{\mathrm{DEP}}$ can be obtained by [34]

$$\langle F_{\mathrm{DEP}}\rangle =2\mathsf{\pi }{a}^3{\mathsf{\epsilon }}_f\mathrm{Re}[K(\mathsf{\omega })]\nabla {\left|\tilde{E}\right|}^2$$

(7)

2.3. The Modified Stokes Formula of a Large Number of Dense Particles

When a particle is suspended in a still fluid, and an electric field is applied to drive the particle by dielectrophretic (DEP) force, a hydrodynamic force is passively generated on the particle due to the particle motion. Generally, the fluid flow is induced by the particle motion under a low Reynolds number ($\mathrm{Re}~{10}^{-3}$), and the traditional Stokes formula has been widely accepted to calculate the hydrodynamic force on a particle, written as

$$F_d=6\mathsf{\pi }a\mathsf{\eta }(u-v)$$

(8)

where η is the viscosity of the fluid, $u$ and $v$ are the velocities of the fluid and the particle, respectively. However, when a large number of particles are densely distributed in the fluid and move close to each other, the near-field hydrodynamic hindrance must be taken into consideration. In this work, a modified Stokes formula is proposed to consider the hydrodynamic aspects.

When a spherical particle moves in a fluid, the stream function and the velocity distribution of the fluid around the particle in local Cartesian coordinates $(\overline{x},\overline{y},\overline{z})$ located in the particle center can be expressed as [37]

$$\begin{array}{l}{\overline{u}}_x=U(-\frac{3}{4}\frac{a}{r}+\frac{3a^3}{4r^3})\frac{\overline{x}\overline{z}}{r^2}\\ {\overline{u}}_y=U(-\frac{3}{4}\frac{a}{r}+\frac{3a^3}{4r^3})\frac{\overline{y}\overline{z}}{r^2}\\ {\overline{u}}_z=U(1-\frac{3}{2}\frac{a}{r}+\frac{a^3}{2r^3}+\frac{3}{4}\frac{a}{r}\frac{{\overline{x}}^2+{\overline{y}}^2}{r^2}-\frac{3a^3}{4r^3}\frac{{\overline{x}}^2+{\overline{y}}^2}{r^2})\end{array}$$

(9)

where $U$ is the magnitude of the relation velocity of the fluid and the particle, and the direction of the relative velocity is the same as $\overline{z}$ axis and $r=\sqrt{{\overline{x}}^2+{\overline{y}}^2+{\overline{z}}^2}$.

Suppose that a large number of particles with global Cartesian coordinates $(x_i,y_i,z_i)$ are randomly suspended in the fluid, and the relative velocity of the fluid and the particles is $u_{ri}=u_i-u_{pi}$, $u_{ri}=\left|u_{ri}\right|$. Two angles are introduced by $\sin \mathsf{\alpha }=\sqrt{u_{rix}^2+u_{riy}^2}/u_{ri}$, $\cos \mathsf{\alpha }=u_{riz}/u_{ri}$, $\sin \mathsf{\beta }=u_{riy}/\sqrt{u_{rix}^2+u_{riy}^2}$, $\cos \mathsf{\beta }=u_{rix}/\sqrt{u_{rix}^2+u_{riy}^2}$. The flow velocity at particle j disturbed by particle i in local Cartesian coordinates can be written as

$$\begin{array}{l}{\overline{u}}_{\overline{x}ij}=u_{ri}(-\frac{3}{4}\frac{a_i}{r_{ij}}+\frac{3a_i^3}{4{r_{ij}}^3})\frac{{\overline{x}}_{ij}{\overline{z}}_{ij}}{{r_{ij}}^2}\\ {\overline{u}}_{\overline{y}ij}=u_{ri}(-\frac{3}{4}\frac{a_i}{r_{ij}}+\frac{3a_i^3}{4{r_{ij}}^3})\frac{{\overline{y}}_{ij}{\overline{z}}_{ij}}{{r_{ij}}^2}\\ {\overline{u}}_{\overline{z}ij}=u_{ri}(1-\frac{3}{2}\frac{a_i}{r_{ij}}+\frac{a_i^3}{2{r_{ij}}^3}+\frac{3}{4}\frac{a_i}{r_{ij}}\frac{{{\overline{x}}_{ij}}^2+{{\overline{y}}_{ij}}^2}{{r_{ij}}^2}-\frac{3a_i^3}{4{r_{ij}}^3}\frac{{{\overline{x}}_{ij}}^2+{{\overline{y}}_{ij}}^2}{{r_{ij}}^2})\end{array}$$

(10)

where $({\overline{x}}_{ij},{\overline{y}}_{ij},{\overline{z}}_{ij})$ is the coordinate of particle j with respect to particle i and $r_{ij}=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+{(z_j-z_i)}^2}$. The net disturbed velocity can be obtained by

$$\begin{array}{l}\Delta {\overline{u}}_{\overline{x}ij}={\overline{u}}_{\overline{x}ij}-u_{rix}=u_{ri}(-\frac{3}{4}\frac{a_i}{r_{ij}}+\frac{3a_i^3}{4{r_{ij}}^3})\frac{{\overline{x}}_{ij}{\overline{z}}_{ij}}{{r_{ij}}^2}\\ \Delta {\overline{u}}_{\overline{y}ij}={\overline{u}}_{\overline{y}ij}-u_{riy}=u_{ri}(-\frac{3}{4}\frac{a_i}{r_{ij}}+\frac{3a_i^3}{4{r_{ij}}^3})\frac{{\overline{y}}_{ij}{\overline{z}}_{ij}}{{r_{ij}}^2}\\ \Delta {\overline{u}}_{\overline{z}ij}={\overline{u}}_{\overline{z}ij}-u_{riz}=u_{ri}(-\frac{3}{2}\frac{a_i}{r_{ij}}+\frac{a_i^3}{2{r_{ij}}^3}+\frac{3}{4}\frac{a_i}{r_{ij}}\frac{{{\overline{x}}_{ij}}^2+{{\overline{y}}_{ij}}^2}{{r_{ij}}^2}-\frac{3a_i^3}{4{r_{ij}}^3}\frac{{{\overline{x}}_{ij}}^2+{{\overline{y}}_{ij}}^2}{{r_{ij}}^2})\end{array}$$

(11)

Then we can get the net disturbed velocity at particle j by particle i in global Cartesian coordinates as

$$\begin{array}{l}\Delta u_{xij}=\cos \mathsf{\phi }(\sqrt{\Delta {\overline{u}}_{\overline{x}ij}^2+\Delta {\overline{u}}_{\overline{y}ij}^2}\cos \mathsf{\alpha }+\Delta {\overline{u}}_{\overline{z}ij}\sin \mathsf{\alpha })\\ \Delta u_{yij}=\sin \mathsf{\phi }(\sqrt{\Delta {\overline{u}}_{\overline{x}ij}^2+\Delta {\overline{u}}_{\overline{y}ij}^2}\cos \mathsf{\alpha }+\Delta {\overline{u}}_{\overline{z}ij}\sin \mathsf{\alpha })\\ \Delta u_{zij}=-\sqrt{\Delta {\overline{u}}_{\overline{x}ij}^2+\Delta {\overline{u}}_{\overline{y}ij}^2}\sin \mathsf{\alpha }+\Delta {\overline{u}}_{\overline{z}ij}\cos \mathsf{\alpha }\end{array}$$

(12)

By adding the net disturbed velocities of the surrounding particles to the flow field at particle j, the flow velocity at particle j can be corrected as

$$\begin{array}{l}{u}_{xj}^1=u_{xj}+{\displaystyle \sum _{i=1,i\ne j}^N\Delta u_{xij}^1}\\ u_{yj}^1=u_{yj}+{\displaystyle \sum _{i=1,i\ne j}^N\Delta u_{yij}^1}\\ u_{zj}^1=u_{zj}+{\displaystyle \sum _{i=1,i\ne j}^N\Delta u_{zij}^1}\end{array}$$

(13)

where $u_{xj},u_{yj},u_{zj}$ are the initial flow field without disturbance.

The net disturbed velocities and flow field will be iteratively corrected, until a converged flow field $u_{xj}^n,u_{yj}^n,u_{zj}^n$ is obtained as

$$\left|\frac{u_{xj}^n-u_{xj}^{n-1}}{u_{xj}^n}\right|\le \mathsf{\tau },\left|\frac{u_{yj}^n-u_{yj}^{n-1}}{u_{yj}^n}\right|\le \mathsf{\tau },\left|\frac{u_{zj}^n-u_{zj}^{n-1}}{u_{zj}^n}\right|\le \mathsf{\tau }$$

(14)

where $\mathsf{\tau }={10}^{-3}$ denotes the convergence criteria of the iteration error; then we can get the corrected fluid drag force by substituting the initial flow field with the corrected one in the Stokes formula as

$$\begin{array}{l}{F}_{dxi}=6\mathsf{\pi }\mathsf{\mu }{a}_i(u_{xi}^n-u_{xpi})\\ F_{dyi}=6\mathsf{\pi }\mathsf{\mu }{a}_i(u_{yi}^n-u_{ypi})\\ F_{dzi}=6\mathsf{\pi }\mathsf{\mu }{a}_i(u_{zi}^n-u_{zpi})\end{array}$$

(15)

2.4. The Governing Equation of the Particles and the Dimensionless Method

The motion of the particles is governed by the kinematic equation as follows

$$m_p\frac{d^2{r}_p}{d{t}^2}=F_{\mathrm{DEP}}+F_d$$

(16)

where $m_p$, $r_p$ are the mass and the position vector of the particles, respectively.

All of the governing equations can be normalized by the characteristic length $a_0=5\mathsf{\mu }\mathrm{m}$, potential ${\mathsf{\phi }}_0=10\mathrm{V}$, velocity $U_0={\mathsf{\epsilon }}_f{\mathsf{\phi }}_0^2/(\mathsf{\eta }{a}_0)$, time $t_0=a_0/U_0$, force $F_0=a_0\mathsf{\eta }{U}_0$, and mass $m_0=a_0^2\mathsf{\eta }/U_0$. The normalized equations are written as:

$$F_{\mathrm{DEP}}^{\ast }=2\mathsf{\pi }{a}^{\ast }\mathrm{Re}[K(\mathsf{\omega })]{\nabla }^{\ast }{\left|{\tilde{E}}^{\ast }\right|}^2$$

(17)

$$F_d^{\ast }=6\mathsf{\pi }{a}^{\ast }(u^{\ast }-v^{\ast })$$

(18)

$$m_p^{\ast }\frac{d^2{r}_p^{\ast }}{d{t}^{\ast 2}}=F_{\mathrm{DEP}}^{\ast }+F_d^{\ast }$$

(19)

where $\ast$ denotes the dimensionless variables.

2.5. The Validation of the Accuracy of the Modified Stokes Formula

To verify the accuracy of the modified Stokes formula in consideration of near-field hydrodynamic hindrance, the dimensionless hydrodynamic forces $F_h^{\ast }$, $F_s^{\ast }$ and $F_m^{\ast }$, respectively calculated by integrating the hydrodynamic stress tensor, the traditional Stokes formula and the modified Stokes formula, are compared. Three spherical particles ($a=5\mathsf{\mu }\mathrm{m}$) are suspended in a channel ($L\times W\times H=300 \mathsf{\mu }\mathrm{m}\times 200 \mathsf{\mu }\mathrm{m}\times 200 \mathsf{\mu }\mathrm{m}$) with a pressure-driven flow ($\Delta p=1\mathrm{Pa}$) along the length direction. The three particles are distributed on the plane parallel to the flow direction, as shown in Figure 2a, where d is the distance between the particle center and the origin. The dependences of the hydrodynamic forces on the ratio of the distance d and the particle radius a of the three particles are shown in Figure 2b–d. As we can see, when the particles get closer to each other, the traditional Stokes formula is not applicable for calculating the hydrodynamic force, while the modified Stokes formula can describe the hydrodynamic hindrance accurately. In this paper, the modified Stokes formula is able to calculate the hydrodynamic forces of dense particles.

3. Numerical Examples and Discussions

The fluid and particles in the present example were water and polystyrene beads, respectively. The densities of the particles were the same as the fluid, $\mathsf{\rho }={10}^3{\mathrm{kg}/\mathrm{m}}^3$; thus, the gravity was ignored. The permittivities and conductivities of the fluid were ${\mathsf{\epsilon }}_f=78{\mathsf{\epsilon }}_0$ and ${\mathsf{\sigma }}_f=2\times {10}^{-3}\mathrm{S}/\mathrm{m}$, respectively, where ${\mathsf{\epsilon }}_0=8.854\times {10}^{-12}\mathrm{F}/\mathrm{m}$ was the permittivity in a vacuum. The viscosity of the fluid was $\mathsf{\eta }={10}^{-3}\mathrm{Pa}\cdot \mathrm{s}$. The properties of the fluid were specified as the same in all examples of this work. Two types of particles with same permittivities but different conductivities ${\mathsf{\epsilon }}_{p1}={\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p1}=2\times {10}^{-4}\mathrm{S}/\mathrm{m},{\mathsf{\sigma }}_{p2}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ were adopted in this work to study the particle interactions. The real part $K_R$ of the C–M factors of the two types of particles versus frequency are shown in Figure 3. As we can see, in the low-frequency region, for particle 1, $K_R=-0.4286<0$, suffering a negative DEP force, called an nDEP particle. For particle 2, $K_R=0.4>0$, suffering a positive DEP force, called a pDEP particle. In this case, they are called dissimilar particles in this work. In the transition region, particle 2 was rapidly transformed from a pDEP particle into an nDEP particle while particle 1 remained as an nDEP particle. In the high-frequency region, both types of particles were nDEP particles as $K_R=-0.4759$. Accordingly, they are called similar particles in the present example.

3.1. The Interaction of Five Particles with Different Conductivities

Five particles with the same size ($a^{\ast }=1$) were positioned on a plane perpendicular to the uniform AC field ${\tilde{E}}_0=0.02$ in the electrolyte. The initial locations of the five particles are listed in

Table 1. The initial locations of the five particles.

Coordinate	Particle 1	Particle 2	Particle 3	Particle 4	Particle 5
$x^{\ast }$	0	2	−2	−2	2
$y^{\ast }$	0	3	3	−3	−3
$z^{\ast }$	0	0	0	0	0

. Four particles were symmetrically distributed about the central particle 1. The permittivities and conductivities of the particles were ${\mathsf{\epsilon }}_{pi}=2.5{\mathsf{\epsilon }}_0,i=1,5$ and ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$, ${\mathsf{\sigma }}_{p2}={\mathsf{\sigma }}_{p3}={\mathsf{\sigma }}_{p4}={\mathsf{\sigma }}_{p5}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$, respectively. When a uniform AC field with $f={10}^3\mathrm{Hz}$ was applied, particle 1 behaved as a pDEP particle while the other four particles behaved as nDEP particles, according to Figure 3. The initial particle locations and the final particle chain are shown in Figure 4, where the blue and red colors denote pDEP and nDEP particles, respectively. As we can see, the four nDEP particles were attracted by the central pDEP particle and moved towards the center. Finally the nDEP particles distributed symmetrically around the pDEP particle as the neighboring nDEP particles kept almost the same distance from each other, and a cross-shaped particle chain was formed.

The particle trajectories are shown in Figure 5. It can be seen that all particles moved only on the initial plane in a perpendicular electric field. The central pDEP particle 1 did not move all the time due to symmetrical forces. The nDEP particles 2 and 3 simultaneously suffered attractive forces from particle 1 and repulsive forces from each other. As a result, they moved towards particle 1 while moving away from each other, as did nDEP particles 4 and 5.

These in-plane motions and DEP behaviors can be connected to a special case observed in previous studies. It is well known that similar particles form particle chains parallel to the electric field and dissimilar particles form particle chains perpendicular to the electric field [2]. However, when the particles are initially positioned in the special situation where the connecting line of the particles is perpendicular to the electric field, the similar particles will repel each other while dissimilar particles attract each other. Both the directions of the repulsive force and the attractive force are parallel to their connecting line. In this present 3D model, any two particles are always in the special situation as the electric field is perpendicular to the particle plane. As a result, the interactive forces of the particles are always parallel to their initial plane and they can only move in the plane.

The variations of the DEP interactive forces of particle 1 and particle 2 with time are shown in Figure 6, where the blue color denotes particle 1 and the red color denotes particle 2. It can be seen that the force of particle 1 was zero all the time due to the symmetrical distribution of the other four particles, which is consistent with Figure 5, while the force of particle 2 increased rapidly with time and finally reached a constant value when particle 2 came into contact with particle 1. The DEP interactive forces of particles 3, 4, and 5 are similar to particle 2 and will not be presented.

When the conductivities of particles 2 and 4 were changed to ${\mathsf{\sigma }}_{p2}={\mathsf{\sigma }}_{p4}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ and the rest of the parameters remained the same as in the previous example, particles 1, 2, and 4 were pDEP particles and particles 3 and 5 were nDEP particles as shown in Figure 3. The initial locations and the final particle chain are shown in Figure 7, where the blue and red colors denote pDEP and nDEP particles, respectively. It can be seen that the two nDEP particles were respectively attracted by the neighboring two pDEP particles, moving towards the regions between the two pDEP particles. The pDEP particles repelled each other and, finally, an alternately arranged particle chain was formed.

The particles’ trajectories are shown in Figure 8, where the blue and red colors denote pDEP and nDEP particles, respectively. It can be seen that the central particle 1 did not move due to symmetrical forces, while the nDEP particle 3 was attracted by two pDEP particles 1 and 2, moving towards the region between them. Similarly, particle 5 moved towards the region between particles 1 and 4. Finally, a line-styled particle chain consisting of alternating arrangements of nDEP and pDEP particles was formed on the initial plane.

3.2. The Interaction of a Large Number of Particles with the Same Size

One hundred particles with the same size $a^{\ast }=1$ were randomly distributed on a plane chip filled with electrolyte perpendicular to a uniform AC electric field ${\tilde{E}}_0=0.02$, as shown in Figure 9a. The electrical properties of the fluid were the same as that of the previous examples, while 50 particles had the properties of ${\mathsf{\epsilon }}_{p1}=2.5{\mathsf{\epsilon }}_0,{\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$, and the other 50 particles had the properties of ${\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0,{\mathsf{\sigma }}_{p2}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$. When $f={10}^3\mathrm{Hz}$, the particles with ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ behaved as pDEP particles while the others were nDEP particles, according to Figure 3. As we can see in Figure 9c, similar particles repelled each other and dissimilar particles attracted each other and, finally, many particle chains with different styles and lengths were formed, for example line-styled, cross-shaped or a combination of both. Each particle chain consisted of alternately arranged nDEP and pDEP particles. The directions and shapes of the particle chains only depended on their initial distributions and the electric properties.

When the particle number increased to 200, half of the particles had the properties of ${\mathsf{\epsilon }}_{p1}=2.5{\epsilon }_0,{\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ and the other half had the properties of ${\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0,{\mathsf{\sigma }}_{p2}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$. The initial positions and the final particle chains when $f={10}^3\mathrm{Hz}$ are shown in Figure 9b,d. It can be seen that the particle chain patterns have similar styles but with more complexity, and more branches in the particles chains are observed than in the case of 100 particles.

3.3. The DEP Interaction of Large Numbers of Particles on a Bounded Circular Plate Chip

The working chambers of microfluidic chips are all confined spaces, and the chip wall’s effects on DEP particle interactions cannot be neglected in the cases of a larger number of dense particles. DEP particle interactions on a bounded circular plane chip are studied in this section, where 140 particles of uniform size $a^{\ast }=1$ were randomly distributed on a circular plane chip with a diameter $d^{\ast }=44$. A uniform AC electric field ${\tilde{E}}_0^{\ast }=0.02$ was applied perpendicularly to the plane chip, as shown in Figure 10a. Half of the particles had the electrical properties of ${\mathsf{\epsilon }}_{p1}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$, and the others had the properties of ${\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p2}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$. Three cases of frequencies, $f={10}^3\mathrm{Hz}$, $f={10}^6\mathrm{Hz}$, and $f={10}^7\mathrm{Hz}$, were studied in the present work. When $f={10}^3\mathrm{Hz}$, the two groups of particles experienced negative and positive DEP forces, respectively, according to Figure 3. Various particle chains were formed with a more complex structure compared with Figure 9. The scattered short particle chains tended to aggregate with a much larger size and more branches due to the restriction of the boundary, as shown in Figure 10b. When $f={10}^6\mathrm{Hz}$, the particles with ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ suffered much weaker positive DEP forces, while the others still suffered strong negative DEP forces. The attractive forces between the pDEP and nDEP particles became much weaker, while the nDEP particles still experienced strong repulsive forces between each other. As a result, the particle chains were more scattered compared to those at low frequency, as we can see in Figure 10c. When $f={10}^7\mathrm{Hz}$, all particles behaved as negative DEP particles. It can be expected that all particles experienced strong repulsive forces and moved away from each other, until all the repulsive forces of each particle reached a balance, resulting in an almost uniform particle distribution, finally, as shown in Figure 10d. It was also found that the final particle distributions were quite similar for different random distributions of particles.

Particle chain behaviors of different particle sizes and conductivities were also investigated in this work. One hundred particles were split into two groups (50, 50) and the particle properties were specified as $a_1^{*}=1$, ${\mathsf{\epsilon }}_{p1}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ and $a_2^{*}=2$, ${\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p2}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$, respectively, as shown in Figure 11a. The frequency of the electric field was $f={10}^3\mathrm{Hz}$. The rest of the parameters remained the same as in the previous examples. The small particles behaved as positive DEP particles and the large particles behaved as negative DEP particles in the present example. It can be seen that the large particles experienced large repulsive forces between each other and small attractive forces from the small particles. The large particles were repelled from the boundary, while the small particles were attracted by the large particles and, finally, long particle chains along the wall were formed, as shown in Figure 11c. In the case of particles with random radii, from $a^{\ast }=1$ to $a^{\ast }=2$, 50 particles with the properties of ${\mathsf{\epsilon }}_{p1}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p1}=6\times {10}^{-3}\mathrm{S}/\mathrm{m}$ behaved as positive DEP particles, and the other particles with the properties of ${\mathsf{\epsilon }}_{p2}=2.5{\mathsf{\epsilon }}_0$, ${\mathsf{\sigma }}_{p2}=2\times {10}^{-4}\mathrm{S}/\mathrm{m}$ behaved as negative DEP particles, as shown in Figure 11b. The particles were more likely to form clusters than long chains, as shown in Figure 11d.

4. Conclusions

The DEP interactions of a large number of dense spherical particles on a plane perpendicular to a uniform AC field were studied by an iterative dipole moment (IDM) model in this work. Numerical examples indicated that the IDM model is an effective method to predict DEP interaction behaviors of a large number of dense spherical particles in an electric field. Some novel phenomena were found in the results and are consistent with the previous studies. When the particles are initially positioned on a plane perpendicular to the electric field, they cannot move out of the plane. Similar DEP spherical particles on a plane perpendicular to an electric field always repel each other, and scatter away. A large number of dense particles in the bounded region eventually form a uniformly distributed particle pattern (equidistant between every two particles) that is almost independent of the initial positions. Dissimilar DEP particles in a plane perpendicular to an electric field always attract each other, forming alternately arranged particle chains, clusters, or groups. The directions and shapes of the final particle chains depend on the particle sizes, number, electrical properties, initial particle distributions, and electric field frequencies. Above all, similar particles and dissimilar DEP particles can easily be transformed into each other via frequency modulation which can be an effective method for particle manipulation in applications of biology, chemistry, and micro-fabrications.