Microparticle Acoustophoresis in Aluminum-Based Acoustofluidic Devices with PDMS Covers

We present a numerical model for the recently introduced simple and inexpensive micromachined aluminum devices with a polydimethylsiloxane (PDMS) cover for microparticle acoustophoresis. We validate the model experimentally for a basic design, where a microchannel is milled into the surface of an aluminum substrate, sealed with a PDMS cover, and driven at MHz frequencies by a piezoelectric lead-zirconate-titanate (PZT) transducer. Both experimentally and numerically we find that the soft PDMS cover suppresses the Rayleigh streaming rolls in the bulk. However, due to the low transverse speed of sound in PDMS, such devices are prone to exhibit acoustic streaming vortices in the corners with a relatively large velocity. We predict numerically that in devices, where the microchannel is milled all the way through the aluminum substrate and sealed with a PDMS cover on both the top and bottom, the Rayleigh streaming is suppressed in the bulk thus enabling focusing of sub-micrometer-sized particles.


Introduction
Acoustofluidic devices based on bulk acoustic waves for microparticle handling are traditionally made by microfabrication in acoustically hard materials, such as silicon or glass, yielding acoustic resonators with relatively high Q values [1]. Such devices can be fabricated controllable and with high accuracy, but may suffer from a high production cost. Several solutions have been proposed and successfully demonstrated to achieve simple and inexpensive acoustofluidic devices. Simple mass-produced glass capillary tubes have been used for trapping of microparticles [2][3][4][5][6][7] and even nanoparticles [8,9]. The first attempts of using polymer-based devices have been published showing applications such as focusing of polymer beads [10][11][12][13], lipids [10], and red blood cells [11,14], as well as blood-bacteria separation [15] and purification of lymphocytes [16]. Although a cheaper material, polymers are difficult to use in acoustofluidics due to their low acoustic contrast relative to water, but recently it was shown how to circumvent this problem by use of the whole-system-resonance principle [17]. According to this principle, the acoustic contrast between the ambient air and the whole device allows the excitation of specific whole-system vibrational resonance modes, which support strong acoustofluidic responses in a given embedded liquid-filled cavity, without the corresponding acoustic pressure being a localized cavity resonance as in conventional devices.
Another solution is to combine acoustically hard materials with the soft rubber polydimethylsiloxane (PDMS). In 2012, Adams et al., presented experiments and simulations of a high-throughput, temperature-controlled microchannel acoustophoresis device made by rapid prototyping. The device was based on a PDMS-gasket defining the side-walls and shape of the acoustofluidic chamber, which then was sealed using standard microscope slides [18]. Also Xu et al., used such a glass-PDMS-glass structure in their recent device for isolation of cells from dilute samples using bead-assisted acoustic trapping [19]. Similarly in 2018, Gautam et al., designed, fabricated, and tested simple and inexpensive micromachined PDMS-covered aluminum-based microfluidic devices for acoustic focusing of particles and cells [20]. These devices appear to be versatile and truly simple to fabricate, as the desired microchannel system is micromilled into the surface of an aluminum base and bonded with a PDMS cover. Since neither theory nor simulation was presented by Gautam et al. [20], we develop in this work a numerical model for such devices and validate them experimentally: In Section 3 for a basic design similar to that in Ref. [20], where a microchannel is micromilled into the surface of an aluminum base and sealed with a PDMS cover as sketched in Figure 1 and with the dimensions given in Table 1; and in Section 4 for a geometrically symmetric, but anti-symmetrically ac-voltage-actuated device.
In the main part of our work, we use the validated numerical model to predict the answer to the following three questions: (1) Does the exceptionally low transverse speed of sound ∼100 m/s in the nearly incompressible PDMS imply a singular behavior of the acoustic streaming near the PDMS-aluminum corners of the device? (2) Does anti-symmetric excitation using a split-top-electrode configuration as in Ref. [17] lead to better acoustophoresis? And finally, (3) does the use of the acoustically soft PDMS cover lead to a suppression of the bulk-streaming rolls generated by the water-PDMS interface, and if so, would PDMS covers sealing both the top and bottom part of devices, where the microchannel is milled all the way through the aluminum base, lead to a suppression of all bulk-streaming rolls in the device? The answers to these questions turn out to be predominantly affirmative.

