4. Discussion

In §§ 2 and 3, we have obtained the velocity field as well as the time-average stream function driven by a body force as was proposed by Hamilton et al. (Reference Hamilton, Ilinskii and Zabolotskaya2003a). The rheological Maxwell upper convective model for the fluid was assumed, which shows the viscoelastic effect through the Deborah number (Bird et al. Reference Bird, Armstrong and Hassager1987). To estimate the values of the dimensionless parameters involved in the analysis, we have considered the following resonator parameters values that have been reported in the specialized literature (Hamilton et al. Reference Hamilton, Ilinskii and Zabolotskaya2003a; Oyama et al. Reference Oyama, Imashiro, Kuriyama, Usui, Ando, Azuma, Morikawa, Kodeki, Takahara and Takemura2021): $H_0=20\delta _{\nu }=21\ \mathrm {\mu }{\rm m}$, $L=789\ \mathrm {\mu }{\rm m}$, $f=1$ MHz, $A=600$ nm and $\lambda _{a}=1.578$ mm. As a particular case, we consider values of the physical properties of blood (Hasgall et al. Reference Hasgall, Gennaro, Baumgartner, Neufeld, Lloyd, Gosselin, Payne, Klingenbock and Kuster2018), hence $\rho _{f}=1060\ {\rm kg}\ {\rm m}^{-3}$, $c_{0}=1578\ {\rm m}\ {\rm s}^{-1}$, $\mu =3.5\times 10^{-3}\ {\rm kg}\ {\rm m}^{-1}\ {\rm s}^{-1}$, $\delta _{\nu }=1.05\ \mathrm {\mu }{\rm m}$. With the above, the dimensionless parameters for the calculations assume the following values (unless otherwise specified for each figure): $\eta =7.6\times 10^{-4}$, $\epsilon =2.7\times 10^{-2}$ and $\alpha =20$.

In figure 3(a) we show the level curves of velocity $U(x=0,y)$, the real part of (2.60), projected onto the $y\text {--}De$ plane for $\alpha =20$. Here, the velocity is evaluated at the coordinate $x=0$, as suggested in the classical acoustic theory (Eckart Reference Eckart1948), since this is the only pressure node ($\hat {P}=0$), as indicated in (2.61). Figure 3 where computed at discrete values of $De$ from $y=-1$ to $y=1$. A fascinating outcome is shown in figures 3(a) and 3(b), where at certain values of $y$ and $De$ the velocity $U(x=0,y)$ acquires particularly large amplitudes. These points are known as resonance modes, useful for establishing acoustic handling of cells and particles in microfluidic systems (Bruus Reference Bruus2008). The resonant values of $De$ must satisfy $\partial U / \partial De =0$, corresponding with local maxima and/or minima of $U(x,y)$. From the figures mentioned above, the value of $\alpha =20$ makes the resonance modes located near the walls ($y=\pm 1$). For increasing values of $De$, they increase in quantity with the tendency to be located near $y=0$. We also observe that these resonance points arise near the zeros of the oscillatory horizontal velocity (marked with dashed white lines), where a change in sign of $U(0,y)$ can be seen. Therefore, to identify these zeros ($U(0,y)=0$) where resonance modes are in the vicinity, we note that (2.60) can be represented as a product of two complex functions $z_1$ and $z_2$, where $z_{j}=z(a_{j},b_{j})$, $a_{j}=\mathrm {Re}(z_{j})$, $b_{j}=\mathrm {Im}(z_{j})$ and $j=1,2$. We can rewrite (2.60) as $U(x,y)=z_{1}\boldsymbol {\cdot } z_{2}$, where

(4.1)\begin{equation} z_{1}= 1 - \frac{\cosh\left( \xi x\right) }{\cosh\left( \xi/2\right) } \end{equation}

and

(4.2)\begin{equation} z_{2}=\frac{\cosh(\zeta y)}{\cosh({\zeta })}-1. \end{equation}

It is of interest to rewrite $z_{1}$ and $z_{2}$ in the form $z_{j}=a_{j}+{\rm i}b_{j}$. Evaluating $z_{1}$ at $x=0$, yields

(4.3)\begin{equation} z_{1}= 1 - \frac{2\cosh\left( a_{\xi}/2\right) \cos\left( b_{\xi}/2\right)}{\cosh\left( a_{\xi}\right) +\cos\left( b_{\xi}\right) }+{\rm i}\frac{2\sinh\left( a_{\xi}/2\right) \sin\left( b_{\xi}/2\right)}{\cosh\left( a_{\xi}\right)+ \cos\left( b_{\xi}\right) } \end{equation}

