Modeling of Microdevices for SAW-Based Acoustophoresis — A Study of Boundary Conditions

Abstract

We present a finite-element method modeling of acoustophoretic devices consisting of a single, long, straight, water-filled microchannel surrounded by an elastic wall of either borosilicate glass (pyrex) or the elastomer polydimethylsiloxane (PDMS) and placed on top of a piezoelectric transducer that actuates the device by surface acoustic waves (SAW). We compare the resulting acoustic fields in these full solid-fluid models with those obtained in reduced fluid models comprising of only a water domain with simplified, approximate boundary conditions representing the surrounding solids. The reduced models are found to only approximate the acoustically hard pyrex systems to a limited degree for large wall thicknesses and but not very well for acoustically soft PDMS systems shorter than the PDMS damping length of 3 mm.

1. Introduction

Separation of particles and cells is important in a wide array of biotechnological applications [1,2,3,4,5,6,7]. This has traditionally been carried out by bulk processes including centrifugation, chromatography, and filtration. However, during the last three decades, microfluidic devices have proven to be a valuable alternative [1,7,8], as they allow for lower sample sizes and decentralized preparations of biological samples, increasing the potential for point-of-care testing. Microfluidic methods for separating particles suspended in a medium include passive methods where particle separation is solely determined by the flow and the size or density of particles [2,9,10,11,12], and active methods where particles migrate due to the application of various external fields each targeting specific properties for particle sorting [1,3,4,6,13,14,15,16]. Acoustophoresis is an active method, where emphasis is on gentle, label-free, precise handling of cells based on their density and compressibility relative to the suspension medium as well as their size [17]. Within biotechnology, acoustophoresis has been used to confine, separate, sort or probe particles such as microvesicles [6,18], cells [1,16,19,20,21,22], bacteria [23,24], and biomolecules [25]. Biomedical applications include early detection of circulating tumor cells in blood [26,27] and diagnosis of bloodstream infections [28].

The acoustic fields used in acoustophoresis are mainly one of the following two kinds: (1) bulk acoustic waves (BAW), which are set up in the entire device and used in systems with acoustically hard walls. BAW depend critically on the high acoustic impedance ratio between the walls and the water. In addition, (2) surface acoustic waves (SAW), which are defined by interdigital electrodes on the piezoelectric transducer and propagate along the transducer surface. SAW are nearly independent of the acoustical impedance ratio of the device walls and the microchannel, and this feature makes the SAW technique versatile. SAW can be used both with hard- and soft-walled acoustophoretic devices, often in the generic setup sketched in Figure 1, where the fluid-filled microchannel is encased by a solid material and is placed directly on top of the piezoelectric substrate to ensure optimal coupling to the SAW induced in the substrate.

Because SAW-based acoustophoretic microdevices are very promising as powerful and versatile tools for manipulation of microparticles and cells, numerical modeling of them is important, both for improved understanding of the acoustofluidic conditions within the devices and to guide proper device design. In the literature, such modeling has been performed in numerous ways. For many common elastic materials, the dynamics of the walls are straightforward to compute fully through the usual Cauchy model of their displacement fields $\mathit{u}$ and stress tensors $\mathsf{\sigma }$. The coupling to the acoustic pressure p and velocity $\mathit{v}$ in the microchannel, described by the Navier–Stokes equation, is handled by the continuity conditions $\mathit{n}·{\mathsf{\sigma }}_{\mathrm{s}}=\mathit{n}·{\mathsf{\sigma }}_{\mathrm{f}}$ and ${\partial }_t\mathit{u}=\mathit{v}$ of the stress and velocity fields. This full model is discussed in detail in Section 4. For acoustically hard walls, such as borosilicate glass (pyrex) with a high impedance ratio ($\tilde{Z}=8.4$) relative to water, the full model is often replaced by a reduced model (exact for $\tilde{Z}=\infty$) with less demanding numerics, where only the fluid domain in the microchannel is treated, and where the elastic walls are replaced by the so-called hard-wall boundary condition demanding zero acoustic velocity at the boundary of the fluid domain [29,30,31]. For rubber-like polymers such as the often used PDMS, the full device modeling is more challenging. For large strains (above 40 %), a representation of the underlying macromolecular network of polymer chains is necessary [32], while for the moderate strains appearing in typical acoustophoretic devices, standard linear elasticity suffices [33,34]. Some authors argue that the low ratio of the transverse to longitudinal speed of sound justifies a fluid-like model of PDMS based on a scalar Helmholtz equation [30,35]. Furthermore, since the acoustic impedance ratio $\tilde{Z}=0.7$ between PDMS and water is nearly unity, the full model has in the literature been replaced by a reduced model, consisting of only the fluid domain with the so-called lossy-wall boundary condition condition representing in an approximate manner the acoustically soft PDMS walls [16,36,37].

