Three-dimensional numerical modeling of surface acoustic wave devices: Acoustophoresis of micro- and nanoparticles including streaming

TL;DR

This paper introduces a comprehensive 3D numerical model for surface acoustic wave devices, capturing electromechanical, acoustic, and streaming effects to better understand and optimize acoustofluidic particle manipulation.

Contribution

It presents the first fully 3D simulation including electromechanical, acoustic, and streaming effects, validated against experimental data for SAW microfluidic devices.

Findings

Accurately predicts acoustic resonances and particle responses.

Explains differences between soft and hard lid device behaviors.

Reproduces observed streaming flow patterns.

Abstract

Surface acoustic wave (SAW) devices form an important class of acoustofluidic devices, in which the acoustic waves are generated and propagate along the surface of a piezoelectric substrate. Despite their wide-spread use, only a few fully three-dimensional (3D) numerical simulations have been presented in the literature. In this paper, we present a 3D numerical simulation taking into account the electromechanical fields of the piezoelectric SAW device, the acoustic displacement field in the attached elastic material, in which the liquid-filled microchannel is embedded, the acoustic fields inside the microchannel, as well as the resulting acoustic radiation force and streaming-induced drag force acting on micro- and nanoparticles suspended in the microchannel. A specific device design is presented, for which the numerical predictions of the acoustic resonances and the acoustophoretic repsonse of suspended microparticles in 3D are successfully compared with experimental observations. The simulation provides a physical explanation of the the observed qualitative difference between devices with an acoustically soft and hard lid in terms of traveling and standing waves, respectively. The simulations also correctly predict the existence and position of the observed in-plane streaming flow rolls. The presented simulation model may be useful in the development of SAW devices optimized for various acoustofluidic tasks.

Tables

Table 1. Table 1: Elasticity constants C α ​ β subscript 𝐶 𝛼 𝛽 C_{\alpha\beta} , mass density ρ sl subscript 𝜌 sl \rho_{\mathrm{sl}} , piezoelectric coupling constants e i ​ α subscript 𝑒 𝑖 𝛼 e_{i\alpha} and relative dielectric constants ε i ​ k subscript 𝜀 𝑖 𝑘 \varepsilon_{ik} of materials used in this work. 128 ∘ YX-cut lithium niobate values are defined in the global system x , y , z 𝑥 𝑦 𝑧 x,y,z , for derivations see Appendix  A . Note that C 12 = C 11 − 2 ​ C 44 subscript 𝐶 12 subscript 𝐶 11 2 subscript 𝐶 44 C_{12}=C_{11}-2C_{44} for isotropic materials (Pyrex and PDMS).

Parameter	Value	Parameter	Value
128∘ YX-cut lithium niobate Weis and Gaylord (1985)
$C_{11}$	202.89 GPa	$C_{12}$	72.33 GPa
$C_{13}$	60.17 GPa	$C_{14}$	10.74 GPa
$C_{22}$	194.23 GPa	$C_{23}$	90.59 GPa
$C_{24}$	8.97 GPa	$C_{33}$	220.29 GPa
$C_{34}$	8.14 GPa	$C_{44}$	74.89 GPa
$C_{55}$	72.79 GPa	$C_{56}$	$-$8.51 GPa
${\rho }_{\mathrm{sl}}$	$4628\text{kg}\text{m}\text{<sup><span>-3</span></sup>}$	$C_{66}$	59.51 GPa
$e_{15}$	1.56 $\text{C}{\text{m}}^{-2}$	$e_{16}$	-4.23 $\text{C}{\text{m}}^{-2}$
$e_{21}$	$-$1.73 $\text{C}{\text{m}}^{-2}$	$e_{22}$	4.48 $\text{C}{\text{m}}^{-2}$
$e_{23}$	$-$1.67 $\text{C}{\text{m}}^{-2}$	$e_{24}$	0.14 $\text{C}{\text{m}}^{-2}$
$e_{31}$	1.64 $\text{C}{\text{m}}^{-2}$	$e_{32}$	$-$2.69 $\text{C}{\text{m}}^{-2}$
$e_{33}$	2.44 $\text{C}{\text{m}}^{-2}$	$e_{34}$	0.55 $\text{C}{\text{m}}^{-2}$
${\epsilon }_{11}$	44.30	${\epsilon }_{22}$	38.08
${\epsilon }_{23}$	$-$7.96	${\epsilon }_{33}$	34.12
Pyrex Narottam P. Bansal (1986)
$C_{11}$	69.73 GPa	$C_{12}$	17.45 GPa
${\rho }_{\mathrm{sl}}$	$2230\text{kg}\text{m}\text{<sup><span>-3</span></sup>}$	$C_{44}$	26.14 GPa
$\epsilon$	4.6	${\Gamma }_{\mathrm{sl}}$	0.0002
PDMS Madsen (1983); Zell et al. (2007); MIT
$C_{11}$	1.13 GPa	$C_{12}$	1.11 GPa
${\rho }_{\mathrm{sl}}$	$1070\text{kg}\text{m}\text{<sup><span>-3</span></sup>}$	$C_{44}$	0.011 GPa
$\epsilon$	2.5	${\Gamma }_{\mathrm{sl}}$	0.0213

Table 2. Table 2: Material parameters of water from Ref.  Muller and Bruus ( 2014 ) .

Parameter	Symbol	Value
Speed of sound	$c_{\mathrm{fl}}$	1497 $\text{m}{\text{s}}^{-1}$
Mass density	${\rho }_{\mathrm{fl}}$	997 $\text{kg}{\text{m}}^{-3}$
Dynamic viscosity	${\eta }_{\mathrm{fl}}$	0.89 mPa s
Bulk viscosity	${\eta }_{\mathrm{fl}}^b$	2.485 mPa s
Compressibility	${\kappa }_{\mathrm{fl}}$	452 ${\text{TPa}}^{-1}$

Table 3. Table 3: Dimensions of the numeric 2D and 3D models.

Parameter	Symbol	2D	3D	Unit
Device depth ($y$)	$L_{\mathrm{sl}}$	-	1200	µm
Solid height ($z$)	$H_{\mathrm{sl}}$	40-1000	500	µm
Solid width ($x$)	$W_{\mathrm{sl}}$	200	80	µm
Channel height	$H_{\mathrm{fl}}$	50-200	50	µm
Channel width	$W_{\mathrm{fl}}$	3500	900	µm
Piezo height	$H_{\mathrm{pz}}$	100-500	300	µm
PML length	$L_{\mathrm{PML}}$	80	80	µm
Electrode depth ($y$)	$L_{\mathrm{el}}$	-	400	µm
Electrode height ($z$)	$H_{\mathrm{el}}$	0.4	0.4	µm
Electrode width ($x$)	$W_{\mathrm{el}}$	20	20	µm
Electrode gap	$G_{\mathrm{el}}$	20	20	µm
SAW wavelength	${\lambda }_{\mathrm{SAW}}$	80	80	µm
No. of electrode pairs	$n_{\mathrm{el}}$	24	4	-
No. of reflectors	$n_{\mathrm{rf}}$	0-6	0	-
Actuation frequency	$f_0$	30-60	50	MHz
Driving voltage	$V_0$	1	1	V
Degrees of freedom	$n_{\mathrm{DOF}}$	$𝒪\left({10}^5\right)$	$𝒪\left({10}^6\right)$	-
Memory requirements	$R$	$𝒪\left(10\right)$	$𝒪\left({10}^3\right)$	GB

Table 4. Table 4: The devices D1 and D2, used in the experimental validation of the numerical model, differs by the choice of lid. The other parameters of D1 and D2 are listed in Table  3 .

Device	Lid material	Lid thickness
D1	PDMS	15 mm
D2, see Fig. 1(a)	Pyrex	0.45 mm

Table 5. Table 5: Measured and simulated values of the frequencies f 𝑓 f near the ideal (unloaded) frequency f SAW = 49.9 subscript 𝑓 SAW 49.9 f_{\mathrm{SAW}}=49.9  MHz, where | Z el ​ ( f ) | superscript 𝑍 el 𝑓 |Z^{\mathrm{el}}(f)| and ψ ​ ( f ) 𝜓 𝑓 \psi(f) have local minima and maxima in Pyrex device D2.

Extremum	$f_{\exp }$	$f_{\mathrm{num}}$	Relative error
	[GHz]	[GHz]	[%]
