Numerical study of acoustophoretic manipulation of particles in microfluidic channels

Abstract

The microfluidic technology based on surface acoustic waves (SAW) has been developing rapidly, as it can precisely manipulate fluid flow and particle motion at microscales. We hereby present a numerical study of the transient motion of suspended particles in a microchannel. In conventional studies, only the microchannel’s bottom surface generates SAW and only the final positions of the particles are analyzed. In our study, the microchannel is sandwiched by two identical SAW transducers at both the bottom and top surfaces while the channel’s sidewalls are made of poly-dimethylsiloxane (PDMS). Based on the perturbation theory, the suspended particles are subject to two types of forces, namely the Acoustic Radiation Force (ARF) and the Stokes Drag Force (SDF), which correspond to the first-order acoustic field and the second-order streaming field, respectively. We use the Finite Element Method (FEM) to compute the fluid responses and particle trajectories. Our numerical model is shown to be accurate by verifying against previous experimental and numerical results. We have determined the threshold particle size that divides the SDF-dominated regime and the ARF-dominated regime. By examining the time scale of the particle movement, we provide guidelines on the device design and operation.

Introduction

As a powerful lab-on-a-chip technology, acoustofluidics has attracted widespread interest because of its ability of the label-free manipulation of particles and cells suspended in liquid samples.^1,2 Devices based on acoustofluidics are easy to operate, cost-effective, non-contact, and cause little damage to cells.^3–5 When Surface Acoustic Waves (SAW) are generated at the substrate surface of a microchannel, part of the acoustic energy will radiate into the liquid inside the channel. ⁶ The scattering of the acoustics waves on the microparticles suspended in the liquid leads to the momentum transfer from the ultrasound field to particles, giving rise to a second-order time-averaged Acoustic Radiation Force (ARF) on solid particles. ⁶ Subject to the ARF, randomly dispersed particles at the inlet of the microchannel move laterally to several equilibrium positions at the channel outlet. Owing to the nonlinear terms in the Navier-Stokes equation, the fluid response to the harmonic acoustic perturbation is not exactly sinusoidal, especially in the viscous boundary layers near solid walls. ⁷ The time-averaged residual flow is referred to as acoustic streaming.^7,8 The effect of streaming on the particles’ motion can be evaluated by the time-averaged Stokes Drag Force (SDF). ⁸ The overall motion of the suspended particles is determined by the combined effect of the ARF due to the acoustic scattering and the SDF due to the acoustic streaming. ⁸ For large particles, the ARF is much greater than the SDF, so effective particle manipulation can be achieved. ⁸ However, as the particle radius is reduced to below a threshold value, the influence of the SDF is expected to be the dominant force, which prevents effective manipulation. ⁸

We neglect the influence of the lift force on particle motion, which may be induced by the velocity nonuniformity (Saffman force and shear-gradient lift force), particle spinning (Magnus effect), or solid boundary (wall-induced lift force). Given that the streaming velocity is small, its first-order and second-order derivatives in the cross-flow direction are expected to be even smaller. Therefore, the drag force is expected to be much larger than the Saffman force and shear-gradient lift force due to the fluid shear and the gradient of shearing, respectively. With regard to the Magnus effect, You et al. ⁹ found that for a small spherical particle such as one with a diameter of about 100 μm, even if the rotational speed reaches 1 million revolutions per minute, the lift force can still be neglected as compared with the drag force. The wall-boundary-induced lift arises from the bounded flow domain because of the adjacent wall boundary, which is only important when there is significant relative motion between the fluid and particle. Hence, the exclusion of the lift force is justified. It is worth noting, though, that the inclusion of the lift force is crucial for analyzing inertial microfluidics, which entirely relies on the channel geometry and flow-induced hydrodynamic effects to manipulate particles. ¹⁰

Due to the difficulty to measure the forces and velocities at microscales, numerical modeling has been extensively used to predict particle movement and interpret experimental findings.^8,11–14 The Finite Element Method (FEM) is one of the most widely-used computational methods to study multiple physical processes. ⁸ Muller et al. ⁸ built a FEM model based on the Helmholtz equation and wave equation. The model estimates the pressure and velocity distributions inside the acoustofluidic channel. ⁸ Using the perturbation theory for the acoustic wave propagation, Nama et al. ¹¹ built a model where only the microchannel’s bottom surface generates the SAW. The impedance boundary condition has been specified at all other boundaries to represent the wave absorption by the PDMS material. ¹¹ Similar studies were carried out by Guo et al., ¹² Mao et al., ¹³ and Sun et al. ¹⁴ who performed numerical and experimental studies to examine the particle trajectories and their final positions in the PDMS microchannel. They indicated the same finding that particles are focused at three positions within the microfluidic channel, corresponding to the locations of the acoustic nodes of the standing SAW imposed at the boundary. These numerical simulations are in qualitative agreement with experimental observations.