The main aim of this paper is to investigate to which extent the numerically less demanding hard- and lossy-wall reduced models compare with the full models for SAW-based acoustofluidic devices. In the full models, we study the two generic cases of acoustically hard pyrex walls and acoustically soft PDMS walls, both treated as linear elastic materials. In the reduced models, the pyrex and PDMS walls are represented by hard-wall and lossy-wall boundary conditions, respectively. In all the models, the fluid (water) is treated as a Newtonian fluid governed by the continuity equation and the Navier–Stokes equation. Our main result is that, for pyrex walls, the reduced model approximates the full model reasonably well for sufficiently thick walls, but fails for thin walls, while for PDMS walls, the lossy-wall boundary condition fails regardless of the wall thickness.

2. Results: Comparing the Full and Reduced 2D Models

In the following, we present our results for the numerical simulations of the acoustic fields in the reduced and full models with SAW actuation, and we compare the two cases. As the microchannels are long and straight along the x-direction, we assume translational invariance along x and restrict the calculational domain to the two-dimensional (2D) cross section in the $yz$ plane. The full model consists of coupled fluid and solid domains, whereas the reduced model consists of a single fluid domain with boundary conditions that in an approximate manner represent the walls. The principle of our model approach is illustrated in Figure 1, while the models are described in detail in Section 4.

2.1. Pyrex Devices: Full Model and Reduced Hard-Wall Model

We consider first the full model of a pyrex microdevice, in which a rectangular water-filled channel of width w and height h is encased by a pyrex wall of height $h+H$ and width $w+H$ (see Figure 1b). We simulate the case of actuating the system both at the horizontal standing half-wave resonance in the water $f_{\mathrm{res}}=c_0/2w=1.24 \mathrm{M}\mathrm{Hz}$ often exploited in experiments, and at the off-resonance frequency $f_{\mathrm{off}}=6.65 \mathrm{M}\mathrm{Hz}$ chosen to facilitate comparisons with the literature [36]. An example of a full-model result for the velocity field $-\mathrm{i}\mathsf{\omega }\mathit{u}$ and relative volume change $|\mathsf{\nabla }·{\mathit{u}}_{\mathrm{s}}|$ in the pyrex as well as ${\mathit{v}}_{\mathrm{f}}$ and $p_{\mathrm{f}}$ in the water, is shown in Figure 2a,c.

We then investigate to which extent the full model can be approximated by the reduced hard-wall model often used in the literature [29,38], where only the water domain is considered, while the pyrex walls are represented by the hard-wall condition. In Figure 3, we show for both off-resonance (left column) and on-resonance (right column) actuation, a qualitative comparison between the reduced and the full model, with wall thickness H ranging from 60 to 1800 $\mathsf{\mu }$m. Considering the resulting amplitudes $|p_{\mathrm{f}}|$ of first-order pressure field $p_{\mathrm{f}}$ in the water domain, we note that, for off-resonance actuation at the frequency $f_{\mathrm{off}}$, the full model with thick walls $H=1500 \mathsf{\mu }\mathrm{m}$ has some features in common with the reduced model.

There are pressure anti-nodes in the corners and an almost horizontal pressure node close to the horizontal centerline. For decreasing wall thickness H in the full model, the pressure field changes qualitatively, as the pressure anti-nodes detach from the side walls and shift towards the center of the fluid domain. When actuated on resonance at the frequency $f_{\mathrm{res}}$, for wall thicknesses as low as $H=120 \mathsf{\mu }\mathrm{m}$, the full-model pressure is nearly indistinguishable from that of the hard-wall reduced model, namely a cosine function with vertical pressure anti-nodal lines along the side walls and a vertical pressure nodal line in the center. For the smallest wall thickness $H=60 \mathsf{\mu }\mathrm{m},$ the iso-bars in the full model tilt relative to vertical. In summary, the correspondence between the full and the reduced model is overall better for on-resonance actuation, but for a large wall thickness, the reduced hard-wall model describes the full pyrex model reasonably well.

Table 1. The ratio $H/\lambda$ for various values of the material thickness H and the different acoustic wavelengths λ present in the system at the two frequencies $f_{\mathrm{res}}$ and $f_{\mathrm{off}}$.

${\mathit{f}}_{\mathbf{res}}=1.24$ MHz					${\mathit{f}}_{\mathbf{off}}=6.65$ MHz
$\mathit{\lambda }\mathbf{(}\mathsf{\mu }$m)		H ($\mathsf{\mu }$m)			$\mathit{\lambda }\mathbf{(}\mathsf{\mu }$m)		H ($\mathsf{\mu }$m)
		60	600	1500			60	600	1500