A	$47.35$	$48.25$	1.9
B	$48.05$	$49.00$	2.0
C	$49.40$	$50.25$	1.7
D	$48.65$	$47.50$	1.7
E	$50.00$	$50.75$	1.7
F	$47.70$	$48.75$	2.2
G	$49.10$	$50.00$	1.8
H	$45.50$	$46.00$	1.1
I	$49.40$	$48.50$	1.8

Equations

$${ϵ}_{ik}^{}$$

$${\sigma }_{ik}$$

$${\mathbf{ϵ}}_{\mathrm{V}}^{}$$

$${\sigma }_{\alpha }=C_{\alpha \beta }{ϵ}_{\beta },$$

$${\nabla }_{\mathrm{V}}=\left(\begin{array}{cccccc}{\partial }_x& 0& 0& 0& {\partial }_z& {\partial }_y\\ 0& {\partial }_y& 0& {\partial }_z& 0& {\partial }_x\\ 0& 0& {\partial }_z& {\partial }_y& {\partial }_x& 0\end{array}\right).$$

$${\nabla }_{\mathrm{V}}\cdot {\mathbf{\sigma }}_{\mathrm{V}}$$

$$\nabla \cdot \mathbf{D}$$

$${\mathbf{\sigma }}_{\mathrm{V}}$$

$$\mathbf{D}$$

$$\left(\!\!\begin{array}[]{c}\sigma_{1}\\ \sigma_{2}\\ \sigma_{3}\\ \hline\cr\sigma_{4}\\ \sigma_{5}\\ \sigma_{6}\\ \hline\cr D_{x}\\ D_{y}\\ D_{z}\\ \end{array}\!\!\right)\!\!=\!\!\left(\!\!\begin{array}[]{c@{\:}c@{\:}c@{}|c@{}c@{}c@{}|c@{}c@{}c}C_{11}\hfil\>&C_{12}\hfil\>&C_{13}&C_{14}&0&0&0&\text{-}e_{21}&\text{-}e_{31}\\ C_{12}\hfil\>&C_{22}\hfil\>&C_{23}&C_{24}&0&0&0&\text{-}e_{22}&\text{-}e_{32}\\ C_{13}\hfil\>&C_{23}\hfil\>&C_{33}&C_{34}&0&0&0&\text{-}e_{23}&\text{-}e_{33}\\ \hline\cr C_{14}\hfil\>&C_{24}\hfil\>&C_{34}&\!C_{44}&0&0&0&\text{-}e_{24}&\text{-}e_{34}\\ 0\hfil\>&0\hfil\>&0&0&C_{55}&C_{56}&\text{-}e_{15}&0&0\\ 0\hfil\>&0\hfil\>&0&0&C_{56}&C_{66}&\text{-}e_{16}&0&0\\ \hline\cr 0\hfil\>&0\hfil\>&0&0&e_{15}&e_{16}&\varepsilon_{11}&0&0\\ e_{21}\hfil\>&\!e_{22}\hfil\>&\!e_{23}&\!e_{24}&0&0&0&\varepsilon_{22}&\varepsilon_{23}\\ e_{31}\hfil\>&e_{32}\hfil\>&e_{33}&\!e_{34}&0&0&0&\varepsilon_{23}&\varepsilon_{33}\\ \end{array}\!\!\right)\!\!\!\left(\!\!\!\!\begin{array}[]{c}\epsilon_{1}\\ \epsilon_{2}\\ \epsilon_{3}\\[-4.2679pt] \rule{22.76219pt}{0.42677pt}\\[-1.42262pt] \epsilon_{4}\\ \epsilon_{5}\\ \epsilon_{6}\\[-8.53581pt] \rule{22.76219pt}{0.42677pt}\\[-1.42262pt] E_{x}\\ E_{y}\\ E_{z}\\ \end{array}\!\!\!\!\right)\!,$$

$$\left(\!\!\begin{array}[]{c}\sigma_{1}\\ \sigma_{2}\\ \sigma_{3}\\ \hline\cr\sigma_{4}\\ \sigma_{5}\\ \sigma_{6}\end{array}\!\!\right)\!\!=\!\!\left(\!\!\begin{array}[]{c@{\:}c@{\:}c@{}|c@{}c@{}c}C_{11}\hfil\>&C_{12}\hfil\>&C_{12}&0&0&0\\ C_{12}\hfil\>&C_{11}\hfil\>&C_{12}&0&0&0\\ C_{12}\hfil\>&C_{12}\hfil\>&C_{11}&0&0&0\\ \hline\cr 0\hfil\>&0\hfil\>&0&C_{44}&0&0\\ 0\hfil\>&0\hfil\>&0&0&C_{44}&0\\ 0\hfil\>&0\hfil\>&0&0&0&C_{44}\\ \end{array}\!\!\right)\!\!\!\left(\!\!\!\!\begin{array}[]{c}\epsilon_{1}\\ \epsilon_{2}\\ \epsilon_{3}\\[-5.69054pt] \rule{22.76219pt}{0.42677pt}\\[-1.42262pt] \epsilon_{4}\\ \epsilon_{5}\\ \epsilon_{6}\end{array}\!\!\!\!\right)\!,$$

$$\displaystyle\bm{\nabla}\cdot\big{(}\bm{\nabla}\!p_{1}\big{)}$$

$${\mathbf{v}}_1$$

$${\varphi }_{\mathrm{ce}}=V_0{\mathrm{e}}^{-\mathrm{i}\omega t},{\varphi }_{\mathrm{ge}}=0,$$

$$(\mathbf{I}-\mathbf{n}\mathbf{n})\cdot \nabla {\varphi }_{\mathrm{fe}}=0,$$

$${\mathbf{\sigma }}_{\mathrm{sl}}\cdot \mathbf{n}$$

$$\mathbf{n}\cdot {\mathbf{v}}_1$$

$$\text{with shear wavenumber }{k}_s=\frac{1+\mathrm{i}}{\delta }.$$

$$\mathbf{\sigma }\cdot \mathbf{n}=0.$$

$$\displaystyle\bm{\nabla}\cdot\bm{\sigma}_{2}-\rho_{0}\bm{\nabla}\cdot\big{\langle}\bm{v}_{1}\bm{v}_{1}\big{\rangle}$$

$$\displaystyle\bm{\nabla}\cdot\big{(}\rho_{0}{\bm{v}_{2}}+\big{\langle}\rho_{1}\bm{v}_{1}\big{\rangle}\big{)}$$

$$\bm{\sigma}_{2}=-{p_{2}}\bm{I}+\eta\,\big{[}\bm{\nabla}\!{\bm{v}_{2}}+\left(\bm{\nabla}\!{\bm{v}_{2}}\right)^{\mathrm{T}}\big{]}+(\beta-1)\,\eta\,(\bm{\nabla}\cdot{\bm{v}_{2}})\,\bm{I}.\\$$

$${\mathbf{v}}_2$$

$$\mathbf{A}$$

$$\displaystyle\hskip 42.67912pt+\bigg{[}\frac{2-\mathrm{i}}{2}\bm{\nabla}\!\!\cdot\!\bm{v}^{\delta 0*}_{1}+\mathrm{i}\Big{(}\bm{\nabla}\!\!\cdot\!\bm{v}_{\mathrm{sl}}^{*}-\partial_{\zeta}v_{1\zeta}^{*}\Big{)}\bigg{]}\bm{v}^{\delta 0}_{1}\bigg{\}},$$

$$\mathbf{B}$$

$$\frac{4\pi }{3}{a}_{\mathrm{pt}}^3$$

$${\mathbf{F}}^{\mathrm{rad}}$$

$${\mathbf{F}}^{\mathrm{drag}}$$

$${\mathbf{f}}^{\mathrm{rad}}=\frac{3}{4\pi a^3}{\mathbf{F}}^{\mathrm{rad}}.$$

$$\partial x_i\to \partial {\tilde{x}}_i$$

$$\mathrm{d}{x}_i\to \mathrm{d}{\tilde{x}}_i$$

$$s(\mathbf{r})$$

$${\partial }_x{u}_x$$

$$u_z$$

$$\varphi$$

$$p_1$$

$$Z^{\mathrm{el}}=\big{|}Z^{\mathrm{el}}\big{|}\>e^{-\mathrm{i}\psi}=\frac{V_{0}}{I},$$

$$\mathbf{P}$$

$${\mathbf{J}}_{\mathrm{pol}}$$