and $z_{2}$ is given by

(4.4) \begin{align} z_{2}&=\frac{\cosh[\alpha\kappa^{{1}/{2}}(y+1)]\cos[\alpha\kappa^{-{1}/{2}}(y-1)]+\cosh[\alpha\kappa^{{1}/{2}}(y-1)] \cos[\alpha\kappa^{-{1}/{2}}(y+1)]}{\cosh[2\alpha\kappa^{{1}/{2}}]+\cos[2\alpha\kappa^{-{1}/{2}}]} \nonumber\\ &\quad \frac{-1\sinh[\alpha\kappa^{{1}/{2}}(y\!+\!1)]\sin[\alpha\kappa^{-{1}/{2}}(y-1)]\!+\! \sinh[\alpha\kappa^{{1}/{2}}(y-1)]\sin[\alpha\kappa^{-{1}/{2}}(y-1)]}{+ {\rm i}\,\cosh[2\alpha\kappa^{{1}/{2}}]+ \cos[2\alpha\kappa^{-{1}/{2}}]}, \end{align}

where

(4.5)\begin{equation} \kappa={-}De \pm \sqrt{De^{2} +1} \end{equation}

and

(4.6)\begin{equation} \xi=a_{\xi}+{\rm i}b_{\xi}. \end{equation}

In (4.4), $\alpha \kappa ^{1/2}$ is the attenuation of the velocity wave, and $\alpha \kappa ^{-1/2}$ is the phase shift. Since the trigonometric functions in (4.4) had periods of $\alpha \kappa ^{-1/2}$, then the angular frequency (through increments of $De=\lambda _{1}\omega$) can periodically attain these values of the phase shift, creating resonance. In the following lines, what interests us is the real part of the product of $z_{1}\cdot z_{2}$ of (2.60), which can be expressed as $a_{1}a_{2}-b_{1}b_{2}$, where $a_{1}$ and $b_{1}$ are the real and imaginary parts of (4.3), respectively. The selection of the real part comes from choosing the time sine function solution of (2.31). As mentioned, the resonance points arise near the zeros of $U(x,y)$ with the tendency to be located near $y=0$, where only $z_{2}$ depends on the $y$ coordinate. Thus, from (4.4) we evaluate $a_{2}$ and $b_{2}$ at $y=0$ to estimate these zeros, where we found the following approximation for $De_{n}$ values corresponding to the number of cycles of the $n{\rm th}$ harmonic in the analysis interval of $b_{2}=0$:

(4.7)\begin{equation} De_{n}=\frac{{\rm \pi}^{4}n^4-\alpha^4}{2\alpha^2{\rm \pi}^{2}n^2}. \end{equation}

Here, $n \in \mathbb {N^{+}}$, $De_{n}\geqslant 0$. Since $b_{2}=0$ at $De_{n}$, we evaluate $a_{1}$ at $De_{n}$ to figure out why $a_{1} a_{2}=0$. First, we numerically calculate the complex values of $\xi$ from (2.58). This computation shows that $\xi$ at $De_{n}$ has a small real part ($\gamma \ll 1$) while its imaginary part approximates ${\rm \pi}$ $( \xi \rvert _{De_{n}}=\gamma +{\rm i}{\rm \pi} )$. If we rewrite (4.3) with the results of $\xi$, and expand using power series of $\gamma$, we get the following equation:

(4.8)\begin{equation} \left.z_{1}\right\vert_{De_{n}}= 1 +{\rm i} \frac{2}{\gamma}+O(\gamma). \end{equation}

From (4.8) we get $a_{1}\neq 0$, but the product $z_{1}\cdot z_{2}$ still makes $Re[U(0,y)]=0$. This suggests that $z_{1}$ rotates $z_{2}$ by ${\rm \pi} /2$ in the complex plane. Therefore, we express the leading order in (4.8) in its polar form in such a manner that we can observe this rotation:

(4.9)\begin{equation} z_{1}= \frac{\sqrt{\gamma^{2}+4}}{\gamma} \exp\left[{\rm i}\phi(\gamma)\right], \end{equation}

where

(4.10)\begin{equation} \phi(\gamma) =\arctan\left(\frac{2}{\gamma}\right). \end{equation}

In the limit as $\gamma \rightarrow 0$, (4.10) reduces to the form

(4.11)\begin{equation} \lim_{\gamma \to 0} \phi(\gamma) =\frac{\rm \pi}{2}. \end{equation}