${\lambda }^{\mathrm{SAW}}$	600	0.100	1.000	2.500	${\lambda }^{\mathrm{SAW}}$	600	0.100	1.000	2.500
${\lambda }^{\mathrm{wa}}$	1200	0.050	0.500	1.250	${\lambda }^{\mathrm{wa}}$	225	0.267	2.667	6.667
${\lambda }_{\mathrm{T}}^{\mathrm{pyrex}}$	2745	0.022	0.219	0.546	${\lambda }_{\mathrm{T}}^{\mathrm{pyrex}}$	515	0.117	1.165	2.913
${\lambda }_{\mathrm{L}}^{\mathrm{pyrex}}$	4483	0.013	0.134	0.335	${\lambda }_{\mathrm{L}}^{\mathrm{pyrex}}$	841	0.071	0.713	1.784
${\lambda }_{\mathrm{T}}^{\mathrm{PDMS}}$	80	0.750	7.500	18.750	${\lambda }_{\mathrm{T}}^{\mathrm{PDMS}}$	15	4.000	40.000	100.000
${\lambda }_{\mathrm{L}}^{\mathrm{PDMS}}$	826	0.073	0.726	1.816	${\lambda }_{\mathrm{L}}^{\mathrm{PDMS}}$	155	0.387	3.871	9.677

shows the values of the thickness-to-wavelength ratios $H/\lambda$.

Finally, in the bottom row of Figure 3, we investigate for the full pyrex model model the displacement at the upper boundary in units of the imposed displacement amplitude $u_0$ at the SAW-actuated lower boundary. If the hard-wall condition of the reduced model is good, this displacement should be very small. However, from the figures it is clearly seen that for the thin wall $H=60 \mathsf{\mu }\mathrm{m}$, the upper-wall displacement is significant, with an amplitude of $4u_0$ at $f_{\mathrm{off}}$ and $2u_0$ at $f_{\mathrm{res}}$. As the wall thickness H increases, the upper-wall displacement amplitudes decreases towards $u_0$. Again, this reflects that the reduced hard-wall model is in fair agreement with the full model for a large wall thickness H, and it is better on resonance, where the specific values at the boundaries are less important as the pressure field is dominated by the pressure eigenmode that does in fact fulfill the hard-wall condition (see Section 3.2).

2.2. PDMS Devices: Full Model and Reduced Lossy-Wall Model

We then move on to show the same comparisons, but where the full model has PDMS walls, and the reduced model has the lossy-wall boundary condition, which takes deformation in the normal direction of the wall into account. The reduced lossy-wall model for PDMS, actuated at the off-resonance frequency $f_{\mathrm{off}}=6.65$ MHz, is exactly the one used by Nama et al. [36]. Given the low impedance ratio $\tilde{Z}=0.7$ between PDMS and water, there is no resonance. Results for the full PDMS model are shown in Figure 2b,d at $f_{\mathrm{res}}=1.24$ MHz, and results at $f_{\mathrm{off}}=6.65$ MHz for the reduced lossy-wall model, and the full PDMS model is shown in Figure 4 with plots similar to the ones in the left column of Figure 3 for the reduced hard-wall model and the full pyrex model.

Initially, we compare in Figure 4a–d the amplitude $|p_{\mathrm{f}}|$ of the first-order pressure field $p_{\mathrm{f}}$ of the reduced lossy-wall model with that of the full PDMS model for the wall thickness H varying from 60 to 1500 $\mathsf{\mu }$m. Due to the lossy-wall boundary condition (Equation (14)), the ellipsoidal pressure anti-nodes in Figure 4a traverse the fluid domain upwards during one oscillation cycle. This is in stark contrast to the pressure structures of the full PDMS model in Figure 4b–d, which are stationary due to the free stress condition (Equation (12)) imposed on the exterior of the PDMS. Moreover, the pressure structure of the reduced lossy-wall model consists of only two pressure antinodes, which is much simpler than the multi-node structure of the full PDMS model. In fact, the only common feature in the pressure fields is the appearance of a well-defined pressure node along the vertical centerline.

The poor qualitative agreement between the pressure field in the reduced lossy-wall model and in the full PDMS model is further supported in Figure 4e, where the upper-wall displacement amplitudes of the models are shown. We introduce the unit $u_{\max }$ as the maximum displacement along the upper-wall in the reduced lossy-wall model, and note that the lossy-wall condition imposes a broad single-node sinusoidal velocity amplitude of unity magnitude, while each of the four full model cases ($H=60$, 600, and 1500 $\mathsf{\mu }\mathrm{m}$) shows a rippled, multi-peaked displacement amplitude of relative magnitudes ranging from 2 to 6. The ripples are caused by the small wavelength (15 $\mathsf{\mu }\mathrm{m}$) of the transverse waves in PDMS at the given frequency (see Table 1).

3. Discussion

3.1. Physical Limitations of the Hard-Wall Condition

As illustrated in Figure 3, there are clear discrepancies between the fields obtained by the reduced hard-wall model and those found using the full pyrex models. This can likely be attributed to two factors in particular: the finite stiffness and density of pyrex, and the non-local SAW actuation imposed along the bottom edge in the model.

The hard-wall condition is physically correct for an infinitely stiff and dense wall, which does not undergo any deformation or motion regardless of the stress exerted by the fluid. A hard wall thus reflects all acoustic energy incident on it back into the fluid. However, pyrex has a finite stiffness and density, and it will thus deform and allow for a partial transmittance of acoustic energy from the fluid. This aspect is part of the full pyrex model, but not of the reduced hard-wall model.

