Efficient modeling of sharp-edge acoustofluidics

Abstract

Sharp-edge structures exposed to acoustic fields are known to produce a strong non-linear response, mainly in the form of acoustic streaming and acoustic radiation force. The two phenomena are useful for various processes at the microscale, such as fluid mixing, pumping, or trapping of microparticles and biological cells. Numerical simulations are essential in order to improve the performance of sharp-edge-based devices. However, simulation of sharp-edge structures in the scope of whole acoustofluidic devices is challenging due to the thin viscous boundary layer that needs to be resolved. Existing efficient modeling techniques that substitute the need for discretization of the thin viscous boundary layer through analytically derived limiting velocity fail due to large curvatures of sharp edges. Here, we combine the Fully Viscous modeling approach that accurately resolves the viscous boundary layer near sharp edges with an existing efficient modeling method in the rest of a device. We validate our Hybrid method on several 2D configurations, revealing its potential to significantly reduce the required degrees of freedom compared to using the Fully Viscous approach for the whole system, while retaining the relevant physics. Furthermore, we demonstrate the ability of the presented modeling approach to model high-frequency 3D acoustofluidic devices featuring sharp edges, which will hopefully facilitate a new generation of sharp-edge-based acoustofluidic devices.

1 Introduction

Acoustofluidic devices commonly exploit two phenomena, namely, acoustic streaming (AS) and acoustic radiation force (ARF), for an increasing variety of tasks due to their ability to precisely manipulate fluids and objects on the microscale and in a contactless manner [1]. A particular type of acoustofluidic device features sharp-edge structures inside an acoustically excited device. In such a configuration, sharp edges produce a strong response in the form of AS [2–5] and ARF [6, 7], as illustrated in Figure 1. The literature indicates that in most 2D cases, the AS around a single sharp edge manifests itself in the form of two counter rotating vortices, one on each side of the tip of the sharp edge, with the flow above the apex directed away from the edge tip [4]. As for the ARF, it has been shown that “heavy” microparticles get attracted to the excited sharp edge, whereas “light” particles get repelled from the edge [7].

FIGURE 1

Acoustically excited sharp edges can be used for micropumping [8, 9], trapping of biological cells [6, 10, 11], and powerful mixing [12–14]. Numerical models capable of predicting the behavior inside such devices and, thus, their performance are a valuable tool for design optimization. Especially, since the production of most acoustofluidic devices requires a cleanroom environment, making prototyping laborious and expensive. Furthermore, numerical simulations can give deeper insight into the underlying physics [2, 4, 5].

The modeling of harmonic acoustic fields at the scale of a typical acoustofluidic device for micromanipulation requires the modeling of one-half to several acoustic wavelengths λ = c₀/f, with the speed of sound in the fluid c₀ and the driving frequency f. Each wavelength is normally discretized with several tens of elements, making such a device possible to simulate as a whole 3D multiphysics problem. However, when analyzing higher-order steady phenomena in viscous fluids, such as AS, another characteristic length needs to be discretized-the viscous boundary layer thickness , with the dynamic fluid viscosity η and the density of the fluid ρ₀. In a typical acoustofluidic device, filled with water and driven at f = 1 MHz, we have λ ≈ 1.5 mm and δ = 0.56 µm. This indicates a large separation of length scales, with the device being , whereas the discretization at boundaries is required to resolve . Muller et al. [15] introduced a Fully Viscous modeling approach based on the perturbation theory that resolves both λ and δ, requiring significant computational effort and is, thus, in most cases limited to simplified 2D models. Over the years, alternative approaches were developed that analytically solve AS in the viscous boundary layer and use the solution as a boundary condition for the flow in the fluid bulk. This so-called limiting velocity method (LVM) was already used in [16], [17, 18], [19], [20, 21], and others. Recently, [22] and [23] improved the LVM by accounting for the Stokes drift and thermal boundary layer, respectively. The LVM approach has already enabled 3D simulations of AS in a whole acoustofluidic device [19, 21, 24]. There is, however, a general limitation of the LVM, in which the radius of the curvature of boundaries in the system needs to be much larger than δ. This makes LVM as such unusable for the simulation of AS in acoustofluidic devices featuring sharp edges, since the behavior in such devices is greatly influenced by the behavior near the tips of sharp edges [5].

Sharp-edge acoustofluidics has been successfully modeled in several studies using the Fully Viscous approach, but typically in 2D [4, 13, 25, 26] and under further simplifications, such as assuming a periodicity across the device [3], using systems at low frequencies with the acoustic wavelength larger than the system (subwavelength system) [2, 5, 27], or exploiting symmetries [28]. Under these simplified conditions, a few studies featured the direct numerical simulations that revealed the threshold amplitude of the oscillatory excitation velocity at which the perturbation-based computation of the acoustic streaming near a sharp edge fails [2, 5, 29].