The previous experimental and numerical researches have been limited to the device configuration where the SAW passage into the microchannel is restricted to the bottom surface.^11–14 Recently, Cardiff University proposed an idea of generating the SAW at both the bottom and top boundaries of the channel. As such a configuration allows more acoustic energy to be introduced into the liquid and then scattered on the suspended particles, it is expected that the particles’ movement can be better handled this way. In this paper, the Model-P and Model-W denote the device configurations where acoustic actuation is imposed at the bottom surface and at both the bottom and top surfaces, respectively. The performance of these two configurations is compared in this study. Another shortcoming of the past research is that only the equilibrium positions of the microparticles are examined at the channel outlet. To the best knowledge of the authors, no investigation has been conducted on the time scale of the particle motion. This paper gives a detailed investigation of the transient motion of particles. The time scale required for particles to move to their designated positions is an important parameter for determining the length and flowrate of the microchannel.

In our studies, the first-order acoustic pressure and velocity fields, the time-averaged second-order velocity field, ARF, SDF, and particle trajectories are simulated using a FEM package. The simulations reveal the detailed characteristics of the acoustic field and the transient movement of microparticles. Prior the aforementioned model application, the paper first introduces the perturbation model and then verifies it against previous numerical and experimental results.

Mathematical model

Governing equations and perturbation solutions

The mass and momentum conservation laws govern the motion of the viscous compressible fluid, which can be described as ¹⁵ :

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{v}\right)=0$$

(1)

$$\rho \frac{\partial \mathbf{v}}{\partial t}+\rho (\mathbf{v}·\nabla )\mathbf{v}=-\nabla p+\eta {\nabla }^2\mathbf{v}+({\eta }_b+\frac{1}{3}\eta )\nabla (\nabla ·\mathbf{v})$$

(2)

where $\rho$ is density, $p$ is pressure, v is the flow velocity vector, $\eta$ is the shear dynamic viscosity coefficient, and ${\eta }_b$ is the bulk viscosity coefficient. ¹⁵ In order to enable a solution of these equations, an equation of state is needed to link the pressure and density. We assume the fluid to be weakly compressible, so there is a linear relationship between pressure $p$ and density $\rho$ . ¹⁵

$$p=c_0^2\rho$$

(3)

Combining equations (1)–(3) with suitable boundary and initial conditions, the problem is well-posed and ready to be solved.

The acoustic waves propagating through a fluid pose small perturbations to the fluid density, pressure, and velocity. In our analyses, we only deal with small perturbations, so the deviations of the velocity, pressure, and density from their static values are small. In our studies of the SAW-induced flow, all first-order quantities follow the harmonic variation, whose frequency is the same as the imposed ultrasonic excitation. Therefore, these independent functions can be expressed into perturbation series. ¹⁵

$$\rho ={\rho }_0+\epsilon {\rho }_1+{\epsilon }^2{\rho }_2$$

(4a)

$$p=p_0+\epsilon p_1+{\epsilon }^2{p}_2$$

(4b)

$$\mathbf{v}={\mathbf{v}}_0+\epsilon {\mathbf{v}}_1+{\epsilon }^2{\mathbf{v}}_2$$

(4c)

where $\epsilon$ is a small non-dimensional parameter. ¹⁵ The subscripts 0, 1, and 2 stand for the static, first-order, and second-order terms, respectively. In our problems, the static pressure ( $p_0$ ) and velocity ( ${\mathbf{v}}_0$ ) are independent of the acoustic actuation and are not considered. The first-order terms are the harmonic response with the same frequency as the input acoustic field. The second-order terms lead to the steady streaming fields as their time-averaged values over an oscillation period are non-zero. Higher-order terms are not considered in this study.

Substituting equations (4a)–(4c) into the governing equations, we can then arrange each term of the equations into the first-order, second-order, and higher-order small components. By retaining only the first-order components on both sides of the equations, we can derive the equations for the first-order pressure and velocity. ¹⁵

$$\frac{\partial {\rho }_1}{\partial t}=-{\rho }_0·(\nabla {\mathbf{v}}_1)$$

(5)

$${\rho }_0\frac{\partial {\mathbf{v}}_1}{\partial t}=-\nabla p_1+\eta {\nabla }^2{\mathbf{v}}_1+({\eta }_b+\frac{1}{3}\eta )\nabla (\nabla ·{\mathbf{v}}_1)$$

(6)

The second-order terms contain higher-frequency oscillations and non-zero mean components. By retaining only the second-order perturbation terms of the Navier-Stokes and continuity equations and then taking the time average, the second-order equations can be derived ¹⁵ :

$${\rho }_0\nabla ·<{\mathbf{v}}_2>=-\nabla ·<{\rho }_1{\mathbf{v}}_1>$$

(7)

$$\begin{array}{c}\multicolumn{1}{c}{-\nabla <p_2>+\eta {\nabla }^2<{\mathbf{v}}_2>+({\eta }_b+\frac{1}{3}\eta )}\\ \multicolumn{1}{c}{\nabla (\nabla ·<{\mathbf{v}}_2>)=<{\rho }_1\frac{\partial {\mathbf{v}}_1}{\partial t}>+{\rho }_0<({\mathbf{v}}_1·\nabla ){\mathbf{v}}_1>}\end{array}$$

(8)