The specific SAW actuation is also different in the full and the reduced model. The microdevice rests on top of the piezoelectric substrate, so in the full model, the standing SAW along the surface of the piezoelectric substrate (typically lithium niobate) will transmit significant amounts of acoustic energy directly into both the pyrex wall and the water, but only the latter is taken into account in the reduced hard-wall model. The coupling between lithium niobate and pyrex is strong since the direction-dependent elastic stiffness coefficients of lithium niobate lies in the range from 53 to 200 $\mathrm{G}\mathrm{Pa}$ [39] and the Young’s modulus of pyrex of 64 $\mathrm{G}\mathrm{Pa}$ lies in the same range [40]. Consequently, the interface between the pyrex wall and the water will move under the combined action of the acoustic fields loaded into the pyrex and the water, respectively.

3.2. Acoustic Eigenmodes

Due to the high impedance ratio $\tilde{Z}=8.4$ for pyrex relative to water (see

Table 2. List of parameters used for geometry, materials, and surface acoustic waves (SAW) in the numerical model. The values for the damping paramters Γ are from Reference [41].

Quantity	Symbol	Unit	Pyrex	Polydimethylsiloxane (PDMS)	Water	SAW
			[40]	[42,43]	[44]	[36]
Width	$\frac{1}{2}H$ or w	$\mathsf{\mu }\mathrm{m}$	30–900	30–750	600	-
Height	H or h	$\mathsf{\mu }\mathrm{m}$	60–1800	60–1500	125	-
Density	${\mathsf{\rho }}_{\mathrm{f}}$ or ${\mathsf{\rho }}_{\mathrm{s}}$	$\mathrm{kg}·{\mathrm{m}}^{-3}$	2230	1070	997	-
Bulk modulus	$K_{\mathrm{f}}$ or $K_{\mathrm{s}}$	$\mathrm{GPa}$	38.46	1.12	2.23	-
Longitudinal sound speed	$c_{\mathrm{L}}$ or $c_0$	$\mathrm{m}·{\mathrm{s}}^{-1}$	5591	1030	1496	-
Transversal sound speed	$c_{\mathrm{T}}$	$\mathrm{m}·{\mathrm{s}}^{-1}$	3424	100	-	-
Damping coefficient	${\mathrm{\Gamma }}_{\mathrm{f}}$ or ${\mathrm{\Gamma }}_{\mathrm{s}}$	1	0.001	0.001	0.002	0
Acoustic impedance ratio	$\tilde{Z}=\frac{{\mathsf{\rho }}_{\mathrm{s}}{c}_{\mathrm{L}}}{{\mathsf{\rho }}_{\mathrm{f}}{c}_0}$	1	8.4	0.7	1	-
SAW wavelength	${\lambda }^{\mathrm{SAW}}$	$\mathsf{\mu }\mathrm{m}$	-	-	-	600
SAW displacement amplitude	$u_0$	$\mathrm{n}\mathrm{m}$	-	-	-	0.1
SAW on-resonance frequency	$f_{\mathrm{res}}=\frac{2w}{c_0}$	MHz	-	-	-	1.24
SAW off-resonance frequency	$f_{\mathrm{off}}$	MHz	-	-	-	6.65

), it is possible in the full pyrex model to excite a resonance in the device at the frequency $f_{\mathrm{res}}=2w/c_0=1.24$ MHz, which is close to the ideal standing half-wave pressure eigenmode of the reduced hard-wall system. At this resonance frequency, the pressure amplitude $\left|p_{\mathrm{f}}\right|$ in the water is several times larger than the pressure amplitude ${\mathsf{\rho }}_{\mathrm{f}}{c}_0\mathsf{\omega }{u}_0$ set by the imposed SAW displacement, and the resonance field mainly depends on the frequency and not significantly on the detailed actuation along the boundary [29]. The full pyrex model and the reduced hard-wall model are therefore expected to be in good agreement at $f_{\mathrm{res}}$, as is verified by the right column in Figure 3.

In contrast, at off-resonance frequencies, such as $f_{\mathrm{off}}=6.65$ MHz in the left column of Figure 3, the detailed actuation does matter. The lower left panel of Figure 3 is an example of this, as it highlights an aspect that restricts the validity of the reduced hard-wall model. For the full model with 60-$\mathsf{\mu }\mathrm{m}$-thick pyrex walls, the maximum displacement along the top boundary of the water domain is approximately four times larger than the displacement amplitude $u_0$ of the imposed SAW boundary condition on the bottom boundary of the water domain. This indicates that the system is actuated close to a structural acoustic eigenmode of the pyrex. An amplification is also seen in the lower right panel of Figure 3 although to a smaller degree. This amplification of boundary displacements brought about by the existence of structural eigenmodes is not taken into account in the reduced hard-wall model.

3.3. Physical Limitations of the Lossy-Wall Condition

