Three-Dimensional Numerical Modeling of Acoustic Trapping in Glass Capillaries

TL;DR

This paper develops a detailed 3D numerical model of acoustic trapping in glass capillaries, predicting resonance modes and forces, and validating results with experiments to improve trap design.

Contribution

It introduces a comprehensive 3D simulation approach for acoustic traps in glass capillaries, considering complex wave interactions and absorption effects.

Findings

Model accurately predicts resonance frequencies and quality factors.

Simulation results agree with experimental data.

Provides insights for designing more stable acoustic traps.

Abstract

Acoustic traps are used to capture and handle suspended microparticles and cells in microfluidic applications. A particular simple and much-used acoustic trap consists of a commercially available, millimeter-sized, liquid-filled glass capillary actuated by a piezoelectric transducer. Here, we present a three-dimensional numerical model of the acoustic pressure field in the liquid coupled to the displacement field of the glass wall, taking into account mixed standing and traveling waves as well as absorption. The model predicts resonance modes well suited for acoustic trapping, their frequencies and quality factors, the magnitude of the acoustic radiation force on a single test particle as a function of position, and the resulting acoustic retention force of the trap. We show that the model predictions are in agreement with published experimental results, and we discuss how improved and more stable acoustic trapping modes might be obtained using the model as a design tool.

Tables

Table 1. Table 1: Material parameters at 25 C ∘ and 4 MHz used in the numerical modeling of the capillary acoustic trap.

Water Muller and Bruus (2014):
Mass density	${\rho }_{\mathrm{fl}}$	997.05	kg m-3
Compressibility	${\kappa }_{\mathrm{fl}}$	452	TPa-1
Speed of sound	$c_{\mathrm{fl}}$	1496.7	m s-1
Dynamic viscosity	${\eta }_{\mathrm{fl}}$	0.890	mPa s
Damping coefficient Hahn and Dual (2015)	${\Gamma }_{\mathrm{fl}}$	$0.004$
Viscous boundary layer	$\delta$	0.27	µm
Pyrex glass Cor :
Mass density	${\rho }_{\mathrm{sl}}$	2230	kg m-3
Young’s modulus	$E$	62.75	GPa
Poisson’s ratio	$\nu$	0.2
Speed of sound, longitudinal	$c_{\mathrm{lo}}$	5592	m s-1
Speed of sound, transverse	$c_{\mathrm{tr}}$	3424	m s-1
Damping coefficient Hahn and Dual (2015)	${\Gamma }_{\mathrm{sl}}$	0.0004
Polystyrene (ps) for 12 $\mu$m-diameter ps particles in water:
Mass density crc (2014)	${\rho }_{\mathrm{ps}}$	1050	kg m-3
Compressibility Muller et al. (2012)	${\kappa }_{\mathrm{ps}}$	249	TPa-1
Poisson’s ratio Mott et al. (2008)	${\nu }_{\mathrm{ps}}$	0.35
Speed of sound at 20 ${}_{}{}^{\circ }\text{C}$ Bergmann (1954)	$c_{\mathrm{ps}}$	2350	m s-1
Monopole coefficient, Eq. (6b)	$f_0$	0.48
Dipole coefficient, Eq. (6c)	$f_1$	0.052-0.003i

Table 2. Table 2: Capillary geometries modeled in this work: C1 from Hammarström  et al. ( 2012 ) , C2 from Lei  et al. ( 2013 ) , C3 from Mishra  et al. ( 2014 ) , C4 from Gralinski  et al. ( 2014 ) , and C5 a capillary design proposed in this work. Besides the symbols defined in Fig.  2 , R cu subscript 𝑅 cu R_{\mathrm{cu}} is the radius of curvature of the fluid channel corners, f res exp superscript subscript 𝑓 res exp f_{\mathrm{res}}^{\mathrm{exp}} and f res subscript 𝑓 res f_{\mathrm{res}} is the experimental and numerical resonance frequency, respectively, and Q 𝑄 Q is the numerically calculated quality factor.

Device	$L$ [mm]	$W$ [mm]	$H$ [mm]	$H_{\mathrm{gl}}$ [mm]	$R_{\mathrm{cu}}\left[\text{µm}\right]$	$L_{p,x}$ [mm]	$L_{p,y}$ [mm]	$f_{\mathrm{res}}^{\exp }$ [MHz]	$f_{\mathrm{res}}$ [MHz]	$Q$	$d_0$ [nm]
C1	2.0	2.0	0.20	0.14	25	1.16	2.0	3.970	3.906	53	0.10
C2	2.0	6.0	0.30	0.30	25	1.0	1.0	2.585	2.495	78	0.10
C3	8.5	0.1	0.10	0.05	19	15	1.0	7.900	6.406	222	0.10
C4	10.0	0.85	0.85	0.225	425	4.0	1.0	1.055	0.981	109	0.10
C5	2.0	0.5	0.20	0.14	25	1.16	0.50	N/A	4.201	53	0.10

Table 3. Table 3: Boundary conditions imposed on the solid and fluid domains in the model shown in Fig.  2 .

Domain $\leftarrow$ boundary	Boundary condition
Solid domain $\leftarrow$ air	${𝝈}_{\mathrm{sl}}\cdot 𝒏=\mathrm{𝟎}$
Solid domain $\leftarrow$ fluid	${𝝈}_{\mathrm{sl}}\cdot 𝒏=-p𝒏$
Solid domain $\leftarrow$ transducer	Eq. (3b) or Eq. (3a)
Solid domain $\leftarrow$ symmetry	$𝒖\cdot 𝒏=0$, $𝒕\cdot {𝝈}_{\mathrm{sl}}\cdot 𝒏=0$
Fluid domain $\leftarrow$ solid	$𝒗\cdot 𝒏=-i\omega 𝒖\cdot 𝒏$
Fluid domain $\leftarrow$ air	$p=0$
Fluid domain $\leftarrow$ symmetry	$𝒏\cdot \nabla p=0$

Table 4. Table 4: List of the mesh parameters used in the model. The maximum mesh element size is set to d mesh = H / 7 subscript 𝑑 mesh 𝐻 7 d_{\mathrm{mesh}}=H/7 with three exceptions: (1) For the reference field g ref subscript 𝑔 ref g_{\mathrm{ref}} of C1 in Eq. ( 7 ) and Fig.  5 , where d mesh ref = 0.9 ​ d mesh superscript subscript 𝑑 mesh ref 0.9 subscript 𝑑 mesh d_{\mathrm{mesh}}^{\mathrm{ref}}=0.9d_{\mathrm{mesh}} , (2) for Fig.  6 , where d mesh = 40.0 ​ µm subscript 𝑑 mesh 40.0 µm d_{\mathrm{mesh}}=40.0~{}\textrm{\textmu{}m} , and (3) For C4 where d mesh = H / 20 subscript 𝑑 mesh 𝐻 20 d_{\mathrm{mesh}}=H/20 .