where the point bracket represents the time average.

Time-averaged forces on suspended microparticles

Typical microfluidic experiments use polystyrene particles. We assume these particles are sparse and thus ignore the disturbance made to the flow by the existence of particles. The particle-particle interactions are not considered either. These microparticles are neutrally buoyant, so their movement in the fluid are subject to the ARF and SDF. The ARF ( ${\mathbf{F}}_{\mathrm{rad}}$ ) originates from the scattering of acoustic waves, while the SDF ( ${\mathbf{F}}_{\mathrm{drag}}$ ) is related to the second-order acoustic streaming. ⁶

The ARF on a single suspend particle can be calculated from the first-order solutions to the Navier-Stokes equations. Considering that the particle radius $r$ is much smaller than the wavelength of the imposed acoustic wave, the time-averaged ARF on a single suspended particle can be calculated as follows. ⁶

$${\mathbf{F}}_{\mathrm{rad}}=-\pi r^3\{\frac{2{\kappa }_0}{3}\mathrm{Re}[f_1^{*}{p}_1^{*}\nabla p_1]-{\rho }_0\mathrm{Re}[f_2^{*}{\mathbf{v}}_1^{*}·\nabla {\mathbf{v}}_1]\},$$

(9)

where ${\kappa }_0$ is the isentropic compressibility of the fluid, the asterisk represents the complex conjugation, the function Re [ ] represents the real part of the complex number in the square bracket, and the compressibility factor $f_1$ and density factor $f_2$ are given by the following formulae. ⁶

$$f_1=1-\frac{{\kappa }_p}{{\kappa }_0},$$

(10)

$$f_2=\frac{2(\frac{{\rho }_p}{{\rho }_0}-1)(1-\gamma )}{2\frac{{\rho }_p}{{\rho }_0}+1-3\gamma },$$

(11)

$$\gamma =-\frac{3}{2}[1+\mathrm{i}(1+\frac{\delta }{r})]\frac{\delta }{r}$$

(12)

where ${\rho }_p$ is density of solid particles, $\mathrm{i}$ is the imaginary unit, and $\delta$ is the thickness of the viscous boundary layer. The compressibility of water ${\kappa }_0$ and the compressibility of solid particles ${\kappa }_p$ are determined by the following formulae. ⁶

$${\kappa }_0=\frac{1}{{\rho }_0{c}_0^2}$$

(13a)

$${\kappa }_p=\frac{3(1-{\sigma }_p)}{[{\rho }_p{c}_p^2(1+{\sigma }_p)]}$$

(13b)

where ${\sigma }_p$ is the Poisson’s ratio of the solid particle. ⁶ The thickness of the viscous boundary layer can be expressed in terms of the dynamic viscosity coefficient η, angular frequency of the acoustic wave ω, and the initial density of the fluid ${\rho }_0$ . ¹¹

$$\delta =\sqrt{\frac{2\eta }{\omega {\rho }_0}}$$

(14)

The time-averaged drag force ${\mathbf{F}}_{\mathrm{drag}}$ on a single suspended particle of radius $r$ and velocity ${\mathbf{v}}_{\mathbf{p}}$ in a fluid with the streaming velocity $<{\mathbf{v}}_2>$ is given by ⁶ :

$${\mathbf{F}}_{\mathrm{drag}}=6\pi \eta r(<{\mathbf{v}}_2>-{\mathbf{v}}_{\mathbf{p}})$$

(15)

The SDF increases linearly with the particle size, as seen in equation (15), while the ARF increases with the particle size cubed, as seen in equation (9). Hence, once the particle size exceeds the threshold value, the ARF quickly becomes dominant with a further increase in the particle size. Conversely, for particles smaller than the threshold value, the ARF rapidly tends to zero with a further decrease of the particle size and then the particles’ motion is exclusively governed by the SDF. The motion of the solid particles is governed by Newton’s second law of motion ⁶ :

$$m_p\frac{d({\mathbf{v}}_{\mathbf{p}})}{\mathrm{dt}}={\mathbf{F}}_{\mathrm{drag}}+{\mathbf{F}}_{\mathrm{rad}}$$

(16)

where, $m_p$ is the mass of the particle.

Boundary conditions

The first-order equations (5) and (6) and second-order equations (7) and (8) allow us to numerically solve the pressure field and velocity field in the microfluidic channel to the second order approximation, but the solution has to be accompanied by appropriate boundary conditions. The time-averaged non-zero deformation at the channel boundary is generally in the sub-nanometer scale. Therefore, it is reasonable to ignore the time-averaged deformation of the channel boundaries. In calculating the time-averaged second-order flow field, all the walls correspond to no-slip boundary conditions. ⁷ The boundary conditions for the first-order acoustic response are more complicated and are the focus of the remaining of this section. The impedance boundary condition ( ${\Gamma }_{\mathrm{p}}$ ) is specified at all passive channel walls. ¹¹ The actuated boundary ( ${\Gamma }_a)$ is implemented by enforcing the Dirichlet actuation boundary at the interface of the substrate and fluid. ¹¹

The impedance boundary condition can be expressed to be ¹¹ :