The comparison between the reduced lossy-wall model and the full PDMS model in Figure 4 shows a clear mismatch. The most important reasons for this are that the lossy-wall model neglects the actuation of both the solid and fluid domain, and that it neglects the transverse motion of the PDMS along the PDMS-water interface.

As for the hard-wall model, the lossy-wall model neglects the strong direct transfer of acoustic energy from the SAW to the PDMS wall, and the implications are the same: the lossy-wall model underestimates the deformation and motion of the PDMS-water boundaries due to this. Moreover, due to the low impedance ratio $\tilde{Z}=0.7$, there are no strong resonances in the water domain like the one at $f_{\mathrm{res}}$ for which the detailed boundary conditions do not matter.

In contrast to the reduced hard-wall model, some aspects of the deformation and motion of the PDMS-water boundaries are taken into account in the reduced lossy-wall model, as it includes the partial reflection and absorption waves from the water domain with perpendicular incidence on the PDMS wall. While this approach would be a good description of a planar or weakly curving interface between two fluids, where all the acoustic excitation takes place in one of the fluids, it is of limited use in the present system, for three reasons: (1) as discussed above, the acoustic energy is injected by the SAW into both the water and the PDMS domain; (2) the PDMS-water boundary is not planar, but consists of three linear segments joined at right angles; and (3) PDMS is not a fluid, but supports shear waves, which are neglected in the reduced lossy-wall model. These three aspects are all part of the full PDMS model, in which PDMS is described as a linear elastic material supporting both longitudinal and transverse waves.

3.4. Modeling PDMS as a Linear Elastic

When modeling large strains above 0.4 in PDMS, non-linear effects are commonly included using hyperelasticity models in the form of a constitutive relation for the stress and strain for which the elastic moduli depends on the stress instead of being constant. For small strains below 0.4, PDMS becomes a usual linear elastic material [45,46,47,48,49]. The magnitude of the strain in terms of the relative volume change $|\mathsf{\nabla }·{\mathit{u}}_{\mathrm{s}}|$ is shown for our system in Figure 2c,d for $H=60 \mathsf{\mu }\mathrm{m}$. The maximum value of $|\mathsf{\nabla }·{\mathit{u}}_{\mathrm{s}}|$ for PDMS is seen to be $5.59\times {10}^{-6}$, which justifies the use of linear elastics as the governing equations of the PDMS walls in our system. Similarly for pyrex, where the maximum value for $|\mathsf{\nabla }·{\mathit{u}}_{\mathrm{s}}|$ is $8.51\times {10}^{-7}$. The use of linear elasticity is further validated in the literature, where linear elastic models of PDMS yield results comparable to those found when using more complex approaches, such as a Mooney–Rivlin constitutive model [33], a neo–Hookian approach [34], and a Maxwell–Wiechert model [50].

Further simplifications based on neglecting the transverse motion of PDMS, such as modeling it as a fluid [30,35] and applying the lossy wall conditions [36], are not advised, since PDMS does have a non-zero transverse bulk modulus and does support transverse sound waves [42,48,49].

As characterization results for PDMS are scarce in the literature, we had to combine the material parameters found in References [42,43] in our simulations.

4. Materials and Methods

Our modeling is based on the generic device design [4,14] illustrated in Figure 1. The device consists of a long, straight, fluid-filled microchannel surrounded by an elastic solid wall on the sides and top. The microchannel and walls rest on a piezoelectric substrate, along which a standing SAW is imposed as a boundary condition. We assume translational invariance along the axial x direction, and only model the transverse $yz$ plane. We implement 2D numerical models in COMSOL Multiphysics 5.2 (COMSOL, Stockholm, Sweden) [51] using the parameters listed in Table 2. All acoustic fields are treated using an Eulerian description, and they have a harmonic time-dependence of the form ${\mathit{u}}_{\mathrm{s}}(y,z){\mathrm{e}}^{-\mathrm{i}\mathsf{\omega }t}$, such that ${\partial }_t$ becomes $-\mathrm{i}\mathsf{\omega }$, where $\mathrm{i}=\sqrt{-1}$, while $\mathsf{\omega }=2\mathsf{\pi }f$ is the angular frequency and f the frequency of the imposed SAW. For simplicity, we often suppress the spatial and temporal variable and write a field simply as ${\mathit{u}}_{\mathrm{s}}$. Finally, following Hahn and Dual [41], we introduce damping in the fluid and the solid using the complex-valued frequency $(1-\mathrm{i}{\mathrm{\Gamma }}_{\mathrm{m}})\mathsf{\omega }$, where ${\mathrm{\Gamma }}_{\mathrm{m}}$ is the damping coefficient in the medium with the values listed in Table 2. For simplicity, we used ${\mathrm{\Gamma }}_{\mathrm{s}}=0.001$ for both pyrex and PDMS, however this implies a damping length for PDMS longer than the 3 mm given in [36], so our model is only valid for PDMS devices with walls thinner than 3 mm.