Location	Linear size	No. of elements	Mesh size
$y$-$z$ plane
Bulk fluid	$H$	7	$d_{\mathrm{mesh}}$
Bulk solid	$H_{\mathrm{gl}}$	6	$d_{\mathrm{mesh}}$
Corner solid-fluid	$\frac{1}{2}\pi R_{\mathrm{cu}}$	10	0.168 $d_{\mathrm{mesh}}$
fluid symmetry edge	$H$	15	$0.5d_{\mathrm{mesh}}$
$x$ direction (extruded)
Bulk fluid	$L+L_{\mathrm{PML}}$	80	$d_{\mathrm{mesh}}$
Bulk solid	$L+L_{\mathrm{PML}}$	80	$d_{\mathrm{mesh}}$

Table 5. Table 5: Simulation results for the quadratic capillary C3: the levitating half-wave resonance frequency f res subscript 𝑓 res f_{\mathrm{res}} , the quality factor Q 𝑄 Q , and the acoustic energy density E ac fl superscript subscript 𝐸 ac fl E_{\mathrm{ac}}^{\mathrm{fl}} . The side length W = H 𝑊 𝐻 W=H and glass thickness H gl subscript 𝐻 gl H_{\mathrm{gl}} are varied within the fabrication specifications, and we use either the average actuation condition ”A” Eq. ( 3a ) or the rigid condition ”R” Eq. ( 3b ). The experimental resonance frequency is f res exp = 7.90 subscript superscript 𝑓 exp res 7.90 f^{\mathrm{exp}}_{\mathrm{res}}=7.90  MHz.

$W=H$	$H_{\mathrm{gl}}$	Actu-	$f_{\mathrm{res}}$	$Q$	$E_{\mathrm{ac}}^{\mathrm{fl}}$
$\left[\text{µm}\right]$	$\left[\text{µm}\right]$	ation	[MHz]	–	[Pa]
90	40	A	8.230	114	67
90	60	A	7.964	60	28
90	60	R	7.950	117	163
100	50	A	7.350	91	40
100	50	R	7.313	157	120
110	40	A	6.670	124	26
110	40	R	6.640	148	84
110	60	A	6.590	79	28

Equations

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

$$\frac{1}{{\rho }_{\mathrm{sl}}}{\mathbf{\sigma }}_{\mathrm{sl}}$$

$$c_{\mathrm{tr}}^2$$

$${\nabla }^{2}p$$

$$\mathbf{v}$$

$$⟨u_z{⟩}_{\text{pz}}$$

$$u_z$$

$$s(x)=k_{\mathrm{PML}}{\left(\frac{x-L}{L_{\mathrm{PML}}}\right)}^2,\text{ for }L\le x\le L+L_{\mathrm{PML}},$$

$${\partial }_x\to {\partial }_{\tilde{x}}$$

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

$$E_{\mathrm{ac}}^{\mathrm{fl}}$$

$$E_{\mathrm{ac}}^{\mathrm{sl}}$$

$$\bm{F}^{\mathrm{rad}}=-\frac{4}{3}\pi a_{\mathrm{pa}}^{3}\Big{[}\kappa_{\mathrm{fl}}\big{\langle}(f_{0}p)\bm{\nabla}p\big{\rangle}-\frac{3}{2}\rho_{\mathrm{fl}}\big{\langle}(f_{1}\bm{v})\cdot\bm{\nabla}\bm{v}\big{\rangle}\Big{]},$$

$$f_0$$

$$f_1$$

$$C(g)=\sqrt{\frac{{\int }_{\Omega }\mathrm{d}V(g-g_{\mathrm{ref}}{)}^2}{{\int }_{\Omega }\mathrm{d}V(g_{\mathrm{ref}}{)}^2}},\text{ for}0\le x\le L.$$

$${\int }_0^{\frac{1}{2}{L}_{\mathrm{p},x}}\frac{2\mathrm{d}x}{L_{\mathrm{p},x}}{\int }_0^{\frac{1}{2}{L}_{\mathrm{p},y}}\frac{2\mathrm{d}y}{L_{\mathrm{p},y}}(u_z-d_0)F(x)=0,$$

$$\displaystyle\text{with }\;F(x)=-\tanh\Bigg{(}\frac{x-\frac{1}{2}L_{\mathrm{p}}}{\Delta L_{\mathrm{p}}}\bigg{)},$$

$$F_{\max }^{\mathrm{rad}}=\alpha \frac{4\pi }{3}{a}^3{E}_{\mathrm{ac}}^{\mathrm{fl}}\frac{\omega }{c_{\mathrm{fl}}},\alpha =0.0139,$$

$$F_{\max }^{\mathrm{rad}}=6\pi {\eta }_{\mathrm{fl}}a\chi v_{\max }^{\mathrm{flow}},\chi ={\chi }_{\mathrm{paral}}^{z=0}=1.064,$$

$$E_{\mathrm{ac}}^{\mathrm{fl}}=\frac{9\beta \chi }{2\alpha }\frac{{\eta }_{\mathrm{fl}}{c}_{\mathrm{fl}}{Q}_{\max }^{\mathrm{flow}}}{\omega a^{2}HW}.$$

Full text

I Introduction

Microscale acoustofluidic devices are used increasingly in biology, environmental and forensic sciences, and clinical diagnostics Bruus et al. (2011); Laurell and Lenshof (2014). Examples include cell synchronization Thevoz et al. (2010), enrichment of prostate cancer cells in blood Augustsson et al. (2012), manipulation of C. elegans Ding et al. (2012), and single-cell patterning Collins et al. (2015). Acoustics can also be used for non-contact microfluidic trapping and particle enrichment Hammarstrom et al. (2010); Hammarström et al. (2012, 2014). Trapping, an important unit operation in sophisticated cell and bioparticle handling systems, can also be obtained using other technologies such as in hydrodynamic Zhang et al. (2016), electro- and dielectro-phoretic Voldman (2006), and optical trapping systems Grier (2003), see the review by Nilsson et al. (2009). However, one of the key parameters for successful commercialization of these technologies is throughput, and here only the hydrodynamic and acoustic systems are competitive Cetin et al. (2014). Hydrodynamic systems offer simple device designs and can support high flow rates, as they rely solely on passive transport from inertial migration and drag induced by secondary flows Zhang et al. (2016). However, they do not have the capabilities of active transport.

In contrast, acoustofluidics does support active migration, basically due to the acoustic contrast between the suspending medium and the particles Karlsen and Bruus (2015), and it is therefore an inherently label-free technique. Moreover, it exhibits high cell-viability, even over several days Wiklund (2012). These aspects and other appealing traits have increased the interest in acoustofluidic systems in recent years Bruus et al. (2011); Laurell and Lenshof (2014).

