3D modeling of acoustofluidics in a liquid-filled cavity including streaming, viscous boundary layers, surrounding solids, and a piezoelectric transducer

We present a full 3D numerical simulation of the acoustic streaming observed in full-image micro-particle velocimetry by Hagsäter et al., Lab Chip 7, 1336 (2007) in a 2 mm by 2 mm by 0.2 mm microcavity embedded in a 49 mm by 15 mm by 2 mm chip excited by 2-MHz ultrasound. The model takes into account the piezo-electric transducer, the silicon base with the water-filled cavity, the viscous boundary layers in the water, and the Pyrex lid. The model predicts well the experimental results.


Introduction and definition of the model system
For the past 15 years, ultrasound-based microscale acoustofluidic devices have successfully and in increasing numbers been used in the fields of biology, environmental and forensic sciences, and clinical diagnostics [1][2][3][4][5]. However, it remains a challenge to model and optimize a given device including all relevant acoustofluidic aspects. Steadily, good progress is being made towards this goal. Examples of recent advances in modeling include work in two dimensions (2D) by Muller and Bruus [6,7] on thermoviscous and transient effects of acoustic pressure, radiation force, and streaming in the fluid domain, and work by Nama et al. [8] on acoustophoresis induced by a given surface acoustic wave in a fluid domain capped by a PDMS lid. Examples of 3D modeling include work by Lei et al. [9,10] on boundary-layer induced streaming in fluid domains with hard wall and outgoing plane-wave boundary conditions, work by Gralinski et al. [11] on the acoustic pressure fields in circular capillaries including the fluid and glass domains and excited by a given wall vibration, a model later extended by Ley and Bruus [12] to take into account absorption and outgoing waves, and work by Hahn and Dual [13] on the acoustic pressure and acoustic radiation force in the fluid domain including the surrounding transducer, silicon and glass domains, as well as bulk, boundary-layer, and thermal dissipation.  Table 1. The length, width, and height L × W × H (in mm) of the six rectangular elements in the acoustofluidic device model of Figure 1(b): The piezoelectric transducer (pz), the silicon base (si), the Pyrex lid (py), the main cavity (ca), and the two inlet channels (c1) and (c2).

Pz26
Silicon Pyrex Cavity Channel 1 Channel 2 L pz ×W pz ×H pz L si ×W si ×H si L py ×W py ×H py L ca ×W ca ×H ca L c1 ×W c1 ×H c1 L c2 ×W c2 ×H c2 49 × 15 × 1.0 49 × 15 × 0.5 49 × 15 × 0.5 2.02 × 2 × 0.2 11.3 × 0.4 × 0.2 12.4 × 0.4 × 0.2 In this paper, we present a 3D model and its implementation in the commercial software COMSOL Multiphysics [14] of a prototypical acoustofluidic silicon-glass-based device that takes into account the following physical aspects: the piezo-electric transducer driving the system, the silicon base that contains the acoustic cavity, the fluid with bulk-and boundary-layer-driven streaming, the Pyrex lid, and a dilute microparticle suspension filling the cavity. This work represents a synthesis of our previous modeling of streaming in 2D [6], acoustic fields in 3D [12], and boundary-layer analysis [15] enabling effective-model computation of streaming in 3D, and it combines and extends the 3D streaming study in the fluid domain by Lei et al. [10] and the 3D study of acoustics in the coupled transducer-sold-fluid system by Hahn and Dual [13]. To test the presented coupled 3D model, we have, as Lei et al. [10], chosen to model the system studied experimentally by Hagsäter et al. in 2007 [16] and shown in Figure 1. It consists of a rectangular 0.5-mm high silicon base, into the surface of which is etched a shallow square-shaped cavity with two inlet channels attached. The cavity is sealed with a 0.5-mm high Pyrex lid that exactly covers the silicon base. At the bottom of the silicon base is attached a 1-mm high rectangular Pz26 piezo-electric transducer. All three solid layers are 49 mm long and 15 mm wide. The nearly-square cavity is 2.02 mm long and 2 mm wide and has attached two inlet channels both 0.4 mm wide, but of unequal lengths 11.3 mm and 12.4 mm, respectively. The channels and cavity are 0.2 mm deep. A sketch of the model device is shown in Figure 1, and its geometrical parameters are summarized in Table 1. The transducer is grounded at the top and driven by an ac voltageφ of amplitude ϕ 0 = 1 V and a frequency around 2.2 MHz applied to its bottom surface.