In total, four models are set up, all with the imposed SAW as a boundary condition representing the actual piezoelectric lithium niobate substrate: (1) the full pyrex model, Figure 5a, where the solid wall is modeled as a linearly elastic material with the parameters of pyrex, while the fluid is modeled as water; (2) the reduced hard-wall model, Figure 5b, where only the water is modeled, while hard-wall boundary conditions replace the pyrex wall; (3) the full PDMS model, Figure 5a, which is the full pyrex model in which the pyrex parameters are replaced by PDMS parameters; and (4) the reduced hard-wall model, Figure 5b, where only the water is modeled, while lossy-wall boundary conditions replace the PDMS wall.

4.1. Governing Equations

The unperturbed fluid at constant temperature $T=298$ K in the fluid domain is characterized by its density ${\mathsf{\rho }}_0$, viscosity ${\mathsf{\eta }}_0$, and speed of sound $c_0$. The governing equations for the acoustic pressure $p_{\mathrm{f}}$, density ${\mathsf{\rho }}_{\mathrm{f}}$, and velocity ${\mathit{v}}_{\mathrm{f}}$ are the usual mass and momentum equations. The constitutive equation between the acoustic pressure $p_{\mathrm{f}}$ and density ${\mathsf{\rho }}_{\mathrm{f}}$ is the usual linear expression, $p_{\mathrm{f}}=c_0^2{\mathsf{\rho }}_{\mathrm{f}}$. Neglecting external body forces on the fluid, while applying perturbation theory [38] and inserting the harmonic time-dependence, the governing equations and the constitutive equation are linearized to the following first-order expressions:

$$\mathrm{i}(1-\mathrm{i}{\mathrm{\Gamma }}_{\mathrm{f}})\mathsf{\omega }{p}_{\mathrm{f}}={\mathsf{\rho }}_0{c}_0^2\mathsf{\nabla }·{\mathit{v}}_{\mathrm{f}}$$

(1)

$$-{\mathsf{\rho }}_0\mathrm{i}(1-\mathrm{i}{\mathrm{\Gamma }}_{\mathrm{f}})\mathsf{\omega }{\mathit{v}}_{\mathrm{f}}=\mathsf{\nabla }·{\mathsf{\sigma }}_{\mathrm{f}}$$

(2)

$${\mathsf{\sigma }}_{\mathrm{f}}=-p_{\mathrm{f}}\mathit{I}+{\mathsf{\eta }}_0[\mathsf{\nabla }{\mathit{v}}_{\mathrm{f}}+\left(\mathsf{\nabla }{\mathit{v}}_{\mathrm{f}}{)}^{\mathrm{T}}\right]+\mathsf{\beta }{\mathsf{\eta }}_0\mathsf{\nabla }\left(\mathsf{\nabla }·{\mathit{v}}_{\mathrm{f}}\right)\mathit{I}$$

(3)

where we have introduced the Cauchy stress tensor ${\mathsf{\sigma }}_{\mathrm{f}}$, and where superscript ”T” denotes tensor transpose, $\mathsf{\beta }$ is the bulk-to-shear viscosity ratio, and $\mathit{I}$ is the unit tensor. With appropriate boundary conditions, the first-order acoustic fields $p_{\mathrm{f}}$, ${\mathsf{\rho }}_{\mathrm{f}}$, and ${\mathit{v}}_{\mathrm{f}}$, can be fully determined by Equations (1)–(3). The specific model-dependent boundary conditions are presented and discussed in Section 1 and Section 4.2.

The dynamics in the solid of unperturbed density ${\mathsf{\rho }}_{\mathrm{s}}$ is described by linear elastics through the momentum equation in terms of the displacement field ${\mathit{u}}_{\mathrm{s}}$ and the solid stress tensor ${\mathsf{\sigma }}_{\mathrm{s}}$. The constitutive equation relating displacement and stress is defined using the longitudinal $c_{\mathrm{L},\mathrm{s}}$ and transverse $c_{\mathrm{T},\mathrm{s}}$ speeds of sound of the given solid:

$$-{\mathsf{\rho }}_{\mathrm{s}}{(1-\mathrm{i}{\mathrm{\Gamma }}_{\mathrm{s}})}^2{\mathsf{\omega }}^2{\mathit{u}}_{\mathrm{s}}=\mathsf{\nabla }·{\mathsf{\sigma }}_{\mathrm{s}}$$

(4)

$${\mathsf{\sigma }}_{\mathrm{s}}={\mathsf{\rho }}_{\mathrm{s}}\left[c_{\mathrm{T},\mathrm{s}}^2(\mathsf{\nabla }{\mathit{u}}_{\mathrm{s}}+\mathsf{\nabla }{\mathit{u}}_{\mathrm{s}}^T)+(c_{\mathrm{L},\mathrm{s}}^2-2c_{\mathrm{T},\mathrm{s}}^2)\left(\mathsf{\nabla }·{\mathit{u}}_{\mathrm{s}}\right)\mathit{I}\right]$$

