Numerical study of thermoviscous effects in ultrasound-induced acoustic streaming in microchannels

We present a numerical study of thermoviscous effects on the acoustic streaming flow generated by an ultrasound standing-wave resonance in a long straight microfluidic channel containing a Newtonian fluid. These effects enter primarily through the temperature and density dependence of the fluid viscosity. The resulting magnitude of the streaming flow is calculated and characterized numerically, and we find that even for thin acoustic boundary layers, the channel height affects the magnitude of the streaming flow. For the special case of a sufficiently large channel height, we have successfully validated our numerics with analytical results from 2011 by Rednikov and Sadhal for a single planar wall. We analyzed the time-averaged energy transport in the system and the time-averaged second-order temperature perturbation of the fluid. Finally, we have made three main changes in our previously published numerical scheme to improve the numerical performance: (i) The time-averaged products of first-order variables in the time-averaged second-order equations have been recast as flux densities instead of as body forces. (ii) The order of the finite-element basis functions has been increased in an optimal manner. (iii) Based on the International Association for the Properties of Water and Steam (IAPWS 1995, 2008, and 2011), we provide accurate polynomial fits in temperature for all relevant thermodynamic and transport parameters of water in the temperature range from 10 to $50 {}^{\circ }\text{C}$.

Figure 1

(a) Sketch of a physical acoustophoresis setup with a silicon microfluidic chip on top of a piezoactuator such as in Ref. [30]. (b) Sketch of the model system for the numerical scheme with the viscous fluid domain surrounded by hard walls. The thick arrows indicate in-phase oscillating displacement of the left and right walls. By default, we set the width $w=380\mu \text{m}$ and the height $h=160\mu \text{m}$ as in Ref. [30].

Figure 2

(a) The used triangular mesh with a gradually increasing element size from $0.5 \mu \text{m}$ at the boundaries to $20 \mu \text{m}$ in the bulk. This mesh, chosen as the default mesh, contains 30 246 elements. (b) Rectangular mesh with thin elongated 0.1-$\mu \text{m}$-by-10-$\mu \text{m}$ elements at the boundaries and gradually changing to nearly square 10-$\mu \text{m}$-by-10-$\mu \text{m}$ elements in the bulk. This mesh contains 3308 elements. (c) The convergence parameter $C$ of Eq. (20) for all first- and second-order fields vs the numerical resolution defined by ${\delta }_{\mathrm{s}}^{{}_{}}/d_{\mathrm{bd}}$, where $d_{\mathrm{bd}}$ is the mesh-element size at the boundary. The fields are solved on triangular meshes with different boundary element sizes but all with fixed bulk element size $d_{\mathrm{bk}}=20\mu \text{m}$ and growth rate $\alpha =1.3$, while the reference solution is calculated for $d_{\mathrm{bd}}=0.15\mu \text{m},d_{\mathrm{bk}}=2\mu \text{m}$, and $\alpha =1.3$. The vertical dotted line indicates the solution for $d_{\mathrm{bd}}=0.5\mu \text{m}$, which is chosen as the default value for the following simulations.

Figure 3

Graph of the acoustic energy density $E_{\mathrm{ac}}^{}$ Eq. (16) as a function of the frequency $f$ of the oscillating boundary condition. $f_{\mathrm{res}}^{}$ is the resonance frequency at the center of the peak, while $f_0^{}$ is the ideal frequency corresponding to matching a half-wavelength with the channel width. The inset shows the magnitude of the resonant oscillating first-order velocity field $v_{y1}^{\mathrm{a}}$ relative to the amplitude of the oscillating velocity boundary condition $v_{\mathrm{bc}}^{{}_{}}$ as a function of the actuation frequency $f$.

Figure 4

(a) Time-averaged second-order fluid velocity field ${\mathit{v}}_2^{}$ (vectors) and its magnitude [color plot ranging from 0 mm/s (black) to 0.13 mm/s (white)] in the vertical channel cross section calculated at $T_0^{}=25{}^{\circ }\text{C}$ and $f_{\mathrm{res}}^{}=1.9669923$ MHz. (b) The horizontal velocity component $v_{y2}^{}$ plotted along $y=w/4$ indicated by the magenta line in (a). The velocity field has been calculated for the five actuation frequencies shown in the inset resonance curve, and normalized to the analytical Rayleigh streaming magnitude $v_{\mathrm{str}}^{\mathrm{R}}=(3/8){\left(v_{y1}^{\mathrm{a}}\right)}^2/c_s^{}$, which is calculated based on the corresponding first-order solutions. The symbols are plotted in selected points illustrating the five numerical solutions (black lines) that coincide. (c) Normalized streaming magnitude $v_{\mathrm{str}}^{}/v_{\mathrm{str}}^{\mathrm{R}}$ (symbols) vs equilibrium temperature $T_0^{}$, calculated for different channel heights $h$. The full curve shows the analytical single-wall result by Ref. [15]. (d) Normalized streaming magnitude $v_{\mathrm{str}}^{}/v_{\mathrm{str}}^{\mathrm{R}}$ (symbols) vs channel height $h$, calculated for different equilibrium temperatures $T_0^{}$. The full curves show the analytical single-wall result by Ref. [15], while the dashed lines show the results of a one-dimensional analytical model with a Poiseuille backflow; see Eq. (21).

Figure 5

Time-averaged second-order temperature $T_2^{}$ calculated for the default 380-$\mu \text{m}$-by-160-$\mu \text{m}$ geometry, actuation frequency $f_{\mathrm{res}}^{}=1.9669923$ MHz, and equilibrium temperature $T_0^{}=25{}^{\circ }\text{C}$. (a) Color plot (black 0 mK to white 0.17 mK) of $T_2^{}$ in the channel cross section. (b) and (c) Line plots of $T_2^{}$ along the horizontal and vertical dashed lines in (a), respectively.

Figure 6

Time-averaged heat current densities in the channel cross section. (a) Sketch of the heat currents (arrows) with indication of the responsible terms in the time-averaged energy conservation Eq. (15). (b) Heat current density (arrows) and its magnitude [color plot ranging from zero (light green) to $7\times {10}^2$ W ${\text{m}}^{-2}$ (dark red)] in the bulk of the channel. The strong currents inside the boundary layers are not shown. (c) Magnitude of the $y$ component of the heat current density [color plot ranging from zero (light green) to $5\times {10}^4$ W ${\text{m}}^{-2}$ (dark red)] inside the boundary layer at the bottom wall. (d) Magnitude of the $z$ component of the heat current density [color plot ranging from $-6\times {10}^2$ W ${\text{m}}^{-2}$ (dark blue) to zero (light green)] inside the boundary layer at the bottom wall.

Figure 7

Color plot of the time-averaged second-order temperature $T_2^{}$ [from zero (black) to maximum (white)] in three cases. (a) $T_2^{}$ calculated from the complete governing equations, identical to Fig. 5. (b) $T_2^{}$ calculated without bulk viscosity ${\eta }^{\mathrm{b}}=0$. (c) $T_2^{}$ calculated with zero viscous stress $\mathit{\tau }=\mathbf{0}$ in the bulk defined by $\left|y\right|<(w/2-4\mu \text{m})$ and $\left|z\right|<(h/2-4\mu \text{m})$.

Figure 8

Fifth-order polynomial fits (dashed light green lines) to IAPWS data [24, 25, 26] (full black lines) of ten thermodynamic and transport parameters of water in the temperature range from 283 K ($10{}^{\circ }\text{C}$) to 323 K ($50{}^{\circ }\text{C}$).