Here, we introduce a Hybrid modeling approach that fully resolves the viscous boundary layers in the vicinity of sharp edges, analogous to the approach of [15], while the rest of the acoustofluidic device is modeled using LVM-based approach derived by [22]. First, we validate the Hybrid approach against the Fully Viscous approach on a simple single-sharp-edge geometry, which is used in several previous studies [4, 5, 26]. Second, we apply the Hybrid approach to model a more complex acoustofluidic device that features a pair of sharp edges and is capable of pumping, as introduced recently in [9]. At last, we apply the Hybrid approach to model high-frequency sharp-edge acoustofluidics in 3D. The developed approach, depending on the complexity of the device, significantly decreases the required degrees of freedom (DOF) in the model compared to fully resolving the viscous boundary layer.

2 Methodology

In the current work, we model only the behavior of fluids within microfluidic channels, whereas the solids typically bounding such channels are replaced by suitable boundary conditions.

We assume a viscous fluid, the motion of which can be described by the compressible Navier–Stokes equationsand the continuity equationwith the velocity v, the dynamic viscosity η, and the bulk viscosity η_B. The fluid is assumed to be barotropic, and the density ρ, therefore, depends only on the pressure p, namely,

The equations are linearized using the regular perturbation approach [30]. Accordingly, the physical fields are expanded in a series □ = □₀ + □₁ + □₂ + …, where □ represents the field and the subscript denotes the respective order. The amplitude of the first-order velocity v₁ is, therefore, assumed to be small with respect to the speed of sound c₀—small Mach number assumption.

2.1 First-order (acoustic) problem

For a fluid quiescent at the zeroth order (v₀ = 0), the substitution of the first-order perturbed fields into governing Eqs. 1, 2 yields the set of first-order equationswith the equilibrium density ρ₀. The equation of state, namely Eq. 3, translates toand connects the first-order density with the first-order pressure. The first-order fields are assumed to have a harmonic time dependency through the factor e^iωt, with the angular frequency ω = 2πf and the imaginary unit i.

2.2 Second-order (streaming) problem

Applying the perturbation theory to the governing equations up to second order and taking the time average over an oscillation period T lead to the equations of acoustic streaming [31]The time average of a product of first-order fields , with the spatially-dependent and time-independent complex amplitude , is computed as with □* denoting a complex conjugate of □ [32]. For in Eq. 7, the tensor product is implied. In Eqs 7, 8, the streaming field is considered to be incompressible, which is in line with related studies [33–36].

At the second order, the no-slip boundary condition is imposed on the Lagrangian velocity of the fluid at the fluid boundary. The Lagrangian velocity is defined as the sum of the Eulerian streaming velocity and the Stokes drift [37, 38]leading to the boundary condition on the Eulerian streaming velocityFor rigid boundaries that we assume in this work, v_SD reduces to zero.

This theoretical framework is only applicable as long as the perturbation theory is valid and as long as the streaming remains laminar. The validity of the perturbation theory approach in the context of sharp-edge streaming is further discussed in [2], [5].

2.3 Limiting velocity method

To compute the acoustic streaming one needs to solve the first- and the second-order problems outlined in the previous sections. Based on the complexity of the formulation, analytical solutions are rare and rely on basic geometries, simplified boundary conditions, and various assumptions regarding the relevant characteristic lengths, as, for example, performed in [16], [39], [34], [40]. However, under additional assumptions, the so-called limiting velocity method has been developed by [18], [38], [22], wherein the aforementioned equations are solved separately in the vicinity of boundaries, where the viscous boundary layer develops. The main assumption of the LVM is the smallness of the viscous boundary layer thickness δ relative to the acoustic wavelength λ = c₀/f and relative to the curvature R of the boundary surface. Additionally, the displacement of a moving boundary surface tangential to the boundary surface needs to be small relative to R, while the displacement normal to the boundary surface needs to be small relative to δ. The derivation of the fields in the boundary layer assumes a solenoidal first-order velocity field, and the resulting (limiting) streaming velocity is then applied to a simplified streaming problem originating from Eqs 7, 8 as a slip velocity at certain small distance from the boundary [18] or at the boundary itself [22], depending on the formulation. It is important to note that the derivations in [22] are developed for e^−iωt time dependence, whereas the interfaces of the Acoustics Module of COMSOL Multiphysics that we use for implementation assume e^iωt. In the continuation, the formulation from [22] is, therefore, adapted to the COMSOL Multiphysics-compatible e^iωt.

The first-order problem, in the scope of the LVM, is simplified by assuming that the velocity field in the fluid bulk is irrotational, which leads to the Helmholtz equationwith p₁ as the only unknown variable and with the viscous compressional wavenumberwith the damping coefficientThe first-order velocity is separated asinto a sum of the irrotational (∇ ×□ = 0) long-range velocity in the bulk of the fluid and the solenoidal (∇ ⋅□ = 0) short-range velocity in the viscous boundary layer . The velocity in the bulk is connected to the pressure throughTo account for the damping of the viscous boundary layer, the following condition needs to be applied at the boundary:with the outward pointing normal of the unit length n, the velocity of the harmonic boundary , and the shear wavenumber

The second-order problem for the streaming in the bulk of the fluid is in LVM governed bywith the bulk Lagrangian densityand the energy-flux density