$$\mathbf{n}·\nabla p_1=\mathrm{i}\frac{\omega {\rho }_0}{{\rho }_{\mathrm{pdms}}{c}_{\mathrm{pdms}}}{p}_1$$

(17)

where $c_{\mathrm{pdms}}$ is the speed of sound in the PDMS material for the channel walls, ${\rho }_{\mathrm{pdms}}$ is the density of PDMS, and $\mathit{n}$ is the outward pointing normal vector of the boundary surface.

Rayleigh waves are the main type of waves produced by the SAW generator. The wave propagates along the y-axis and decays exponentially, which penetrates through the microchannel wall to actuate the fluid inside. Under the acoustic actuation of a single wave, the displacement function at the substrate boundary is ¹¹ :

$$u_y(t,y)=0.6u_0{e}^{-C_{\mathrm{d}}y}\sin [\frac{-2\pi (y-\frac{w}{2})}{\lambda }+\omega t]$$

(18a)

$$u_z(t,y)=-u_0{e}^{-C_{\mathrm{d}}y}\cos [\frac{-2\pi (y-\frac{w}{2})}{\lambda }+\omega t]$$

(18b)

where ${\mathit{u}}_y$ and $u_z$ denote the SAW displacements in the y and z directions, respectively, $u_0$ is amplitude of the displacement, $\lambda$ is the wavelength, w is the width of the channel, and $C_{\mathrm{d}}$ is the decay coefficient. ¹¹ With the acoustic frequency of 6.65 MHz as considered in this paper, the decay coefficient $C_{\mathrm{d}}$ takes the value of 116  ${\mathrm{m}}^{-1}$ , which has been obtained by Nama et al. ¹¹ for a SAW propagating beneath an infinitely thick layer of water. Superposing the two SAWs traveling in opposite directions with a phase difference $\pi$ , the boundary displacement for the standing wave actuation is as follows.

$$\begin{array}{c}\multicolumn{1}{c}{u_y(t,y)=0.6u_0{e}^{-C_{\mathrm{d}}y}\{\sin [\frac{-2\pi (y-\frac{w}{2})}{\lambda }+\omega t-\pi ]+\sin [\frac{-2\pi (\frac{w}{2}-y)}{\lambda }+\omega t]\}}\end{array}$$

(19a)

$$u_z(t,y)=-u_0{e}^{-C_{\mathrm{d}}y}\{\cos [\frac{-2\pi (y-\frac{w}{2})}{\lambda }+\omega t-\pi ]+\cos [\frac{-2\pi (\frac{w}{2}-y)}{\lambda }+\omega t]\}$$

(19b)

We can then obtain the substrate velocity by taking the partial derivative of the displacement function with respect to time to impose over ${\Gamma }_a$ . ¹¹

$${\mathbf{v}}_1(t,y)=\frac{d\mathit{u}(t,y)}{\mathrm{dt}}$$

(20)

Numerical model

Model setup

Acoustophoresis chips are often manufactured by bonding the PDMS channel onto a piezoelectric substrate. ² The PDMS material is a silicon-based polymer and consists of Sylgard-184-silicone elastomer base and Sylgard-184-silicone elastomer curing agent at a mixing ratio of 10:1 by weight. ¹⁶ The material properties of this PDMS (10:1) have been specified in our numerical simulations. The piezoelectric substrate is composed of Lithium Niobate ( $\mathrm{LiNb}{\mathrm{O}}_3$ ). ² A pair of metallic interdigitated transducers (IDTs) are placed at the two ends of the piezoelectric substrate. ² Triggered by a harmonic electric signal, the IDTs generate surface acoustic waves, which radiate upwards into the microchannel. ² Stable acoustic pressure gradients are then formed in the fluid inside the microchannel. In the numerical simulation, we neglect any streamwise (x direction) variations caused by the acoustic waves, as the flow is dominated by the uniform flow along a long and straight channel along the x axis. Hence, we restrict the analysis to be two-dimensional over a rectangular cross-section in the y-z plane, as shown in Figure 1. The rectangular cross-section has height $h=125\mathrm{\mu }\mathrm{m}$ and width $w=600\mathrm{\mu }\mathrm{m}$ . The resonance frequency $f$ of SAWs is $6.65\mathrm{MHz}$ , and the wavelength $\lambda$ of SAWs is $600\mathrm{\mu }\mathrm{m}$ . All key parameters are listed in Table 1.

Table 1.

Key parameters used in the study.