From (4.9) and (4.11) we conclude that $z_{1}\rvert _{De_{n}}=i$. Therefore, $U(0,y)$ can be expressed as $U(0,y)\rvert _{De_{n}}=-b_{2}+{\rm i} a_{2}$, where $b_{2}\rvert _{De_{n}}=0$, showing that its real part is zero.

Figure 3.

(a) Level curves of the real part of $U(0,y)$ at $x=0$ as a function of $y$ and $De$ number at $\alpha =20$. (b) Three-dimensional plot of the real part of $U(0,y)=0$ at $x=0$ as a function of $y$ and $De$ number at $\alpha =20$. (c) Three-dimensional plot of the real part of $\hat {P}(x)$ (from (2.61)) as a function of $x$ and $De$ number at $\alpha =20$. (d) Three-dimensional plot of the maximum value of the real part of $\hat {P}(x)$ as a function of $\alpha$ and $De$ number at $x=-0.5$.

With this expression for $De_{n}$, we find that resonance modes are in the vicinity of the values obtained from (4.7). Furthermore, the number $n$ defines the number of resonance modes near $De_{n}$, e.g. $De_{28}=9.6$ at $\alpha =20$ corresponds to 28 resonance modes (see figure 3a). From this statement we observe that for increasing values of $De$, the resonance modes of $U(x,y)$ increase in quantity along the coordinate $y$, where $De_n$ depends notably on $\alpha$, as can be observed from (4.7). Another feature that is worth noting corresponds to the real part of the pressure $\hat {P}(x)$, which has local maxima and minima at these values of $De_{n}$ along the coordinate $x$; unfortunately, at $\alpha =20$, the $De_{n}$ are so close to each other that $\hat {P}(x)$ looks like a smooth surface, as seen in figure 3(c). This smoothness can be observed in figure 3(d), showing that with increasing values of $\alpha$ the acoustic pressure becomes flatter. We have found that the maximum acoustic pressure varies according to $\alpha ^{6/5}De^{1/2}$, as shown in figure 3(d). Equation (4.7) can be rewritten as follows:

(4.12)\begin{equation} De\left(q\right)=\frac{1}{2}\left(q^{2}-\frac{1}{q^{2}}\right), \end{equation}

where $q=n {\rm \pi}/\alpha$. In figure 4, we plot $De(q)$ versus $q$, which can be used in either of two forms: for a given $q$, the corresponding $De_n$ is obtained; or for a given $De_n$, $q$ is determined. For this latter case, from (4.12), we obtain

(4.13)\begin{equation} q=\sqrt{De+\sqrt{De^{2}+1}}. \end{equation}

For a fluid with a given $De$ number, we can obtain the corresponding value of $q$ from (4.13), implying that we can estimate the parameter $\alpha$, fixing $n$, or vice versa; in this case, care must be taken in considering that the resonances modes $n$ must be a natural number. From figure 4, as a particular case, in the limit of $De \rightarrow 0$, $q=1$, meaning that $\alpha = n {\rm \pi}$.

Figure 4.

$De(q)$ evaluated at different values of $q$, where $q=n{\rm \pi} /\alpha$.

To understand the time-average stream function, (3.19), we first need to analyse the time-average vorticity function, $\langle \varOmega _{1} \rangle$. The time-average vorticity function is shown in (3.20), where the left-hand side term is of $O(1)$, while the right-hand terms are of $O(\alpha ^{3}De^{1/2})$, neglecting terms of $O(\epsilon )$. From (3.20) the previous order of magnitude for the right-hand side terms is not easily observed, that is why we present how the order of magnitude for the first term appearing on the right-hand side, $2\alpha ^2 \langle u_0 ({\partial \varOmega _0}/{\partial x}\rangle )$, is obtained. From (2.60) we know that $u_0$ is $O(1)$. Taking the derivative with respect to $x$ of (3.3), and considering the leading order, we have ${\partial \varOmega _{0}}/{\partial x}=-{\partial ^2 u_0 }/{\partial x \partial y} + O(\epsilon ^2)$. Doing the cross derivative of (2.60) we obtain ${\partial ^2 u_0}/{\partial x \partial y} = \zeta \xi ({\sinh (\xi x)}/{\cosh ({\xi }/{2})}) {\sinh (\zeta y)}/{\cosh (\zeta )}$. The prefactors of this latter equation, $\zeta$ and $\xi$, are of $O(\alpha De^{1/2})$ and $O(1)$, respectively. Multiplying all factors we obtain that $2\alpha ^2 \langle u_0 ({\partial \varOmega _0}/{\partial x}\rangle )$ is approximately of $O(\alpha ^3 De ^{1/2})$. Conceptually, the first two right-hand terms are the convective derivative of Cauchy equation. The third, fourth and fifth right-hand terms come from the shear stresses of the Maxwell's model, where the third and fourth are convective terms while the fifth is related to the deformation of the material. Due to these orders of magnitudes, ($\alpha ^{3} De^{{1}/{2}}<\eta ^{-1}$), so that the perturbation method is valid.