PDMS Pz26
Water Aluminum (a) (b) inlet outlet Figure 1. (a) 3D sketch drawn to scale of the single-PDMS-cover (beige) aluminum-based (light gray) device driven by a piezoelectric PZT-Pz26 transducer (dark gray, placed on a protrusion next to the PDMS) with silver electrodes (not shown). The inlet and outlet channels are marked by small circles. (b) Cross-sectional view of the device at the vertical center plane, used in the 2D model presented in Section 3, where the 9-µm-thick electrodes of the Pz26 transducer are connected to ground (purple) and to the driving ac-voltage (red). The gap between the Pz26 and the PDMS cover is 1 8 W fl .

Theory: The Governing Equations and Boundary Conditions
We follow the work by Skov et al. [21] in our theoretical and numerical modeling of the acoustofluidic system sketched in Figure 1. The continuum fields of the model are the following: An electric potential scalar fieldφ(r, t) is present in the piezoelectric domain (Pz26), an elastic displacement vector fieldũ(r, t) is present in all four solid domains (Pz26, aluminum, silver electrodes, and PDMS), and an acoustic pressure scalar fieldp 1 (r, t) as well as a steady acoustic streaming velocity vector field v 2 (r) and pressure p 2 (r) are present in the water-filled microchannel. As in typical experiments, the system is driven by a time-harmonic ac-voltage V 0 e −iωt of amplitude V 0 and angular frequency ω = 2π f , where f is the ultrasound frequency, applied to the silver electrodes of the Pz26 transducer. Consequently, to first order in the applied voltage, all fields acquire a time-harmonic phase factor e −iωt and a complex-valued spatially varying amplitude, In the following, the common phase factor e −iωt is left out of the linear equations. Only the space-dependent amplitudes are computed, but multiplying them by e −iωt recovers the time dependence.

The Piezoelectric Transducer
We model a standard lead-zirconate-titanate (PZT) transducer of type Pz26 polarized in the z direction. The mechanical stress tensor σ and the electric displacement field D are given in terms of the gradients of the elastic displacement field u and the electric potential ϕ by the electromechanical coupling matrix. Using the compact Voigt notation, this constitutive relation is given by Here, C ik are the elastic coefficients, ε ik are the electric permittivities, and e ik are the piezoelectric coupling constants. The remaining components of σ are given by the symmetry relations σ ik = σ ki .
The governing equations in the piezoelectric Pz26 material of mass density ρ sl are the Cauchy equation for the elastic displacement field u and, for frequencies less than 100 MHz in systems smaller than 1 m without free charges, the quasi-static Gauss law for the electric potential ϕ, The boundary conditions for the Pz26 domain with surface normal vector n are continuity of ϕ, u and σ · n across the elastic solid interfaces to the silver electrodes, and zero stress as well as zero electric charge on surfaces exposed to air, Pz26 domain ← charged electrode: ϕ = V 0 , σ · n and u continuous, Pz26 domain ← grounded electrode: ϕ = 0, σ · n and u continuous, Pz26 domain ← air: n · D = 0 and σ · n = 0. (4c) Here, the notation A ← B refers to the influence of B on domain A with outward surface normal n.

The Elastic Aluminum Base, PDMS Cover, and Silver Electrodes
The aluminum base, the PDMS cover, and the silver electrodes on the Pz26 transducer can all be described as isotropic, linear elastic materials. Again using the compact Voigt notation, the corresponding constitutive equation relating the strain and the stress are given by For isotropic materials only C 11 and C 44 are independent, since the constraint C 12 = C 11 − 2C 44 applies. The transverse and longitudinal speed of sound becomes c lo = C 11 /ρ sl and c tr = C 44 /ρ sl , respectively. The governing equation for the displacement field u in an elastic solid with mass density ρ sl is the Cauchy equation The boundary conditions for the elastic solid domains with surface normal vector n are continuity of u and σ · n across the elastic solid interfaces, and zero stress on surfaces exposed to air, solid domain ← adjacent solid: σ · n and u continuous, solid domain ← air: The boundary conditions for the fluid-solid interfaces are given below.