Polystyrene
Density	${\rho }_p$	1050 $\mathrm{kg}/{\mathrm{m}}^3$
Speed of sound	$c_p$	2350  $\mathrm{m}/\mathrm{s}$
Poisson’s ratio 17	${\sigma }_p$	0.35
Compressibility	${\kappa }_p$	249 $1/\mathrm{Pa}$
Poly-dimethylsiloxane (PDMS 10:1)
Density	${\rho }_{\mathrm{pdms}}$	920  $\mathrm{kg}/{\mathrm{m}}^3$
Speed of sound	$c_{\mathrm{pdms}}$	1076.5  $m/\mathrm{s}$
Water
Density 18	${\rho }_0$	997  $\mathrm{kg}/{\mathrm{m}}^3$
Speed of sound 18	$c_0$	1497  $\mathrm{m}/\mathrm{s}$
Shear viscosity 18	$\eta$	0.89 mPa s
Bulk viscosity 18	${\eta }_b$	2.47 mPa s
Compressibility	${\kappa }_0$	448 $1/\mathrm{Pa}$
Acoustic actuation parameters
Wavelength (set by IDTs)	$\lambda$	600 $\mathrm{\mu }\mathrm{m}$
Forcing frequency	$f$	6.65 MHz

Figure 1.

Cross-sectional view of the Model-P (left) and Model-W (right) configurations: (a) positions of the Lithium Niobate substrate and water-filled PDMS channel and (b) locations of the impedance boundary ( ${\mathrm{\Gamma }}_{\mathrm{p}}$ ) and velocity boundary ( ${\mathrm{\Gamma }}_a)$ .

Numerical procedure

The governing equations are numerically solved according to the following procedure in COMSOL Multiphysics 5.4. First of all, the first-order acoustic field is obtained by solving equations (5) and (6). Then, the second-order acoustic field is acquired by solving equations (7) and (8). Finally, by combining information from the first-order and second-order fields, we estimate the forces on and trajectories of particles in the microfluid channel.

Mesh convergence analysis

We use triangular mesh for spatial discretization of the computational domain. The mesh convergence analysis is needed to find the suitable mesh resolution. Therefore, we continually reduce element size to decide the threshold at which the result becomes unaffected by further decreasing of the minimum element length $d_{\mathrm{mesh}}$ along the domain boundaries.

The mesh convergence function $C(g)$ is defined as follows. ⁸

$$C(g)=\sqrt{\frac{{\int }^{}{(g-g_{ref})}^2\mathrm{dydz}}{{\int }^{}{(g_{ref})}^2\mathrm{dydz}}}$$

(21)

where $g$ stands for the computed value of the first-order pressure $p_1$ , the first-order y-component velocity ${\mathrm{v}}_{1y}$ , the first-order z-component velocity ${\mathrm{v}}_{1z}$ , the second-order y-component velocity ${\mathrm{v}}_{2y}$ , or the second-order z-component velocity ${\mathrm{v}}_{2z}$ at a given mesh resolution $d_{\mathrm{mesh}}$ , $g_{\mathrm{ref}}$ represents the reference values, and the integration is over the whole computational domain. A non-dimensional parameter $\delta /$ $d_{\mathrm{mesh}}$ is introduced to quantify the mesh resolution, which is the ratio of the oscillatory boundary layer thickness over the element size. In order to sufficiently resolve the small-scale physics inside the boundary layer, this ratio needs to be greater than a critical value. As illustrated in Figure 2, large elements are arranged in the middle of domain, while small elements are arranged near the domain boundaries. To obtain the reference values $g_{\mathrm{ref}}$ , we adopt an extremely fine mesh with $\delta /$ $d_{\mathrm{mesh}}=5$ , which results in $5.2\times {10}^5$ elements in total. The coarsest mesh considered corresponds to $\delta /$ $d_{\mathrm{mesh}}=0.03$ , which results in 2478 triangular elements only. Figure 3 demonstrates a typical mesh convergence study result. The five curves in each plot represent the mesh dependence of the five computed variables, that is, $p_1$ , ${\mathrm{v}}_{1y}$ , ${\mathrm{v}}_{1z}$ , ${\mathrm{v}}_{2y}$ , and ${\mathrm{v}}_{2z}$ . The dashed line indicates that the convergence function reaches a value of 0.002 when $\delta /$ $d_{\mathrm{mesh}}=2$ , which is consistent with previous studies. ⁸ When $\delta /$ $d_{\mathrm{mesh}}=3.3$ , all the five computed variables have reached sufficient convergence. Consequently, we use the mesh size $\delta /$ $d_{\mathrm{mesh}}=3.3$ in the rest of the study. As can be seen from equations (9) and (15), the ARF and SDF are determined by the flow field. Therefore, once the flow field computation has converged, the force computation should become independent of the computational mesh as well.

Figure 2.

Typical computational mesh with element size $\delta /$ $d_{\mathrm{mesh}}=3.3$ .

Figure 3.

Variations of the mesh convergence functions $C(g)$ with the mesh resolution: (a) Model-P and (b) Model-W.

Model verification

In the first verification case, we set up the model to be exactly the same as that reported in Nama et al. ¹¹ The polystyrene microparticles are suspended in a water-filled PDMS channel (Model-P). Figure 4(a) demonstrates the comparison of the first-order pressure distribution. An apparent standing wave pattern is formed along the y-axis. The upwards-pointing magenta arrows indicate wave radiation from the bottom wall to the top wall. The first order pressure oscillates with a maximum amplitude of 12.9 kPa according to both studies. Figure 4(b) compares the first-order acoustic velocity distribution ( ${\mathbf{v}}_1$ ), and the maximum magnitude is 5.3 mm/s according to both studies. Figure 4(c) shows the distribution of the time-averaged second-order velocity $<{\mathbf{v}}_2>$ , with the maximum magnitude of 1.47  $\mathrm{\mu }\mathrm{m}/\mathrm{s}$ . The value of $<{\mathbf{v}}_2>$ is negligibly small and decreases rapidly away from the bottom surface. The present numerical results match well with previous numerical results, which verifies that the present numerical model is reliable.