The boundary condition that constrains the second-order problem is given in a form of a slip velocitywith orthogonal unit vectors t₁ and t₂ that are tangential to the boundary surface and withwith ; follows from the no-slip boundary condition at the first order asThe superscript □⁰ indicates that □ is evaluated at the boundary.

2.4 Numerical model

The numerical models implementing the first- and the second-order problems are developed in the finite-element method framework of COMSOL Multiphysics 5.6 [41]. In the analysis that follows, we use three kinds of numerical models, namely, (I) the Fully Viscous model that numerically solves Eqs. 4–8, 10 in line with the approach introduced in [15] and used in many studies since [5, 9, 28, 42], (II) the Full Limiting Velocity Method (FLVM) model that is based on the LVM of [22] and incorporates Eqs. 11, 16, 18, 19, and 22, and (III) the Hybrid model that uses the formulation of the Fully Viscous model in the vicinity of sharp edges and the LVM formulation in the rest of the fluid domain, as illustrated in Figure 2.

FIGURE 2

2.4.1 Fully viscous model

The Fully Viscous model is implemented in COMSOL Multiphysics by using the adiabatic form of the Thermoviscous Acoustics interface for the first-order problem, which is then solved with a Frequency Domain study. The second-order problem is implemented using the Creeping Flow interface, with the streaming source term from Eq. 7, which is the spatial variation of the Reynolds stress, and added to the governing equations of the Creeping Flow interface as a volume force. The second-order problem is solved with the Stationary study.

The boundary conditions in the Fully Viscous model at the first order are the no-slip condition applied on all the non-moving (rigid) walls; specifically,while at the moving walls, the applied wall velocity is imposed aswith the velocity v_wall defined on the case-by-case basis. At the second order, the no-slip boundary condition is applied to all the walls, but on the Lagrangian velocity, as in Eq. 10. However, the Stokes drift that appears in Eq. 10 only at the moving walls.

2.4.2 Full Limiting Velocity Method model

The first-order problem of the FLVM model is implemented through the Pressure Acoustics interface, with the boundary condition on the pressure gradient from Eq. 16 that accounts for viscous boundary layer damping imposed as an inward velocity at all boundaries. The first-order problem is solved with a Frequency Domain study. The second-order problem is implemented via the Creeping Flow interface, same as for the Fully Viscous model. The source term from Eq. 19 is implemented as a volume force.

The boundary conditions of the FLVM model at the first order are the slip condition with the additional artificial velocity in the direction normal to the wall that imposes viscous damping resulting from the viscous boundary layer; specifically,where is defined in Eq. 16. The case-specific oscillatory velocity of the moving walls is imposed directly in the expression from Eq. 16, by setting , with the case-specific v_wall. At the second order, the slip velocity defined in Eq. 22 is imposed to the fluid at the walls. Specifically,For the 2D cases, the boundary condition in Eq. 29 simplifies as t₂ = 0.

2.4.3 Hybrid model

In the Hybrid model, the Fully Viscous domain is implemented analogously to the implementation of the Fully Viscous model, whereas the LVM domain is implemented analogously to the FLVM model. The Thermoviscous Acoustics interface and the Pressure Acoustics interface are coupled through the Acoustic-Thermoviscous Acoustic Boundary interface. In the Creeping Flow interface, the gradient of the Lagrangian density is omitted from the LVM domain for numerical stability, since it does not contribute to the streaming [43]. Consequently, the gradient of the Lagrangian density needs to be subtracted from the streaming source term in the Fully Viscous domain, to avoid a non-physical source of streaming at the interface of the LVM and the Fully Viscous domains. In the LVM domain, we also omit the streaming source from Eq. 19 for numerical stability, as it contributes to bulk-driven streaming, which is insignificant in the cases analyzed here, but might be important in some other configurations [44].

The boundary conditions in the Fully Viscous domain of the Hybrid model resemble the boundary conditions applied in the Fully Viscous model, while the boundary conditions imposed in the LVM domain of the Hybrid model resemble those in the FLVM model. There are two exceptions to this at the first order: The walls in the Fully Viscous domain are coupled to the LVM domain, where the slip velocity is applied (indicated in Figure 2) through the following conditions:withwith the identity tensor I; the second exception is at the coupling boundary between the Fully Viscous domain and the LVM domain, where the following coupling conditions are applied:with the left-hand side evaluated in the LVM domain and the right-hand side in the Fully Viscous domain, at the coupling boundary. At the second order, there is no special coupling between the two domains, except for the volume force accounting for the spatial variation of the Reynolds stress on the right-hand side of Eq. 7, which is applied only in the Fully Viscous domain and not in the LVM domain.

The size of the LVM domain in the Hybrid model, here defined through the length b in Figure 2, should be chosen such that it surrounds the region of the high curvature of the boundary, satisfying the condition b ≫ δ under the assumption that δ ≳ 1/R.

For all the models, the Creeping Flow interface contains a pressure constraint in the form of a zero average across the whole fluid domain. All the studied 2D cases were solved with the default solvers, suggested by COMSOL Multiphysics. For the 3D cases featured here, the solver had to be manually switched to a Direct solver.

