Incorporation of SWAT & WEAP models for analysis of water demand deficits in the Kala Oya River Basin in Sri Lanka: perspective for climate and land change

1.   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.

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 $\tilde{\varphi}$ of amplitude $\varphi_0 = 1$ V and a frequency around 2.2 MHz applied to its bottom surface.

Figure 1.  (a) Top-view photograph of the original transducer-silicon-glass device studied in 2007 by Hagsäter et al. ^[16]. (b) A cut-open 3D sketch of the device in the red-dashed area of panel (a) showing the Pz26 piezo-electric transducer (green), the silicon base (gray), the water-filled cavity (blue) in the top of the silicon base, and the Pyrex lid (orange).

DownLoad: Full-Size Img PowerPoint

Table 1.  The length, width, and height $L \times W \times 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_{\mathrm{pz}}\!\times\!W_{\mathrm{pz}}\!\times\! H_{\mathrm{pz}}$	$L_{\mathrm{si}}\!\times\!W_{\mathrm{si}}\!\times\! H_{\mathrm{si}}$	$L_{\mathrm{py}}\!\times\!W_{\mathrm{py}}\!\times\! H_{\mathrm{py}}$	$L_{\mathrm{ca}}\!\times\!W_{\mathrm{ca}}\!\times\! H_{\mathrm{ca}}$	$L_{\mathrm{c1}}\!\times\!W_{\mathrm{c1}}\!\times\! H_{\mathrm{c1}}$	$L_{\mathrm{c2}}\!\times\!W_{\mathrm{c2}}\!\times\! H_{\mathrm{c2}}$
$\!49 \times 15 \times 1.0\!$	$\!49 \times 15 \times 0.5\!$	$\!49 \times 15 \times 0.5\!$	$\!2.02 \times 2 \times 0.2\!$	$\!11.3 \times 0.4 \times 0.2\!$	$\!12.4 \times 0.4 \times 0.2\!$

---

 | Show Table

DownLoad: CSV

2.   Theoretical background

We summarize the coupled equations of motion for a system driven by a time-harmonic electric potential, $\tilde{\varphi} = \varphi_0\: \mathrm{e}^{- {{\rm{i}}\omega t}}$ applied to selected boundaries of a piezo-electric Pz26 ceramic. Here, tilde denotes a field with harmonic time dependency, $\omega$ is the angular frequency in the low MHz range, and "$\mathrm{i}$" is the imaginary unit. This harmonic boundary condition excites the time-harmonic fields: the electric potential $\tilde{\varphi}(\mathit{\boldsymbol{r}}, t)$ in the Pz26 ceramic, the displacement $\tilde{ \mathit{\boldsymbol{u}}}(\mathit{\boldsymbol{r}}, t)$ in the solids, and the acoustic pressure $\tilde{p}^{{}}_1(\mathit{\boldsymbol{r}}, t)$ in the water,

In our simulation, we first solve the linear equations of the amplitude fields $\varphi(\mathit{\boldsymbol{r}})$, $\mathit{\boldsymbol{u}}(\mathit{\boldsymbol{r}})$, and $p^{{}}_1(\mathit{\boldsymbol{r}})$. Then, based on time-averaged products (over one oscillation period) of these fields, we compute the nonlinear acoustic radiation force $\mathit{\boldsymbol{F}}^{\mathrm{rad}}$ and the steady-state acoustic streaming velocity $\mathit{\boldsymbol{v}}_2(\mathit{\boldsymbol{r}})$.

2.1.   Linear acoustics in the fluid

In the fluid (water) of density $\rho^{{}}_{\mathrm{fl}}$, sound speed $c^{{}}_{\mathrm{fl}}$, dynamic viscosity $\eta^{{}}_{\mathrm{fl}}$, and bulk viscosity $\eta^{{{\mathrm{b}}}}_{\mathrm{fl}}$, we model the acoustic pressure $p^{{}}_1$ as in Ref. ^[12],

Here, $\boldsymbol{v}^{{}}_1$ is the acoustic velocity which is proportional to the pressure gradient $\boldsymbol{\nabla} p^{{}}_1$, while $\Gamma^{{}}_\mathrm{fl} \ll 1$ is a weak absorption coefficient, and $\kappa^{{}}_{\mathrm{fl}} = (\rho^{{}}_{\mathrm{fl}} c^{{2}}_{\mathrm{fl}})^{-1}$ is the isentropic compressibility of the fluid, see Table 2 for parameter values. The time-averaged acoustic energy density $E_{{\mathrm{ac}}_{}}^{{\mathrm{fl}}}$ 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 $^\circ$C for isotropic Pyrex borosilicate glass ^[17], cubic-symmetric silicon ^[18], and water ^[6]. Note that $c^{{}}_{12} = c^{{}}_{11} - 2 c^{{}}_{44}$ for isotropic solids.