Figure 4.

Comparison between the present results (right) and those in Nama et al. ¹¹ (left): (a) first-order pressure field, (b) first-order velocity field, and (c) time-averaged second-order velocity field.

In the second verification case, we set up the model to be exactly the same as that reported in Sun et al. ¹⁴ The $10\mathrm{\mu }\mathrm{m}$ -diameter polystyrene microparticles are suspended in a PDMS channel (Model-P). The resonance frequency of SAWs is $19.40$  MHz with a wavelength of $280\mathrm{\mu }\mathrm{m}$ . The width and height of the microchannel is 280  $\mathrm{\mu }\mathrm{m}$ and 60  $\mathrm{\mu }\mathrm{m}$ , respectively. Figure 5(a) demonstrates a snapshot of the particles’ positions observed in the experiment from the top of the channel. ¹⁴ Randomly-dispersed particles exposed to the SAW field are eventually focused into three parallel lines along the microfluidic channel. Figure 5(b) and (c) show the numerical simulations of the particles’ trajectories, indicated with blue curves, and the particles’ final positions, indicated with big red dots. It is seen that the present simulations agree well with the simulations in Sun et al. As the channel width is equal to the SAW wavelength, there are three pressure nodes (PN) over the width of the channel, as labeled in Figure 5(b). The microparticles are initially uniformly distributed over the cross section. Subject to the SAWs, they finally accumulate to eight positions over the cross section, four above the middle PN and two above each of the two sidewall PNs. When looking from the top of the channel, microparticles are aligned into three parallel queues along the channel, with two immediately by the sidewalls and one coinciding the channel centerline. Therefore, the numerical simulations are consistent with the experimental observations.

Figure 5.

Comparison between the present and Sun et al. ¹⁴ : (a) plan view of particles’ trajectories in experiments of Sun et al., ¹⁴ (b) cross-sectional view of particles’ trajectories and final positions in simulation of Sun et al., ¹⁴ and (c) cross-sectional view of particles’ trajectories and final positions in present simulation.

In our subsequent study, the channel dimensions and acoustic wave properties are taken to be the same as those used in Nama et al. ¹¹ In a separate research by the authors, we varied the channel dimensions and showed that the flow fields over adjacent half wavelength distances in the transverse direction (y-direction) demonstrate an antisymmetric relationship, as can also been seen in Figure 4. In the vertical direction (z-direction), the flow field roughly repeat itself for every wavelength of the acoustic propagation in liquid. Hence, we can easily infer the flow fields and the particle trajectories in microchannels of other dimensions and with acoustic waves of different frequencies. The flow field and particle motion also positively correlate with the amplitude of the acoustic waves specified at the boundary and the correlation is roughly linear.

Comparison between Model-P andModel-W configurations

Flow field

In the conventional Model-P setup, acoustic waves are introduced into the fluid from the bottom wall only. As seen in Figure 4(a), this setup results in the maximum amplitude of the first-order pressure oscillation of 12.9 kPa. The imposed standing wave has three pressure nodes at the bottom surface, located by the sidewalls and in the middle. The acoustic pressure induced by the bottom excitation attains the maximum amplitude some distance above the two antinodes at the bottom surface. In the Model-W setup, the acoustic waves are introduced from both the bottom and top walls to the fluid. As expected, Figure 6(a) shows a significant increase in pressure strength, with the maximum amplitude of the first-order pressure amplitude increased to 228 kPa. In the Model-W setup, the pressure distribution is antisymmetric about the horizontal line through z = 62.5  $\mathrm{\mu }\mathrm{m}$ – half height of the channel. Such an antisymmetric pattern is caused by the two SAWs of the same magnitude and 180° phase difference propagating from the bottom and top surfaces, respectively. In the Model-P setup, the first-order velocity ( ${\mathbf{v}}_1$ ) has a maximum magnitude of 5.3 mm/s and the velocity distribution is antisymmetric about the vertical line through y = 0, as indicated in Figure 4(b). Figure 6(b) shows that the Model-W setup produces the greatly increased magnitude of the first-order velocity (70.4 mm/s) and the velocity distribution is antisymmetric about both the vertical line through y = 0 and the horizontal line through z = 62.5  $\mathrm{\mu }\mathrm{m}$ . Hence, the Model-W setup produces much stronger wave resonance in the fluid than the conventional Model-P setup.

Figure 6.

Flow field in a Model-W setup: (a) first-order pressure field, (b) first-order velocity field, and (c) time-averaged second-order velocity field.