A sub-class of acoustofluidic devices is based on simple, cheap and commercially available glass capillaries. Early studies from the 1990s include improving agglutination between erethrocytes Grundy et al. (1993), simple non-contact particle manipulation in a circular wave tube by acoustic streaming and radiation forces Vainshtein et al. (1996), and highly sensitive size-selective particle separation Wiklund et al. (2001). However, the ultrasound capillary systems remains an active research topic, as cell-handling applications have become more numerous and refined over the years. Recent examples include cardiac myoblast viability in ultrasonic fields Ankrett et al. (2013), sonoporation Carugo et al. (2011), squeezing of red blood cells Mishra et al. (2014), and trapping and retaining particles against an external flow Gralinski et al. (2014).

In this paper, we study in particular a mm-sized glass capillary system used as a versatile acoustic trap in many experimental studies Carugo et al. (2011); Mishra et al. (2014); Lei et al. (2013); Hammarström et al. (2012, 2014) and clinical applications, such as isolation of cell-secreted membrane vesicles Evander et al. (2015) and capture and enrichment of bacteria from blood samples for rapid sepsis diagnostics Ohlsson et al. (2016). A sketch of the generic capillary system is shown in Fig. 1, while typical material parameters are listed in Table 1. The capillary is typically actuated locally by a piezoelectric transducer, which is coupled to the capillary either by epoxy glue or by a small drop of glycerol, the latter allowing for removal and reuse of both the capillary and the transducer. The device is driven at a resonance frequency to obtain the largest possible acoustic field in the fluid. In the optimal case, the main component of the resonance field is a vertically standing half-wave, but with an unavoidable admixture of a traveling wave in the axial direction away from the transducer. The resulting acoustic radiation forces trap particles above or near the transducer in all spatial directions and are often strong enough to retain a single microparticle against the drag from a fluid flow speed of the order of 10 µm/s.

In contemporary acoustofludics it is a challenge to model and optimize the design of a given device. Examples of recent advances in modeling include Lei et al. Lei et al. (2013), who modeled the three-dimensional (3D) fluid domain without taking the solid domain into account; Muller and Bruus Muller and Bruus (2014, 2015), who made detailed models in 2D of the thermoviscous and transient effects in the fluid domain; Gralinski et al. Gralinski et al. (2014), who modeled circular capillaries in 3D with fluid and glass domains without taking absorption and outgoing waves into account; Hahn and Dual Hahn and Dual (2015), who calculated the acoustic field in a 3D model for a glass-silicon device (not a capillary system) and characterized the various loss mechanisms; and Garofalo et al., who studied a coupled transducer-silicon-glass-water system in 2D Garofalo et al. (2017).

In this work we present present a 3D numerical model of the capillary acoustic trap of Fig. 1. We model the acoustic pressure field in the liquid coupled to the displacement field of the glass wall, taking into account mixed standing and traveling waves as well as absorption. We model the outlets, which in practise are connected to long tubes, either as being terminated by reflecting no-stress surfaces or by absorbers that absorbs all outgoing acoustic waves from the capillary. We compare prior experimental results from the four devices listed in Table 2 with predictions of our model, in particular the frequency response, the levitating resonance modes, and the acoustic trapping forces. For one of the devices we perform two convergence analyses to show to which degree numerical convergence is obtained. Lastly, we demonstrate how to apply the model as a design tool, by studying the effects of narrowing the width of a given capillary.

II Theory and numerical model

We model single-frequency harmonic actuation at frequency $f$ and angular frequency $\omega=2\pi f$, such that any first-order acoustic field $g$ in the complex-phase representation has the time dependence $g(\bm{r},t)=g(\bm{r})\mathrm{e}^{-\mathrm{i}\omega t}$. We introduce effective absorption in our equations by modifying the time derivative as $\partial_{t}\rightarrow-\mathrm{i}\omega(1+\mathrm{i}\Gamma)$, where $\Gamma\ll 1$ is an effective absorption parameter. The tubings at the ends of the capillary are modeled as ideal acoustic absorbers using perfectly matched layers (PML). Due to mirror symmetry about the $y$-$z$ and $x$-$z$ planes, we only model a quarter of the system as shown in Fig. 2.

II.1 Governing equations

The dynamics of the solid (Pyrex) is modeled by the elastic displacement $\bm{u}$, stress $\bm{\sigma}_{\mathrm{sl}}$, density $\rho_{\mathrm{sl}}$, and the transverse and longitudinal speed of sound, $c_{\mathrm{tr}}$ and $c_{\mathrm{lo}}$,

[TABLE]

Here, I is the unit tensor, $E$ is Young’s modulus, and $\nu$ is Poisson’s ratio. The fluid (water) with its acoustic pressure $p$, velocity $\bm{v}$, density $\rho_{\mathrm{fl}}$, and sound speed $c_{\mathrm{fl}}$ is modeled as pressure acoustics with absorption Hahn and Dual (2015),

[TABLE]

II.2 Boundary conditions

The applied boundary conditions are summarized in Fig. 2 and Table 3. In short, the stress is zero on all outer boundaries facing the air. The stress and the velocity fields are continuous across all internal boundaries. The piezoelectric actuation on the interface $\partial\Omega_{\mathrm{pz}}$ is applied by imposing either an area-averaged displacement $\langle u_{z}\rangle_{\text{pz}}$ representing a transducer attached via a glycerol layer, or a rigid displacement on $u_{z}$ representing a transducer fixed by glue to the solid,

[TABLE]

Symmetry conditions are used on the two symmetry planes, while outlets with perfect absorption, which effectively removes all reflections, are modeled using perfectly matched layers as described in the following.

II.3 Perfectly matched layers

In most setups, the capillary is connected to the microfluidic circuit through tubings at the ends. If a large section of the capillary is inserted into a connecting elastic tube Lei et al. (2013), this might cause significant absorption. Furthermore, the water domain continues uninterrupted from the capillary into the connecting tube, so an extreme case is that none of the acoustic waves going out from the actuation region are reflected at the ends of the capillary. This situation we model using perfectly matched layers (PML). A PML is a domain acting as an artificial perfect absorber of all outgoing waves, here waves traveling along the $x$ axis. In our model the region $x<L$ is normal, while the PML of length $L_{\mathrm{PML}}$ is at $L<x<L+L_{\mathrm{PML}}$, see Fig. 2. Following Collino and Monk Collino and Monk (1998), ideal absorption along the $x$ direction for $x>L$ can be obtained by a complex-valued coordinate stretching of $x$ in the PML domain. This stretching is based on a real-valued function $s(x)$, which is zero in the normal region $x<L$ and is smoothly increasing from zero in the PML as

[TABLE]

where the latter appears in the weak-form implementation of the governing equations and boundary conditions, see Section II.5. The values of $L_{\mathrm{PML}}$ and $k_{\mathrm{PML}}$ are also discussed in Section II.5. The coordinate stretching in the PML introduces a positive imaginary part in the wave numbers in the $x$ direction, and thus dampens the outgoing waves exponentially along the $x$ axis.