It is important to note that the first-order problem implemented using COMSOL Multiphysics through the native Thermoviscous Acoustics or Pressure Acoustics interfaces needs to be consistent with the time dependency with a factor of e^iωt, whereas the acoustofluidic theories often assume e^−iωt, as, for example, assumed in [22], [23].

For all the studies herein, we assume the fluid to be water with ρ₀ = 1000 kg m⁻³, c₀ = 1481 m s⁻¹, η = 1.02 mPa s, and η_B = 3.09 mPa s [45].

Exemplary 2D meshes for the Fully Viscous model and the Hybrid model are compared in Figure 3. The inherently larger bulk element size of the Hybrid model and the absence of mesh refinements at the boundaries promise a substantial reduction in the computational effort for the Hybrid model compared to the Fully Viscous model, especially for larger problems. The larger bulk element size admissible in the LVM domain, compared to the Fully Viscous domain, stems from the lack of the divergence of the Reynolds stress in Eq. 19 that requires a refined mesh [22]. The difference in the total DOF is further increased due to the Fully Viscous model solving the velocity field and the pressure field in the fluid bulk in the first-order problem, whereas the hybrid model only solves for the pressure field in the fluid bulk.

FIGURE 3

The convergence of representative models with respect to the mesh element size is discussed in Supplementary Material.

In the Fully Viscous model, the first-order pressure and velocity were discretized using cubic Lagrange and quartic Lagrange elements, respectively. The second-order pressure and velocity were discretized with quadratic and cubic elements, respectively. In the Hybrid model, cubic Lagrange elements were used for the first-order pressure in the LVM domain, while the pressure and velocity in the Fully Viscous domain were discretized with quadratic and cubic Lagrange elements, respectively. At the second order, the discretization was quadratic and cubic for pressure and velocity, respectively.

The results of numerical simulations are evaluated through the base variables (p₁, v₁, and ) as well as through the time- and space-averaged acoustic energy density [46]within a spatial domain Ω that corresponds to the volume in a 3D or the area in a 2D configuration. The first term within in Eq. 35 represents the potential energy density, while the second term represents the kinetic energy density. In the LVM subdomains of Ω, v₁ is computed from Eq. 15.

3 Results

We study several cases using the described numerical models, which validate the Hybrid model on various 2D geometries in combination with fluidic channel resonances. In the last part, the benefit of the Hybrid model is demonstrated on the analysis of an exemplary 3D case that would have been very difficult if not impossible to model using preexisting approaches.

3.1 Rectangular channel in 2D

We first analyze a simple case of a 2D rectangular fluidic channel that is geometrically equivalent to the case presented in Figure 2, but without the sharp edge. The analysis, presented in Figure 4, compares the FLVM model and the Hybrid model against the Fully Viscous model, in order to validate the implementation in COMSOL Multiphysics and to estimate the baseline error caused by the coupling of Fully Viscous and LVM domains within the Hybrid model. In the Hybrid model, the Fully Viscous domain is defined with length b, while the LVM formulation is applied to the rest of the rectangular domain defined with length a, as specified in Figure 2.

FIGURE 4

The 2D geometry of the rectangular fluidic channel in Figure 4 is frequently used for studies in acoustofluidics [15, 25, 39, 46, 47]. The geometry is defined through the lengths a = 1 mm and b = 200 µm. Predictions of the three models are compared around the first three resonance modes across the width of the rectangular channel in the x-direction. We apply the excitation at x = ±a/2 with v_a = 1 mm s⁻¹ imposed in the same direction at both sides of the channel (in phase) to support the build-up of the resonance modes 1 (Figure 4A) and 3 (Figure 4C) across the channel width. For the resonance mode 2 in Figure 4B, the velocity v_a = 1 mm s⁻¹ is imposed in the opposing directions (π phase shift) at the two sides of the channel.

Qualitatively, the first- and second-order fields look the same for all models and modes; therefore, only reference normalized fields are shown in Figures 4A-F. The frequency sweeps around the three resonance frequencies in Figures 4G-L, which further validate our COMSOL Multiphysics implementation of the LVM, as particularly the FLVM model matches the Fully Viscous model perfectly. While the Hybrid model matches the Fully Viscous model well for all three modes in terms of the corresponding resonance frequency, it overestimates the acoustic energy density E_ac and maximal streaming velocity . The error for the first mode is 7% and 6% for E_ac and , respectively. In the case of the two higher modes, the errors reach up to 76%. The error comes from the coupling of the LVM with Fully Viscous domains and can be minimized by reducing the length of the boundary section that assumes a Fully viscous formulation and by positioning the Fully Viscous domain close to the pressure node (see Supplementary Material for further analysis). However, for frequencies away from the center of the resonance peak, the errors are minimized and all three models agree well.

3.2 Rectangular channel with a sharp-edge structure in 2D

In Figure 5, we look at the same 2D rectangular fluidic channel as in Figure 4, whereas a sharp edge has been added. A comparable geometry is featured in several prior experimental [2, 4, 7, 14] and numerical [2, 5, 25, 26] studies of sharp-edge acoustofluidics.