$$I={\int }_{\partial {\Omega }_{\mathrm{ce}}}{\mathbf{J}}_{\mathrm{pol}}\cdot \mathbf{n}\mathrm{d}A.$$

$$\big{|}Z^{\mathrm{el}}\big{|}=\left|\frac{V_{0}}{I}\right|,\qquad\psi=\arg\left(\frac{V_{0}}{I}\right).$$

Full text

I Introduction

During the past decade, surface acoustic wave (SAW) devices have been developed for a multitude of different types of acoustofluidic handling of micrometer-sized particles inside closed microchannels. Examples include acoustic mixing Sritharan et al. (2006), continuous particle or droplet focusing Shi et al. (2008); Franke et al. (2009) and separation Tan et al. (2009); Shi et al. (2009a), single-particle handling Ding et al. (2012a); Tran et al. (2012), acoustic tweezing Shi et al. (2009b); Collins et al. (2016a); Riaud et al. (2017), two-dimensional single patterning Collins et al. (2015, 2018), on-chip studies of microbial organisms Zhou et al. (2017); Zhang et al. (2019), and non-trivial electrode shapes to generate chirped, focused, and rotating acoustic waves Ding et al. (2012b); Riaud et al. (2015); Collins et al. (2016b); Riaud et al. (2017).

The development of effective handling of submicrometer-sized particles has been less successful. It remains a challenge to handle this in biotechnology highly important class of particles including small bacteria, exosomes, and viruses. Could these particles be handled in a controlled way, it would be of particular interest for developing new and more efficient diagnostics Liga et al. (2015). The first steps towards acoustofluidics handling of nanometer-sized particles have been taken relying on acoustic streaming effects with both bulk acoustic waves (BAW) Antfolk et al. (2014) and SAW Mao et al. (2017), or using seed particles to enhance acoustic trapping in BAW devices Hammarström et al. (2012). However, these methods have a low selectivity. However, recently SAW devices have been developed to focusing nanoparticles Collins et al. (2017) and separation of nanoparticles Sehgal and Kirby (2017); Wu et al. (2017). In particular Sehgal and Kirby Sehgal and Kirby (2017) demonstrated separation between 100- and 300-nm-diameter particles on the proof-of-concept stage. To fully utilize the potential of this and similar devices, further development is necessary to increase the efficiency and sorting flow rates. Here, numerical simulations may play a crucial role, both in improving the understanding of the underlying physical acoustofluidic processes, and to ease the cumbersome development cycle consisting of an iterative series of creating, fabricating, and testing device designs.

An increasing amount of numerical studies include piezoelectric dynamics in two-dimensional (2D) models Johansson et al. (2012); Garofalo et al. (2017); Darinskii et al. (2016); Tan et al. (2010), but mostly the piezoelectric transducers are introduced in numeric models in the form of analytic approximations Köster (2007); Zhang et al. (2018); Lei et al. (2013); Ley and Bruus (2017); Nama et al. (2015); Vanneste and Bühler (2011), and designs are often based on a priori knowledge of the piezoelectric effect in the unloaded substrates typically applied in telecommunication. In acoustofluidic devices, the acoustic impedance of the contacting fluid is much closer to that of the substrate causing waves to behave much differently from those in telecommunications devices. It is thus prudent to include the piezoelectric effect and the coupling between the fluid and substrate in numeric models to accurately describe the device behavior. Additionally, three-dimensional (3D) simulations in the literature are scarce, but they are essential for making full-device acoustophoresis predictions as many actual acoustofluidic devices do exhibit non-trivial features in 3D due to asymmetric and intricate shapes of electrodes and channels.

In this paper, we present 3D numerical simulations taking into account the electromechanical fields of the piezoelectric SAW device, the acoustic displacement field in the attached elastic material, in which the liquid-filled microchannel is embedded, the acoustic fields inside the microchannel, as well as the resulting acoustic radiation force and streaming-induced drag force acting on microparticles suspended in the microchannel. The model is validated experimentally with devices based on the SAW device described by Sehgal and Kirby Sehgal and Kirby (2017). In Section II we describe the physical model system representing the SAW device and state the governing equations, and in Section III we treat the implementation of the model system in a weak-form, finite-element model. The results of the model in reduced 2D and in full 3D are presented in Sections V and VI, and finally in Sections VII and VIII we discuss our findings and summarize our conclusions.

II The model SAW system and the governing equations

The model SAW system is shown in Fig. 1a. Essentially, it consists of a piezoelectric lithium niobate substrate with a specific interdigitated transducer (IDT) metal-electrode configuration on the surface. On top of the substrate a microfluidic channel is defined in an elastic material, either the acoustically soft rubber polydimethylsiloxane polymer (PDMS) or the acoustically hard borosilicate glass (Pyrex).

We follow Sehgal and Kirby Sehgal and Kirby (2017) and place the IDT electrodes directly underneath the microchannel and choose the periodicity of the electrode pattern to result in a SAW wavelength $\lambda_{\mathrm{SAW}}=80~{}\textrm{\textmu{}m}$ and a (unloaded) resonance frequency $f_{\mathrm{SAW}}=c_{\mathrm{SAW}}/\lambda_{\mathrm{SAW}}=(3995~{}\textrm{m}/\textrm{s})/(80~{}\textrm{\textmu{}m})=49.9$ MHz. The driving electrodes are flanked by Bragg-reflector electrodes to (partially) reflect the outgoing SAWs traveling along the surface from the driving electrodes. As described in more detail in Appendix A, the lattice coordinate system $X,Y,Z$ of the 128∘ YX-cut lithium niobate wafer is rotated the usual $38^{\circ}=128^{\circ}-90^{\circ}$ about the $x$-axis to obtain an optimal SAW configuration.

To facilitate separation of nanoparticles, the axis of the microchannel is tilted $10^{\circ}$ angle relative to the IDT electrodes. At both ends, the microchannel branches out in a number of side channels with vertical openings for inlet and outlet tubing. In the numerical model, this inlet/outlet structure is represented by ideally absorbing boundary conditions.

The SAW device is actuated by a time-harmonic voltage difference at frequency $f$ applied to the IDT electrodes. The corresponding angular frequency is $\omega=2\pi f$.

The following formulation of the governing equations, is a further development of our previous work presented in Refs. Skov and Bruus (2016); Ley and Bruus (2017); Skov et al. (2019) to take into account SAW in 3D models of lithium-niobate-driven ultrasound acoustics in liquid-filled microchannels.

II.1 The Voigt notation for elastic solids

In linear elastodynamics with the elasticity tensor $C_{iklm}$, the stress $\sigma_{ik}$ and strain $\epsilon^{\mathrm{}}_{ik}$ tensors with $i,k=1,2,3$ (or $x$, $y$, $z$) are defined in index notation as

[TABLE]

In the Voigt notation (subscript V) Auld (1990), the symmetric stress and strain double-index tensor components $\sigma_{ik}=\sigma_{ki}$ and $\epsilon^{\mathrm{}}_{ik}=\epsilon_{ki}$ are organized in single-index vectors $\sigma_{\alpha}$ and $\epsilon_{\alpha}$ with $\alpha=1,2,\ldots,6$, as,

[TABLE]

and the stress-strain relation is written,

[TABLE]

where $C_{\alpha\beta}$ is the 6$\times$6 Voigt elasticity matrix. We also introduce the 3$\times$6 Voigt matrix gradient operator $\bm{\nabla}_{\mathrm{V}}$,

[TABLE]

The equations governing the device are divided into three sets. One set is the first-order time-harmonic equations for the acoustic fields, the second set contains the steady time-averaged second-order fields, and the third set are the time-dependent equations describing the acoustophoretic motion of suspended particles.

II.2 The time-harmonic first-order fields

By construction, all first-order fields are proportional to the time-harmonic electric potential actuating the SAW device at angular frequency $\omega$. Consequently, all first-order fields are time-harmonic acoustic fields of the form $\hat{g}(\bm{r},t)=g(\bm{r})\>\mathrm{e}^{-\mathrm{i}\omega t}$, where $g(\bm{r})$ is the complex-valued field amplitude. The corresponding physical field is the real part \operatorname{Re}\!\big{[}\hat{g}(\bm{r},t)\big{]}. All terms thus have the same explicit time dependence $\mathrm{e}^{-\mathrm{i}\omega t}$, so this factor is divided out, leaving us with the governing equations for the amplitude $g$, where we for brevity suppress the spatial argument $\bm{r}$.