Theoretical background
We summarize the coupled equations of motion for a system driven by a time-harmonic electric potential,φ = ϕ 0 e −iωt applied to selected boundaries of a piezo-electric Pz26 ceramic. Here, tilde denotes a field with harmonic time dependency, ω is the angular frequency in the low MHz range, and "i" is the imaginary unit. This harmonic boundary condition excites the time-harmonic fields: the electric potentialφ(r, t) in the Pz26 ceramic, the displacementũ(r, t) in the solids, and the acoustic pressurep 1 (r, t) in the water, (2.1) In our simulation, we first solve the linear equations of the amplitude fields ϕ(r), u(r), and p 1 (r). Then, based on time-averaged products (over one oscillation period) of these fields, we compute the nonlinear acoustic radiation force F rad and the steady-state acoustic streaming velocity v 2 (r).

Linear acoustics in the fluid
In the fluid (water) of density ρ fl , sound speed c fl , dynamic viscosity η fl , and bulk viscosity η b fl , we model the acoustic pressure p 1 as in Ref. [12], Here, v 1 is the acoustic velocity which is proportional to the pressure gradient ∇p 1 , while Γ fl 1 is a weak absorption coefficient, and κ fl = (ρ fl c 2 fl ) −1 is the isentropic compressibility of the fluid, see Table 2 for parameter values. The time-averaged acoustic energy density E fl ac in the fluid domain is the sum of the time-averaged (over one oscillation period) kinetic and compressional energy densities, Table 2. Material parameters at 25 • C for isotropic Pyrex borosilicate glass [17], cubicsymmetric silicon [18], and water [6]. Note that c 12 = c 11 − 2c 44 for isotropic solids.

Linear elastic motion of the solids
In the solid materials, each with a given density ρ sl , we model the displacement field u using the equation of motion given by [12] − ρ sl ω 2 (1 + iΓ sl ) u = ∇ · σ, (2.4) where Γ sl 1 is a weak damping coefficient. Here, σ is the stress tensor, which is coupled to u through a stress-strain relation depending on the material-dependent elastic moduli. The time-averaged acoustic energy density in the solids is given by the sum of kinetic and elastic contributions, where "Re" denotes the real value and "*" the complex conjugate of a complex number, respectively.

Stress-strain coupling in elastic solids
For a crystal with either cubic or isotropic symmetry, the relation between the stress tensor σ i j and strain components 1 2 (∂ i u j + ∂ j u i ) is given in the compact Voigt representation as [19] , for Pyrex and silicon. (2.6) Here, c i j are the elastic moduli which are listed for Pyrex and silicon in Table 2.

Stress-strain coupling in piezoelectric ceramics
Lead-zirconate-titanate (PZT) ceramics are piezoelectric below their Curie temperature, which typically is 200 − 400 • C. Using Cartesian coordinates and the Voigt notation for a PZT ceramic, the mechanical stress tensor σ i j and electric displacement field D i are coupled to the mechanical strain components 1 2 (∂ i u j + ∂ j u i ) and the electrical potential ϕ through the relation [19] The values of the material parameters for the PZT ceramic Pz26 are listed in Table 3. Due to the high electric permittivity of Pz26, we only model the electric potential ϕ in the transducer, and since we assume no free charges here and only low-MHz frequencies, ϕ must satisfy the quasi-static equation, Table 3. Material parameters of Ferroperm Ceramic Pz26 from Meggitt A/S [20]. Isotropy in the x-y plane implies c 66 = 1 2 (c 11 − c 12 ). The damping coefficient is Γ sl = 0.02 [13].