FIGURE 5

We apply the excitation at x = ±a/2 with v_a = 1 mm s⁻¹ imposed in the same direction at both sides of the channel (in phase) to facilitate the build-up of the first resonance mode across the channel width. The first resonance mode was chosen for the analysis as it induces a pressure node and, thus, a velocity antinode at the location of the sharp-edge tip, producing strong sharp-edge phenomena.

The presence of the sharp edge shifts the resonance frequency of the first mode in Figure 5 from ∼ 740.1 kHz down to ∼ 720.8 kHz, when compared to a channel without the sharp edge in Figure 4. The lower resonance frequency in the rectangular fluidic channel featuring a sharp edge, compared to the channel without the sharp edge, was already observed in simulations in [26]. This frequency shift might be due to the increased path of waves, as they have to circumvent the sharp edge, analogously to the presence of other objects dispersed in an acoustic field that decrease the resonance frequency, for instance, bubbles [30].

Figure 5 reveals a good agreement between the Hybrid model and the Fully Viscous model, whereas the FLVM model quantitatively and qualitatively disagrees with the Fully Viscous model. The latter is expected based on the underlying assumptions of the LVM, and the FLVM model is, therefore, omitted from further analysis. At the resonance peak, the Hybrid model slightly underestimates the magnitude of the first-order fields, reflected in ∼16% lower E_ac maximum compared to the Fully Viscous model. Similarly, the maximal magnitude of the Eulerian streaming velocity is underestimated by ∼16% at the resonance peak. However, the agreement between the Hybrid model and the Fully Viscous model is very good away from the center of the resonance peak. Specifically, the difference in E_ac and is 5% for frequencies 2 kHz | from the center frequency of the peak.

In order to understand how the newly developed Hybrid model is compared to the Fully Viscous model in terms of the computational demand relative to an increasing size of the domain, we keep track of the DOF as the channel width a is increased. The results shown in Figure 6 present E_ac, , and the DOF for the rectangular channel with a sharp edge from Figure 2 at f = 10 kHz and across the channel width range of 0.3 mm ≤ a ≤ 50 mm, making this a subwavelength problem (λ ≈ 150 mm). The results indicate a good agreement between the two models in terms of E_ac and across the whole range of a. Interestingly, the behavior below a ≈ 1 mm depends on a, while it remains independent of a for larger domain sizes. The DOF in the Fully Viscous model, however, significantly exceeds the DOF in the Hybrid model, with the DOF reaching -fold reduction as the domain size is increased, indicating the benefit of the Hybrid model.

FIGURE 6

To compare the scaling of the Hybrid model and the Fully Viscous model at higher frequencies, we vary a from 0.3 mm to 3.5 mm at f = 720 kHz. The results in Figure 7 reveal three resonances in the channel across the investigated range of a, with the first peak at a ≈ 1 mm corresponding to the resonance peak previously analyzed in Figure 5. The two models match very well in terms of E_ac and . The DOF generally exceed those from Figure 6 due to the higher frequency, and the DOF in the Hybrid model again show a significant up to -fold decrease relative to the DOF in the Fully Viscous model.

FIGURE 7

3.3 Sharp-edge-based micropump in 2D

The main advantage of the developed Hybrid model is the ability to model larger systems that feature sharp-edge structures. We analyze the performance of the Hybrid model for a larger system in Figure 8, in the case of a 2D acoustofluidic micropump inspired by [9] that is capable of fluid pumping and mixing due to the sharp-edge streaming, as well as microparticle/cell focusing and trapping by the acoustic radiation force due to the fluidic channel resonances.

FIGURE 8

Similar to the previous results, the behavior predicted by the Fully Viscous model and the Hybrid model matches qualitatively and quantitatively, except for the resonance peak, where the Hybrid model underestimates the magnitude of acoustic fields.

3.4 Rectangular channel with a sharp-edge structure in 3D

Finally, we demonstrate the capability of our efficient modeling approach by solving a high-frequency acoustic streaming problem in 3D using a Hybrid model with the geometry from Figure 2 extruded to a depth of d. The exemplary mesh and results are shown in Figure 9.

FIGURE 9

We apply the 3D model to study the influence of the depth d on the Eulerian streaming pattern on the xy-plane at z = d/2 in Figure 10. The results show that the pattern is visibly affected by the presence of the front (z = d) and back (z = 0) covers when d = 10 µm in Figure 10A, but the pattern gradually converges to the streaming pattern in the 2D model from Figure 5 as d is increased to 100 µm in Figure 10D. However, as shown in Figure 10E, the maximal Eulerian streaming velocity magnitude remains impaired by the 3D effects even at d = 100 μm, when compared to the magnitudes from a comparable 2D configuration at f = 720 kHz from Figure 5C.

FIGURE 10

4 Discussion

We have demonstrated that the novel Hybrid model, combining the Fully Viscous modeling of the acoustic phenomena [15] near sharp edges with the efficient modeling of the viscous boundary layers according to [22] in the rest of the device, provides up to -fold reduction in the required DOF in 2D validation studies while capturing the relevant physics. Furthermore, this approach enables the modeling of 3D sharp-edge phenomena at high frequencies, where the viscous boundary layer thickness δ reduces to the rounding radius of sharp edges or below.