(5)

4.2. Boundary Conditions

For simplicity, the full dynamics of the piezoelectric substrate is not modeled. Instead, the standing SAW is implemented by prescribing displacements ${\mathit{u}}_{\mathrm{pz}}=\left(u_{y,\mathrm{pz}},u_{z,\mathrm{pz}}\right)$ in the y- and z-directions, respectively, on the bottom boundary of our domain using the following analytical expression from the literature [36,52], where the damping coefficient of $116 \mathrm{m}$⁻¹ has been neglected given the small dimensions (<0.002 m) of the microfluidic device [36]:

$$u_{y,\mathrm{pz}}=0.6u_0\left\{sin\left[k\left(\frac{1}{2}w-y\right)+\mathsf{\omega }t\right]+sin\left[k\left(y-\frac{1}{2}w\right)+\mathsf{\omega }t\right]\right\}$$

(6)

$$u_{z,\mathrm{pz}}=-u_0\left\{cos\left[k\left(\frac{1}{2}w-y\right)+\mathsf{\omega }t\right]+cos\left[k(y-\frac{1}{2}w)+\mathsf{\omega }t\right]\right\}$$

(7)

$${\mathit{v}}_{\mathrm{f}}=-\mathrm{i}\mathsf{\omega }{\mathit{u}}_{\mathrm{pz}}\mathrm{imposed}\mathrm{on}\mathrm{the}\mathrm{fluid}\mathrm{at}\mathrm{the}\mathrm{fluid}-\mathrm{SAW}\mathrm{interface},$$

(8)

$${\mathit{u}}_{\mathrm{s}}={\mathit{u}}_{\mathrm{pz}}\mathrm{imposed}\mathrm{on}\mathrm{the}\mathrm{solid}\mathrm{at}\mathrm{the}\mathrm{solid}-\mathrm{SAW}\mathrm{interface},$$

(9)

where $k=2\mathsf{\pi }/\lambda$ is the wavenumber and $u_0$ the displacement amplitude of the SAW.

In the full models, a no-stress condition for $\mathsf{\sigma }$ is applied along the exterior boundary of the solid. On the interior fluid-solid boundaries, continuity of the stress is implemented as a boundary condition on ${\mathsf{\sigma }}_{\mathrm{s}}$ in the solid domain imposed by the fluid stress ${\mathsf{\sigma }}_{\mathrm{f}}$, while continuity of the velocity is implemented as a boundary condition on ${\mathit{v}}_{\mathrm{f}}$ in the fluid domain imposed by the solid velocity $-\mathrm{i}\mathsf{\omega }{\mathit{u}}_{\mathrm{s}}$. Along the free surfaces of the solid, a no-stress condition is applied:

$${\mathit{n}}_{\mathrm{s}}·{\mathsf{\sigma }}_{\mathrm{s}}={\mathit{n}}_{\mathrm{s}}·{\mathsf{\sigma }}_{\mathrm{f}}\mathrm{imposed}\mathrm{on}\mathrm{the}\mathrm{solid}\mathrm{at}\mathrm{the}\mathrm{fluid}-\mathrm{solid}\mathrm{interface},$$

(10)

$${\mathit{v}}_{\mathrm{f}}=-\mathrm{i}\mathsf{\omega }{\mathit{u}}_{\mathrm{s}}\mathrm{imposed}\mathrm{on}\mathrm{the}\mathrm{fluid}\mathrm{at}\mathrm{the}\mathrm{fluid}-\mathrm{solid}\mathrm{interface},$$

(11)

$${\mathit{n}}_{\mathrm{s}}·{\mathsf{\sigma }}_{\mathrm{s}}=\mathbf{0}\mathrm{imposed}\mathrm{on}\mathrm{the}\mathrm{solid}\mathrm{at}\mathrm{exterior}\mathrm{boundaries}.$$

(12)

We have also performed simulations, where the stress-free condition Equation (12) on the exterior boundaries is changed into a lossy-wall conditions involving high acoustic impedance ratios (PDMS/air 3600 and pyrex/air 41000). As expected, the resulting fields are almost unchanged: we observe the same morphology, and, quantitatively, the average value pressure field in the water domain exhibits relative changes less than $4\times {10}^{-5}$, hence we employ the simpler Equation (12).

In the reduced models, boundary conditions are imposed on the fluid to represent the surrounding material. Stiff and heavy materials such as pyrex are represented by the hard-wall (no motion) condition at the boundary of the fluid domain. Soft and less heavy materials such as PDMS are represented by the lossy-wall condition for partial acoustic transmittance perpendicular to the boundary of the fluid domain. For both conditions, a no-slip condition is applied on the tangential velocity component. The specific expression implemented in COMSOL are:

$${\mathit{v}}_{\mathrm{f}}=\mathbf{0}\mathrm{boundary}\mathrm{condition}\mathrm{representing}\mathrm{hard}\mathrm{walls},$$

(13)