In (3.19) the time-average stream function depends on two terms, the time-average vorticity and $\langle \rho _{1} \partial u_{0}/\partial y\rangle$. Nyborg (Reference Nyborg and Mason1965) proposed using the time-average mass-transport velocity, which takes the form $\boldsymbol {u}^{M}=\langle \boldsymbol {u}_{1}\rangle +\langle \rho _{1}\boldsymbol {u}_{0}\rangle$, to simplify the latter term. The average mass-transport velocity relies upon the following points: (1) $\langle \rho _{1}\boldsymbol {u}_{0}\rangle =0$ which occur if $\rho _{1}$ and $\boldsymbol {u}_{0}$ differ in phase by $90^{\circ }$ (Hamilton et al. Reference Hamilton, Ilinskii and Zabolotskaya2003a), or if $\rho$ is aligned with the acoustic pressure $P$ ($\rho =\rho (P)$) then its nodes appear at the same time as $\boldsymbol {u}$ antinodes and both get cancelled; (2) $\langle \rho _{1} \partial u_{0}/\partial y\rangle =0$, following the same logic as point (1). Further calculations show that, due to the nature of the complex standing waves in two dimensions, $\langle \rho _{1} \partial u_{0}/\partial y\rangle$ do not differ in phase. Nevertheless, taking into account this term only modifies the magnitude of the streaming cells, while their size and direction stay the same even at different values of $\alpha$ or $De$. Therefore, we neglected it in the following discussion.

We begin this analysis with the known phenomenon that occurs in Newtonian fluids for the parameters presented in this work. Then, we observe how its behaviour is modified for the case of using a Maxwell fluid. In figures 5(a) and 5(b), the Newtonian fluid case ($De=0$) of the time-average mass-transport is shown for different values of $\alpha$. Figure 5(a) shows both Stokes’ and Rayleigh's streaming cells at $\alpha =4.22$, while in figure 5(b), with an increase of the microchannel height ($\alpha =20$), the outer vortices overcome the inner streaming cells. This transition where the Rayleigh's vortices start overlapping the Stokes streaming cells in narrow channels has already been reported (Rayleigh Reference Rayleigh1884; Hamilton et al. Reference Hamilton, Ilinskii and Zabolotskaya2003a). Furthermore, four rotating streaming cells appear in the domain, which correspond to the physical lengths selected (half the wavelength along the $x$ axis, and $2H_{0}$ on the $y$ axis).

Figure 5.

Acoustic streaming streamlines for the real part of the dimensionless time-average mass-transport field for (a) $\alpha =4.22$ with $De=0$, and $\alpha =20$ with (b) $De=0$, (c) $De=1.4$, (d) $De=1.5$, (e) $De=1.7$ and (f) $De=5$.

For a Maxwell fluid, its behaviour is shown in figures 5(c), 5(d), 5(e) and 5(f). In figure 5(c) the resurgence of the Stokes streaming cells is visible, while the Rayleigh vortices decrease in size (in comparison with its Newtonian counterpart shown in figure 5b). With a small increment of $De$, as shown in figure 5(d), the inner vortices dominate over the outer streaming cells, where the latter decrease in size and are pushed to the lateral walls. In figure 5(e) the reverse flow appears, where the inner streaming cells occupy the whole channel. This second transition between inner and outer vortices due to viscoelastic effects occurs in the interval $1< De<2$, and is induced by the fourth term on the right-hand side of (3.20) since the third and fifth terms cancel each other. That fourth term (convective term) is responsible for the growth of inner vortices and the reverse flow arises. The above assertion is based on the numeric experiments that we carried out. In figure 5(f) we observe that the magnitude of the time-average mass-transport streamlines is larger than the Newtonian case (see figure 5b), showing that elasticity can raise notably the velocity magnitude in acoustic streaming.