In a linear piezoelectric material with a mass density $\rho_{\mathrm{sl}}$ and no free charges, the solid displacement field $\bm{u}$ and the electric potential field $\phi$ are governed by the Cauchy equation and Gauss’s law,

[TABLE]

Here, $\varepsilon_{0}$ is the vacuum permittivity, $\bm{\varepsilon}_{\mathrm{r}}$ is the relative permittivity tensor of the material, and superscript “T” denotes the transpose of a matrix, see Table 1.

For anisotropic lithium niobate, Eqs. (5a) and (5b) are turned into equations for $\bm{u}$ and $\phi$ by using the explicit form of Eqs. (5c) and (5d) written as the coupling-matrix,

[TABLE]

For isotropic elastic solids with no charges and no piezoelectric coupling $\bm{e}=\bm{0}$, only Eq. (5a) is relevant, and it becomes an equation for $\bm{u}$, as Eq. (5c) reduces to

[TABLE]

with only two independent elastic constants, $C_{11}$ and $C_{44}$, because $C_{12}=C_{11}-2C_{44}$ for isotropic material.

In a fluid with speed of sound $c_{\mathrm{fl}}$, mass density $\rho_{\mathrm{fl}}$, dynamic viscosity $\eta_{\mathrm{fl}}$, viscous boundary layer thickness $\delta=\sqrt{\frac{2\eta_{\mathrm{fl}}}{\rho_{\mathrm{fl}}\omega}}$, viscosity ratio $\beta=\frac{\eta_{\mathrm{fl}}^{\mathrm{b}}}{\eta_{\mathrm{fl}}}+\frac{1}{3}$, and effective damping coefficient $\Gamma_{\mathrm{fl}}=\frac{1+\beta}{2}(k_{0}\delta)^{2}$, the first-order pressure field ${p_{1}}$ is governed by the Helmholtz equation, and the acoustic velocity field $\bm{v}_{1}$ is given by the pressure gradient,

[TABLE]

where $k_{c}$ is the weakly damped compressional wave-numberBach and Bruus (2018). See Table 2 for parameter values.

Turning to the boundary conditions, we introduce $\bm{n}$ as the normal vector for a given surface. The SAW device in Fig. 1 is actuated by a time-harmonic potential of amplitude $V_{0}$ on the surfaces of the charged electrodes (ce) and 0 V on the grounded electrodes (ge), respectively,

[TABLE]

At a given fluid-solid interface we impose the usual continuity conditions Ley and Bruus (2017) with the recently developed boundary-layer corrections included Bach and Bruus (2018): the solid stress $\bm{\sigma}_{\mathrm{sl}}$ is given by the acoustic pressure $p_{1}$ with the addition of the boundary-layer stress, and the fluid velocity $\bm{v}_{1}$ is given by the solid-wall velocity $\bm{v}_{\mathrm{sl}}=-\mathrm{i}\omega\bm{u}$ with the addition of the boundary-layer velocity $\bm{v}_{\mathrm{sl}}-\bm{v}_{1}$,

[TABLE]

The terms containing the shear wavenumber $k_{s}$ represent the corrections arising from taking the 400-nm wide, viscous boundary layer into account analytically Bach and Bruus (2018).

All exterior solid surfaces facing the air have a stress-free boundary condition prescribed,

[TABLE]

This is a good approximation because the surrounding air has an acoustic impedance 3 to 4 orders of magnitude lower than that of the solids causing 99.99 % of incident acoustic waves from the solid to be reflected. Moreover the shear stress from the air is negligible.

II.3 The time-averaged second-order fields

The slow timescale or steady fields in the fluid are the time-averaged second-order velocity ${\bm{v}_{2}}$ and pressure ${p_{2}}$ field. These are governed by the time-averaged momentum and mass-conservation equations,

[TABLE]

Along a fluid-solid interface with tangential vectors $\bm{e}_{\xi}$ and $\bm{e}_{\eta}$ and the normal vector $\bm{e}_{\zeta}=\bm{n}$, we use for ${\bm{v}_{2}}$ the effective boundary condition derived in Ref. Bach and Bruus (2018). Here, the viscous boundary layer is taken into account analytically by introducing the boundary-layer velocity field $\bm{v}^{\delta 0}=\bm{v}_{\mathrm{sl}}-\bm{v}_{1}$ in the fluid along the fluid-solid interface,

[TABLE]

where the asterisk denotes complex conjugation.

II.4 Acoustophoresis of suspended particles

To predict the acoustophoretic motion of a dilute suspension of spherical micro- and submicrometer-sized particles in the fluid of density $\rho_{\mathrm{fl}}$, compressibility $\kappa_{\mathrm{fl}}$, and viscosity $\eta_{\mathrm{fl}}$, we implement a particle tracing routine in the model. We consider Newton’s second law for a single spherical particle of radius $a_{\mathrm{pt}}$ and density $\rho_{\mathrm{pt}}$ moving with velocity $\bm{v}_{\mathrm{pt}}$ under the influence of gravity $\bm{g}$, the acoustic radiation force $\bm{F}^{\mathrm{rad}}$ Settnes and Bruus (2012), and the Stokes drag force $\bm{F}^{\mathrm{drag}}$ Bruus (2011) induced by acoustic streaming of the fluid,

[TABLE]

By direct time integration of Eq. (13a) applied to a set of particles initially placed on a square grid, the acoustophoretic motion of the particles can be predicted and compared to the experimentally observed one. We note that gravity effects are negligible as \rho_{\mathrm{pt}}\bm{g}~{}\ll~{}\kappa_{\mathrm{fl}}\big{\langle}(f_{0}p_{1})\bm{\nabla}\!p_{1}\big{\rangle}.

III Numerical implementation

Inspired by our previous experimental work, Sehgal and Kirby Sehgal and Kirby (2017), we study the SAW test system shown in Fig. 1 with actuating electrodes and Bragg-reflector electrodes placed directly underneath the microchannel. The parameter values used in the numerical simulation are listed in Table 3, and a sketch of the vertical cross section of the test system is shown in Fig. 2. Note that the SAW wavelength $\lambda_{\mathrm{SAW}}$ is set by the IDT electrode geometry as $\lambda_{\mathrm{SAW}}=2(W_{\mathrm{el}}+G_{\mathrm{el}})$. We study microcavities defined in either acoustically soft PDMS, see Fig. 2(a), or the acoustically hard borosilicate glass (Pyrex), see Fig. 2(b), and we perform numerical simulation in both 2D and 3D.

Following the procedure of our previous numerical simulations Ley and Bruus (2017); Skov et al. (2019), the coupled governing equations from Sections II.2-II.4 are implemented in the finite-element-method software COMSOL Multiphysics 5.3a COM , using the weak-form partial differential equation interface “PDE Weak Form” in the mathematics module. For a given driving voltage $V_{0}$, actuation frequency $f$, and angular frequency $\omega=2\pi f$ specified in the actuation boundary condition (8a), the numeric model is solved in three sequential steps: (1) the first-order equations (5) and (7a) presented in Section II.2 for the pressure $p_{1}$, displacement $\bm{u}$, and electric potential $\phi$, together with the corresponding boundary conditions (8)-(10); (2) the steady second-order streaming velocity $\bm{v}_{2}$ in Section II.3 governed by (11) and (12), where time-averaged products of the first-order fields appear as source terms; and (3) the acoustophoretic motion of suspended test particles in Section II.4 found by time integration of Eq. (13). As in previous works Ley and Bruus (2017), we have performed convergence analyses of the model to verify that the model converges towards a single solution as the mesh size decreases.

Simulations of the full 3D model are time and computer-memory consuming. Therefore, part of the analysis has been performed on 2D models to study the resonance behavior of the device and the acoustic radiation force in the vertical $y$-$z$ plane normal to the electrodes in the horizontal $x$-$y$ plane. In these simulations, presented in Section V, it is possible to model a cross-section of the device to scale. To investigate effects that have non-trivial behavior in full 3D, such as the acoustic streaming and the acoustophoretic motion of suspended particles presented in Section V, we must perform full 3D modeling. However, in this case, the extended computer memory requirements has necessitated a scale down of the model. The parameters for the 2D and 3D simulations are listed in Table 3.

III.1 Perfectly matched layers