Parameter		Pyrex	Si	Unit	Parameter		Water	Unit
Mass density	${\rho _{{\rm{sl}}}}$	2230	2329	kg m$^{-3}$	Mass density	${\rho _{{\rm{fl}}}}$	997.05	kg m$^{-3}$
Elastic modulus	$c^{{}}_{11}$	69.72	165.7	GPa	Sound speed	${c _{{\rm{fl}}}}$	1496.7	m s$^{-1}$
Elastic modulus	$c^{{}}_{44}$	26.15	79.6	GPa	Dyn. viscosity	${\eta _{{\rm{fl}}}}$	2.485	mPa s
Elastic modulus	$c^{{}}_{12}$	17.43	63.9	GPa	Bulk viscosity	$\eta _{{\rm{fl}}}^{\rm{b}}$	0.890	mPa s
Damping coeff.	${\Gamma _{{\rm{sl}}}}$	0.0004	0.0000	1	Damping coeff.	${\Gamma _{{\rm{fl}}}}$	0.00002	1
–	–	–	–	–	Compressibility	${\kappa _{{\rm{fl}}}}$	452	TPa$^{-1}$

---

 | Show Table

DownLoad: CSV

2.2.   Linear elastic motion of the solids

In the solid materials, each with a given density $\rho^{{}}_{\mathrm{sl}}$, we model the displacement field $\mathit{\boldsymbol{u}}$ using the equation of motion given by ^[12]