With this analysis we can stipulate the following. What is observed in the average mass-transport streamlines is the growth in size of the Stokes vortices and the diminishing of outer vortices due to increments of $De$. This phenomenon is related to the attenuation of the velocity wave, as $De$ increases. The above can be inferred from (2.48), which can be rewritten in the following form:

(4.14)\begin{equation} \zeta =\frac{\alpha}{\sqrt{-De + \sqrt{De^{2} +1}}} \left({-}De + \sqrt{De^{2} +1}+{i}\right). \end{equation}

The above equation, (4.14), suggests that the Stokes vortices expand to the entire channel. The above can be explained from the growth of the thickness of the boundary layer ($\delta _{\nu }$), which increases with increments of $De$ (Bird et al. Reference Bird, Armstrong and Hassager1987). When $De>1$, a second transition period, where the inner vortex dominates over the outer streaming cells, due to the convective terms of the constitutive Maxwell equation, appears. Then, the reverse flow arises when the Stokes vortices cover the whole microchannel (overlapping the Rayleigh vortices). Our hypothesis for the above sentence is that the maximum magnitude of the horizontal velocity, obtained for the Newtonian case, must be recovered for the viscoelastic fluid case at a distance from the wall equal to $6\delta _{\nu }/5$. In addition, while there are resonance modes for $De$ of $O(1)$ at $\alpha \gg 1$, they do not modify the time-average mass-transport streamlines because they are located close to the walls (see figures 3a and 3b as a reference).

Finally, it is of interest to know the maximum value of the time-average mass-transport speed $\langle s^{M}_{1} \rangle$ in terms of $\alpha$ and $De$. To calculate $\langle s^{M}_{1} \rangle$, it is necessary to evaluate the time average of the first order of (3.1a,b) as follows:

(4.15)\begin{equation} \langle s^{M}_{1} \rangle =\sqrt{\frac{\partial \langle \psi^{M}_{1}\rangle}{\partial x}^{2}+ \frac{\partial \langle \psi^{M}_{1}\rangle}{\partial y}^{2}}. \end{equation}

Through multiple numerical iterations of (4.15) the maximum value of $\langle s_1^{M}\rangle$ is obtained, finding that $\langle s^{M}_{1} \rangle \sim \alpha ^{2.65}De^{1.9}$. In figures 6(a) and 6(b) the time-average mass-transport speed is shown for different values of $\alpha$ and $De$. One can observe that at high Deborah numbers the effects of the resonance modes modify $\langle s_1^{M}\rangle$, and when $\alpha \gg 1$ and $De>8$ the maximum magnitude of the mass-transport speed will show higher affinity for the resonance modes. In this latter case, $\alpha ^{2.65}De^{1.9}$ is no longer valid.

Figure 6.

Maximum value of the time-average mass-transport speed $\langle s^{M}_{1} \rangle$ as a function of (a) $\alpha$ for representative values of $De$ ($2\leqslant De \leqslant 10$), and as a function of (b) $De$ for representative values of $\alpha$ ($5\leqslant \alpha \leqslant 20$).

From a biological standpoint of view, blood is mainly composed of red blood cells (RBCs) and platelets, where its viscoelastic properties come from the RBCs (Qi & Shaqfeh Reference Qi and Shaqfeh2018). If blood is under a shear flow, the platelets are redistributed uniformly near the walls due to enhanced RBC mixing in the central region of the channel, where the latter migrate towards the centre due to their own deformation (Wohl & Rubinow Reference Wohl and Rubinow1974; Crowl & Fogelson Reference Crowl and Fogelson2011). As described by Antfolk & Laurell (Reference Antfolk and Laurell2019), sometimes RBCs are mixed with bacteria. In the case of medium switching, the higher acoustic impedance ($Z_{j}\equiv \rho _{j}c_{j}$, for a material $j$) can induce that the bacteria end up together with the RBCs rather than staying at the side, ensuring that they would transit to the central stream. While this occurs in a Newtonian fluid, for a viscoelastic fluid, RBCs will migrate to an intermediate equilibrium position between the centreline and the wall, where the deformation and viscoelastic effects cancel each other (Chan & Leal Reference Chan and Leal1979). With the aforementioned we can compare our solutions in figure 5. If it follows increments of $De$, the acoustic streaming will show a transition between inner and outer vortices, where a reverse flow will appear at $De>1$ and the RBCs will remain at the centre. This behaviour doesn't correspond to the one obtained by Chan & Leal (Reference Chan and Leal1979), which suggests that an analysis where we consider $De\gg 1$ (fluid as an elastic solid) could show this behaviour.