We reduce the numeric footprint of the model by implementing perfectly matched layers (PMLs) in the model as described by Ley and Bruus Ley and Bruus (2017): Large passive domains surrounding the acoustically active region are replaced by much smaller domains, in which PMLs act as ideal absorbers of out-going acoustic waves thus completely removing reflections. In contrast to Ref. Ley and Bruus (2017), the PMLs in the present model are functions of all three spatial coordinates.

In the small surrounding domains, the PMLs are implemented in the weak-form governing equations by a complex-valued coordinate transformation of the spatial derivatives $\partial x_{i}$ and integral measures $\mathrm{d}x_{i}$ appearing,

[TABLE]

where $s(\bm{r})$ is a real-valued function of position. Here, $s(\bm{r})$ is given for the specific case shown in Fig. 2 with a PML of width $L_{\mathrm{PML},i}$ in the three coordinate directions $i=x,y,z$ placed outside the region $x<x_{0}$, $y<y_{0}$, and $z<z_{0}$, $\Theta(x)$ is the Heaviside step function ($=1$ for $x>0$, and 0 otherwise), and $k_{\mathrm{PML}}$ is an adjustable parameter for the strength of the PML absorbtion. The bottom PML in the niobate substrate is used because SAWs decay exponentially in the depth on the scale of the wavelength, whereas the top and side PMLs are used to mimic attenuating in respective materials over large distances.

III.2 Symmetry planes

As in previous numeric works Muller and Bruus (2015); Ley and Bruus (2017), we use an antisymmetry line to reduce the numerical cost of our 2D models. The antisymmetry line is realized by boundary conditions on the solid displacement, the electric potential, and the fluid pressure along the line,

[TABLE]

We check these conditions against the values along the device centerline in a 2D simulation for a fully symmetric device and observe that they are in good agreement.

In 3D we cannot use symmetry planes, as the device is manifestly asymmetric due to the $10^{\circ}$ angle between the IDT and the walls of the microchannel.

III.3 Embedded electrodes

In the actual device, the $400$-nm thick electrodes protrude into the fluid domain. In our numeric model we simplify the device by submerging them into the substrate to form a planar solid-fluid interface as shown in Fig. 2. Thereby the fluid-solid interface has no sharp corners, at which singularities appear in the numeric gradients. Furthermore, the planar interface mitigates the need for an enormous number of mesh elements ranging from nm to µm in the fluid domain, which would either lower the element quality greatly or add massive computational costs. This reduction in model complexity is justified by the height of the electrodes being less than 1 % of the channel height and having no influence on the pressure acoustics of the system. On the other hand, we cannot completely neglect the electrodes, because jumps in acoustic impedance between the metal electrodes and the niobate substrate cause partial reflections of SAWs running along the substrate. Thus we choose to keep but submerge the electrodes.

IV Experimental methods

To validate the numerical models, we have performed experiments on two type of devices listed in table IV, namely microchannels defined in slabs of either PDMS (D1) or Pyrex (D2) bonded on top of the lithium niobate substrate equipped with the IDT and Bragg reflectors. The PDMS device (D1) is fabricated by standard photolithography techniques listed in our previous work Sehgal and Kirby (2017). The Pyrex device (D2) is fabricated by glass microfabrication techniques, briefly described in the following. A microchannel of desired dimensions is wet-etched in a borosilicate glass wafer by 49% hydrofluoric (HF) acid using a multilayered mask of chrome, gold, and SPR220 photoresist. The input and output ports of the microchannel are obtained from the laser cutting of glass. The bonding between glass microchannel and lithium niobate substrate is achieved by coating a 5 µm layer of SU-8 epoxy on the surface of lithium niobate. The microchannel is gently placed on the uncured SU-8 and the epoxy is baked following standard steps. The SU-8 outside the microchannel region is selectively crosslinked to achieve bonding and the SU-8 inside the microchannel region is dissolved away with a developer, thus obtaining a Pyrex lid microchannel on top of the lithium niobate substrate (D2). The devices are tested with 1.7-µm-diameter fluorescent polystyrene particles (Polysciences, Inc.) that are suspended in de-ionized water (18.2 M$\Omega$/cm, Labconco WaterPro PS) containing 0.7% (w/v) Pluronic F-127 to prevent particle aggregation. The particle solution is injected into the microchannel after priming the devices with 70% ethanol solution to avoid the formation of air bubbles. An ultrasound field is set up in the devices by applying an RF signal at desired frequency to the IDT with a HP 8643A signal generator and an ENI 350L RF power amplifier. The acoustophoretic motion of the tracer particles are visualized on a fixed-stage, upright fluorescent microscope (Olympus BX51WI) with a digital CCD camera (Retiga 1300, Q Imaging). The images are acquired with Q-Capture Pro 7 software and post processed in ImageJ. The electrical impedance of the devices is measured directly from an impedance analyzer (Agilent 4395A).

V Results of the 2D modeling

In the following, we compare the results of the 2D modeling in the vertical $x$-$z$ with experiments carried out on the two devices D1 and D2 listed in Table 4. Such a comparison is reasonable because the low channel height of 50 µm implies an approximate translation invariance along the $y$-axis spanning the length (aperture) 2400 µm of the IDT electrodes, as seen in the 3D geometry of Fig. 1. Also the variation along the $x$ axis given by the width $20~{}\textrm{\textmu{}m}$ of the individual electrodes, and the periodicity $\lambda_{\mathrm{SAW}}=80~{}\textrm{\textmu{}m}$ the IDT, are much smaller than IDT aperture along $y$ axis. We can therefore obtain a reasonable estimate of the electrical and acoustical response of the device, by just considering the 2D domain in the vertical $x$-$z$ plane shown in Fig. 2.

V.1 Electrical response

As a first validation of the model, we study the electrical impedance

[TABLE]

in terms of the driving voltage $V_{0}$ and the complex-valued current $I$ through the device, because this quantity is relatively easy to obtain both in simulation and in experiment. We compare model predictions of the magnitude \big{|}Z^{\mathrm{el}}\big{|} and phase $\psi$ of the impedance with the experimentally measured counterparts.

In the model, we compute $Z^{\mathrm{el}}$ from the time-harmonic dielectric polarization density $\bm{P}$ and the corresponding polarization current $\bm{J}_{\mathrm{pol}}$ in the lithium niobate substrate, which we treat as an ideal dielectric without free charges,

[TABLE]

In Fig. 3(a) and (b), we compare the values of $|Z^{\mathrm{el}}|$ computed by Eq. (17d) for our 2D model with those measured on Pyrex device D2 of Fig. 1(a) and Table 4 for microchannels with air or with DI water. The numerical simulation predicts correctly the value of the resonance observed near 48 MHz in the experiments. As shown in Table 5, the relative difference between computed and measured values of the frequencies $f$, where $|Z^{\mathrm{el}}(f)|$ and $\psi(f)$ have local minima or maxima, is about 2 % or less. We also see that simulation also predicts the monotonically decreasing background signal for $|Z^{\mathrm{el}}(f)|$ before and after the resonance relatively well for both an air- and water-filled microchannel. However, the simulation fails to predict the correct ratio of the resonance peak heights.

For the phase $\psi$ shown in Fig. 3(c) and (d), the simulation predicts the resonance frequencies correctly, but fails to predict the monotonically increasing background signal. By adding external stray impedances to our 2D model to simulate the surrounding 3D system, it is however possible to generate a slant in the phase curves by fitting the values of these stray impedances. We do not show these results as they are descriptive and not predictive in nature.

V.2 Wall material: hard pyrex versus soft PDMS

Our previous device Sehgal and Kirby (2017) features a soft PDMS polymer lid, as is commonly used due to the ease of fabrication and handling. However, the acoustic properties of PDMS are far from ideal: its impedance is nearly equal to that of water (20 % lower) and the attenuation is about two orders of magnitude larger than that of the boundary layer in water. In the following, we therefore simulate the acoustic properties of the device D1 with a PDMS lid and contrast them with those of device D2 with a much stiffer Pyrex lid, using the two models shown in Fig. 2 and Table 4. Compared to water, the acoustic impedance of Pyrex is 8.3 times larger and its attenuation 10 times smaller.