Boundary conditions and boundary layers in the fluid at the fluid-solid interfaces
The applied boundary conditions are the usual ones, namely that (1) the stress and the velocity fields are continuous across all fluid-solid and solid-solid interfaces, (2) the stress is zero on all outer boundaries facing the air, (3) the piezoelectric ceramic is driven by a given electric potential at specified surfaces that represent the presence of infinitely thin, massless electrodes, and (4) there are no free charges on the surface of the ceramic. The influence (A ← B) on domain A from domain B with the surface normal n pointing away from A, is given by While the overall structure of these boundary conditions is the usual continuity in stress and velocity, the details of Eqs. (2.9d) and (2.9e) are not conventional. They are the boundary conditions for the surface stress σ · n of Eq. (2.4) and the acoustic velocity v 1 of Eq. (2.2) (proportional to the gradient of the acoustic pressure p 1 ) derived by Bach and Bruus using their recent effective pressure-acoustics theory [15]. In this theory, the viscous boundary layer of thickness δ = 2η fl /(ρ fl ω) (≈ 0.35 µm at 2.3 MHz) has been taken into account analytically. As a result, terms appear in Eqs. (2.9d) and (2.9e) that involve the shear-wave number k s = (1 + i) δ −1 as well as the tangential divergence of the tangential component of the difference between the solid-wall velocity v sl = −iωu and the acoustic velocity v 1 at the fluid-solid interface. This boundary condition also takes into account the large dissipation in the boundary layers, which leads to an effective damping coefficient Γ eff fl ≈ δ H ≈ 0.002, the ratio of the boundary layer width δ to the device height H [6,13,15]. Remarkably, this boundary-layer dissipation dominates dissipation in the fluid domain, because Γ fl Γ eff fl 1.

The acoustic streaming
The acoustic streaming is the time-averaged (over one oscillation period), steady fluid velocity v 2 that is induced by the acosutic fields. In our recent analysis [15], we have shown that the governing equation of v 2 corresponds to a steady-state, incompressible Stokes flow with a body force in the bulk due to the time-averaged acoustic dissipation proportional to Γ fl . Further, at fluid-solid interfaces, the slip velocity v bc 2 takes into account both the motion of the surrounding elastic solid and the Reynolds Here, we have used a special case of the slip velocity v bc 2 , which is only valid near acoustic resonance, where the magnitude |v 1 | of the acoustic velocity in the bulk is much larger than ω |u bc sl | of the walls.

The acoustic radiation force and streaming drag force on suspended microparticles
The response of primary interest in acoustofluidic applications, is the acoustic radiation force F rad and the Stokes drag from the acoustic streaming v 2 acting on suspended microparticles. In this work, we consider 1-and 5-µm-diameter spherical polystyrene "Styron 666" (ps) particles with density ρ ps and compressibility κ ps . For such large microparticle suspended in water of density ρ fl and compressibility κ fl , thermoviscous boundary layers can be neglected, and the monopole and dipole acoustic scattering coefficients f 0 and f 1 are real numbers given by [21], Given an acoustic pressure p 1 and velocity v 1 , a single suspended microparticle of radius a, experience an acoustic radiation force F rad , which, since f 0 and f 1 are real, is given by the potential U rad [22], The microparticle is also influenced by a Stokes drag force F drag = 6πη fl a v 2 −v ps , where v 2 and v ps is the streaming velocity and the polystyrene particle velocity at the particle position r ps (t), respectively. In the experiments, the streaming and particle velocities are smaller than v 0 = 1 mm/s, which for a 5-µm-diameter particle corresponds to a small particle-Reynolds number 1 ρ fl η fl av 0 = 0.6. Consequently, we can ignore the inertial effects and express the particle velocity for a particle at position r from the force balance F rad + F drag = 0, between the acoustic radiation force and streaming drag force, v ps (r) = v 2 (r) + 1 6πη fl a F rad (r). (2.12) The particle trajectory r ps (t) is then determined by straightforward time integration of d dt r ps = v ps (r ps ).

Numerical implementation
Following the procedure described in Ref. [12], including mesh convergence tests, the coupled field equations (2.2) and (2.4) for the fluid pressure p 1 and elastic-solid displacement u are implemented directly in the finite-element-method software Comsol Multiphysics 5.3a [14] using the weak form interface "PDE Weak Form". A COMSOL script with a PDE-weak-form implementation of acoustofluidics is available as supplemental material in Ref. [7]. Here, we extend the model of Ref. [12] by including the transducer with the piezoelectric stress-strain coupling Eq. (2.6) and implementing the governing equation (2.8) for the electric potential ϕ in weak form. Similarly, the boundary conditions Eq. (2.9) are implemented in weak form. Specifically, the effective-model boundary conditions are implemented as "Weak Contributions" as follows.
We optimize the mesh to obtain higher resolution in the water-filled cavity, where we need to calculate numerical derivatives of the resulting fields to compute the streaming and radiation forces, and less in the surrounding solids and in the transducer. We ensure having at least six nodal points per