$${\mathit{v}}_{\mathrm{f}}=\frac{p_{\mathrm{f}}}{c_s{\mathsf{\rho }}_s}\mathit{n}\mathrm{boundary}\mathrm{condition}\mathrm{representing}\mathrm{lossy}\mathrm{walls}.$$

(14)

4.3. Numerical Implementation and Validation

We follow our previous work [29,53], and implement the governing equations in weak form in the commercial software COMSOL Multiphysics 5.2 [51]. To fully resolve the thin acoustic boundary layer of width $\mathsf{\delta }$,

$$\mathsf{\delta }=\sqrt{\frac{2{\mathsf{\eta }}_0}{{\mathsf{\rho }}_0\mathsf{\omega }}}=0.21 \mathsf{\mu }\mathrm{m}\mathrm{at}\mathsf{\omega }=2\mathsf{\pi }\times 6.5\mathrm{MHz}$$

(15)

in the water domain near its edges, the maximum mesh size $h_{\mathrm{edge}}$ at the solid-fluid boundary is much smaller than that in the bulk called $h_{\mathrm{bulk}}$. Both of these are controlled by the mesh parameter $k_{\mathrm{p}}$,

$$h_{\mathrm{edge}}=\frac{1}{k_{\mathrm{p}}}\mathsf{\delta }{h}_{\mathrm{bulk}}=50h_{\mathrm{edge}}\mathrm{with}{k}_{\mathrm{p}}=1\mathrm{in}\mathrm{the}\mathrm{main}\mathrm{runs}.$$

(16)

The coarse mesh with $k_{\mathrm{p}}=0.1$ is shown in Figure 6a. In our largest (full) models using $k_{\mathrm{p}}=1$, the implementation resulted in $8.1\times {10}^6$ degrees of freedom and a computational time of 30 min on a standard PC work station. The implementation of the model in the fluid domain has been validated both numerically and experimentally in our previous work [29,53]. The solid domain implementation was validated by calculating resonance modes for a long rectangular cantilever, clamped at one end and free at the other, and comparing them successfully against analytically known results. Finally, for both the full and the reduced models, we performed a mesh convergence analysis using the relative mesh convergence parameter $C\left(g\right)$ for a given field $g(y,z)$ as introduced in Reference [29]:

$$C\left(g\right)=\sqrt{\frac{{\int }_{\mathsf{\Omega }}(g-g_{\mathrm{ref}}{)}^2\mathrm{d}y\mathrm{d}z}{{\int }_{\mathsf{\Omega }}(g_{\mathrm{ref}}{)}^2\mathrm{d}y\mathrm{d}z}}$$

(17)

Here, $g_{\mathrm{ref}}$ is the solution obtained with the finest possible mesh resolution, in our case the one with mesh parameter $k_{\mathrm{p}}=5$. For all fields, our mesh analysis revealed that satisfactory convergence was obtained with the mesh parameter set to $k_{\mathrm{p}}=1$. For this value, the relative mesh convergence parameter was both small, $C\approx 0.01$, and exhibited an exponential asymptotic behavior, $C\simeq {\mathrm{e}}^{-k_{\mathrm{p}}}$, as a function of the mesh parameter $k_{\mathrm{p}}$ (for two examples, see Figure 6b,c).

5. Conclusions

A numerical method has been presented for 2D full modeling of a generic SAW microdevice consisting of a long, straight, fluid-filled microchannel encased in a elastic wall and resting on a piezoelectric substrate in which a low-MHz-frequency standing SAW is imposed. We have also presented reduced models consisting only of the fluid domain, where boundary conditions are used as simplified representations of the elastic wall. An acoustically hard wall, such as pyrex, is represented by a hard-wall boundary condition, while an acoustically soft wall, such as PDMS, is represented by a lossy-wall boundary condition. Our results show that the full pyrex model is approximated fairly well for thick pyrex walls using the hard-wall model, when the SAW is actuated on a frequency corresponding to a resonance frequency of the water domain, but less well for thinner walls at resonance and for any wall thickness off resonance. The reduced lossy-wall model was found to poorly approximate the full PDMS model for walls thinner than the 3-mm PDMS damping length, especially regarding the resulting running pressure waves in the reduced lossy-wall model in contrast the standing waves in the full PDMS model.

Modeling of acoustofluidic devices should thus be performed in full to take into account all effects relating to the elastic walls defining the microchannel. At higher frequencies or higher acoustic power levels, even the full model presented here must be extended to take into account thermoviscous effects in the form of increased heating and temperature-depending effects [44,54]. Finally, to obtain quantitatively better results for the pressure fields driving acoustophoresis in the water domain, the piezoelectric substrate should be included in future simulations. Hopefully, such an analysis will be of interest to experimentalists, who in turn may provide improved experimental data to validate the model. Moreover, with such an extended model including the dynamics of the piezoelectric substrate, a study could be carried out for other actuation conditions than the ones studied here, such as bulk acoustic wave actuation, for which more experimental results exist.