The opposite limit of zero absorption is handled by dropping the PML and imposing the no-stress condition at $x=L$.

II.4 Acoustic energy density and radiation force

Good acoustic trapping in a capillary requires a high acoustic energy density $E_{\mathrm{ac}}^{\mathrm{fl}}$ in the fluid Pierce (1991) and $E_{\mathrm{ac}}^{\mathrm{sl}}$ in the solid Landau and Lifshitz (1986), both given by the sum of the space- and time-averaged kinetic and potential energy density, in the respective volumes $V_{\mathrm{act}}^{\mathrm{fl}}$ and $V_{\mathrm{act}}^{\mathrm{sl}}$ above the actuator,

[TABLE]

where the angled bracket denotes the time average over one acoustic oscillation period.

The trapping force acting on an elastic particle of radius $a_{\mathrm{pa}}$, density $\rho_{\mathrm{pa}}$, and compressibility $\kappa_{\mathrm{pa}}$ suspended in the fluid, is the acoustic radiation force $\bm{F}^{\mathrm{rad}}$, which for mixed standing and traveling waves is given by Settnes and Bruus (2012)

[TABLE]

with $\tilde{\delta}=\frac{1}{a_{\mathrm{pa}}}\sqrt{\frac{2\eta_{\mathrm{fl}}}{\omega\rho_{\mathrm{fl}}}}$ being the viscous boundary-layer thickness normalized by the particle radius $a_{\mathrm{pa}}$.

II.5 Numerical implementation in weak form

Following Gregersen et al. (2009), we implement the governing equations and boundary conditions in the weak-form PDE module of the finite-element software COMSOL Multiphysics 5.2 COMSOL Multiphysics 5.2a . We use Lagrangian test functions of second order for $p$, $u_{x}$, $u_{y}$, and $u_{z}$, and solve for given geometry, materials, actuation frequency $\omega$, and actuation displacement $d_{0}$.

We test the PML as follows. The length $L_{\mathrm{PML}}$ of the PML is comparable to the longest wavelength of the system Oskooi et al. (2008), here the longitudinal wave length $\lambda_{\mathrm{lo}}=2\pi c_{\mathrm{lo}}/\omega=1.4$ mm in Pyrex at 4 MHz. By running parametric sweeps in $L_{\mathrm{PML}}\approx\lambda_{\mathrm{lo}}$ and $k_{\mathrm{PML}}>0$, we find that for $L_{\mathrm{PML}}$ as short as $0.25\lambda_{\mathrm{lo}}$, the results are independent of the PML strength in the range $10<k_{\mathrm{PML}}<10^{4}$, showing that for this broad parameter range, the outgoing waves are absorbed before reaching the PML back edge at $x=L+L_{\mathrm{PML}}$, and without traces of back-scattering from the PML front edge at $x=L$. To further show the independence of the results on the PML, we set $L_{\mathrm{PML}}=0.25\lambda_{\mathrm{lo}}$ and $k_{\mathrm{PML}}=10^{3}$, and vary the length $L_{\mathrm{ch}}$ of the capillary between the piezoelectric transducer and the PML front edge. In Fig. 3, $p$ and $u_{z}$ are plotted versus $x$ for four fixed $(y,z)$ coordinates, and it is seen that for $L_{\mathrm{ch}}\gtrsim\lambda_{\mathrm{lo}}$, the resulting fields coincide everywhere outside the PML domain. For smaller values $L_{\mathrm{ch}}\lesssim\lambda_{\mathrm{lo}}$ we do observe deviations in the response, indicating that the outgoing waves are not fully established before they enter the PML and are absorbed. Consequently, in the following, we choose the fixed parameter values $L_{\mathrm{ch}}=1$ mm, $L_{\mathrm{PML}}=0.25\lambda_{\mathrm{lo}}=0.35$ mm, and $k_{\mathrm{PML}}=10^{3}$.

The mesh is generated by first defining a mesh in the $y$-$z$ plane, see Fig. 4(a), and then extruding it equidistantly along the $x$ axis using the ’Swept Mesh’ function in COMSOL Multiphysics, see Fig. 4(b). The size of the mesh elements is controlled by the maximum mesh size $d_{\mathrm{mesh}}$, except elements on the $z$ axis (the intersection of the two symmetry planes) and on the corner edges of the fluid channel (where the curvature is large), both of which are assigned a smaller maximum mesh size. All mesh sizes used in the model are listed in Table 4. Unless stated otherwise, we have set $d_{\mathrm{mesh}}=H/7$.

For a given field variable $g$, we perform a mesh convergence test based on the convergence parameter $C(g)$, which is defined in Ref. Muller et al. (2012) as

[TABLE]

Here, the integration volume $\Omega$ is the domain in which $g$ is defined, but excluding the PML domain. The field $g$ is calculated with the above-mentioned characteristic mesh size $d_{\mathrm{mesh}}$, and $g_{\mathrm{ref}}$ is the reference field calculated with the finer mesh $d_{\mathrm{mesh}}^{\mathrm{ref}}=0.9d_{\mathrm{mesh}}$. We cannot use a smaller value of $d_{\mathrm{mesh}}^{\mathrm{ref}}$, because this value combined with the length $L=2.5~{}\textrm{mm}$ results in a memory consumption of 85 GB RAM out of the 128 GB available on our workstation, see Section III. With these values, our simulation involve $2\times 10^{6}$ degrees of freedom and a computation time around 15 minutes per parameter set. In Fig. 5, semilog plots of $C$ versus $d_{\mathrm{mesh}}^{\mathrm{ref}}/d_{\mathrm{mesh}}$ for all four fields are shown. The plots exhibit an exponential decrease of $C$, which indicates good numerical mesh convergence. For the chosen mesh sizes we obtain $C\approx 0.002$, which is an acceptable level for the present study.

Finally, the area-averaged actuation condition (3a) is imposed on the displacement field $u_{z}$ on the actuation interface using an envelope function $F(x)$,

[TABLE]

where the transition length is set to $\Delta L_{\mathrm{p}}=100~{}\textrm{\textmu{}m}$.

III Results

We characterize numerically the four different glass capillaries C1, C2, C3, and C4 listed in Table 2, for which we find experimental results in the literature: C1 from Hammarström et al. (2012) (and also employed in Refs. Hammarstrom et al. (2010) and Hammarström et al. (2014)), C2 from Lei et al. (2013), C3 from Mishra et al. (2014), and C4 from Gralinski et al. (2014). We also characterize a low-aspect ratio version of C1 denoted C5. We simulate the acoustic fields for these five capillary systems, each with the two different actuation boundary conditions (3a) and (3b), and each both with and without the absorbing PML region. Supplementary to these studies, we simulate geometrical variations to study the effects of the capillaries manufacturing tolerances and the length-dependent reflections from the capillary end. This amounts to more than 20 different capillary configurations, each swept in frequency using 50 values or more, resulting in over 1000 calculations, all of which show good numerical convergence without spurious effects. We remark that a simulation of a given frequency and geometry takes between 2 and 15 minutes on our workstation, a Dell Inc Precision T3610 Intel Xeon CPU E5-1650 v2 at 3.50 GHz with 128 GB RAM and 6 CPU cores.