Results for the transducer-glass-silicon acoustofluidic device
We apply the 3D model of Section 2 to the transducer-glass-silicon acoustofluidic device by Hagsäter et al. [16], shown in Figure 1 and using the parameter values listed in Tables 1, 2, and 3. In Figure 2 we compare the experimental results from Ref. [16] with our model simulations. (b1) Numerical 3D COMSOL modeling with actuation voltage ϕ 0 = 1 V of the acoustic potential U rad from 0 fJ (black) to 7 fJ (orange) and the velocity (yellow arrows, maximum 170 µm/s) after 1 ms of 5-µm-diameter polystyrene particles in the horizontal center plane of the water-filled cavity at the resonance f = 2.166 MHz. (b2) Numerical modeling at the same conditions as in panel (b1), but at the slightly lower frequency 2.163 MHz, of the particle velocity v ps (magenta vectors) and its magnitude v ps from 0 (black) to 200 µm/s (white) of 1-µm-diameter polystyrene particles.
In Figure 2(a1) we show the measured micro-particle image velocimetry (micro-PIV) results obtained on a large number of 5-µm-diameter tracer particles at an excitation frequency of 2.17 MHz. The yellow arrows indicate the velocity of the tracer particles 1 ms after the ultrasound has been turned on, and the black bands are the tracer particles focused at the minimum of the acoustic potential U rad after a couple of seconds of ultrasound actuation. A clear pattern of 3 wavelengths in each direction is observed. Similarly, in Figure 2(a2) is shown the micro-PIV results for the smaller 1-µm-diameter tracer particles. It is seen that these particles, in contrast to the larger particles, are not focused but keep moving in a 6-by-6 flow-roll pattern. This result from Ref. [16] is remarkable, as the conventional Rayleigh streaming pattern [6,7,23] has four streaming rolls per wavelength oriented in the vertical plane, but here is only seen two rolls per wavelength, and they are oriented in the horizontal plane.
In Figure 2(b1) and (b2) we see that our model predicts the observed acoustofluidics response qualitatively for both the larger and the smaller tracer particles at a resonance frequency slightly below 2.17 MHz. Even the uneven local amplitudes of the particle velocity v ps in the 6-by-6 flow-roll pattern, which shifts around as the frequency is changed a few kHz, is in accordance with the observations. In Ref. [16] it is mentioned that "If the frequency is shifted slightly in the vicinity of 2.17 MHz, the same vortex pattern will still be visible, but the strength distribution between the vortices will be altered.". We have chosen the 3-kHz lower frequency in Figure 2(b2) compared to (b1) to obtain a streaming pattern similar to the observed one for the small 5-µm-diameter particles.
Quantitatively, we find the following. The acoustic resonance is located at 2.166 MHz, only 0.2 % lower than the experimental value of 2.17 MHz. This good agreement should not be over emphasized, as we had to assume a certain length and width of the Pz26 transducer, because its actual size was not reported in Ref. [16]. Another source of error is that we have not modeled the coupling gel used in the experiment between the Pz26 transducer and the silicon base. The actual actuation voltage in the experiment has not been reported, so we have chosen ϕ 0 = 1 V, well within the range of the 20 V peak-to-peak function generator mentioned in Ref. [16], as it results in velocities v ps ≈ 170 µm/s for the large 5-µm-diameter, in agreement with the 200 µm/s reported in the experiment.
In Figure 3 we show another result that is in agreement with the experimental observations, namely the particle trajectories r ps (t) for suspensions of tracer particles of different size. The larger 5-µm-diameter particles are focused along the bottom of the troughs in the acoustic potential U rad , shown in Figure 2(b1), after a short time 1 12 (2 mm)/(170 µm/s) ≈ 1 s, forming the red wavy bands in Figure 3(a) very similar to the observed black bands in Figure 2(a1). In contrast, the smaller 1-µm-diameter particles are caught by the 6-by-6 streaming vortex pattern and swirl around without being focused, at least within the first 1.5 s as shown in Figure 3(b), in full agreement with the experimental observation shown in Figure 2(a2).