Pressure Acoustics in the Water-Filled Microchannel
The governing equation for the acoustic pressure p 1 inside a fluid with mass density ρ fl , speed of sound c fl , dynamic viscosity η fl , and bulk viscosity η b fl is the Helmholtz equation The acoustic velocity v 1 in the bulk of the fluid is a potential flow given by the acoustic pressure p 1 , The boundary conditions for the fluid-solid interface are continuity of the stress and the velocity. In Ref. [21] the explicit form of these boundary conditions were derived taking into account the viscous boundary layer in the fluid. By introducing the shear wavenumber k s = (1 + i)/δ, where δ = 2η fl /(ωρ fl ) is the thickness of the viscous boundary layer, and expressing the displacement velocity as v sl = −iωu, the fluid-solid boundary conditions become, Here, the viscous loss in the fluid is taken into account in the boundary conditions, and it appears through the slip velocity v δ = (v sl − v 1 ).

Acoustic Streaming in the Water-Filled Microchannel
Steady acoustic streaming arises in the fluid due to the inherent nonlinear fluid dynamics, here the time-average of products of time-harmonic fields. Following Bach and Bruus [22], the streaming velocity v 2 is governed by a real-valued incompressible Navier-Stokes equation with a body force due to the real part of the acoustic energy-flux density, The no-slip condition on the fluid-solid interface requires that v 2 = − (u · ∇) v 1 , which leads to the boundary condition for the acoustic streaming v 2 at the fluid-solid interface, where the superscript "0" indicates a field evaluated along the fluid-solid interface. The subscripts and ⊥ and the unit vectors e and e ⊥ refer to the parallel and normal direction to the solid wall respectively. The steady pressure p 2 accompanying v 2 is governed by the continuity for v 2 , and since it only appears as a gradient in Equation (11), we must fix the level of p 2 by the constraint Ω fl p 2 dydz = 0.

The Acoustic Radiation and Drag Force on Suspended Microparticles
Suspended spherical particles of radius a and mass density ρ pa will experience the radiation force F rad due to acoustic scattering, which is given by the potential U rad with monopole and dipole scattering coefficients f 0 and f 1 , respectively [23], The velocity v pa of a suspended particle at position r(t) in the water is determined by a balance between the Stokes drag force F drag = 6πη fl a(v 2 − v pa ), the acoustic radiation force F rad , and the buoyancy force The particle trajectory is then given by integration in time as r pa (t) = t 0 v pa r(t ) dt . Inertia is neglected here, because the largest particle with radius a = 2.4 µm moves slower than v pa < 1 cm/s yielding a small particle Reynolds number, Re pa = 1 η fl ρ fl av pa ≈ 0.03 1. We define the horizontal y axis such that the channel is centered around y = 0. To determine at which frequencies f acoustic resonance modes appear that leads to good particle focusing towards the vertical center plane located at y = 0, we introduce the focusing figure of merit F ( f ). This is a modified version of the figure of merit R defined in Ref. [17], as follows: F should be large, when at the same time the average acoustic energy density E ac is large, and the acoustic radiation force F rad has the property that its average horizontal y component is large and points toward the center, whereas its average vertical z component has a small magnitude, Here, sgn(−y) designates minus the sign of the y coordinate.

Numerical Implementation and Experimental Validation of the Single-PDMS-Cover Model
To illustrate the numerical implementation of the model, we choose the geometry as sketched in Figure 1, with the dimensions listed in Table 1, and the material parameters given in Appendix A. This device geometry is similar to the one in Ref. [20], namely an aluminum-based device with a straight channel and a single PDMS cover. Our general numerical modeling is then validated experimentally for this device and further supported by the results presented in Section 4, for the anti-symmetric actuated design.

Model Implementation in COMSOL Multiphysics
The model system of Figure 1 as well as its governing equations and boundary conditions given in Section 2, are implemented in the commercial finite-element method software COMSOL Multiphysics 5.4 [24] using the weak form PDE module, closely following the method presented in Ref. [21]. We use quartic Lagrange shape functions for p 1 , cubic for u, ϕ, and v 2 , and quadratic for p 2 . For simplicity, we approximate the device by an infinitely long straight channel, and thus restrict the computation to the 2D cross section. The material parameters used in the model are given in Appendix A. The simulations were performed on a workstation with a 3.5-GHz Intel Xeon CPU E5-1650 v2 dual-core processor and with 128 GB RAM.