To gain an overall understanding of these systems, we focus in the following subsections on different acoustofluidic aspects for each capillary.

III.1 Analysis of capillary C1

For capillary C1 Hammarström et al. (2012), we mimic the experimental use of glycerol to attach the transducer by applying the area-average condition (3a) for the actuation. Moreover, to avoid the complications arising from reflections at the end of the capillary, we first assume perfect absorption of outgoing waves and thus use the PML.

First, we study the frequency response. In Fig. 6 we show a semilog plot of the volume-averaged acoustic energy density in the fluid volume above the actuator ($0<x<\frac{1}{2}L_{\mathrm{p},x}$), $E_{\mathrm{ac}}^{\mathrm{fl}}$ for the fluid volume $V^{\mathrm{fl}}_{\mathrm{act}}$, and $E_{\mathrm{ac}}^{\mathrm{sl}}$ for the solid volume $V^{\mathrm{sl}}_{\mathrm{act}}$, as a function of actuation frequency $f$. A levitating resonance, strongly dominated by the energy density in the fluid, is identified at $f_{\mathrm{res}}=3.906~{}\textrm{MHz}$, only 2% lower then the experimental value $f^{\mathrm{exp}}_{\mathrm{res}}=3.970$ MHz and 4% from the simple half-wave value $c_{\mathrm{fl}}/(2H)=3.74$ MHz. The resonance peak is fitted well by a Lorentzian line shape centered around $f_{\mathrm{res}}=3.906~{}\textrm{MHz}$ with a quality factor $Q=53.0$.

We successfully check that this $Q$ factor is consistent with the relation $E_{\mathrm{ac}}^{\mathrm{tot}}=QP_{\mathrm{in}}/(\omega V^{\mathrm{fl}}_{\mathrm{act}})$, where P_{\mathrm{in}}=\int_{A_{\mathrm{act}}}\big{\langle}(-\mathrm{i}\omega\bm{u})\cdot\bm{\sigma}_{\mathrm{sl}}\cdot\bm{n}\big{\rangle}\>\mathrm{d}a is the time-averaged power delivered by the actuator, and $E_{\mathrm{ac}}^{\mathrm{tot}}=E_{\mathrm{ac}}^{\mathrm{fl}}+E_{\mathrm{ac}}^{\mathrm{sl}}$ is the total stored acoustic energy density. Here, we use that in steady state the input power equals the dissipated power. Besides the main levitating resonance peak, minor resonances without specific structure are also present, reflecting the many complex modes in the coupled fluid-solid system. None of these have particular good trapping properties, and they would probably not be observed in acoustophoretic experiments.

Next, we study the acoustic field. Upon inspection of the pressure and displacement fields at the resonance $f_{\mathrm{res}}$, we find them to be mainly a standing wave above the actuation plane combined with an outgoing traveling wave away from the actuator for $x>\frac{1}{2}L_{\mathrm{p},x}$, see Figs. 7 and 8 (VIDEO 1 and 2). The amplitude of the acoustic fields $p$ and $\bm{u}$ are largest above the actuator. The fluid pressure $p$ is a nearly-perfect vertical half-wave with a horizontal nodal plane near the channel center, which enables microparticle levitation given its maximum of 0.19 MPa for an actuation amplitude $d_{0}$ = 0.1 nm. Because of this feature, we call this type of resonance the levitating half-wave resonance in the following. At the edge of the transducer, the pressure amplitude drops an order of magnitude, which according to Eq. (6a) gives rise to the lateral forces that constitute the acoustic trap.

For the displacement field in the solid, we find its maximum value to be 1.2 nm, which is an order of magnitude larger than the average actuation amplitude $d_{0}=0.1$ nm. Besides the outgoing axial displacement waves leaving the channel, we also note the presence of circumferential displacement waves with short wave lengths around $\frac{1}{4}W$. The glycerol-like, area-averaged actuation (3a) fixes only the mean vertical displacement $u_{z}$ while allowing for fluctuations. At resonance $f_{\mathrm{res}}$ we find $u_{z}=(1.0\pm 2.8)d_{0}$ on the actuation plane. We note that the standard deviation $2.8d_{0}=0.28$ nm is very small compared to a typical coupling-layer thickness of around 0.1 mm.

To characterize the trapping capabilities of the device, we study the radiation force $\bm{F}^{\mathrm{rad}}$ Eq. (6a) acting on a 12-µm-diameter polystyrene (ps) test particle suspended in the fluid above the actuator. The $z$ component $F^{\mathrm{rad}}_{z}$ of $\bm{F}^{\mathrm{rad}}$, derived from the nearly perfect standing half-wave in the vertical direction, levitates the test particle and holds it close to the pressure nodal plane. In the actuation region, this plane is nearly horizontal with a vertical position between $-0.1H$ and 0. Away from the actuation region, it gradually changes its vertical position and fades away, see Fig. 9 (VIDEO 3). The maximum levitating force on the test particle is $F^{\mathrm{rad}}_{z}=22$ pN, which is nearly 50 times larger than the buoyancy-corrected gravitational force $F^{\mathrm{grv}}=0.44$ pN.

The lateral trapping force $\bm{F}^{\mathrm{rad}}_{\parallel}$ is given by the $x$ and $y$ components $F^{\mathrm{rad}}_{x}$ and $F^{\mathrm{rad}}_{y}$ of Eq. (6a), and it acts as a retention force when a fluid flow is imposed. In Fig. 10(a) for the C1 capillary at resonance $f_{\mathrm{res}}$, $\bm{F}^{\mathrm{rad}}_{\parallel}$ is plotted in the plane $z=-0.05H$ (vector plot) together with its magnitude $F^{\mathrm{rad}}_{\parallel}$ (color plot). We see two clear trapping points within the actuation region, which are surrounded by strong gradients in the $x$ and $y$ directions, thus enabling lateral particle trapping. The maximum value $F^{\mathrm{rad}}_{x}$ within the actuation region is of particular interest, because it serves as a measure for the overall trapping strength of the device. We find $F^{\mathrm{rad}}_{x}=0.44$ pN at the position $(x,y,z)=(0.3L_{\mathrm{p},x},0,-0.05H)$. Outside the actuation region, a number of weaker, secondary trapping points are observed as shown in Fig. 10(b). Their presence is emphasized by re-scaling the color range, and we note their positions appear to be random. This finding is in qualitative agreement with experimental observations under stop-flow conditions Ley et al. (2016).

We also verify the capability of the device to trap particles by simulating the trajectories of test particles initially placed on a regular grid throughout the trapping region and ending up in the trapping points (not shown).