We study by numerical simulation the acoustic fields of device D1 and D2 near the ideal (unloaded) frequency $f_{\mathrm{SAW}}=49.9$ MHz. By locating the maximum of the average acoustic energy in the water-filled channel plotted versus the actuation frequency $f$ (not shown), we determine the (loaded) resonance frequency $f_{\mathrm{res}}$ of the two devices to be $f_{\mathrm{res}}^{\mathrm{D1}}=47.75$ MHz and $f_{\mathrm{res}}^{\mathrm{D2}}=46.50$ MHz, respectively. In Fig. 4 we show line plots along the height ($z$ direction) and across the width ($x$ direction) of numerically simulated acoustic fields for these two devices.

In Fig. 4(a) and (b) is shown the magnitude $|u_{z}|$ of the $z$ component of the acoustic displacement $\bm{u}$, which in water is defined through acoustic velocity Eq. (7b) as $\bm{v}_{1}=-\mathrm{i}\omega\>\bm{u}$, along a vertical cut-line through the entire device. In D1, $|u_{z}|$ has the characteristics of a traveling wave emitted from the SAW substrate (maximum amplitude), traversing the water with little reflection (a small oscillation amplitude), and being absorbed in the PDMS lid (decaying amplitude). In contrast, $|u_{z}|$ in D2 has the characteristics of a standing wave localized in the water channel with reflections from the surrounding solids: huge oscillations in the water domain with minima close to zero and an amplitude exceeding that in the emitting substrate and the receiving lid. We also notice that in the stiff Pyrex the attenuation is weak, and that the wave is reminiscent of a standing wave between the water interface below the lid and the air interface above. The corresponding acoustic energy flux density \bm{S}_{\mathrm{ac}}=\big{\langle}p_{1}\bm{v}_{1}\big{\rangle} in both systems is non-zero and predominantly vertical, but with a much larger amplitude in D1 compared to D2.

In Fig. 4(c) is shown the the magnitude $|\bm{u}|$ of the acoustic displacement $\bm{u}$ along horizontal cut-lines following the top ($z=H_{\mathrm{fl}}$) and the bottom ($z=0$) of the water channel across the region containing the IDT. In both devices the periodicity of the IDT electrodes is clearly seen, but the amplitude in the nearly-standing wave case of D2 is 2-3 times larger than in the traveling wave case of D1. Moreover, it is seen that the acoustic waves dies out faster in D1 than in D2 away from the IDT region. The tiny oscillations in the PDMS lid (green curve) for $x<-12\lambda_{\mathrm{SAW}}$ stems from the minute transverse wavelength $\sim 11~{}\textrm{\textmu{}m}=0.13~{}\lambda_{\mathrm{SAW}}$ in PDMS.

In Fig. 4(d) is shown the the magnitude $|p_{1}|$ of the acoustic pressure $p_{1}$ in the water along the horizontal cut-lines $z/H_{\mathrm{fl}}=\frac{1}{6},\frac{2}{6},\frac{3}{6}$. Here the traveling versus standing wave nature of the two devices mentioned above, is prominent: In D1, $|p_{1}|$ is nearly independent of the height, and its envelope amplitude is steadily decaying from 90 to 55 kPa from the center to the edge of the IDT region. In contrast, $|p_{1}|$ has large amplitude fluctuations as a function of the horizontal position $x$ and for the three vertical $z$ positions. Moreover, $|p_{1}|$ does not decay away from the the IDT. Clear, $p_{1}$ in the water channel of D2 is dominated by reflections between the solid-water interfaces. This observation can be quantified by the the standing wave ratio, $\mathrm{SWR}={\max(|{p_{1}}|)}/{\min(|{p_{1}}|)}$ that describes the ratio of standing to traveling waves in a given field. In an ideal resonator and an ideally transmitting system, $\mathrm{SWR}=\infty$ and 1, respectively. Here, we find $\mathrm{SWR(D2)}=12.7$ and $\mathrm{SWR(D1)}=1.3$. These numbers underlines the good acoustic properties of the water-Pyrex systems compared to the bad one of the PDMS system. The ratio of the SWR numbers is 9.8, almost equal to the impedance ratio 10.5, which emphasizes the nearly perfect vertical energy flux density $\bm{S}_{\mathrm{ac}}$ discussed above, as the impedance extracted from the properties of a plane wave with a vertical incident on a planar surface.

V.3 Acoustophoresis

Whereas we have not made experimental validation of the above simulation results for the acoustic fields $p_{1}$ and $\bm{u}$, we compare in the following the experimentally observed acoustophoretic motion at the SAW resonance frequency $f_{\mathrm{SAW}}$ of microparticle suspensions in the water-filled microchannel, with that obtained by numerical simulation in our 2D model. The central experimental and numerical results are shown in Fig. 5, in the left column for the PDMS-lid device D1 and in the right column for the Pyrex-lid device D2. In the Supplemental Material 111The Supplemental Material at [URL] contains four animations of the numerically simulated acoustophoresis of 0.1- and 1.7-µm-diameter particles in device D1 and D2, respectively, corresponding to Fig. 5(c) and (g). are shown four animations of the acoustophoresis in Fig. 5(c) and (g) of 0.1- and 1.7-µm-diameter particles in device D1 and D2.

In Fig. 5(a) and (e) we observe that the suspended 1.7-$\upmu\textrm{m}$-diameter particles in D1 focus on the edges of the electrodes, whereas in D2 they mainly focus along the center line of each electrode. This difference in acoustophoretic focusing is caused solely by choice of lid material and its thickness. Already in Fig. 4, we saw how the change from the PDMS lid to the Pyrex lid led to a change from a predominantly traveling wave, to a nearly standing wave in the $z$ direction. As a consequence, both the pressure and its gradients in device D1 are smaller than those in D2, and from Eq. (13b) follows that the acoustic radiation force $\bm{F}^{\mathrm{rad}}$ changes significantly.

This change in $\bm{F}^{\mathrm{rad}}$ per particle volume, named $\bm{f}^{\mathrm{rad}}$ in Eq. (13d), is shown as the vector and gray-scale plots for device D1 and D2 in the right half of Fig. 5(b) and (f), respectively. Compared to D2 having $|\bm{f}^{\mathrm{rad}}|=7.4~{}\textrm{pN}/\upmu\textrm{m}^{3}$, the magnitude $|\bm{f}^{\mathrm{rad}}|=0.4~{}\textrm{pN}/\upmu\textrm{m}^{3}$ is 18 times smaller in D1, and $|\bm{f}^{\mathrm{rad}}|$ is more smeared out (even smaller gradients). Both force fields have a three-period structure along the vertical $z$ axis, reflecting that $H_{\mathrm{fl}}\approx\frac{3}{2}c_{\mathrm{fl}}/f_{\mathrm{SAW}}$. In D1, the center of the force-field structure is displayed relative to the center of the electrode, whereas in D2 it is above the electrode center. Moreover, whereas $\bm{f}^{\mathrm{rad}}$ has four less-marked, unstable nodal planes in D1 at $z/H_{\mathrm{fl}}=0,\frac{1}{3},\frac{2}{3},1$, it has three well-defined, stable ones in D2 at $z/H_{\mathrm{fl}}=\frac{1}{6},\frac{3}{6},\frac{3}{6}$.

The corresponding streaming velocity field $\bm{v}_{2}$ in D1 and D2 is shown as the vector and color plots in the left half of Fig. 5(b) and (f). The streaming appears strikingly equal both in magnitude ($66~{}\textrm{\textmu{}m}/\textrm{s}$ for D1 and $76~{}\textrm{\textmu{}m}/\textrm{s}$ for D2), shape and topology, but again with the center of the pattern in D1 shifted slightly away from the electrode center. The reason for this resemblance in $\bm{v}_{2}$ stems from the energy flux density $\bm{S}_{\mathrm{ac}}$, which in both devices points (nearly) vertically up along the $z$ axis above the electrodes, and is weak in between. As the (Eckart) streaming is proportional to $\bm{S}_{\mathrm{ac}}$ Eckart (1948), even in microcavities Skov et al. (2019), the streaming moves upward due to $\bm{S}_{\mathrm{ac}}$ above the electrodes and downward by recirculation between the electrodes. $\bm{S}_{\mathrm{ac}}$ has nearly the same amplitude in D1 and D2 because, although the acoustic field in D2 is much larger than in D1, it is mostly a standing wave with zero energy flux density, and the little part that is a traveling wave in D2 that carries the energy flux density, is nearly of the same magnitude as the traveling wave that constitutes the main part of the weaker acoustic field in D1.