Manufacturing Method of the Chip for the Experimental Validation
A sketch of the device fabrication method is shown in Figure 2. A more detailed description of the process steps is given in the following.

Plasma activation
Thin film PDMS cover (1a) The micromachined substrate was cleaned with acetone, ethanol, and Milli-Q (Millipore Corporation, Burlington, MA, USA), and dried on a 140 • C hotplate for 2-3 min prior to bonding with the PDMS film covering the channel. This constitutes the base of the device.
The PDMS covers. Sylgard 184 silicone elastomer (Dow Corning, Ellsworth Adhesives, Germantown, WI, USA) was mixed with the curing agent at the commonly used weight ratio of 10:1 and degassed. Two types of covers were made: (a) Thin PDMS-film covers, by deposition of 1 mL of elastomer on a 100-µm-thick 100-mm-by-100-mm plastic transparency sheet (Mylar), followed by spin-coating and curing at 65 • C for 15-30 min, and (b) 1.5-mm-thick PDMS covers by conventional mould casting.
Device assembly. Afterwards, the cured PDMS film and cleaned aluminum substrates were treated with air plasma in a Zepto plasma cleaner (Diener electronic GmbH + Co. KG, Ebhausen, Germany) for 60 s. PDMS and aluminum were subsequently bonded together and cured at 80 • C for 4 min. After curing, the Mylar sheet was removed from the PDMS-aluminum assembly. For flow connections, silicone tubes with inner diameters that match 1/16" Teflon tubings, were glued to the base of the device. A PZT ceramic transducer (Pz26, Meggitt A/S, Kvistgaard, Denmark) designed for 2 MHz actuation was superglued to the final device.

Experimental Validation of the Numerical Model
The electrodes of the Pz26 transducer were coupled to an ac-voltage generator operating at 20 V peak-to-peak (V 0 = 10 V) at frequencies ranging from 1.5 to 2.5 MHz. After stopping the particle-loading flow, the position and velocity v pa of fluorescently-marked polystyrene particles (see Table A3 in Appendix A) was measured using the single-camera general defocusing particle tracking (GDPT) technique [25,26] with fluorescent polystyrene beads. We use a 10 µm × 5 µm grid size and a recorded image frame rate of 10 Hz, and the motion of 2a = 4.8, 2.0, and 1.0 µm-diameter tracer particles is tracked for 30, 60, and 120 s, respectively. During the data processing, the outliers were filtered out by limiting the displacement deviation of all the particles, by limiting the velocity magnitude, and by restricting the particle count to 2 in each grid.
To validate our 2D model, we compare in Figure 3 simulated and measured particle behavior. The top row show top-view micrographs of the microchannel under flow-through condition, where 4.8-µm-diameter particles in bright-field are seen to focus in the center, Figure 3a. This resonance mode is called "S" for "side actuated" and was located at f = 2.048 MHz. This focusing is confirmed by the fluorescent image Figure 3b, which however also reveals particles accumulating near the corners of the device. In Figure 3c,d, we show measured and simulated particle velocities v pa in the vertical cross section for different particle diameters. It is seen that the 2D model captures the following five main features in the measured particle velocity field, even though the 2D model geometry of Figure 1b assumes a translational invariant cross section along the x direction, which the experimental 3D geometry of Figure 1a clearly does not have. (1) The numerically predicted resonance is located at f = 1.803 MHz, only 12% below the experimental value. (2) As expected, for large particles 2a = 4.8 µm the motion is dominated by the radiation force, which is partly focusing in the vertical nodal plane at y = 0 and partly pointing towards the soft PDMS cover and the top corners. As the particle size is reduced to 2a = 2.0 µm and then further to 2a = 1.0 µm, acoustic streaming becomes more dominant for the particle motion. This cross-over is clearly seen both in model and measurement of Figure 3c

Modeling of Single-PDMS-Cover Devices with Anti-Symmetric Voltage Actuation
Following the results of Ref. [17], we investigate numerically, if better acoustophoresis, quantified by the focusing figure of merit F defined in Equation (15), is obtained by exciting the half-wave resonance mode (which is anti-symmetric around the nodal plane at y = 0). This is achieved by splitting the top electrode and applying an anti-symmetric ac voltage, as sketched in Figure 4 for a device with a split-gap of size 50 µm × 40 µm cut into the top surface of the Pz26 transducer. Here, the PDMS, the channel, and the aluminum base are translational invariant along the x direction.