where $\Gamma^{{}}_\mathrm{sl}\ll1$ is a weak damping coefficient. Here, $\mathit{\boldsymbol{\sigma}}$ is the stress tensor, which is coupled to $\mathit{\boldsymbol{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 "$\mathrm{Re}$" denotes the real value and "*" the complex conjugate of a complex number, respectively.

2.3.   Stress-strain coupling in elastic solids

For a crystal with either cubic or isotropic symmetry, the relation between the stress tensor $\sigma^{{}}_{ij}$ and strain components $\frac12(\partial^{{}}_i u^{{}}_j+ \partial^{{}}_j u^{{}}_i)$ is given in the compact Voigt representation as ^[19]

Here, $c^{{}}_{ij}$ are the elastic moduli which are listed for Pyrex and silicon in Table 2.

2.4.   Stress-strain coupling in piezoelectric ceramics

Lead-zirconate-titanate (PZT) ceramics are piezoelectric below their Curie temperature, which typically is $200 - 400 ^\circ \rm{C}$. Using Cartesian coordinates and the Voigt notation for a PZT ceramic, the mechanical stress tensor $\sigma^{{}}_{ij}$ and electric displacement field $D^{{}}_i$ are coupled to the mechanical strain components $\frac12(\partial^{{}}_i u^{{}}_j+ \partial^{{}}_j u^{{}}_i)$ and the electrical potential $\varphi$ 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 $\varphi$ in the transducer, and since we assume no free charges here and only low-MHz frequencies, $\varphi$ 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} = \frac12(c^{{}}_{11} - c^{{}}_{12})$. The damping coefficient is $\Gamma^{{}}_\mathrm{sl} = 0.02$ ^[13].

Parameter	Value	Parameter	Value	Parameter	Value
${\rho _{{\rm{sl}}}}$	7700 kg/m$^3$	${\varepsilon_{11}}$	828 ${\varepsilon _0}$	${\varepsilon _{33}}$	700 ${\varepsilon _0}$
$c^{{}}_{11}$	168 GPa	$c^{{}}_{33}$	123 GPa	$e^{{}}_{31}$	$-$2.8 C/m$^2$
$c^{{}}_{12}$	110 GPa	$c^{{}}_{44}$	30.1 GPa	$e^{{}}_{33}$	14.7 C/m$^2$
$c^{{}}_{13}$	99.9	GPa $c^{{}}_{66}$	29.0 GPa	$e^{{}}_{15}$	9.86 C/m$^2$

---

 | Show Table

DownLoad: CSV

2.5.   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 \leftarrow B$) on domain A from domain B with the surface normal $\mathit{\boldsymbol{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 $\mathit{\boldsymbol{\sigma}}\cdot \mathit{\boldsymbol{n}}$ of Eq. (2.4) and the acoustic velocity $\boldsymbol{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 $\delta = \sqrt{2 \eta^{{}}_{\mathrm{fl}}/(\rho^{{}}_{\mathrm{fl}}\omega)}$ ($\approx 0.35 \;{\rm{ \mathsf{ μ} }m}$ at 2.3 MHz) has been taken into account analytically. As a result, terms appear in (2.9d) and (2.9e) that involve the shear-wave number $k^{{}}_\mathrm{s} = (1+ \mathrm{i})\:\delta^{-1}$ as well as the tangential divergence of the tangential component of the difference between the solid-wall velocity $\boldsymbol{v}^{{}}_{\mathrm{sl}} = - \mathrm{i} \omega \mathit{\boldsymbol{u}}$ and the acoustic velocity $\boldsymbol{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 $\Gamma_{\mathrm{fl}}^{\mathrm{eff}} \approx \frac{\delta}{H} \approx 0.002$, the ratio of the boundary layer width $\delta$ to the device height $H$ ^[6,13,15]. Remarkably, this boundary-layer dissipation dominates dissipation in the fluid domain, because $\Gamma^{{}}_\mathrm{fl} \ll \Gamma_{\mathrm{fl}}^{\mathrm{eff}} \ll 1$.

2.6.   The acoustic streaming

The acoustic streaming is the time-averaged (over one oscillation period), steady fluid velocity $\boldsymbol{v}^{{}}_2$ that is induced by the acosutic fields. In our recent analysis ^[15], we have shown that the governing equation of $\boldsymbol{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 $\Gamma^{{}}_\mathrm{fl}$. Further, at fluid-solid interfaces, the slip velocity $\boldsymbol{v}^{{{\mathrm{bc}}}}_2$ takes into account both the motion of the surrounding elastic solid and the Reynolds stress induced in viscous boundary layer in the fluid,

Here, we have used a special case of the slip velocity $\mathit{\boldsymbol{v}}^{\mathrm{bc}}_{2}$, which is only valid near acoustic resonance, where the magnitude $| \boldsymbol{v}^{{}}_1|$ of the acoustic velocity in the bulk is much larger than $\omega\:| \mathit{\boldsymbol{u}}^{\mathrm{bc}}_{\mathrm{sl}}|$ of the walls.

2.7.   The acoustic radiation force and streaming drag force on suspended microparticles

The response of primary interest in acoustofluidic applications, is the acoustic radiation force $\mathit{\boldsymbol{F}}^\mathrm{rad}$ and the Stokes drag from the acoustic streaming $\boldsymbol{v}^{{}}_2$ acting on suspended microparticles. In this work, we consider 1- and 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter spherical polystyrene "Styron 666" (ps) particles with density $\rho^{{}}_\mathrm{ps}$ and compressibility $\kappa_\mathrm{ps}$. For such large microparticle suspended in water of density $\rho^{{}}_{\mathrm{fl}}$ and compressibility $\kappa^{{}}_{\mathrm{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 $\mathit{\boldsymbol{v}}_1$, a single suspended microparticle of radius $a$, experience an acoustic radiation force $\mathit{\boldsymbol{F}}^\mathrm{rad}$, which, since $f_0$ and $f_1$ are real, is given by the potential $U^{\mathrm{rad}_{}}$ ^[22],

The microparticle is also influenced by a Stokes drag force $\mathit{\boldsymbol{F}}^\mathrm{drag} = 6\pi \eta^{{}}_{\mathrm{fl}} a \big(\boldsymbol{v}^{{}}_2- \boldsymbol{v}^{{}}_{\mathrm{ps}}\big)$, where $\boldsymbol{v}^{{}}_2$ and $\boldsymbol{v}^{{}}_{\mathrm{ps}}$ is the streaming velocity and the polystyrene particle velocity at the particle position $\mathit{\boldsymbol{r}}_{\mathrm{ps}}(t)$, respectively. In the experiments, the streaming and particle velocities are smaller than $v_0 = 1 {\mathrm{mm}}/{\mathrm{s}}$, which for a 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter particle corresponds to a small particle-Reynolds number $\frac{1}{ \rho^{{}}_{\mathrm{fl}}} \eta^{{}}_{\mathrm{fl}} a v_0 = 0.6$. Consequently, we can ignore the inertial effects and express the particle velocity for a particle at position $\mathit{\boldsymbol{r}}$ from the force balance $\mathit{\boldsymbol{F}}^\mathrm{rad} + \mathit{\boldsymbol{F}}^\mathrm{drag} = \boldsymbol{0}$, between the acoustic radiation force and streaming drag force,

The particle trajectory $\mathit{\boldsymbol{r}}_{\mathrm{ps}}(t)$ is then determined by straightforward time integration of $\frac{ \mathrm{d}}{ \mathrm{d} t} \mathit{\boldsymbol{r}}_{\mathrm{ps}} = \boldsymbol{v}^{{}}_{\mathrm{ps}}(\mathit{\boldsymbol{r}}_{\mathrm{ps}})$.

2.8.   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 $\mathit{\boldsymbol{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 $\varphi$ 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. The stress condition Eq. (2.9d) is given by the weak contribution

where $\mathit{\boldsymbol{n}} = (\mathtt{nX}, \mathtt{nY}, \mathtt{nZ})$ is the normal vector away from the solid domain, and $\mathtt{test(uX)}$ is the finite-element test function corresponding to the $x$-component $u_x$ of the solid displacement field $\mathit{\boldsymbol{u}}$, and similar for $y$ and $z$. The velocity condition Eq. (2.12) is given by the weak contribution

where $\mathit{\boldsymbol{n}} = (\mathtt{nX}, \mathtt{nY}, \mathtt{nZ})$ now is the normal vector away from the fluid, $\mathtt{test(p1)}$ is the test function for $p_1$, and $\mathtt{dtang}$ is the tangent-plane derivative operator available in COMSOL, see Ref. ^[15].

In a second step, we implement Eq. (2.10) for the acoustic streaming $\mathit{\boldsymbol{v}}_2$ in weak form. Specifically, the effective-model slip velocity condition are implemented as a "Dirichlet Boundary Condition" as follows. We use the outward normal vector $(\mathtt{nX}, \mathtt{nY}, \mathtt{nZ})$ as before and also the two perpendicular tangent vectors $(\mathtt{t1X}, \mathtt{t1Y}, \mathtt{t1Z})$ and $(\mathtt{t2X}, \mathtt{t2Y}, \mathtt{t2Z})$, and write the $x$-component $\mathtt{v2bcX}$ of $\mathit{\boldsymbol{v}}_2^{\mathrm{bc}}$ as,

and similarly for the $y$ and $z$ components. Here, $(\mathtt{AX}, \mathtt{AY}, \mathtt{AZ})$ is a vector defined in terms of the tangent-plane derivative $\boldsymbol{\nabla}_\parallel$ and the parallel velocity $\mathit{\boldsymbol{v}}_{1\parallel} = (\mathtt{v1parX}, \mathtt{v1parY}, \mathtt{v1parZ})$ with the $x$-component $\mathtt{v1parX} = (\mathit{\boldsymbol{v}}_1\cdot \mathit{\boldsymbol{t}}_1)\:t_{1x} + (\boldsymbol{v}^{{}}_1\cdot \mathit{\boldsymbol{t}}_2)\:t_{2x}$, as follows,

Finally, the acoustic radiation force $\mathit{\boldsymbol{F}}^\mathrm{rad}$ acting on the particles is calculated from Eq. (2.11) using the acoustic pressure $p^{{}}_1$ and velocity $\mathit{\boldsymbol{v}}_1$, and subsequently in a third step, following Ref. ^[23], we compute the particle trajectories $\mathit{\boldsymbol{r}}_{\mathrm{ps}}(t)$ from the time-integration of Eq. (2.12).

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 wave length in all domains, which for the second-order test function we use, corresponds to maximum mesh sizes of 0.52 mm, 0.59 mm, 0.50 mm, and 0.22 mm in the domains of Pz26, silicon, Pyrex, and water, respectively. The final implementation of the model contains 1.1 and 0.4 million degrees of freedom for the first- and second-order fields, repsectively. On our workstation, a Dell Inc Precision T7500 Intel Xeon CPU X5690 at 3.47 GHz with 128 GB RAM and 2 CPU cores, the model requires 45 GB RAM and takes 18 min per frequency. When running frequency sweeps of up to 70 frequency values, we used the DTU high-performance computer cluster requiring 464 GB RAM and 11 min per frequency.

3.   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.

Figure 2.  (a1) Micro-PIV measurements adapted from Ref. ^[16] of the particle velocity $\boldsymbol{v}^{{}}_{\mathrm{ps}}$ after 1 ms (yellow arrows, maximum $200 \;{\rm{ \mathsf{ μ} }m}/ \rm{s}$) superimposed on a micrograph of the final positions (black curved bands) of 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter polystyrene particles in water with a standing ultrasound wave at 2.17 MHz. (a2) Same as panel (a1), but for 1-$\;{\rm{ \mathsf{ μ} }m}$-diameter polystyrene particles moving in a 6-by-6 flow-roll pattern without specific final positions. (b1) Numerical 3D COMSOL modeling with actuation voltage $\varphi_0 = 1$ V of the acoustic potential $U^{\mathrm{rad}_{}}$ from $0$ fJ (black) to 7 fJ (orange) and the velocity (yellow arrows, maximum $170 \;{\rm{ \mathsf{ μ} }m}/ \rm{s}$) after 1 ms of 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter polystyrene particles in the horizontal center plane of the water-filled cavity at the resonance $f = 2.166 \rm{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 $\boldsymbol{v}^{{}}_{\mathrm{ps}}$ (magenta vectors) and its magnitude $v_{\mathrm{ps}}$ from 0 (black) to $200 \;{\rm{ \mathsf{ μ} }m}/ \rm{s}$ (white) of 1-$\;{\rm{ \mathsf{ μ} }m}$-diameter polystyrene particles.

DownLoad: Full-Size Img PowerPoint

In Figure 2(a1) we show the measured micro-particle image velocimetry (micro-PIV) results obtained on a large number of 5-$\;{\rm{ \mathsf{ μ} }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^{\mathrm{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-$\;{\rm{ \mathsf{ μ} }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 $\mathit{\boldsymbol{v}}_{\mathrm{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-$\;{\rm{ \mathsf{ μ} }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 $\varphi_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_{\mathrm{ps}} \approx 170 \;{\rm{ \mathsf{ μ} }m}/ \rm{s}$ for the large 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter, in agreement with the 200 $\;{\rm{ \mathsf{ μ} }m}/ \rm{s}$ reported in the experiment.

In Figure 3 we show another result that is in agreement with the experimental observations, namely the particle trajectories $\mathit{\boldsymbol{r}}_{\mathrm{ps}}(t)$ for suspensions of tracer particles of different size. The larger 5-$\;{\rm{ \mathsf{ μ} }m}$-diameter particles are focused along the bottom of the troughs in the acoustic potential $U^{\mathrm{rad}_{}}$, shown in Figure 2(b1), after a short time $\frac{1}{12}(2 \rm{mm})/(170 \;{\rm{ \mathsf{ μ} }m}/ \rm{s}) \approx 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-$\;{\rm{ \mathsf{ μ} }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).

Figure 3.  Numerical 3D COMSOL modeling of the trajectories $\mathit{\boldsymbol{r}}_{\mathrm{ps}}(t)$ (blue tracks) of 3600 polystyrene particles of radius $a$ corresponding to the cases shown in Figure 2(b1) and (b2). The particles start from $60\times 60$ regular quadratic grid points in the horizontal center plane of the cavity at $t = 0 \rm{s}$ when the ultrasound field is turned on, and their positions after 1.5 s are represented by red points. (a) $a = 2.5 \;{\rm{ \mathsf{ μ} }m}$ at $f = 2.166 \rm{MHz}$. (b) $a = 0.5 \;{\rm{ \mathsf{ μ} }m}$ at $f = 2.163 \rm{MHz}$ with the Eckart bulk force in Eq. (2.10a) increased by a factor 4. (c) Same as panel (b) but without the Eckart bulk force in Eq. (2.10a).

DownLoad: Full-Size Img PowerPoint

4.   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-$\;{\rm{ \mathsf{ μ} }m}$-diameter particles is caused by the so-called Eckart bulk force, the term in (2.10a) proportional to the acoustic energy flux density or intensity $\mathit{\boldsymbol{S}}_{\mathrm{ac}} = \frac12 \mathrm{Re}\big[p^*_1 \boldsymbol{v}^{{}}_1\big]$. 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_{\mathrm{ca}}$ and looks as the one in the midplane at $0.5\: H_{\mathrm{ca}}$ shown in Figs. 2(b2) and 3(b). Lei et al. ^[10] also pointed out that $\mathit{\boldsymbol{S}}_{\mathrm{ac}}$ could lead to the horizontal 6-by-6 flow-roll pattern in their 3D-fluid-domain model with hard-wall 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 $\mathit{\boldsymbol{S}}_{\mathrm{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_{\mathrm{ca}}$ and absent in the center plane at $0.5\:H_{\mathrm{ca}}$, the focal plane in the experimental studies. As our slip-velocity condition (2.10b) also contains $\mathit{\boldsymbol{S}}_{\mathrm{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^{\mathrm{rad}_{}}$ is fairly robust, the computation of the streaming velocity $\boldsymbol{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_{\mathrm{ca}}/W_{\mathrm{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_{\mathrm{ca}}/W_{\mathrm{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.

5.   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.

Acknowledgments

H. Bruus was supported by the BioWings project funded by the European Union's Horizon 2020 Future and Emerging Technologies (FET) programme, grant No. 801267.

Conflict of interest

All authors declare no conflicts of interest in this paper.