While the novel Hybrid model introduces a quantitative error on the estimated magnitude of fields at the resonance peaks, it correctly predicts the qualitative behavior throughout the modeled device. This is of no concern for problems that do not involve resonances in the fluidic channels, e.g., the configurations from [2, 5, 14]. Sharp-edge-based devices operating at higher frequencies, such as the multifunctional chip reported by [9], can still be effectively modeled, as demonstrated in Figure 8. Shortcomings in terms of erroneous magnitudes at resonances are not too important, since the underlying total damping defining the quality of a given resonance is generally difficult to model anyway and can vary due to a plethora of experimental conditions [19], such as the temperature, material damping, damping in the connecting glue layers, and many more.

The relative error of the Hybrid model is for a given resonance similar for both analyzed quantities, E_ac and . This implies that the error originates from the computation of the first-order fields, since both quantities rely on the square of a first-order quantity—E_ac directly, while depends the square of the first-order velocity as a source term in Eq. 7 and as a boundary condition through Eq. 9.

Our Hybrid model allows the study of 3D sharp-edge effects, as demonstrated in Figure 10, where we showed that the 3D effects in the quasi-2D geometry are limited to a certain depth, as previously hypothesized in [5]. This provides the justification for the use of 2D models when modeling quasi-2D channels of lab-on-a-chip devices that typically exceed the depth of 100 µm.

The herein introduced Hybrid modeling approach can be used to aid the emerging field of sharp-edge acoustofluidics. In the latter, sharp-edge structures such as needles or wedges in microfluidic devices are excited with ultrasound for a wide range of applications, ranging from propulsion of microrobots [48–50] and enhanced drug delivery [27] to advanced lab-on-a-chip technology involving sharp-edge-based micropumps and micromixers [3, 8, 9, 12–14, 51].

References

• 1.

  BruusHDualJHawkesJHillMLaurellTNilssonJet alForthcoming lab on a chip tutorial series on acoustofluidics: Acoustofluidics—Exploiting ultrasonic standing wave forces and acoustic streaming in microfluidic systems for cell and particle manipulation. Lab Chip (2011) 11:3579–80. 10.1039/c1lc90058g

  • CrossRef
  • Google Scholar

• 2.

  OvchinnikovMZhouJYalamanchiliS. Acoustic streaming of a sharp edge. The J Acoust Soc America (2014) 136:22–9. 10.1121/1.4881919

  • CrossRef
  • Google Scholar

• 3.

  NamaNHuangPHHuangTJCostanzoF. Investigation of acoustic streaming patterns around oscillating sharp edges. Lab Chip (2014) 14:2824–36. 10.1039/c4lc00191e

  • CrossRef
  • Google Scholar

• 4.

  DoinikovAAGerltMSPavlicADualJ. Acoustic streaming produced by sharp-edge structures in microfluidic devices. Microfluidics and Nanofluidics (2020) 24:32. 10.1007/s10404-020-02335-5

  • CrossRef
  • Google Scholar

• 5.

  ZhangCGuoXRoyonLBrunetP. Unveiling of the mechanisms of acoustic streaming induced by sharp edges. Phys Rev E (2020) 102:043110. 10.1103/physreve.102.043110

  • CrossRef
  • Google Scholar

• 6.

  LeibacherIHahnPDualJ. Acoustophoretic cell and particle trapping on microfluidic sharp edges. Microfluidics and Nanofluidics (2015) 19:923–33. 10.1007/s10404-015-1621-1

  • CrossRef
  • Google Scholar

• 7.

  DoinikovAAGerltMSDualJ. Acoustic radiation forces produced by sharp-edge structures in microfluidic systems. Phys Rev Lett (2020) 124:154501. 10.1103/physrevlett.124.154501

  • CrossRef
  • Google Scholar

• 8.

  HuangPHNamaNMaoZLiPRufoJChenYet alA reliable and programmable acoustofluidic pump powered by oscillating sharp-edge structures. Lab Chip (2014) 14:4319–23. 10.1039/c4lc00806e

  • CrossRef
  • Google Scholar

• 9.

  PavlicAHarshbargerCLRosenthalerLSnedekerJGDualJ. Sharp-edge-based acoustofluidic chip capable of programmable pumping, mixing, cell focusing and trapping. Phys Fluids (2023) 35:022006. 10.1063/5.0133992

  • CrossRef
  • Google Scholar

• 10.

  NiveditaNGargNLeeAPPapautskyI. A high throughput microfluidic platform for size-selective enrichment of cell populations in tissue and blood samples. Analyst (2017) 142:2558–69. 10.1039/c7an00290d

  • CrossRef
  • Google Scholar

• 11.

  GargNWesterhofTMLiuVLiuRNelsonELLeeAP. Whole-blood sorting, enrichment and in situ immunolabeling of cellular subsets using acoustic microstreaming. Microsystems and Nanoengineering (2018) 4:17085. 10.1038/micronano.2017.85

  • CrossRef
  • Google Scholar