Numerical Optimization of the Thickness of the PDMS Cover
To illustrate how the thickness H pd of the PDMS cover affects the resonances, we vary H pd , and for each value, we sweep the frequency f from 1.5 MHz to 2.5 MHz in steps of 5 kHz to locate the resonances in terms of peaks in the focusing figure of merit F versus f . Such a spectrum is plotted in Figure 5a, where the area of each point is proportional to F . For the same cover thickness H pd = 1.5 mm as in the side-actuated mode "S" of Figure 3 with E ac = 8.9 J/m 3 , we find a resonance mode "A " at f = 2.095 MHz with E ac = 4.6 J/m 3 . An even better resonance mode can be found. In Figure 5a a strong resonance, which depends on the PDMS-cover thickness, is indicated by the blue curve. The optimal cover thickness along this curve is found to be H pd = 80 µm at f = 2.070 MHz, identified as resonance mode "A", which yields an acoustic energy density of 96 J/m 3 , twenty times larger than "A ".  Figure 4, as a function of the actuation frequency f and the thickness H pd of the PDMS cover. The area of the points is proportional to the focusing figure of merit F defined in Equation (15). The "A " marks the resonance mode for H pd = 1.5 mm as in mode "S" of Figure 3. The "A" marks the resonance for the optimal cover thickness H pd = 80 µm. The blue curve indicates resonances sensitive to H pd . (b-d) Color and vector plot of the measured (left column) and simulated (right column) acoustophoretic particle velocity v pa of mode "A ", for particle diameters 4.8, 2.0, and 1.0 µm, respectively. The measured resonance frequency was f = 2.052 MHz and the simulated was found at f = 2.095 MHz (2% higher). All color plots range from 0 µm/s (white) up to v max pa (red) with values 38 or 19 µm/s as indicated in each panel.
In Figure 5b-d are shown experimental and simulated acoustophoretic particle velocities v pa for three different particle diameters 2a = 4.8, 2.0 and 1.0 µm in mode "A ". The magnitude v pa = 38 µm/s for a = 2.4 µm in Figure 5b is similar to that for mode "S", v pa = 40 µm/s in Figure 3c. As for mode "S" in Figure 3, also Figure 5 shows the well-known cross-over from radiation-force-dominated focusing motion of the largest particles to the streaming-roll-dominated motion of the smallest particles. Particular to the single-PDMS-cover device is, firstly, that in Figures 3c and 5d there is only one pair of counter-rotating vortices near the bottom aluminum wall, the other pair near the top-PDMS-cover is suppressed as anticipated, and secondly, very strong localized vortices appear in the top corners of the channel where the PDMS-cover joins the aluminum base.
For both mode "S" in Figure 3 and the anti-symmetric mode "A " in Figure 5b-d there is good qualitative agreement between particle velocities v pa and resonance frequencies f measured in 3D and simulated in 2D assuming translational invariance. Quantitatively, the agreement in f is better for mode "A " (+2%) than for mode "S" (−12%), presumably because the translational invariance is more severely broken in the latter case. In both cases the velocity magnitudes agrees within 30-50%. The skew (white) zero-velocity band observed in Figure 5b (left) indicates that the experimental device was not completely symmetric as in the model Figure 5b (right). The agreement between measured and simulated particle velocities v pa implies that our model can provide reliable estimates of the acoustic energy densities E ac in the devices.