Figure 4(c) shows that there are four streaming eddies over the channel width in the Model-P setup. The streaming flow is mirror-symmetric about the vertical line through y = 0, thus two eddies are clockwise and the other two are anti-clockwise. The time-averaged second-order velocity has the maximum magnitude of 1.47  $\mathrm{\mu }\mathrm{m}/\mathrm{s}$ , which occurs at the bottom wall. In the Model-W setup. Figure 6(c) shows that the velocity field is mirror-symmetric about both the vertical line through y = 0 and the horizontal line through z = 62.5  $\mathrm{\mu }\mathrm{m}$ . There are a total of eight streaming eddies, with four clockwise ones and four anti-clockwise ones. The maximum magnitude of the time-averaged second-order velocity is increased to 30.6  $\mathrm{\mu }\mathrm{m}/\mathrm{s}$ , which occurs at the bottom and top walls. The streaming flow features shown in Figures 4(c) and 6(c) are similar to Rayleigh’s analytical solution, but the current simulations fully consider the effects of the channel sidewalls and the channel height.

Hydrodynamic forces on particles

The ARF and SDF acting on a solid particle can be calculated once the first-order flow field and the time-averaged second-order flow field have been determined. Their magnitude depends on the particle size. In this section, we consider the situation when the particles are stationary with a radius of $5\mathrm{\mu }\mathrm{m}$ as an example to illustrate these forces.

Figure 7(a) shows that the maximum value of the acoustic radiation force is $8.5\times {10}^{-14}$  N in the Model-P setup and similar to the first-order flow field the ARF distribution is mirror-symmetric about the vertical line through y = 0. In the Model-W setup, as shown in Figure 7(b), the maximum value of the ARF is $7.1\times {10}^{-12}$  N, which is almost two orders of magnitude greater than that of the Model-P setup. Consistent with the flow field distribution, the ARF distribution in the Model-W setup is mirror-symmetric about both the vertical line through y = 0 and the horizontal line through z = 62.5  $\mathrm{\mu }\mathrm{m}$ .

Figure 7.

ARF distribution in the y-z plane, $r=5\mathrm{\mu }\mathrm{m}$ : (a) Model-P and (b) Model-W.

Because we assume that particles’ velocity is zero to demonstrate the forces on them, the direction of the SDF remains the same as that of the time-averaged second-order velocity. As seen in Figure 8, the maximum SDF values for the Model-P and Model-W setups are $2.2\times {10}^{-16}$  N and $1.2\times {10}^{-14}$  N, respectively. The maximum SDF values occur close to the boundary where the acoustic waves are generated, corresponding to the locations of the maximum streaming velocities. Hence, when particle radius is $5\mathrm{\mu }\mathrm{m}$ , the SDF is about two orders of magnitude smaller than the ARF in either Model-P or Model-W setup and thus the influence of the SDF on the particles’ migration can be neglected.

Figure 8.

SDF distribution in the y-z plane, $r=5\mathrm{\mu }\mathrm{m}$ : (a) Model-P and (b) Model-W.

Microparticle trajectories

Knowing the forces acting on the microparticles suspended in the water-filled microchannel, we can plot the particle trajectories using the particle tracing module. For the time-domain simulations, we select the generalized alpha solver in COMSOL, setting the alpha parameter to be 0.5 and specifying a fixed time step Δt = 0.01 s. As can be seen in Figure 9, particles are uniformly distributed over the cross-section in the y-z plane at time zero. In Figure 10, every line represents the trajectory made by a single particle, whereas the color of the line indicates the instantaneous particle velocities. In our study, we assume that the acoustic waves propagate perpendicular to the sidewalls of the microchannel. Hence, the acoustic perturbation is normal to the main flow in x direction, which is the reason for the two-dimensional cross-sectional simulation of the problem in the y-z plane. In the longitudinal direction, that is, x direction, the flow is nearly uniform and the suspended microparticles can be assumed to have the same velocity as the fluid. We neglect the acoustic influence on the movement of the fluid and particles in the longitudinal direction. The three-dimensional motion of a microparticle can be inferred from the particle’s two-dimensional trajectory in the y-z plane. Along the x direction, the particle’s velocity can be estimated to be the liquid flowrate divided by the cross-sectional area of the microchannel, while the particle’s displacement is equal to the product of the velocity and time. Hence, the particle follows a 3-D spiral trajectory to reach its equilibrium position and then maintains a straight trajectory further downstream.

Figure 9.

Initial positions of the 297 particles over the cross-section.

Figure 10.

Comparison of the particle in the y-z plane between Model-P and Model-W: (a) $r=0.5\mathrm{\mu }\mathrm{m}$ , t = 25 s, (b) $r=1.0\mathrm{\mu }\mathrm{m}$ , t = 25 s, (c) $r=2.0\mathrm{\mu }\mathrm{m}$ , t = 25 s, (d) $r=3.0\mathrm{\mu }\mathrm{m}$ , t = 25 s and (e) final positions for particles greater than the threshold size.