• 12.

  HuangPHXieYAhmedDRufoJNamaNChenYet alAn acoustofluidic micromixer based on oscillating sidewall sharp-edges. Lab Chip (2013) 13:3847–52. 10.1039/c3lc50568e

  • CrossRef
  • Google Scholar

• 13.

  NamaNHuangPHHuangTJCostanzoF. Investigation of micromixing by acoustically oscillated sharp-edges. Biomicrofluidics (2016) 10:024124. 10.1063/1.4946875

  • CrossRef
  • Google Scholar

• 14.

  ZhangCGuoXBrunetPCostalongaMRoyonL. Acoustic streaming near a sharp structure and its mixing performance characterization. Microfluidics and Nanofluidics (2019) 23:104. 10.1007/s10404-019-2271-5

  • CrossRef
  • Google Scholar

• 15.

  MullerPBBarnkobRJensenMJHBruusH. A numerical study of microparticle acoustophoresis driven by acoustic radiation forces and streaming-induced drag forces. Lab Chip (2012) 12:4617–27. 10.1039/c2lc40612h

  • CrossRef
  • Google Scholar

• 16.

  RayleighL. On the circulation of air observed in Kundt’s tubes, and on some allied acoustical problems. Phil Trans R Soc Lond (1884) 175:1–21. 10.1098/rstl.1884.0002

  • CrossRef
  • Google Scholar

• 17.

  NyborgWL. Acoustic streaming due to attenuated plane waves. The J Acoust Soc America (1953) 25:68–75. 10.1121/1.1907010

  • CrossRef
  • Google Scholar

• 18.

  NyborgWL. Acoustic streaming near a boundary. J Acoust Soc America (1958) 30:329–39. 10.1121/1.1909587

  • CrossRef
  • Google Scholar

• 19.

  HahnPLeibacherIBaaschTDualJ. Numerical simulation of acoustofluidic manipulation by radiation forces and acoustic streaming for complex particles. Lab Chip (2015) 15:4302–13. 10.1039/c5lc00866b

  • CrossRef
  • Google Scholar

• 20.

  LeiJGlynne-JonesPHillM. Acoustic streaming in the transducer plane in ultrasonic particle manipulation devices. Lab Chip (2013) 13:2133–43. 10.1039/c3lc00010a

  • CrossRef
  • Google Scholar

• 21.

  LeiJHillMGlynne-JonesP. Numerical simulation of 3d boundary-driven acoustic streaming in microfluidic devices. Lab Chip (2014) 14:532–41. 10.1039/c3lc50985k

  • CrossRef
  • Google Scholar

• 22.

  BachJSBruusH. Theory of pressure acoustics with viscous boundary layers and streaming in curved elastic cavities. J Acoust Soc America (2018) 144:766–84. 10.1121/1.5049579

  • CrossRef
  • Google Scholar

• 23.

  JoergensenJHBruusH. Theory of pressure acoustics with thermoviscous boundary layers and streaming in elastic cavities. J Acoust Soc America (2021) 149:3599–610. 10.1121/10.0005005

  • CrossRef
  • Google Scholar

• 24.

  SkovNRBachJSWinckelmannBGBruusH. 3d modeling of acoustofluidics in a liquid-filled cavity including streaming, viscous boundary layers, surrounding solids, and a piezoelectric transducer. Aims Math (2019) 4:99–111. 10.3934/math.2019.1.99

  • CrossRef
  • Google Scholar

• 25.

  LeiJHillMde León AlbarránCPGlynne-JonesP. Effects of micron scale surface profiles on acoustic streaming. Microfluidics and Nanofluidics (2018) 22:140–14. 10.1007/s10404-018-2161-2

  • CrossRef
  • Google Scholar

• 26.

  ZhouY. Effect of microchannel protrusion on the bulk acoustic wave-induced acoustofluidics: Numerical investigation. Biomed Microdevices (2022) 24:7–12. 10.1007/s10544-021-00608-6

  • CrossRef
  • Google Scholar

• 27.

  PerraEHaywardNPritzkerKPNieminenHJ. An ultrasonically actuated needle promotes the transport of nanoparticles and fluids. J Acoust Soc America (2022) 152:251–65. 10.1121/10.0012190

  • CrossRef
  • Google Scholar

• 28.

  PavlicANagpurePErmanniLDualJ. Influence of particle shape and material on the acoustic radiation force and microstreaming in a standing wave. Phys Rev E (2022) 106:015105. 10.1103/physreve.106.015105

  • CrossRef
  • Google Scholar

• 29.

  OroscoJFriendJ. Modeling fast acoustic streaming: Steady-state and transient flow solutions. Phys Rev E (2022) 106:045101. 10.1103/physreve.106.045101

  • CrossRef
  • Google Scholar

• 30.

  HamiltonMFBlackstockDT. Nonlinear acoustics. San Diego: Academic press San Diego (1998).

  • Google Scholar