Finally, we comment on calibration of the acoustic energy density $E_{\mathrm{ac}}^{\mathrm{fl}}$ in the fluid, which controls the magnitude of the trapping force in the capillary device. Because in practise it is not possible to calculate $E_{\mathrm{ac}}^{\mathrm{fl}}$ in a capillary device for a known ac voltage applied to the piezoelectric transducer, it is of importance to measure it in any given experiment. The first method is the drop voltage method Lei et al. (2013), where at first $E_{\mathrm{ac}}^{\mathrm{fl}}$ is set high enough to levitate a test particle. Then the ac voltage and thus $E_{\mathrm{ac}}^{\mathrm{fl}}$ is lowered until the particle drops out of the trap, at which point $F^{\mathrm{rad}}_{z}(E_{\mathrm{ac}}^{\mathrm{fl}})=\frac{4}{3}\pi a^{3}(\rho_{\mathrm{pa}}-\rho_{\mathrm{fl}})g$, and $E_{\mathrm{ac}}^{\mathrm{fl}}$ can be found.

The second method is to increase the flow rate, and thus the flow speed $v_{\mathrm{flow}}$, through the device until it reaches critical magnitude, where the maximum retention force $F^{\mathrm{rad}}_{\mathrm{max}}$ cannot balance the Stokes drag force on the particle any longer. From a scaling argument, and introducing a coefficient $\alpha$ dependent on the device composition and geometry, we can write

[TABLE]

where the value for $\alpha$ is found by simulation on the C1 capillary. The flow speed at which the Stokes drag exactly balances $F^{\mathrm{rad}}_{\mathrm{max}}$ is denoted $v^{\mathrm{flow}}_{\mathrm{max}}$, so we can write

[TABLE]

where $\chi$ is the wall-induced drag enhancement in the center plane of a parallel-plate channel Happel and Brenner (1983), and the value is for the C1 capillary and the given test particle. Approximating the fluid channel to be rectangular, the maximum flow velocity can be related to the maximum flow rate $Q^{\mathrm{flow}}_{\mathrm{max}}$ as $v^{\mathrm{flow}}_{\mathrm{max}}=\frac{\beta}{HW}Q^{\mathrm{flow}}_{\mathrm{max}}$, where $\beta$ is an aspect-ratio dependent constant Bruus (2008), which for C1 is $\beta=1.74$. Combining Eqs. (9) and (10) leads to

[TABLE]

With the assumed actuation amplitude $d_{0}=0.1~{}$nm, we find the maximum acoustic energy density to be 1.60 J/m3, which yields the maximum flow speed and rate of $v^{\mathrm{flow}}_{\mathrm{max}}=3.1~{}\textrm{\textmu{}m}/\textrm{s}$ and $Q^{\mathrm{flow}}_{\mathrm{max}}=0.043~{}\textrm{\textmu{}L}/$min. We emphasize that this calibration method only works if the in-plane streaming velocity is small compared to the flow velocity.

III.2 Analysis of capillary C2

In the experiment, capillary C2 is held in place by two rubber-sleaves Lei et al. (2013), so in this case the use of the absorbing PML region is probably justified. However, capillary C2 differs from capillary C1 by being more flat, having a cross-section aspect ratio of 1:20 in contrast to that of 1:10 for the latter Hammarström et al. (2012). A more significant difference is that the transducer of C2 was attached by glue and not by glycerol as for C1. We therefore study the different responses arising from using these two actuation conditions.

Using the glycerol-like area-averaged actuation condition (3a), a strong levitating half-wave resonance resembling the one in C1 and having a quality factor $Q=78$, is found at $f_{\mathrm{res}}^{\mathrm{avg}}=2.495~{}\textrm{MHz}$ equal to the simple half-wave value $c_{\mathrm{fl}}/(2H)$. The maximum levitating force is $F^{\mathrm{rad}}_{z}=18$ pN.

Changing to the glue-like rigid actuation condition (3b), the resonance frequency drops slightly to $f_{\mathrm{res}}^{\mathrm{rgd}}=2.345~{}\textrm{MHz}$, fairly close to the experimentally observed resonance at 2.585 MHz, while the quality factor decreases to $Q=61$, the acoustic energy density increases by 9%, and the maximum levitating force is nearly doubled to become $F^{\mathrm{rad}}_{z}=33$ pN. Qualitatively, the frequency spectrum changes from having an additional resonance peak close to the levitating resonance of $f_{\mathrm{res}}=2.495~{}\textrm{MHz}$ using the glycerol-like condition, to having no extra resonances using the glue-like condition. The structure of the resonance field in the fluid at $f_{\mathrm{res}}$ largely remains the same for the two different boundary conditions.

III.3 Analysis of capillary C3