Discussion
Our full 3D numerical model, which takes into account the piezo-electric transducer, the silicon base with the water-filled cavity, the viscous boundary layers in the water, and the Pyrex lid, has been tested qualitatively and quantitatively by comparing the results for the acoustic radiation force, for the streaming velocity, and for the trajectories of tracer particles of two different sizes with the decade-old experimental results presented by Hagsäter et al. [16]. Remarkably, as predicted by Bach and Bruus [15], we find that the characteristic horizontal 6-by-6 flow-roll pattern of the small 1-µmdiameter particles is caused by the so-called Eckart bulk force, the term in (2.10a) proportional to the acoustic energy flux density or intensity S ac = 1 2 Re p * 1 v 1 . In our simulations this pattern occupies 80 % of the cavity volume stretching from 0.1 to 0.9 in units of the channel height H ca and looks as the one in the midplane at 0.5 H ca shown in Figs. 2(b2) and 3(b). Lei et al. [10] also pointed out that S ac could lead to the horizontal 6-by-6 flow-roll pattern in their 3D-fluid-domain model with hardwall and outgoing-plane-wave boundary conditions of the same device. In their model, the Eckart bulk force was neglected, and the horizontal-flow-roll producing term S ac appears only as part of their limiting-velocity boundary condition. As the remaining curl-free part of the boundary condition is dominating, they found the horizontal 6-by-6 flow-roll pattern to be confined to narrow regions around the two horizontal planes at 0.2 and 0.8 H ca and absent in the center plane at 0.5 H ca , the focal plane in the experimental studies. As our slip-velocity condition (2.10b) also contains S ac , see Eq. (62a) in Ref. [15], we do reproduce their findings, when we suppress the Eckart bulk force in Eq. (2.10b). This is illustrated in Figure 3(c), where we show that the flow-roll behavior is suppressed in the center plane and replaced by a clear divergent behavior. In agreement with Lei et al. [10], we find that although the determination of the first-order pressure p 1 and the acoustic potential U rad is fairly robust, the computation of the streaming velocity v 2 from the Stokes equation (2.10a) is sensitive to the exact value of the frequency and of the detailed shape of the fluid solid interface. In Ref. [24] we have shown in a simplified 3D-rectangular-fluid-domain model that the rotation of the acoustic intensity changes an order of magnitude when the aspect ratio L ca /W ca changes 1 %. In this study we have increased the Eckart bulk force in Eq. (2.10a) by a factor of 4 in order to make the rotating 6x6 pattern dominate clearly over the Rayleigh streaming in the center plane. This amplification may reflect that the chosen aspect ratio L ca /W ca = 1.01 was not
exactly the one realized in the experiment, an effect which should be studied further in experiments and simulations.
Our numerical study indicate that although the cavity in the Hagsäter device has a size of only three acoustic wavelengths, the existence of in-plane flow rolls may be controlled by the Eckart bulk force. This conclusion runs contrary to the conventional wisdom that the Eckart bulk force is only important in systems of a size, which greatly exceeds the acoustic wave length. This phenomenon deserves a much closer study in future work.
While our model takes many of the central aspects of acoustofluidics into account, it can still be improved. One possible improvement would be to include the influence of heating on the material parameters as in Ref. [6]. One big challenge in this respect is to determine the material parameters of the solids, which may be temperature and frequency dependent. Another difficult task is to model the coupling between the transducer and the chip, which in experiments typically are coupled using coupling gels or other ill-characterized adhesives. The last point we would like to raise is use of the simple Stokes drag law on the suspended particles in the cavity. Clearly, this model may be improved by including particle-wall effects and particle-particle interactions. However, as direct simulations of both of these effects are very memory consuming their implementation would require effective models.

Conclusion
We have described the implementation of a full 3D modeling of an acoustofluidic device taking into account the viscous boundary layers and acoustic streaming in the fluid, the vibrations of the solid material, and the piezoelectricity in the transducer. As such, our simulation is in many ways close to a realistic device, which is also reflected in the agreement between the simulation and the experiment shown in Figs. 2 and 3. Our model has correctly predicted the unusual streaming pattern observed in the device at the 2.17-Mz resonance: a horizontal 6-by-6 flow-roll pattern in 80 % of the cavity volume, a pattern much different form the conventional 12-by-2 Rayleigh streaming pattern in the vertical plane. Moreover, our model has revealed the surprising importance of the Eckart bulk force in an acoustic cavity with a size comparable to the acoustic wavelength. In future work, we must analyze the sensitivity of the streaming velocity and improve our understanding of the amplitude of the Eckart bulk force. By introducing the model, we have demonstrated that simulations can be used to obtain detailed information about the performance of an acoustofluidic device in 3D. Such simulations are likely to be useful for studies of the basic physics of acoustofluidics as well as for engineering purposes, such as improving existing microscale acoustofluidic devices. However, To fully exploit such modeling, more accurate determination is needed of the acoustic parameters of the actual transducers, elastic walls, and particle suspensions employed in a given experiment.