• 31.

  LighthillJ. Acoustic streaming. J Sound Vibration (1978) 61:391–418. 10.1016/0022-460x(78)90388-7

  • CrossRef
  • Google Scholar

• 32.

  BruusH. Theoretical Microfluidics, 18. Oxford University Press (2007).

  • Google Scholar

• 33.

  NyborgWLM. Acoustic streaming. In: Physical acoustics, 2. Elsevier (1965). p. 265–331.

  • Google Scholar

• 34.

  DoinikovAAThibaultPMarmottantP. Acoustic streaming in a microfluidic channel with a reflector: Case of a standing wave generated by two counterpropagating leaky surface waves. Phys Rev E (2017) 96:013101. 10.1103/physreve.96.013101

  • CrossRef
  • Google Scholar

• 35.

  DoinikovAAThibaultPMarmottantP. Acoustic streaming induced by two orthogonal ultrasound standing waves in a microfluidic channel. Ultrasonics (2018) 87:7–19. 10.1016/j.ultras.2018.02.002

  • CrossRef
  • Google Scholar

• 36.

  BaaschTDoinikovAADualJ. Acoustic streaming outside and inside a fluid particle undergoing monopole and dipole oscillations. Phys Rev E (2020) 101:013108. 10.1103/physreve.101.013108

  • CrossRef
  • Google Scholar

• 37.

  AndrewsDGMcIntyreM. An exact theory of nonlinear waves on a Lagrangian-mean flow. J Fluid Mech (1978) 89:609–46. 10.1017/s0022112078002773

  • CrossRef
  • Google Scholar

• 38.

  VannesteJBühlerO. Streaming by leaky surface acoustic waves. Proc R Soc A: Math Phys Eng Sci (2011) 467:1779–800. 10.1098/rspa.2010.0457

  • CrossRef
  • Google Scholar

• 39.

  HamiltonMFIlinskiiYAZabolotskayaEA. Acoustic streaming generated by standing waves in two-dimensional channels of arbitrary width. J Acoust Soc America (2003) 113:153–60. 10.1121/1.1528928

  • CrossRef
  • Google Scholar

• 40.

  PavlicADualJ. On the streaming in a microfluidic Kundt’s tube. J Fluid Mech (2021) 911:A28. 10.1017/jfm.2020.1046

  • CrossRef
  • Google Scholar

• 41.

  ComsolAB. COMSOL Multiphysics 5.6 (2020). Available at: www.comsol.com (March 08, 2023).

  • Google Scholar

• 42.

  BaaschTPavlicADualJ. Acoustic radiation force acting on a heavy particle in a standing wave can be dominated by the acoustic microstreaming. Phys Rev E (2019) 100:061102. 10.1103/physreve.100.061102

  • CrossRef
  • Google Scholar

• 43.

  RiaudABaudoinMMatarOBThomasJLBrunetP. On the influence of viscosity and caustics on acoustic streaming in sessile droplets: An experimental and a numerical study with a cost-effective method. J Fluid Mech (2017) 821:384–420. 10.1017/jfm.2017.178

  • CrossRef
  • Google Scholar

• 44.

  BachJSBruusH. Bulk-driven acoustic streaming at resonance in closed microcavities. Phys Rev E (2019) 100:023104. 10.1103/physreve.100.023104

  • CrossRef
  • Google Scholar

• 45.

  HolmesMParkerNPoveyM. Temperature dependence of bulk viscosity in water using acoustic spectroscopy. J Phys Conf Ser (2011) 269:012011. 10.1088/1742-6596/269/1/012011

  • CrossRef
  • Google Scholar

• 46.

  MullerPBBruusH. Theoretical study of time-dependent, ultrasound-induced acoustic streaming in microchannels. Phys Rev E (2015) 92:063018. 10.1103/physreve.92.063018

  • CrossRef
  • Google Scholar

• 47.

  BruusH. Acoustofluidics 2: Perturbation theory and ultrasound resonance modes. Lab Chip (2012) 12:20–8. 10.1039/c1lc20770a

  • CrossRef
  • Google Scholar

• 48.

  KaynakMOzcelikANamaNNourhaniALammertPECrespiVHet alAcoustofluidic actuation of in situ fabricated microrotors. Lab Chip (2016) 16:3532–7. 10.1039/c6lc00443a

  • CrossRef
  • Google Scholar

• 49.

  KaynakMOzcelikANourhaniALammertPECrespiVHHuangTJ. Acoustic actuation of bioinspired microswimmers. Lab Chip (2017) 17:395–400. 10.1039/c6lc01272h

  • CrossRef
  • Google Scholar

• 50.

  DillingerCNamaNAhmedD. Ultrasound-activated ciliary bands for microrobotic systems inspired by starfish. Nat Commun (2021) 12:6455. 10.1038/s41467-021-26607-y

  • CrossRef
  • Google Scholar

• 51.

  ChenZLiuPZhaoXHuangLXiaoYZhangYet alSharp-edge acoustic microfluidics: Principles, structures, and applications. Appl Mater Today (2021) 25:101239. 10.1016/j.apmt.2021.101239

  • CrossRef
  • Google Scholar