According to Newton’s second law (13), the above-mentioned properties of the acoustic radiation force density $\bm{f}^{\mathrm{rad}}$ and streaming velocity field $\bm{v}_{2}$ governs the observable acoustophoretic motion of suspended particles. In Fig. 5(c) and (g), as well as in the Supplemental Material Note (1), is shown the results of simulating such motion for 0.1- and 1.7-µm-diameter polystyrene beads in both D1 and D2, 0.5 s after starting from an initial homogeneous distribution (blue points and percentage numbers). The motion of the large particles is dominated by the radiation force Muller et al. (2012), so the different focusing of these particles seen in the right half of Fig. 5(c) and (g) is explained in terms of $\bm{f}^{\mathrm{rad}}$: Because $\bm{f}^{\mathrm{rad}}$ has no stable nodal planes in D1, all particles accumulate the floor or the ceiling of the channel, and most of them (98 %) are pushed to the regions above the electrode gaps as indicated by the vector plot in the right half of Fig. 5(c). In contrast, the stable nodal planes of $\bm{f}^{\mathrm{rad}}$ in D2, Fig. 5(g) right half, guides 96 % of the particles into the three stable points above the electrode center, with 41 %, 27 %, and 28 % at $z/H_{\mathrm{fl}}=\frac{1}{6}=0.17$, $\frac{3}{6}=0.50$, and $\frac{5}{6}=0.83$, respectively. If we instead, as in the experiments described below, allow for a sedimentation time of 3 min before turning on the acoustics, the distribution of the focused particles changes to 58 %, 40 %, and 2 % at $z/H_{\mathrm{fl}}=\frac{1}{6}=0.17$, $\frac{3}{6}=0.50$, and $\frac{5}{6}=0.83$, respectively.

The acoustophoretic motion of the small 0.1-$\upmu\textrm{m}$-diameter particles are dominated by the Stokes drag from the streaming field $\bm{v}_{2}$, see the left side of Fig. 5(c) and (g) and the videos in the Supplemental Material Note (1). The simulation shows that the particles do not settle in fixed positions but follow oblong paths in the vertical plane similar in shape to the large streaming rolls spanning the entire height of the channel with an upwards motion over the electrodes and downwards in between electrodes, see Fig. 5(b) and (f). In D1, $\bm{F}^{\mathrm{rad}}$ is so small that it plays essentially no role. In D2, however, $\bm{F}^{\mathrm{rad}}$ is stronger and superposes with $\bm{F}^{\mathrm{drag}}$ to govern the acoustophoretic motion. This superposition of forces is similar to the analysis presented by Antfolk et al. Antfolk et al. (2014), but whereas in their system the nanoparticles spirals towards the point at the center of a single flow roll, the nanoparticles above a single electrode in D2 are focused into the center line of each of the two flow rolls shown in Fig. 5(f). The location of these center lines are defined by the vertical and horizontal nodal lines of $\bm{f}^{\mathrm{rad}}$ represented by the black regions at the electrode gaps $x/\lambda_{\mathrm{SAW}}=\frac{n}{2}$ and at the stable nodal planes $z/H_{\mathrm{fl}}=\frac{1}{6},\frac{5}{6}$, respectively, in the gray-scale plot of Fig. 5(f) and (g).

Most of these theoretical predictions are validated by experiments. After loading the particle suspension in to the device, it takes about 3 min for the fluid to come to rest, during which time the 1.7-$\upmu\textrm{m}$-diameter particles sediment slowly. This partial sedimentation shifts the homogeneous particle distribution downwards, so that the particle distribution is inhomogeneous when the acoustic field is turned on. In the experiments on PDMS-lid device D1, the large 1.7-$\upmu\textrm{m}$-diameter particles are observed to accumulate at the floor and the ceiling in the regions between the electrodes, and the small 0.1-$\upmu\textrm{m}$-diameter particles are observed to circulate in broad streaming rolls. In contrast, in the experiments on the Pyrex-lid device D2, the large particle are seen to accumulate above the center of the electrodes near two planes, 36 % of them at $z=(15\pm 5)~{}\textrm{\textmu{}m}=(0.3\pm 0.1)H_{\mathrm{fl}}$ and 64 % of them at $z=(30\pm 5)~{}\textrm{\textmu{}m}=(0.6\pm 0.1)H_{\mathrm{fl}}$. Here, the uncertainty is estimated from the optical focal depth in the setup. These numbers are in fair agreement with the simulation results mentioned above and shown in Fig. 5(g) (purple numbers). Finally, the observed acoustophoretic focusing time of 0.1 s matches the theoretical predictions.

VI Results of the 3D modeling

In this section we address the more realistic, but also more cumbersome simulations in 3D for the Pyrex-lid device D2. Even given our access to the High Performance Computing clusters at the DTU Computing Center (HPC-DTU) HPC , we cannot simulate the entire chip shown in Fig. 1(a). Whereas we keep the correct dimensions in the height, we scale down the width and length to both be around 1 mm. The 3D model geometry is shown in Fig. 6 with the detailed parameter values listed in Table 3. In this reduced geometry, the IDT contains only 4 electrode pairs and no Bragg reflectors. Although the model is down-sized in two of the three dimensions, it still contains all the main components of a acoustofluidic SAW device: A first step, in which the piezo-electric device, the IDT electrodes, the elastic lid, and the microchannel with the fluid and its viscous boundary layer, are combined in the calculation of the electrically induced acoustic fields. A second step, in which the acoustic radiation force and the acoustic streaming velocity are computed, and used in the governing equation predict the acoustophoretic motion of suspended spherical particles.

VI.1 The acoustic fields and radiation force

The 3D model shown in Fig. 6 contains 4.6 million degrees of freedom. The calculation was distributed across 80 nodes on the HPC-DTU cluster and took 14 hours to compute. The first result is that the computed pressure and displacement fields $p_{1}$ and $\bm{u}$ in 3D are both qualitatively and quantitatively similar to the ones computed in the 2D model. For vertical slice planes parallel to the $x$-$z$ plane and place near the center of the IDT at $y=\frac{1}{2}L_{\mathrm{sl}}$, the agreement is of course better than for those near the edge of the IDT near $y=\frac{1}{2}(L_{\mathrm{sl}}\pm L_{\mathrm{el}})$, but in all cases we find the period-3 structure of $|p_{1}|$ along the $z$ direction seen in Fig. 4(b). Likewise, for the acoustic radiation force density, we recover the period-3 structure in $|\bm{f}^{\mathrm{rad}}|$ seen in Fig. 5(f) and for the particle focusing points in Fig. 5(g). The experimental observation of this vertical focusing is thus validating this point in our 3D model.

VI.2 The acoustic streaming rolls

The streaming-dominated, in-plane acoustophoretic motion of 0.75-µm-diameter particles suspended in the device is used in Fig. 7 to compare our model predictions to observed particle motion. As shown in Fig. 7(a), the experimentally observed particle motion in Pyrex-lid device D2 at the edges of the IDT electrodes is dominated by streaming rolls in the horizontal $x$-$y$ plane. We compare this motion with the streaming velocity field ${\bm{v}_{2}}$ calculated using the 3D model and shown in Fig. 7(b). Although the model only includes a mm-sized sub-region of the experimental device, the same streaming pattern is evident in both the model device and in the experimental device. The agreement in terms of direction, position, and magnitude is good, albeit with small differences. In both the simulation and in the experiment, the centers of the streaming rolls are located at the edges of the electrodes, with clockwise-circulating flows. Similar to the 2D streaming pattern in Fig. 5, the observed horizontal streaming rolls are a combination of a recirculating flow and an energy flux density, here perpendicular to and away from the IDT array. The streaming velocity in D2 near the right edge of the blue rectangular region shown in Fig. 7(a) and (b) is measured in the 24-electrode-pair device to be $\sim 200~{}\textrm{\textmu{}m}/\textrm{s}$ and in the simulated 4-electrode-pair device to be $\sim 20~{}\textrm{\textmu{}m}/\textrm{s}$, or $\sim 120~{}\textrm{\textmu{}m}/\textrm{s}$ if multiplied by the ratio of the number of electrode pairs, 24/4.

VII Discussion