For the quadratic capillary C3 Mishra et al. (2014), we use the nominal geometry parameters $W=H=100~{}\textrm{\textmu{}m}$ and $H_{\mathrm{gl}}=50~{}\textrm{\textmu{}m}$ with PML absorption, and we predict numerically the frequency of the levitating half-wave resonance to be $f_{\mathrm{res}}=7.35$ MHz, only 2% from the simple half-wave value $c_{\mathrm{fl}}/(2H)=7.49$ MHz, and 7% lower than the experimental value, $f_{\mathrm{res}}^{\mathrm{exp}}=7.90$ MHz. This discrepancy leads us to study changes in predicted acoustic response as a function of geometrical variations within the fabrication uncertainties listed by the producer (VitroCom #8510, Mountain Lakes, NJ): $\pm 10\%$ for $W$ and $H$, and $\pm 20\%$ for $H_{\mathrm{gl}}$. Some of the resulting numerical predictions for the resonance frequency $f_{\mathrm{res}}$ and $Q$ values, using both the averaged and rigid actuation condition, are listed in Table 5.

For these thin-walled capillaries, the levitating half-wave resonance depends strongly on the channel height $H$ and width $W$, and much less on the glass thickness $H_{\mathrm{gl}}$. The 10% variation $H=W=(100\pm 10)~{}\textrm{\textmu{}m}$ leads to a 10% variation in $f_{\mathrm{res}}\approx(7.4\pm 0.7)$ MHz for fixed $H_{\mathrm{gl}}$, while the 20% variation in $H_{\mathrm{gl}}=(50\pm 10)~{}\textrm{\textmu{}m}$ leads to variations less than 2% in $f_{\mathrm{res}}$ for fixed $H$ and $W$. Moreover, in agreement with the results for C2, $E_{\mathrm{ac}}^{\mathrm{fl}}$ and the $Q$ factor for the glue-like rigid actuation condition R are increased by a factor 3 - 6 and 1.2 - 3, respectively, compared to those of the glycerol-like averaged condition A. The maximum levitating force with $H=W=100~{}\textrm{\textmu{}m}$ for condition R is $F^{\mathrm{rad}}_{z}=3690$ pN, while it drops to 1380 pN for condition A. The two actuation types lead to the same levitating resonance frequency within 0.5%.

III.4 Analysis of capillary C4

In the original work on the circular capillary C4 Gralinski et al. (2014), effects due to reflections at the end of the capillary were discussed. We therefore study numerically the frequency response with PML absorption (perfect absorption of all out-going waves) and without it (reflection of all out-going waves at the end-wall without absorption).

For the case with PML absorption, we find the levitating half-wave resonance to have the following characteristics: $f_{\mathrm{res}}=0.9753$ MHz, $Q=109$, $E_{\mathrm{ac}}^{\mathrm{fl}}=0.96$ Pa, and maximum levitating force $F^{\mathrm{rad}}_{z}=1.6$ pN using the glycerol-like condition (3a), and $f_{\mathrm{res}}=0.9818$ MHz, $Q=79$, $E_{\mathrm{ac}}^{\mathrm{fl}}=0.57$ Pa, and maximum levitating force $F^{\mathrm{rad}}_{z}=2.6$ pN using the glue-like condition (3b). These values for $f_{\mathrm{res}}$ are close to the experimental value of 0.981 MHz and the theoretical value $\gamma_{11}c_{\mathrm{fl}}/(\pi H)=1.03$ MHz based on the simple hard-wall cylinder geometry, where $\gamma_{11}=1.841$ is the first zero of the derivative $J^{\prime}_{1}(x)$ of the first Bessel function. As for the other capillaries, the choice of actuation condition leads to no qualitative, and only minor quantitative differences in the frequency spectrum.

When the PML is removed by imposing the zero-stress condition at $x=L$, the out-going waves from the actuator region are reflected by the end wall, and complex interference effects arise in the capillary. For fixed $L_{\mathrm{p},x}$, the channel length $L$ is varied between 1 mm and 10 mm, and relative to the PML case, we find that the levitating resonance frequency only changes about 1%, the $Q$ value generally increases and becomes $Q\approx 200-300$, and $E_{\mathrm{ac}}^{\mathrm{fl}}$ increases by a factor 4 to 6.

For the longest capillary $L=10~{}\textrm{mm}$, we plot in Fig. 11 $E_{\mathrm{ac}}^{\mathrm{fl}}$ as a function of frequency $f$ near the levitating half-wave resonance $f_{\mathrm{res}}$ with and without PML absorption. Firstly, we note that $f_{\mathrm{res}}$ in the two cases differs less than 0.4%. Secondly, without PML absorption, the energy is not lost by the out-going waves, and consequently, the line shape becomes higher and more narrow. Indeed, the resonance amplitude $E_{\mathrm{ac}}^{\mathrm{fl}}$ increases from 0.58 Pa (PML) to 2.69 Pa (no PML), while the $Q$ factor increases from 109 (PML) to 252 (no PML). Thirdly, without PML, the reflections from the end wall give rise to a number of additional, nearly equidistant, smaller resonance peaks, reminiscent of a Fabry-Pérot-like interference condition for the acoustic wave in the $x$ direction. Finally, without PML, the reflections lead to a more irregular spatial pattern of the levitating resonance mode in the actuation region (not shown), as compared to that of the mode with PML absorption and no reflections.

Due to the complex nature of the coupling between pressure waves in the fluid and the compressional and shear waves in the solid, we have only been able to interpret the levitating resonance mode in simple terms of the three sound speeds of the system, listed in Table 1, combined with the geometrical length scales of capillary C4, listed in Table 2.

III.5 Analysis of capillary C5

We end our numerical modeling of the capillaries by a study of capillary C5, which is a geometry defined by us and not studied in the literature. The focus here is to obtain a more regular spatial pattern of the mode with spatial variations on the longest possible length scale.

The study in Section III.4 of the circular capillary C4 reveals that spatial irregularities in the levitating resonance mode are reduced when the reflections from the end wall are absorbed by a PML. This tendency we also find in the study in Section III.1 of the wide flat capillary C1, where, even in the presence of a PML, the levitating resonance mode has significant spatial variations, see Figs. 8 and 10(b). These spatial variations are caused by the relative high mode numbers in the solid, and therefore we design capillary C5 as a version of capillary C1 with its width reduced by a factor of 4 to $W=0.5~{}\textrm{mm}$, thus reducing the cross-section aspect ratio $\frac{W}{H}$ from 10 to 2.5.

The separation of resonance frequencies increases in the narrow capillary C5. In the range between 3.5 and 4.4 MHz it contains only two resonance peaks (at 3.572 MHz and 4.201 MHz), whereas capillary C1 in the same interval contains three pronounced peaks and four smaller ones. The resonance at $f_{\mathrm{res}}=4.201$ MHz, with $Q=53$, is a levitating mode well suited for acoustic trapping, and in Figs. 12 and 13 (VIDEO 4 and 5) we plot the corresponding pressure field $p$ and displacement field $\bm{u}$. The spatial structure of the resonance is more regular, especially in the solid, where the displacement in the $y$-$z$ plane now resembles a simple fundamental mode with a pivoting motion about a single stationary point away from the symmetry plane (the black point on the $y$ axis in the the inset of Fig. 13). Compared to that of capillary C1, the levitating mode of capillary C5 stretches far outside the actuation region, it is more regular in space, but its amplitude is reduced approximately by a factor of 1.8, and the maximum levitating force drops a factor of 3 to $F^{\mathrm{rad}}_{z}=7$ pN. The pressure nodal surface of the levitating mode is weakly wobbling around the horizontal plane located at $z=-0.03H$, see Fig. 14. This wobbling is less pronounced than the one we find in capillary C1, which moreover is displaced farther downward, namely to $z=-0.05H$.

Compared to capillary C1, the reduced pressure amplitude in capillary C5 leads to a corresponding reduction in the lateral radiation force by a factor of $1.8^{2}=3.2$, plotted for $z=-0.03H$ in Fig. 15, and it has a more regular spatial structure compared to that shown in Fig. 10. The main lateral trapping points are now confined to the intersection of the nodal plane and the $x$-$z$ plane inside the actuation region, although the levitating mode extends far beyond the actuation region, see Fig. 14. There do exist secondary trapping points outside the actuation region, and compared to those of capillary C1 they are relatively strong, but they are fewer in number and appear in a more regular pattern. This finding is in qualitative agreement with experimental observations on the circular capillary C4 Gralinski et al. (2014).

IV Discussion

The main outcome of our numerical modeling of capillary trapping devices is that for all five capillaries C1 - C5 in our study, we find a pronounced levitating mode with good trapping characteristics: a nearly horizontal pressure nodal plane near the center of the channel, see Figs. 9 and 14, with sufficiently strong lateral acoustic forces, see Figs. 10 and 15.

The values of the levitating resonance frequencies $f_{\mathrm{res}}$ are in good agreement with the experimental values obtained for capillary C1 - C4 Hammarström et al. (2012); Lei et al. (2013); Mishra et al. (2014); Gralinski et al. (2014) as listed in Table 2. They vary less than 1% as we change the boundary condition at the outlet from an ideally absorbing PML region to the no-stress condition of a capillary facing air, and as we change the actuation from glycerol-like area-averaged displacement to glue-like rigid displacement. However, we do find that varying these boundary conditions results in significant changes of the acoustic energy density $E_{\mathrm{ac}}^{\mathrm{fl}}$, the acoustic force $\bm{F}^{\mathrm{rad}}$ on a test particle, as well as the resonance quality factor $Q$. In general we find two opposing trends: On the one hand the absorption of the PML region reduces both the acoustic energy and the trapping forces by a factor of 2 - 6, which is a drawback for achieving good trapping. On the other hand, the presence of the PML region has three beneficial aspects: it removes the sensitivity to reflections of acoustic waves at the outlet wall, it results in more regular spatial variations, and it makes the trapping strength independent of the position of the actuator relative to the capillary ends.

The sub-millimeter dimensions of the capillary cross section implies that the levitating resonance frequency $f_{\mathrm{res}}$ is sensitive to the exact geometry. We find that when changing the geometry within the tolerances listed by the manufacturer of the capillaries, the value of the levitating resonance can vary as much as 10% or up to 0.7 MHz for a 7-MHz resonance. Because the levitating resonance is close to an ideal standing half-wave resonance in the channel of height $H$, we have $f_{\mathrm{res}}\propto H^{-1}$, and the most critical length parameter therefore becomes $H$.

Turning to the acoustic radiation force acting on a 12-µm-diameter polystyrene test particle in the device, our model correctly predicts that this force is present in the actuation region with a magnitude large enough to support one or more stable trapping points there. For the assumed actuation amplitude $d_{0}=0.1$ nm, the maximum levitating force $F^{\mathrm{rad}}$ is in the range from 2 pN (the circular capillary C4) to 3690 pN (the square capillary C3), which is between 5 and 8000 times the buoyancy-corrected gravitational force $F^{\mathrm{grv}}$, of course depending on the transducer length $L_{\mathrm{p},x}$. Our model also correctly predicts the existence of weaker, secondary trapping points in the region outside the actuation region, see Figs. 10(b) and 15(b), in agreement with published experiments Gralinski et al. (2014); Ley et al. (2016). However, even without the PML absorption, the magnitude of the lateral acoustic force $\bm{F}^{\mathrm{rad}}_{\parallel}$ is small for most geometries, $F^{\mathrm{rad}}_{\parallel}/F^{\mathrm{grv}}\approx 1-10$, the exception being capillary C3 with $F^{\mathrm{rad}}_{\parallel}/F^{\mathrm{grv}}>100$. This indicates that the performance of actual acoustic capillary traps may depend on additional forces such as streaming-induced drag forces and particle-particle interactions.

While a wide capillary is good for enhanced throughput, it suffers from the existence of high-mode, short-wavelength elastic waves propagating around the perimeter of the capillary. We find that reducing the aspect ratio removes these perimeter waves, while maintaining the good characteristics of the levitating mode. Other irregular spatial patterns in the levitating mode and the acoustic trapping force are induced by reflections of acoustic waves at the end of the capillary. We find that by implementing an ideally absorbing PML region, these reflections are removed, and a more regular spatial behavior of the trap is obtained. While the reduction in aspect ratio reduces the through-put and the introduction of absorption reduces the trapping force by a factor 2 - 6, the resulting increased regularity in spatial behavior, the increased separation of modes in frequency implying reduced mixing of unwanted modes, and the decreased sensitivity to the exact location of the transducer relative to the capillary ends, may be beneficial design considerations worth taking into account. The re-designed devices may be more robust to small perturbations in the device geometry and mechanical actuation.

Our numerical modeling of the capillary systems for acoustic trapping is in fair agreement with existing experimental data in the literature. However, to fully establish its capability as a design tool for optimizing acoustic trapping, more experimental data is needed to characterize the acoustic properties more fully. Here, stop flow experiments as in Ref. Ley et al. (2016) would be helpful, as would more systematic reporting of the acoustic energy density present in the devices using, say, the drop-voltage method Lei et al. (2013) or the flow-versus-retention-force briefly mentioned in Section III.1. Also systematic reporting of the magnitude of the acoustic streaming would be helpful, because there might be situations in the capillary systems, where the acoustic streaming is so large that it is comparable to the flow velocities in the channel. If this were the case, then our analysis of the lateral acoustic forces must be supplemented by a discussion of the drag forces from acoustic streaming.

V Conclusion

The increased use of acoustic-trapping devices in microfluidic handling of particles allows for the development of new functionalities in contemporary biotechnology. In this paper we develop a three-dimensional numerical model for water-filled glass capillaries, which may be a useful design tool in future developments of acoustic microparticle traps. We validate the model and demonstrate, how it captures many of the experimental observations reported in the literature for four different glass capillary devices, in particular regarding the frequency, the $Q$ factor, and the spatial structure and trapping capability of the levitating resonance modes.

We demonstrate the potential of the model as a design tool through an analysis of the sensitivity of the system to changes in the geometry, the specific actuation condition, and the acoustic absorption. We find that while the existence of a main levitating resonance mode and the value of the resonance frequency are relatively insensitive to the these changes, the opposite is true for the magnitude of the trapping force and the lateral spatial structure: The trapping force is strongly increased for rectangular cross-sections with a low aspect ratio, the best being the square shape, and the lateral spatial structure of the levitating mode becomes more uniform. For the case without acoustic absorption, reflections at the outlet of the capillary make these responses very sensitive to the length of the system. The introduction of deliberate absorption of acoustic waves at the outlet would lower this sensitivity and probably result in more robust devices. However, this benefit in device operation must be weighed against a decreased magnitude of the acoustic forces.

Our numerical model takes many, but certainly not all, physical relevant aspects into account. The obvious next step, which we are currently working on, is to include the drag force from acoustic streaming as well as both acoustic and hydrodynamic particle-particle interaction effects. These more elaborate phenomena all involve lateral forces comparable to, or possibly even larger than, the lateral single-particle acoustic radiation force $\bm{F}^{\mathrm{rad}}$ calculated in our model. The simulation of these phenomena relies on numerically well-characterized acoustic fields, and this we believe is provided by the model we present here. Possible model strategies to pursue in such calculations are outlined by Muller et al. (2012); Muller and Bruus (2014) and Lei et al. (2017) for the case of acoustic streaming, by Silva and Bruus (2014) for acoustic particle-particle interactions, and by Ley and Bruus (2016) for hydrodynamic particle-particle interactions.

Our numerical model reveals that the acoustic trapping devices based on simple glass capillaries are perhaps not as simple as they look at a first glance.

Acknowledgements

We thank our collaborators Carl Johannsson, Mikael Evander and Thomas Laurell from Lund University, Sweden, for valuable discussions and experimental data and inputs. This work was supported by Lund University and the Knut and Alice Wallenberg Foundation (Grant No. KAW 2012.0023).