The Role of Variations in the PDMS Material Properties
We suspect that the strong top-corner vortices appearing in Figures 3 and 5 are due to the nearly incompressible nature of PDMS, or equivalently, the very low transverse speed of sound in PDMS. To investigate this hypothesis, we define an artificial polymer alloy of PDMS and Poly(Methyl MethAcrylate), PDMS x PMMA 1−x with mixing ratios x. In Figure 6 we study the four mixing ratios x = 1.0, 0.99, 0.95, and 0 in covers 80 µm thick. For each case, the resonance peak with the largest figure of merit F was found, and in close-up views of the channel are shown the corresponding particle velocity v pa in the water and the displacement u in the surrounding solid.
In the case of a pure PDMS cover Figure 6a, a large displacement is narrowly confined to the top corners, this was also evident in Figures 3 and 5, where the maximum particle velocity was seen near the top corners. Already at 99% PDMS and 1% PMMA, Figure 6b, the transverse speed of sound in the polymer is two times greater than in pure PDMS, and the displacement field at the polymer-water interface is now more evenly distributed. This is even more pronounced for the case of 95% PDMS and 5% PMMA, Figure 6c, similarly for pure PMMA, Figure 6d. Note that in the latter case, the displacement in the cover is also resonating close to its transverse half-wave frequency.
In conclusion, the low transverse speed of sound in PDMS seems to imply the appearance of strong streaming vortices localized near the PDMS-aluminum corners of the device. the resonance frequency f res is noted together with the average acoustic energy density E ac , and the transverse wavelength λ tr in the artificial PDMS-PMMA polymer cover. The color plots indicates the particle velocity magnitude v pa in the water ranging from 0 (white) to the maximum value v max pa = 258 µm/s (red), and of the displacement u in surrounding solids ranging from 0 (dark blue) to 70 nm (yellow). The deformation is scaled 500 times to be visible.

Modeling Dual-PDMS-Cover Aluminum Devices with Anti-symmetric Voltage Actuation
In Figures 3c,d and 5d showing the acoustophoretic velocity fields for small microparticles suspended in the single-PDMS-cover device, we notice that whereas the usual Rayleigh flow rolls are present near the bottom aluminum wall, they seem to be suppressed near the PDMS cover. We therefore study the dual-PDMS-cover device sketched in Figure 7a  In Figure 8, we study if indeed the streaming flow rolls are suppressed in the dual-PDMS-cover device compared to the single-PDMS-cover device. We compute the position of 400 1-µm-diameter polystyrene particles at discrete time steps (black dots) starting at a uniform distribution at time t = 0, where the acoustic field is turned on, and ending at time t focus = 3 2Φ c 2 fl ω 2 a 2 η fl E ac (red dots), the so-called focusing time determined by the acoustic contrast factor Φ given in Table A3 in Appendix A. Inertial effects can be ignored since the Reynolds number Re pa = ρ fl av pa η fl < 0.03 for the given particles and acoustic fields.
Comparing the particle trajectories of mode "A" in the single-PDMS-cover device Figure 8a to those of mode "B" in the dual-PDMS-cover device Figure 8b, we clearly see that in the latter case the streaming flow rolls are suppressed. The dual-PDMS-cover device is therefore predicted to be a good candidate system for controlled focusing of sub-micrometer-sized particles. For this to work, it is of course essential that the focusing time is sufficiently short. For the given device actuated at V 0 = 10 V, the focusing time of the 1-µm-diameter particles is prohibitively long, namely t focus = 23.6 s. Raising the driving voltage by a factor √ 20 ≈ 4.5 to V 0 = 45 V would lower the focusing time to t focus ≈ 1 s, because t focus ∝ V

Conclusions
We have developed a model for analyzing the single-PDMS-cover aluminum-base device with side actuation, recently introduced by Gautam et al. [20]. The model, currently restricted to the case of a constant 2D cross section in a translational invariant device, is validated experimentally with fair qualitative and quantitative agreement by fabricating and characterizing two types of single-PDMS-cover aluminum-base devices: One which is actuated with a symmetric ac-voltage on a Pz26 transducer placed at the side of the channel, and another with an anti-symmetric ac-voltage on a transducer placed right under the channel. Both numerical simulations and experiments support our hypothesis that using a soft PDMS cover of the acoustophoresis channel, the boundary driven acoustic streaming is suppressed in the bulk. The developed model can thus predict the streaming patterns in such devices, and we subsequently used it to show three aspects: (1) The incompressible nature of the soft PDMS cover introduces strong streaming rolls confined near the corners where the PDMS cover joins the aluminum base, while maintaining the conventional large Rayleigh streaming rolls extending from the aluminum-water interface; (2) An optimal thickness of the PDMS cover can be determined by simulation; (3) In devices with a dual-PDMS cover, the model predicts that the conventional Rayleigh streaming flow rolls should be suppressed and changed into vortices confined near the corners of the channel. Experimental work is in progress to verify these predictions.

Appendix A. Material Parameters
The following three tables contain the parameter values used in the numerical modeling of the PDMS-cover aluminum-based devices.    Water [32] Polystyrene Particles [