By comparing our model simulations to measurable quantities, we find that the model can predict the overall electrical and acoustophoretic behavior of the two types of SAW-devices D1 (PDMS lid) and D2 (Pyrex lid) fairly well. For the electrical response of the device we see a good agreement between the trends near resonance of the predicted and measured values of the electrical impedance, although the predicted values are obtained in an ideal 2D model neglecting stray impedances. The predicted acoustophoretic focusing of the 1.7-µm-diameter polystyrene particles at the ceiling and floor above the edges of the electrodes in D1, and at 1/6 and 3/6 of the channel height above the center of the electrodes in D2, agrees well with experimental observations.

An interesting feature of the model is the three half-wave resonance excited vertically in the Pyrex-lid device D2. It highlights the importance of careful consideration of the material selection for acoustofluidic devices, to fit with the desired purpose of the device. Because a PDMS lid is an acoustically soft material with an acoustic impedance $Z^{\mathrm{ac}}$ similar to that of water ($Z^{\mathrm{ac}}_{\mathrm{PDMS}}=1.19~{}\textrm{MRayl},Z^{\mathrm{ac}}_{\mathrm{H_{2}O}}=1.49~{}\textrm{MRayl}$), most of the energy in an acoustic wave in water impinging on the water-PDMS interface is transmitted into the PDMS, where it dissipates into heat. Only a small fraction of the energy is reflected back into the fluid. As illustrated in Fig. 4, by replacing the PDMS lid of the device in Ref. Sehgal and Kirby (2017) with an acoustically hard ($Z^{\mathrm{ac}}_{\mathrm{Py}}=12.47~{}\textrm{MRayl}$) Pyrex lid, 78.6% of the wave energy is theoretically reflected back into the fluid domain at the channel lid, compared to the 10.9% in a PDMS lid. The resonance build-up in the microchannel is further enhanced, as the height of the channel sustains three half-waves at the resonant frequency of the IDT, $f_{\mathrm{res}}=\frac{c_{\mathrm{SAW}}}{\lambda_{\mathrm{SAW}}}=\frac{c_{\mathrm{fl}}}{3\lambda_{\mathrm{Wa}}}$. This resonance behavior is very similar to the integer-half-wave resonances common in BAW devices, whereas the beneficial energy localization at the surface of the SAW is still retained. Thus, the energy loss and heat generation occurring in the piezo-electric substrate in BAW devices is mitigated in this device whereas strong microchannel resonances can be achieved, when using a Pyrex lid in an IDT-inside SAW design. Considering this the terms ’BAW’ and ’SAW’ seem inadequate when describing acoustofluidic devices, as the actuation scheme of the piezo-electric transducer alone does not suffice to describe the resonance behavior of a device. A more descriptive feature of a device is the nature of the wave field in the fluid, because we show the main factor determining acoustophoresis in the SAW is the difference between traveling and standing wave fields in the fluid.

In acoustofluidic focusing devices, a strong streaming flow is often detrimental to the desired application, as they tend to counteract the radiation force by pulling small particles away from the nodes. In the Pyrex-lid device, however, the vertical part of the streaming enhances particle focusing, as it pulls particles from areas with weak radiation force into the lower node of the acoustic radiation force, increasing the focusing efficiency.

VIII Conclusion

We have presented a 3D model, and implemented it in the finite-element software COMSOL Multiphysics, for numerical simulation of SAW-devices taking into account the piezo-electric substrate, the IDT metal electrodes, the elastic solid defining the microchannel, the water in the microchannel as the viscous boundary layer of the water. With such simulations, we are able to decrease the gap between the systems that we can model and those used in actual experiments. This work thus brings us closer to the point, where numerical simulation can guide rational design of acoustofluidic devices.

To push acoustofluidic devices closer to medical application, the development of novel device designs beyond the proof-of-concept stage is vital. We have presented a close-to-scale numeric model of an acoustofluidic SAW device by expanding on previous model experiences Ley and Bruus (2017); Skov et al. (2019) and the recently developed effective-boundary-layer theory Bach and Bruus (2018). With this we have captured the inner workings of a non-trivial device. The model includes the linear elasticity of the defining material, the scalar pressure field of the microchannel fluid and the piezoelectricity of the lithium niobate substrate.

Using the numeric model, we illustrate the impact that the material selection in acoustofluidic chips has on acoustophoretic performance. Based on the numerically predicted acoustic fields, we propose design improvements over the previous design Sehgal and Kirby (2017), consisting primarily of substituting the original PMDS lid with a Pyrex lid. According to our model, the new lid leads to higher energy densities and more uniform particle focusing. This causes the chip to build up strong resonances in a standing wave field, similar to those in a BAW device. Furthermore, we have used our model to predict the electrical response of the a 2D model of the system, the acoustophoretic focusing of particles suspended over the IDT area of the device, and the streaming motion within devices. For each of these comparison parameters we have found an agreement between predictions and experiments.

Despite our focus on a specific device design in this manuscript, the model can handle a much wider class of acoustofluidic devices. We have a developed a model that can be reshaped to simulate any BAW or SAW device design of well-characterized piezoelectric transducers, Newtonian fluids, and isotropic and anisotropic linear solids.

In future work it would be prudent to improve the model accuracy by including the temperature field to account for the thermal dependence of material parameters, particularly the fluid bulk and dynamic viscosities. To implement the temperature field one must account for the various sources of thermal generation in terms of mechanical losses and viscous dissipation described in Hahn and Dual (2015), which requires a good knowledge of the damping properties of each component of the device.

IX Acknowledgements

This work is partially supported by NSF CBET-1605574, NSF CBET-1804963, and NIH PSOC- 1U54CA210184-01. The device fabrication is performed in part at the Cornell Nanoscale Facility (CNF), which is supported by the National Science Foundation (Grant ECCS-1542081).

Appendix A Bond and rotation matrices

The elasticity, coupling and permittivity properties of mono-crystalline lithium niobate are listed in Ref. Weis and Gaylord (1985) for a Cartesian material coordinate system $X,Y,Z$ defined as shown in Fig. 8(a). The $Z$-axis is oriented in the growth direction, the $X$-axis is the normal to one of the three mirror planes, and the $Y$-axis follows from the right-hand rule, placing it within the mirror plane the $X$-axis is normal to. The device in this manuscript, however, is manufactured on a wafer of the more commonly used 128∘ YX-cut lithium niobate. These are wafers of lithium niobate cut from a single crystal so that the positive surface normal forms a 128∘ angle with the material $Y$-axis. In our model, we define a coordinate system $x,y,z$ with the $x$-axis coinciding with the material $X$-axis, the $z$-axis normal to the wafer surface and the $y$-axis determined by the right-hand rule. This global coordinate system coincides with the material coordinate system rotated an angle $\theta=128^{\circ}-90^{\circ}=38^{\circ}$ counter-clockwise around the $X$-axis, as shown in Fig. 8(b).

In the following, we define the matrix operations necessary to determine the material parameters in the global system $x,y,z$ from the values known in $X,Y,Z$.

In the usual Cartesian notation exists a matrix $\bm{R}$ transforming a 3$\times$1 vector $\bm{P}^{\mathrm{mt}}$ expressed in material coordinates $X,Y,Z$ to the vector $\bm{P}^{\mathrm{gl}}$ expressed in terms of a global coordinate system $x,y,z$.

[TABLE]

whereas 3$\times$3 matrices are transformed as

[TABLE]

For 6$\times$1 vectors in Voigt notation similar matrices called Bond matrices $\bm{M}_{s}$ transform stress vectors $\bm{\sigma}^{\mathrm{mt}}_{\mathrm{V}}$ expressed in material coordinates into the same stress in terms of the global coordinate system $\bm{\sigma}^{\mathrm{gl}}_{\mathrm{V}}$

[TABLE]

and similarly to Eq. (19) 6$\times$6 matrices are transformed as

[TABLE]

It is important to note that Voigt notation stress and strain vectors do not transform alike and two transformation matrices exist in Voigt notation $\bm{M}_{\sigma}\neq\bm{M}_{\epsilon}$. Hence, the transformation rules deviate slightly from those in 3$\times$3 matrices.

Finally, 3$\times$6 matrices such as the coupling tensor, $\bm{e}$ can be transformed using a rotation matrix and Bond matrix.

[TABLE]

Mathematically, a positive rotation $\theta$ degrees about the material $X$-axis is obtained by the rotation ${\bm{R}_{x}(\theta)}$ and Bond $\bm{M}_{\sigma,x}(\theta)$ matrices

[TABLE]

using $\mathcal{C}$ and ${\mathcal{S}}$ as shorthand for $\cos(\theta)$ and $\sin(\theta)$ respectively.