For small particles in the test (r = 0.5 μm and r = 1 μm), as can be seen in Figure 10(a) and (b), both the ARF and SDF acting on the particles are the small, and the SDF arising from the acoustic streaming is the dominant force to determine the particles’ motion. Hence, the particles’ motion is slow and the particles follow the circulation eddies whose characteristics are visualized in Figures 4(c) and 6(c). Corresponding to the stronger streaming flow, the particles’ motion in the Model-W setup is more apparent, as seen in the right graphs in Figure 10(a) and (b). However, the SDF drives particles to loop around streaming vortices rather than concentrate at certain locations. In the Model-P setup, the particles’ displacements is hardly noticeable in 25 s, except in the region very close to the bottom surface.

For the intermediate particle size (r = 2 μm), as seen in Figure 10(c), only the Model-W setup can lead to the gradual concentration of particles in certain regions. In the Model-P setup, the ARF and SDF have similar magnitude and their relative strength varies from region to region. In the upper part of the channel, the ARF is stronger, so particles demonstrate a tendency to slowly move toward the top surface. In the lower part of the channel, the SDF is stronger, so there is no sign of particle concentration. Instead, the large SDF near the bottom is responsible for the circulation of particles inside streaming eddies.

For the larger particles, Figure 10(d) shows a clear trend of particle concentration driven by the dominant ARF in the channel in both setups. In the Model-P setup, particles are pushed toward the middle and four corners of the microchannel. In the Model-W setup, particles are pushed to the middle and the mid-height region in the channel. The final locations where particles will come to rely on the time averaged ARF field. Figure 10(e) highlights the final locations of the particles with red boxes under sufficient duration of the standing SAW exposure. These final locations are the same for all particle sizes larger than the threshold value, but the time required for particles to move to these locations depends on the particle size, which will be discussed in the next section.

Time scale of particle motion

In practical applications, it is extremely important to know the time required for particles to move to their designated positions, as this time determines the minimum length of the channel for successfully microparticle manipulation. Because of the greater magnitude of the ARF acting on large particles, large particles reach their designated positions in shorter time than small particles. To facilitate analyses, we set a target box centered at each of the final positions, as displayed in Figure 10(e). When a particle has moved into a target box, then it is regarded to have arrived at its designated position. The side length of the target box is varied from 30 to 90  $\mathrm{\mu }\mathrm{m}$ in the current analysis. A large range of particle sizes have been examined. Figure 11 plots the increase in the numbers of particles arriving at the designated positions versus time for Model-P and Model-W configurations.

Figure 11.

Variation of the number of particles arriving at designated positions with SAW exposure time and its sensitivity to particle radius, SAW excitation method, and target size: (a) side length of target box 30 $\mathrm{\mu }\mathrm{m}$ , (b) side length of target box 60 $\mathrm{\mu }\mathrm{m}$ and (c) side length of target box 90 $\mathrm{\mu }\mathrm{m}$ .

Our study shows that the critical particle radius is about 2.4 and 1.9  $\mathrm{\mu }\mathrm{m}$ for Model-P and Model-W setups, respectively. For very small particles, they experience negligibly small ARF and SDF and thus remains almost stationary during the computation. For intermediate particles whose radii are still smaller than the critical value, some particles slowly circulate about the streaming vortices. Therefore, some particles move into and then out of the target boxes. However, they will not rest in the target boxes. For particles larger than the threshold size, more and more particles fall in the target boxes and never come out. Eventually, all the 297 particles arrive at their target positions. The speed of the particle accumulation is clearly seen to be positively correlated with the particle size and the target window size.

Figure 12 quantitatively shows how the time, required for complete particle concentration, varies with the particle radius, target size, and model setup. The most important finding is that the Model-W setup outperforms the traditional Model-P setup by slightly reducing the threshold particle radius for effective separation and significantly shortening the time scale of the transient process.

Figure 12.

Variations of the minimum required SAW exposure time with particle radius, SAW excitation method, and target size.

Conclusions

In this paper, we examine the detailed flow and hydrodynamic force fields in acoustophoretic channels, including the first-order acoustic pressure ( $p_1$ ), first-order velocity ${\mathbf{v}}_1$ , time-averaged second-order velocity $<{\mathbf{v}}_2>$ , ARF, and SDF. The study emphasizes the transient motion of particles and the results filled the gap in our understanding of the time scale of the acoustophoretic manipulation. The comparisons between the Model-P and Model-W highlight the importance of the way in which acoustic actuation is exerted. Our results indicate that Model-W observes significant enhancement of the acoustophoretic forces acting on microparticles, including both the ARF and SDF. As a result, the Model-W configuration exhibits a better separation efficiency as it is capable of separating smaller particles and leads to much faster particle movement. For the channel geometry, material property and SAW parameters considered in this study, the threshold particle radius is reduced from 2.4 to 1.9  $\mathrm{\mu }\mathrm{m}$ and particles require only one third or one half of the time to reach their designated positions, as compared with the situation in Model-P. As particles can move to their equilibrium positions in a shortened time, the length of the microchannel can be reduced and the flow rate through the channel can be increased. According to the dimensional analyses of acoustophoretic manipulation, the main conclusions drawn from this study regarding the relevant performance of Model-P and Model-W are expected to remain the same even if the channel dimensions, acoustic wave magnitude and frequency are different.