2.1 Experimental set-up

Surface acoustic waves were generated at the surface of an $X$ -cut lithium niobate piezoelectric substrate in the crystallographic $Z$ -direction by interdigitated electrodes, with a spatial period of 175  $\unicode[STIX]{x03BC}$ m corresponding to a resonant frequency of 19.9 MHz (the sound speed in this direction is $3484~\text{m}~\text{s}^{-1}$ , see Campbell & Jones (Reference Campbell and Jones1970)). In practice, the best actuation efficiency was obtained at 20.37 MHz which was used as the driving frequency for all experiments. A water–glycerol droplet of 12.5  $\unicode[STIX]{x03BC}$ l was placed on the substrate initially treated with OTS (octadecyltrichlorosilane) self-assembled monolayer (SAM) to obtain hydrophobic wetting properties (see figure 1). The fraction of glycerol in the mixture modifies the shear and bulk viscosities with relatively weak variations of the other relevant parameters such as surface tension and contact angle (physical data are shown in table 1). Spherical latex beads of 10  $\unicode[STIX]{x03BC}$ m diameter coated with fluorescent dyes (ThermoScientific) were dispersed in the droplet prior to experiments to visualize its inner flow. We minimized the droplet evaporation by deporting the light source with an optical fibre, using the cold part of the optical spectrum and restraining the experiment duration below 2 min. The images were acquired via an Hamamatsu camera (4 megapixels) and quantitative velocity magnitude was measured using the particle image velocimetry (PIV) module of ImageJ (https://sites.google.com/site/qingzongtseng/piv). Additional details on the PIV experiment are provided in the supplemental material available at https://doi.org/10.1017/jfm.2017.178. We restricted the power of the SAW to a few tens of picometres to minimize the droplet deformation, which are observed at much higher amplitude (Schindler et al. Reference Schindler, Talkner and Hänggi2006; Brunet et al. Reference Brunet, Baudoin, Bou Matar and Zoueshtiagh2010; Alghane et al. Reference Alghane, Chen, Fu, Li, Desmulliez, Mohammed and Walton2012; Baudoin et al. Reference Baudoin, Brunet, Bou Matar and Herth2012). The substrate vertical amplitude of vibration on the central finger of the IDT was calibrated for the range of actuation power used in the present experiments with a laser Doppler vibrometer (SH130, B.M. Industries).

Finally, the droplets inner flow was visualized from below (to avoid optical aberrations by the drop surface) with a Hamamatsu high resolution camera mounted on an inverted microscope (Olympus IX71). The depth of field is estimated to be 16  $\unicode[STIX]{x03BC}$ m for objects of the size of a pixel when using the $4\times$ magnification objective with this inverted microscope, allowing the visualization of a droplet cross-section. The vertical position of focus plane was adjusted by eye as close as possible to half of the drop height, although this condition is achieved within a few per cent accuracy. After turning on the SAW generator, we waited for the droplet flow pattern to reach a steady state. This duration varied widely with viscosity, from seconds for water droplets up to minute for the most viscous mixtures (as expected from theoretical analysis). Figure 2 reproduces a few examples of such particles trajectories viewed from below, obtained at increasing viscosities from (a) to (f). The streamlines are obtained by a simple superposition of successive images.

Figure 2.

Flow visualization of a 12.5  $\unicode[STIX]{x03BC}$ l droplet from below, at various glycerol concentration. The SAW propagates from left to right with an amplitude $u_{0}=62\pm 1$  pm. As the viscosity increases, one remarks the progressive transition from a four-vortex to a two-vortex flow structure. (a) Pure water; (b) 30 wt.% glycerol; (c) 40 wt.% glycerol; (d) 60 wt.% glycerol; (e) 80 wt.% glycerol; (f) 90 wt.% glycerol.

Figure 3.

Experimental velocity field from below, at various viscosities. The SAWs propagate from left to right with an amplitude $u_{0}\simeq 62$  pm. $V_{droplet}=12.5~\unicode[STIX]{x03BC}\text{l}$ . As the viscosity increases, one remarks the progressive transition from a four-vortex to a two-vortex flow structure. (a) 10 wt.% glycerol ( ${\mathcal{U}}_{max}\simeq 180~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ ); (b) 30 wt.% glycerol ( ${\mathcal{U}}_{max}\simeq 100~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ ); (c) 40 wt.% glycerol ( ${\mathcal{U}}_{max}\simeq 70~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ ); (d) 90 wt.% glycerol ( ${\mathcal{U}}_{max}\simeq 10~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ ). The arrow length is indicative of the velocity magnitude for each experiment.

The velocity field in the representative cut was extracted from the pictures of the flow streams presented above by using the PIV ImageJ plugin (see figure 3). The analysis was further refined by discarding the 5 % least reliable velocity vectors. (The reliability criteria was the magnitude of the Laplacian of the velocity field.) We made sure that the system had reached steady state by waiting until the space-averaged magnitude of the velocity field did not vary by more than 10 % between two different time intervals. Each couple of images provides a flow map, and each plot in figure 3 is the average of three different flow maps obtained with the same droplet at steady state. Then the average velocity was determined by averaging spatially the flow field in the drop. The streaming pattern being inhomogeneous in space, this space-averaged flow velocity over one cross-section only gives an order of magnitude of the volume-averaged velocity. The resulting trend is presented in the discussion section together with numerical results (see figure 13).

3.1 Field decomposition

As stated in the introduction, we can resolve each physical quantity $f$ into three contributions: hydrostatic $f_{0}$ , acoustic $\tilde{f_{1}}$ and hydrodynamic $\bar{f_{2}}$ . They represent respectively the system at rest (without acoustic field), the oscillating part of the perturbation induced by sound waves and the time-averaged part of the perturbation over an acoustic period. In our experiments, acoustic and hydrodynamic Mach numbers are small. Moreover, solid displacements hardly exceed 0.5 nm, which restricts acoustic perturbation velocity magnitude to below $10~\text{mm}~\text{s}^{-1}$ and consequently streaming velocities below $1~\text{mm}~\text{s}^{-1}$ . According to this slow streaming approximation (Lighthill Reference Lighthill1978), the fluid density $\unicode[STIX]{x1D70C}$ the pressure $p$ and the Eulerian velocity $v$ , can be expressed as follows:

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}+\tilde{\unicode[STIX]{x1D70C}}_{1}+\bar{\unicode[STIX]{x1D70C}}_{2}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle p=p_{0}+\tilde{p}_{1}+\bar{p}_{2}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle v_{i}=\tilde{v}_{1,i}+\bar{v}_{2,i}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle v_{0,i}=0, & \displaystyle\end{eqnarray}$$

with $i\in \{x,y,z\}$ , $\bar{f_{2}}=\langle f-f_{0}\rangle$ , $\langle \rangle$ the time average, $\tilde{f_{1}}=f-f_{0}-\bar{f_{2}}$ , $\langle \tilde{f_{1}}\rangle =0$ , and $\bar{f_{2}}\ll \tilde{f_{1}}\ll f_{0}$ . The low Mach numbers assumption gives $\tilde{v}_{1},\bar{v}_{2}\ll c_{0}$ , with $c_{0}$ the sound speed in the fluid at hand. To simplify the notations, the indices $1$ and $2$ will be omitted in the following.

3.2 Fundamental equations

The starting point of the derivation is the compressible isentropic Navier–Stokes equations. These equations are relevant to compute acoustic streaming in liquids since, in this case, thermal effects (wave thermal damping, fluid heating) can be neglected compared to their viscous counterparts (viscous damping, acoustic streaming). Indeed, thermal effects are proportional to $(\unicode[STIX]{x1D6FE}-1)$ , with $\unicode[STIX]{x1D6FE}$ the adiabatic index and thus are very weak in liquids (see e.g. Coulouvrat Reference Coulouvrat1992).

In this case, the mass conservation equation for a fluid reads:

(3.5) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}+\unicode[STIX]{x2202}_{i}\unicode[STIX]{x1D70C}v_{i}=0,\end{eqnarray}$$

and the momentum conservation equation:

(3.6) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}v_{i}+\unicode[STIX]{x2202}_{j}(\unicode[STIX]{x1D70C}v_{i}v_{j})=-\unicode[STIX]{x2202}_{i}p+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}v_{i}+\left(\frac{\unicode[STIX]{x1D707}}{3}+\unicode[STIX]{x1D709}\right)\unicode[STIX]{x2202}_{ij}^{2}v_{j}.\end{eqnarray}$$

In these equations, $\unicode[STIX]{x1D707}$ stands for the dynamic viscosity, $\unicode[STIX]{x1D709}$ for the bulk viscosity, $t$ for the time and the indices $i$ and $j$ follow the Einstein summation convention. The second-order isotropic Taylor expansion of the equation of state reads:

(3.7) $$\begin{eqnarray}\text{d}p=c_{0}^{2}\,\text{d}\unicode[STIX]{x1D70C}+\frac{1}{2}\unicode[STIX]{x1D6E4}\,\text{d}\unicode[STIX]{x1D70C}^{2},\end{eqnarray}$$

with $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x2202}^{2}p/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{2}|_{s}=Bc_{0}^{2}/A\unicode[STIX]{x1D70C}_{0}$ . $A$ and $B$ are the first- and second-order coefficients of the Taylor expansion relating the material pressure to its density. These nonlinearity coefficients are classically introduced in nonlinear acoustics.

3.3 Time-averaged equations at second order: acoustic steady streaming

If we take the time average of the mass and momentum conservation equations (3.5) and (3.6) up to second order, and introduce the Poynting vector (also called intensity vector in the field of acoustics) $\unicode[STIX]{x1D6F1}_{i}=\tilde{p}\tilde{v}_{i}$ , we obtain:

(3.8) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\bar{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{i}\bar{v}_{i}+\frac{1}{c_{0}^{2}}\unicode[STIX]{x2202}_{i}\langle \unicode[STIX]{x1D6F1}_{i}\rangle =0,\end{eqnarray}$$

and:

(3.9) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}_{0}\bar{v_{i}}+1/c_{0}^{2}\langle \unicode[STIX]{x1D6F1}_{i}\rangle )+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{j}\langle \tilde{v}_{i}\tilde{v}_{j}\rangle =-\unicode[STIX]{x2202}_{i}\bar{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\bar{v}_{i}+\left(\frac{\unicode[STIX]{x1D707}}{3}+\unicode[STIX]{x1D709}\right)\unicode[STIX]{x2202}_{ij}^{2}\bar{v}_{j},\end{eqnarray}$$

since $\unicode[STIX]{x1D6F1}_{i}=c_{0}^{2}\tilde{\unicode[STIX]{x1D70C}}\tilde{v}_{i}$ at leading order.

These two equations can be simplified for a large span of fluids and systems using the following assumptions. First, if we consider the acoustic streaming produced by bulk acoustic wave (away from boundaries), the third term of the mass conservation equation (3.8) is proportional at leading order to the bulk viscous dissipation of the wave energy, which remains weak in most media. This is quantified by the acoustical Reynolds number $Re_{ac}$ which compares the viscous dissipation to inertia or equivalently the wave attenuation length $L_{a}=\unicode[STIX]{x1D70C}c_{0}^{3}/\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D707}(4/3+\unicode[STIX]{x1D709}/\unicode[STIX]{x1D707})$ to the wavelength $\unicode[STIX]{x1D706}$ :

$$\begin{eqnarray}Re_{ac}=\frac{L_{a}}{\unicode[STIX]{x1D706}}=\frac{\unicode[STIX]{x1D70C}_{0}c_{0}^{2}}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D707}\left(4/3+\displaystyle \frac{\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D707}}\right)}.\end{eqnarray}$$

Except at very high frequency ( ${>}1$  GHz) or for extremely viscous fluids and high driving frequencies, the acoustical Reynolds number is generally high ( $Re_{ac}\gg 1$ ). It is worth to mention that a quite different Reynolds number could be defined for the solenoidal part of the velocity field that results from the high shear rate in the Stokes boundary layer (Rayleigh Reference Rayleigh1915; Vanneste & Bühler Reference Vanneste and Bühler2011). In this case, the attenuation length should be the boundary layer thickness. Although this region represents a tiny fraction of the droplet size at the acoustical frequencies considered here (thickness $\unicode[STIX]{x1D6FF}_{S}=\sqrt{2\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}}\simeq 0.12~\unicode[STIX]{x03BC}\text{m}$ for water at 20 MHz up to 1.8  $\unicode[STIX]{x03BC}$ m for 90 wt.% glycerol mixture), it may be argued that the boundary layer streaming extends much further than the Stokes boundary layer itself. However, according to Vanneste & Bühler (Reference Vanneste and Bühler2011) comprehensive calculations in a simplified 2-D geometry, the boundary layer streaming is negligible compared to the bulk streaming as long as the fluid vessel size and the SAW attenuation length both greatly exceed the SAW wavelength, which is the case here. We will therefore only consider the bulk waves which yields a high acoustical Reynolds number as estimated for the frequency and liquids used in the present experiments and simulations in table 1. Finally, since the boundary layer streaming is unimportant here, we also neglect the tangential displacement associated with the surface acoustic wave (see Vanneste & Bühler (Reference Vanneste and Bühler2011) for a comprehensive analysis).

Moreover, if we consider only steady streaming (stationary flow produced by acoustic waves), the time derivatives in (3.8) and (3.9) can be cancelled out. We obtain in this case:

(3.10) $$\begin{eqnarray}\unicode[STIX]{x2202}_{i}\bar{v}_{i}=0,\end{eqnarray}$$

which amounts to saying that the steady streaming flow is incompressible. Then the time-averaged momentum conservation equation becomes:

(3.11) $$\begin{eqnarray}-\!\unicode[STIX]{x2202}_{i}\bar{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\bar{v}_{i}+{\mathcal{F}}_{i}=0,\end{eqnarray}$$

with ${\mathcal{F}}_{i}$ the Reynolds stress imbalance of the sound wave:

(3.12) $$\begin{eqnarray}{\mathcal{F}}_{i}=-\langle \unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{j}(\tilde{v}_{i}\tilde{v}_{j})\rangle =-\unicode[STIX]{x1D70C}_{0}\langle \tilde{v}_{j}\unicode[STIX]{x2202}_{j}\tilde{v}_{i}+\tilde{v}_{i}\unicode[STIX]{x2202}_{j}\tilde{v}_{j}\rangle .\end{eqnarray}$$

This equation is simply the steady-state Stokes equation driven by a forcing term ${\mathcal{F}}_{i}$ resulting from average nonlinear interactions of the acoustic field.

It is worth noting that the derivation of acoustic streaming constitutive equations follows a similar procedure as the one used for the derivation of the Reynolds-averaged Navier–Stokes equation in the field of turbulence. It describes how some strong fluctuating nonlinear terms influence the steady flow (Bühler Reference Bühler2009; Vanneste & Bühler Reference Vanneste and Bühler2011). Nevertheless, the fundamental differences between the derivation of the constitutive equations of acoustic streaming and turbulence are (i) that owing to the weak amplitude of the acoustic field, a perturbation analysis is possible and (ii) that the source term in the average equations emanates in the former case from the first-order compressible field, namely the acoustic wave.

3.4 Periodic fluctuations up to second order: nonlinear acoustics

The mass and momentum equations for the periodic fluctuations $\tilde{f}$ up to second order can be obtained by subtracting the time-averaged equations (3.8) and (3.9) from the initial Navier–Stokes isentropic equations (3.5) and (3.6):

(3.13) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{i}\tilde{v}_{i}=-\unicode[STIX]{x2202}_{i}\langle \!\langle \tilde{\unicode[STIX]{x1D70C}}\tilde{v_{i}}\rangle \!\rangle ,\end{eqnarray}$$

and:

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{t}\tilde{v}_{i}+\unicode[STIX]{x2202}_{i}\tilde{p}-\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\tilde{v}_{i}-\left(\frac{\unicode[STIX]{x1D707}}{3}+\unicode[STIX]{x1D709}\right)\unicode[STIX]{x2202}_{ij}^{2}\tilde{v}_{j}=-\unicode[STIX]{x2202}_{t}\langle \!\langle \tilde{\unicode[STIX]{x1D70C}}\tilde{v_{i}}\rangle \!\rangle -\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{j}\langle \!\langle \,\tilde{v_{i}}\tilde{v_{j}}\rangle \!\rangle ,\end{eqnarray}$$

with $\langle \!\langle \tilde{f}\tilde{g}\rangle \!\rangle =\tilde{f}\tilde{g}-\langle \,\tilde{f}\tilde{g}\rangle$ .

The left-hand sides of (3.13) and (3.14), along with the equation of state (3.7) at first order, constitute the linear equations of damped acoustic waves. The right-hand sides of these equations correspond to nonlinear terms, which modify the propagation of acoustic waves through energy transfers to harmonic frequencies ( $2\unicode[STIX]{x1D714},3\unicode[STIX]{x1D714},\ldots$ ).

If we assume (see previous section) that the acoustical Reynolds number is high and the acoustical Mach number is small, these equations become at leading order:

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{i}\tilde{v}_{i}=0, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{t}\tilde{v}_{i}+\unicode[STIX]{x2202}_{i}\tilde{p}=0, & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \text{with }\tilde{p}=c_{0}^{2}\tilde{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

which amounts to discarding all nonlinear and dissipative effects. From (3.16), we can infer that the oscillating flow is potential at leading order ( $\tilde{v}_{i}=-\unicode[STIX]{x2202}_{i}(\tilde{\unicode[STIX]{x1D719}})$ , with $\tilde{\unicode[STIX]{x1D719}}$ the velocity potential). A simple combination of (3.15) and (3.17): $[\!c_{0}^{2}\;\unicode[STIX]{x2202}_{i}$ (3.15)  $-\,\unicode[STIX]{x2202}_{t}$ (3.16)] with (3.17) yields the celebrated d’Alembert equation:

(3.18) $$\begin{eqnarray}\unicode[STIX]{x2202}_{tt}^{2}\tilde{\unicode[STIX]{x1D719}}-c_{0}^{2}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}=0,\end{eqnarray}$$

with $\tilde{p}=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D719}}$ and $\tilde{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{0}/c_{0}^{2}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D719}}$ .

Now, if we do the same combination of equations (3.13), (3.14) and (3.7) but up to next order in $M$ and $1/Re_{ac}$ , we obtain the Kuznetsov equation (Kuznetsov Reference Kuznetsov1970) (see e.g. Coulouvrat (Reference Coulouvrat1992) for a detailed demonstration with asymptotic analysis):

(3.19) $$\begin{eqnarray}\unicode[STIX]{x2202}_{tt}^{2}\tilde{\unicode[STIX]{x1D719}}-c_{0}^{2}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}-\frac{\unicode[STIX]{x1D707}b}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x2202}_{jj}^{2}\tilde{\unicode[STIX]{x1D719}}=\unicode[STIX]{x2202}_{t}\left(\frac{B}{2Ac_{0}^{2}}\langle \!\langle (\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D719}})^{2}\rangle \!\rangle +\langle \!\langle (\unicode[STIX]{x2202}_{i}\tilde{\unicode[STIX]{x1D719}})^{2}\rangle \!\rangle \right),\end{eqnarray}$$

with $b=4/3+\unicode[STIX]{x1D709}/\unicode[STIX]{x1D707}$ . Finally in the paraxial approximation (weak diffraction of the beam), we have the Lighthill–Westervelt equation $(\unicode[STIX]{x2202}_{i}\tilde{\unicode[STIX]{x1D719}})^{2}=(1/c_{0}^{2})(\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D719}})^{2}$ , leading to:

(3.20) $$\begin{eqnarray}\unicode[STIX]{x2202}_{tt}^{2}\tilde{\unicode[STIX]{x1D719}}-c_{0}^{2}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}-\frac{\unicode[STIX]{x1D707}b}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x2202}_{jj}^{2}\tilde{\unicode[STIX]{x1D719}}=\frac{\unicode[STIX]{x1D6FD}}{c_{0}^{2}}\unicode[STIX]{x2202}_{t}\langle \!\langle (\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D719}})^{2}\rangle \!\rangle ,\end{eqnarray}$$

with $\unicode[STIX]{x1D6FD}=1+B/2A$ the so-called nonlinear parameter. This equation allows us to compute the damped nonlinear propagation of acoustic waves.

The question then arises as to whether the nonlinear propagation of the acoustic wave must be considered to compute the acoustic streaming sources in (3.11)? An elementary analysis solely based on the order of the nonlinear terms might lead to the misleading premature conclusion that since acoustic nonlinear terms on the right-hand side of (3.19) are of second order, their quadratic combination in (3.12) is of fourth order and thus could be safely neglected when computing the acoustic streaming. In fact these nonlinear terms are weak but nevertheless cumulative, so they can play a significant role over a distance called the ‘shock distance’ $L_{s}=c_{0}^{2}/\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FD}U_{ac}$ , which depends on the first-order velocity magnitude $U_{ac}$ . To answer this question correctly, it must thus be recalled that acoustic streaming is a consequence of the attenuation of the acoustic wave. This attenuation is proportional to the square of the acoustic wave frequency $\unicode[STIX]{x1D714}^{2}$ . Since nonlinear terms in (3.20) induce energy transfers from the driving frequency to higher harmonics, they promote the dissipation and thus the acoustic streaming. In an unbounded medium, the streaming enhancement by harmonics generation can be quantified by the ratio of the wave attenuation length $L_{a}=\unicode[STIX]{x1D70C}_{0}c_{0}^{3}/\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D707}b$ to the shock distance $L_{s}$ :

(3.21) $$\begin{eqnarray}\frac{L_{a}}{L_{s}}=\frac{\unicode[STIX]{x1D70C}_{0}c_{0}\unicode[STIX]{x1D6FD}U_{ac}}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D707}b}.\end{eqnarray}$$

Nonlinear terms in (3.20) can thus be neglected when $L_{a}/L_{s}\ll 1$ . In the present situation, since $U_{ac}<10~\text{mm}~\text{s}^{-1}$ , the maximum value of this ratio is $10^{-1}$ for water and goes down to $10^{-3}$ for water–glycerol mixtures. Moreover, in cavities with water/air interfaces such as drops, the shock distance must also be compared to the size of the cavity $L_{c}$ . Indeed, nonlinear effects are only significant when they are cumulative. Since it was shown by Tanter et al. (Reference Tanter, Thomas, Coulouvrat and Fink2001) that each wave reflection at an air–water interface results in the deconstruction of these nonlinear effects, these latter can only be significant if the characteristic size of the cavity $L_{c}$ is larger than the shock distance. Here $L_{c}\sim 1$  mm while $L_{s}\sim 1$  m. Consequently, nonlinear terms can be safely discarded in (3.20) for the analysis and simulation of the present experiments, leading to the equation of damped acoustic waves:

(3.22) $$\begin{eqnarray}\unicode[STIX]{x2202}_{tt}^{2}\tilde{\unicode[STIX]{x1D719}}-c_{0}^{2}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}-\frac{\unicode[STIX]{x1D707}b}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x2202}_{jj}^{2}\tilde{\unicode[STIX]{x1D719}}=0.\end{eqnarray}$$

3.5 Streaming source term: the hydrodynamic Reynolds stress tensor

In this section we follow the Eckart (Reference Eckart1948) and Vanneste & Bühler (Reference Vanneste and Bühler2011) guidelines to analyse the different contributions of the streaming source (3.12). We discard all the gradient terms that are balanced by a change of static pressure and therefore to do not contribute to acoustic streaming. This simplification is essential since the magnitude of the neglected terms is much larger than that of the relevant terms and can lead to significant numerical error when calculating the flow produced by the acoustic wave.

Since the wave perturbation is irrotational, the first term of the force in (3.12) is easily integrated and identified as the kinetic energy ${\mathcal{K}}=(1/2)\unicode[STIX]{x1D70C}_{0}\tilde{v}^{2}$ :

(3.23) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{0}\langle \tilde{v}_{j}\unicode[STIX]{x2202}_{j}\tilde{v}_{i}\rangle =\unicode[STIX]{x2202}_{i}\langle {\mathcal{K}}\rangle .\end{eqnarray}$$

The second term is computed from the mass conservation equation (3.13) at leading order:

(3.24) $$\begin{eqnarray}\langle \unicode[STIX]{x1D70C}_{0}\tilde{v}_{i}\unicode[STIX]{x2202}_{j}\tilde{v}_{j}\rangle =-\langle \tilde{v}_{i}\unicode[STIX]{x2202}_{t}\tilde{\unicode[STIX]{x1D70C}}\rangle =\langle \tilde{\unicode[STIX]{x1D70C}}\unicode[STIX]{x2202}_{t}\tilde{v}_{i}\rangle ,\end{eqnarray}$$

where we used integration by parts to move from the second expression to the third. Then, if we use the classical vector identity $\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=\unicode[STIX]{x1D6FB}^{2}\tilde{\boldsymbol{u}}+\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times \tilde{\boldsymbol{u}}$ , the wave momentum equation (3.14) yields:

(3.25) $$\begin{eqnarray}\langle \unicode[STIX]{x1D70C}_{0}\tilde{v}_{i}\unicode[STIX]{x2202}_{j}\tilde{v}_{j}\rangle =-\left\langle \frac{\tilde{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{i}\tilde{p}\right\rangle +\left(\frac{4}{3}\unicode[STIX]{x1D707}+\unicode[STIX]{x1D709}\right)\left\langle \frac{\tilde{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{jj}^{2}\tilde{v}_{i}\right\rangle .\end{eqnarray}$$

Finally, the equation of state (3.7) yields: $\langle (\tilde{\unicode[STIX]{x1D70C}}/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}_{i}\tilde{p}\rangle =\langle (c^{2}/2\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}_{i}\tilde{\unicode[STIX]{x1D70C}}^{2}\rangle$ . Consequently, this term is the gradient of the potential energy of the wave in a linear medium:

(3.26) $$\begin{eqnarray}\left\langle \frac{\tilde{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{i}\tilde{p}\right\rangle =\unicode[STIX]{x2202}_{i}\langle {\mathcal{V}}\rangle .\end{eqnarray}$$

If we combine (3.23), (3.25) and (3.26), we obtain the following expression of the force ${\mathcal{F}}_{i}$ (see e.g. Gusev & Rudenko Reference Gusev and Rudenko1979):

(3.27) $$\begin{eqnarray}{\mathcal{F}}_{i}=-\unicode[STIX]{x2202}_{i}\langle {\mathcal{L}}\rangle -\left(\frac{4}{3}\unicode[STIX]{x1D707}+\unicode[STIX]{x1D709}\right)\left\langle \frac{\tilde{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{jj}^{2}\tilde{v}_{i}\right\rangle ,\end{eqnarray}$$

with the acoustic Lagrangian ${\mathcal{L}}={\mathcal{K}}-{\mathcal{V}}$ . We can work out a more practical equation by substituting the linear undamped wave equation in the viscous term and assuming an harmonic wave motion:

(3.28) $$\begin{eqnarray}{\mathcal{F}}_{i}=-\unicode[STIX]{x2202}_{i}\langle {\mathcal{L}}\rangle +\frac{\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D708}b}{c^{4}}\langle \unicode[STIX]{x1D6F1}_{i}\rangle ,\end{eqnarray}$$

with $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}_{0}$ . The first term derives from a potential and is independent of the bulk and shear viscosities, and thus of the wave attenuation. Since acoustic streaming relies on the momentum transfer from the wave mode (irrotational, compressible) to the viscous mode (solenoidal), through the wave attenuation, this term does not contribute to the steady flow; it is simply balanced by a hydrostatic pressure gradient. This can be easily verified by recasting equation (3.11) in the form:

(3.29) $$\begin{eqnarray}-\!\unicode[STIX]{x2202}_{i}\bar{p^{\ast }}+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\bar{v}_{i}+{\mathcal{F}}_{i}^{\ast }=0,\end{eqnarray}$$

with $\bar{p}^{\ast }=\bar{p}+\langle {\mathcal{L}}\rangle$ and ${\mathcal{F}}_{i}^{\ast }=(\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D708}b/c^{4})\langle \unicode[STIX]{x1D6F1}_{i}\rangle$ . ${\mathcal{F}}_{i}^{\ast }$ is related to the wave dissipation and has a much smaller magnitude than its counterpart $\langle {\mathcal{L}}\rangle$ . Nevertheless, it is not potential and hence is the sole source term of acoustic bulk streaming, as emphasized in Lighthill (Reference Lighthill1978) and Eckart (Reference Eckart1948) studies.

3.6 Final simplified system of equations

We recall briefly here the optimized expressions to simulate the experimentally observed acoustic streaming in droplets:

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle \text{acoustic wave: }\unicode[STIX]{x2202}_{tt}^{2}\tilde{\unicode[STIX]{x1D719}}-c_{0}^{2}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}-\frac{\unicode[STIX]{x1D707}b}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x2202}_{jj}^{2}\tilde{\unicode[STIX]{x1D719}}=0, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\displaystyle & \displaystyle \text{streaming Stokes flow: }\left\{\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{i}\bar{v}_{i}=0,\\ -\unicode[STIX]{x2202}_{i}\bar{p^{\ast }}+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\bar{v}_{i}+{\mathcal{F}}_{i}^{\ast }=0,\end{array}\right. & \displaystyle\end{eqnarray}$$

(3.32) $$\begin{eqnarray}\displaystyle & \displaystyle \text{acoustic streaming source term: }{\mathcal{F}}_{i}^{\ast }=\frac{\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D708}b}{c^{4}}\langle \unicode[STIX]{x1D6F1}_{i}\rangle . & \displaystyle\end{eqnarray}$$

4.1 Computation of the acoustic field

The acoustic field is computed in the frequency domain. In this case, equation (3.30) becomes:

(4.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{ii}^{2}\tilde{\unicode[STIX]{x1D719}}+k^{2}\tilde{\unicode[STIX]{x1D719}}=0,\end{eqnarray}$$

with:

(4.2) $$\begin{eqnarray}k^{2}=\frac{k_{0}^{2}}{1+\text{i}/Re_{ac}}.\end{eqnarray}$$

Here, $k_{0}=\unicode[STIX]{x1D714}/c_{0}$ is the wavenumber of the unattenuated wave and $Re_{ac}$ is the acoustical Reynolds number ( $Re_{ac}\gg 1$ ).

The large discrepancy between the droplet size and the acoustic wavelength yields very large and intensive simulations. For instance, direct 3-D simulation of the acoustic field on a 32 GB RAM computer with the finite element method only allows us to simulate 2 mm diameter droplets up to 8 MHz. As shown in the supplemental material, RAM requirements sharply increase with increasing frequency, and extrapolation to 20 MHz culminates at 1.0 TB, thus preventing any direct computation of the acoustic field.

To minimize memory requirements, we took advantage of the droplet rotational symmetry to reduce the problem to a two-dimensional one. The protocol described in the following uses Fourier transform to resolve the incident field as a sum of circular harmonics, solves each of them separately and then reconstructs the field thanks to the superposition principle. In this way, the complete problem is decomposed into subproblems with low memory requirements which can be computed in parallel. This method thus ensures an optimal matching with the capacity of the computer (number of cores, memory).

4.1.1 Method: spatial Fourier transform

Working in cylindrical coordinates, Fourier transform allows resolving of any function into a convenient weighted sum of complex exponentials:

(4.3) $$\begin{eqnarray}f(r,\unicode[STIX]{x1D703})=\mathop{\sum }_{l=-\infty }^{+\infty }f_{l}(r)\text{e}^{\text{i}l\unicode[STIX]{x1D703}},\end{eqnarray}$$

with

(4.4) $$\begin{eqnarray}f_{l}(r)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{+\unicode[STIX]{x03C0}}f(r,\unicode[STIX]{x1D703})\text{e}^{-\text{i}l\unicode[STIX]{x1D703}}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

Here, the only non-axisymmetric boundary condition is the normal displacement $\tilde{u}$ of the substrate due to the incident SAW. It is projected into Fourier harmonics:

(4.5) $$\begin{eqnarray}\tilde{u} _{l}(r)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{+\unicode[STIX]{x03C0}}\tilde{u} (r,\unicode[STIX]{x1D703})\text{e}^{-\text{i}l\unicode[STIX]{x1D703}}\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$

In practice, the value of $l$ can be restricted. Indeed, $l_{max}\simeq \unicode[STIX]{x03C0}D/\unicode[STIX]{x1D706}_{s}$ corresponds to the maximum number of wavelengths $\unicode[STIX]{x1D706}_{s}$ the input surface acoustic wave can travel along the perimeter of the droplet, where $D$ is the droplet diameter. We computed this integral for $l\in [0,2l_{max}]$ , $l$ being an integer. The value of $\tilde{u}$ depends on the incident wave. Neglecting diffraction, the SAW magnitude decreases exponentially as soon as it meets the droplet interface at a given point $x_{0}$ . The attenuation rate $\unicode[STIX]{x1D6FC}$ is provided for instance by Campbell & Jones (Reference Campbell and Jones1970). For a given point $M(x,y)$ , the propagation length beneath the droplet is given by $x-x_{0}(y)$ (see figure 5). The vertical displacement field $\tilde{u}$ at the droplet base is then given by Shiokawa et al. (Reference Shiokawa, Matsui and Ueda1990):

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} =u_{0}\exp (-\text{i}k_{s}x)\exp (-\unicode[STIX]{x1D6FC}(x-x_{0}(y))), & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle x_{0}(y)=-\sqrt{R^{2}-y^{2}}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{0}\ln (10)F/20, & \displaystyle\end{eqnarray}$$

where $F$ is the SAW frequency in Hz, $u_{0}$ is the magnitude of the acoustic perturbation displacement and $\unicode[STIX]{x1D6FC}_{0}$ the attenuation coefficient (in $\text{s}~\text{m}^{-1}$ ) $\unicode[STIX]{x1D6FC}_{0}\simeq 2.0\times 10^{-7}\times \unicode[STIX]{x1D70C}_{0}$ . We computed $\unicode[STIX]{x1D6FC}_{0}$ from the value given by Campbell & Jones (Reference Campbell and Jones1970) in the case of water loading ( $0.2~\text{dB}~\text{MHz}^{-1}~\text{cm}^{-1}$ ).

Figure 5.

Model of the incident leaky SAW. Colours are indicative of the SAW magnitude. At 20 MHz and for a $12.5~\unicode[STIX]{x03BC}\text{l}$ water droplet ( $2R=3.7$  mm), the incident SAW vertical displacement $\tilde{u}$ drops by 80 % as it propagates beneath the droplet.

We then solve all variables in the form: $\tilde{p}(r,\unicode[STIX]{x1D703},z)=\tilde{p}_{l}(r,z)\text{e}^{\text{i}l\unicode[STIX]{x1D703}}$ , $\tilde{v}_{j}(r,\unicode[STIX]{x1D703},z)=\tilde{v}_{j,l}(r,z)\text{e}^{\text{i}l\unicode[STIX]{x1D703}}$ with $j$ either $r,\unicode[STIX]{x1D703},z$ . Equation (4.1) becomes:

(4.9) $$\begin{eqnarray}k^{2}\tilde{p}=-\frac{1}{r}\unicode[STIX]{x2202}_{r}(r\unicode[STIX]{x2202}_{r}\tilde{p})-\frac{1}{r^{2}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}^{2}\tilde{p}-\unicode[STIX]{x2202}_{zz}^{2}\tilde{p},\end{eqnarray}$$

which can be re-cast using the axisymmetric variables:

(4.10) $$\begin{eqnarray}\left(k^{2}-\frac{l^{2}}{r^{2}}\right)\tilde{p}_{l}=-\unicode[STIX]{x2202}_{rr}^{2}\tilde{p}_{l}-\unicode[STIX]{x2202}_{zz}^{2}\tilde{p}_{l}-\frac{1}{r}\unicode[STIX]{x2202}_{r}\tilde{p}_{l}.\end{eqnarray}$$

This equation is solved with a finite element method by the commercial solver COMSOL 4.3b. (The default partial differential equation interface for axisymmetric systems does not include the last right-hand term of (4.10).) The boundary condition at the liquid–air interface reads:

(4.11) $$\begin{eqnarray}\tilde{p}_{l}=0.\end{eqnarray}$$

At the solid liquid interface, we enforce an impedance boundary condition with the source term computed from (4.5):

(4.12) $$\begin{eqnarray}\unicode[STIX]{x2202}_{z}\tilde{p}_{l}=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}^{2}\tilde{u} _{l}-\text{i}\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}{{\mathcal{Z}}_{s}}\tilde{p}_{l},\end{eqnarray}$$

where ${\mathcal{Z}}_{s}=\unicode[STIX]{x1D70C}_{s}c_{s}$ is the acoustic impedance of the solid. (This expression is exact only for plane waves with normal incidence falling on an isotropic unbounded solid medium. Nevertheless, the agreement of our simulations with experimental results was equally good at low attenuation (where this boundary condition might have played a role) and at high enough viscosity (where this boundary condition is irrelevant since the wave is attenuated even before bouncing back to the solid).) The solution $\tilde{p}$ can then be reconstructed thanks to the linearity of the equations:

(4.13) $$\begin{eqnarray}\tilde{p}(r,\unicode[STIX]{x1D703},z)=\mathop{\sum }_{-\infty }^{+\infty }\tilde{p}_{l}\text{e}^{\text{i}l\unicode[STIX]{x1D703}},\end{eqnarray}$$

and the velocity field can be reconstructed in a similar fashion described in appendix A. Thus, the acoustic field calculation can be summed up as follows: we start by computing the incident wave using (4.6)–(4.8), then its Fourier transform as given by (4.5). We combine these data and boundary conditions (4.11), (4.12) with the wave equation (4.10) to obtain the acoustic field for each individual harmonic. Finally, the total field is reconstructed using (A 7), (A 8). The algorithms ensuring the azimuthal Fourier transform were validated carefully by comparing the acoustic field in droplets exposed to 6 MHz SAWs as computed by a direct finite element model and by the Fourier method (figures available in supplemental material).

4.1.2 Resulting acoustical field in the droplet

To the best of our knowledge, the current computation of the acoustic field in a 3-D sessile droplet involves frequencies an order of magnitude above the only other work published so far by Quintero & Simonetti (Reference Quintero and Simonetti2013). As a result, it is significantly different and we will dedicate a few lines to detail the key features of this field.

Figure 6.

Meridian cross-section of the acoustic field in a water droplet excited by a 20 MHz acoustic field. (a) Acoustic pressure $\tilde{p}$ . $\tilde{p}_{max}=40$  kPa. (b) Average Lagrangian density $\langle {\mathcal{L}}\rangle$ . $-0.27~\text{J}<\langle {\mathcal{L}}\rangle <0.020~\text{J}~\text{m}^{-3}$ . (c) Average energy density $\langle {\mathcal{E}}\rangle =\langle {\mathcal{K}}+{\mathcal{V}}\rangle$ . $\langle {\mathcal{E}}\rangle _{max}=0.31~\text{J}~\text{m}^{-3}$ . (d) Poynting vector. $\Vert \langle \unicode[STIX]{x1D72B}\rangle \Vert _{max}=300~\text{W}~\text{m}^{-2}$ . Droplet volume is 12.5  $\unicode[STIX]{x03BC}$ l, base diameter is 3.7 mm. The incident wave comes from the left with a vertical displacement of 10 pm.

In figure 6, we show the acoustic field in a sessile water droplet excited by a 20 MHz SAW. The acoustic pressure (6 a) appears with two caustics superimposed on a quasi-random background field. The incident wave is overwhelmed by the numerous reflections on the droplet surface. The two caustics are much more pronounced than what is found in the two-dimensional analogue shown in Brunet et al. (Reference Brunet, Baudoin, Bou Matar and Zoueshtiagh2010), probably due to the increased ray convergence in three dimensions. These caustics show no symmetry axis within the $x$ – $z$ plane and are therefore structurally stable under the assumption that the droplet contact line is not pinned (Berry Reference Berry1976). This legitimates a posteriori the validity of using hemispherical geometry to model droplets slightly flattened by gravity (De Gennes, Brochard-Wyart & Quéré Reference De Gennes, Brochard-Wyart and Quéré2013). The Lagrangian of the acoustic field (figure 6 b) is mostly focused along the two caustics. Its large magnitude indicates that the acoustic field is not even remotely related to plane waves. The energy of the acoustic field (figure 6 c) clearly shows the predominance of the caustics. The angle of this concentration of energy matches the Rayleigh refraction angle and the symmetry shows the stability of this particular acoustic ray. Finally, the Poynting vector (figure 6 d) is similar to the energy density and the Lagrangian, except that its divergence is almost zero in weakly attenuating media, and it gives some orientation to the flow forcing that is consistent with the experiment.

Figure 7.

Meridian cross-section of the acoustic field in a 90 wt.% glycerol droplet excited by a 20 MHz acoustic field. (a) Acoustic pressure $\tilde{p}_{max}=16.5$  kPa. (b) Poynting vector. $\Vert \langle \unicode[STIX]{x1D72B}\rangle \Vert _{max}=25~\text{W}~\text{m}^{-2}$ Droplet volume is 12.5  $\unicode[STIX]{x03BC}$ l, base diameter is 4.0 mm. The incident wave comes from the left with an amplitude of 10 pm.

Figure 8.

Poynting vector in a 90 wt.% glycerol droplet excited by a 20 MHz acoustic field. $\Vert \langle \unicode[STIX]{x1D72B}\rangle \Vert _{max}=25~\text{W}~\text{m}^{-2}$ Droplet volume is 12.5  $\unicode[STIX]{x03BC}$ l, base diameter is 4.0 mm. The incident wave propagates along the $x$ axis (from left to right) with an amplitude of 10 pm.

In figure 7, we display the radiation patterns obtained in 90 wt.% glycerol droplets. At such high viscosity, the pressure field becomes less symmetrical. Indeed, the incident wave is attenuated faster and hence undergoes less reflections at the droplet surface. Similarly to the low viscosity case, the acoustic power flow is concentrated in the meridian $x{-}z$ plane as shown in figure 8. The ray diameter is approximately the diffraction limit $\unicode[STIX]{x1D706}/2$ (figure 7). The Rayleigh radiation angle appears more clearly at higher dissipation, and the wave pattern loses its front–back symmetry. Remarkably, the Poynting vector becomes completely asymmetrical and forces the flow only on the right side of the droplet. We will detail the consequences of this change on the resulting flow pattern in the next section.

Figure 9.

Surface (a,c,e,g) and volume (b,d,e,h) caustics deduced from geometrical acoustics. Surface caustics focus tangentially to the droplet surface whereas volume caustics focus orthogonally to the surface. The snapshots are taken after different acoustic path length: (a,b) 1.6 droplet radii, (c–h) 40 droplet radii; (g,h) depict the beam density. One unit volume is $r^{3}/10^{6}$ . The small mean free path caustics (a,c,e,g) follow periodical orbits along the droplet surface whereas large mean free path caustics (b,d,e,h) quickly show a more erratic behaviour. Animations of these trajectories are available in supplemental material.

Since the caustics play a major role in driving the droplet internal flow, we used geometrical acoustics to trace their origin. We simulated the propagation of multiple rays originating from the solid surface and propagating along the Rayleigh angle. These rays travel freely in the liquid volume but are reflected at the air–liquid and solid–liquid interfaces. We can then distinguish two types of rays depending on their mean free path. The rays with the shortest mean free path (figure 9 a,c,e) are guided by the droplet surface, whereas the rays with the longest mean free path (b,d,f) propagate in the droplet volume. The snapshots in figure 9(a,b) represent the wave after traveling approximately 1.6 droplet radii. At this stage, surface rays form a surface caustic tangential to the air–liquid interface, whereas volume rays form a volume caustic normal to the interface. As suggested by one referee, we allowed the wave to travel up to 40 times the droplet radius. Unexpectedly, the surface rays follow stable paths along the droplet surface (see figure 9 c,d) and maintain the same two caustics over a large number of periodic trajectories (orbits). Contrasting with the clean and periodic pattern of the surface rays, the volume rays yield a more chaotic arrangement. Interestingly, periodic echoes from sessile droplets have been observed before by Quintero & Simonetti (Reference Quintero and Simonetti2013), but were attributed to pseudo-Scholte waves propagating along the droplet contact line. The orbits observed here were not considered in Quintero et al. investigation and have a period compatible with the experimental echo. In figure 9(g,h), the density of rays faithfully reproduces the patterns observed in low-viscosity droplets. In the high-viscosity case, the acoustic rays attenuate before reaching the caustics at the back, which leads to the observed loss of symmetry.

4.2 Computation of the resulting flow

4.2.2 Streaming source spatial filtering

In the field of turbulence, energy is accumulated and transported by large flow patterns but is dissipated only by the smallest eddies. The gigantic difference of length scale all along this energy cascade is a major issue. A well-established method called LES allows computation of turbulence on relatively coarse grids that account for subgrid-scale (SGS) dynamics through a SGS model (see e.g. Deardorff Reference Deardorff1970; Pope Reference Pope2004; Bou-Zeid Reference Bou-Zeid2015). The case of acoustic streaming appears as a reverse situation wherein a large-scale flow (the acoustic streaming) emerges from small-scale fluctuations (the acoustic wave). The SSSF method presented in this section relies on the fact that the small-scale variations of the streaming source term ${\mathcal{F}}^{\ast }$ do not contribute to the flow since they are filtered by the fluid viscosity. Indeed, the acoustic streaming under investigation is slow and laminar, yielding a linear equation with a momentum source term. It is then interesting to consider the Fourier transform of the velocity field. We can match each wavenumber of the flow field with a (possibly null) forcing term in order to solve each equation independently. It is then straightforward to notice that higher wavenumbers are filtered out by the Laplacian operator of the viscosity (decreasing in $1/k^{2}$ ). In this regard, this SSSF model differs significantly from LES: in the former, smaller scales are sources of momentum and dissipation happens at larger scales, whereas in the latter, smaller scales act as momentum sinks because large-scale flows dissipate little energy.

Such filtering enables us to use grid cell sizes for the flow computation larger than the acoustic wavelength and thus considerably reduces the computational requirements for the resolution of the flow problem. It is important to note that since we work in the small Reynolds number regime, the characteristic length scale of the flow is entirely dictated by the streaming source term and the boundary conditions (no additional scale emerges from the flow itself like in the case of turbulent flows). The filtered source term $\text{}\underline{{\mathcal{F}}^{\ast }}$ is obtained in the real space from the convolution product with the filtering function ${\mathcal{H}}(x,y,z)$ :

(4.14) $$\begin{eqnarray}\text{}\underline{{\mathcal{F}}^{\ast }}=f\ast {\mathcal{H}},\end{eqnarray}$$

where the filtering function ${\mathcal{H}}(x,y,z)$ is defined from the filter transfer function $H(k_{x},k_{y},k_{z})$ according to the formula:

(4.15) $$\begin{eqnarray}{\mathcal{H}}(x,y,z)=\int _{{\mathcal{S}}}H(k_{x},k_{y},k_{z})\text{e}^{\text{i}k_{x}x+\text{i}k_{y}y+\text{i}k_{z}z}\,\text{d}{\mathcal{S}},\end{eqnarray}$$

with ${\mathcal{S}}$ the reciprocal space, $H(k_{x},k_{y},k_{z})=1$ when $k_{x}^{2}+k_{y}^{2}+k_{z}^{2}<k_{c}^{2}$ and zero otherwise, and $k_{c}$ is the critical wavenumber of the filtered flow structures. We choose the critical wavenumber $k_{c}$ as half the acoustical wavenumber in the fluid at working frequency. Indeed, the acoustic forcing term is the product of two acoustic quantities, which halves the spatial period. The exact choice of $k_{c}=k/2$ is somewhat arbitrary provided $k_{c}$ is below $2k$ and larger than $2\unicode[STIX]{x03C0}/L$ , where $L$ is the typical scale of the feature to be observed. This assertion was validated for various $0<L<D/40$ , and results are provided in supplemental material.

Consequently, the equations solved with the SSSF method are simply:

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{i}\bar{v}_{i}=0, & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x2202}_{i}\bar{p^{\ast }}+\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{jj}^{2}\bar{v}_{i}+\text{}\underline{{\mathcal{F}}_{i}^{\ast }}=0. & \displaystyle\end{eqnarray}$$

In the next section, we will show that the flow patterns computed from the DNS and the SSSF method agree qualitatively and quantitatively. Hence, in the remaining part of the paper all the simulations at 20 MHz will be performed with the SSSF method to overcome hardware limitations.

5.1 Acoustic forcing term and flow pattern

The acoustic forcing term, given by (3.32) is proportional to the Poynting vector. This means that the knowledge of the acoustic power flow is tantamount to the knowledge of the forcing term. The Poynting vector in water and glycerol droplet is shown in figures 6 and 7. As stated in § 4.1.2, it is focused on small regions of the droplets corresponding to the caustics, illustrated in figure 9. As we compare the resulting flow pattern in figure 12 to the forcing term shown in figures 6 and 7, we notice that these caustics act as momentum source points to generate the flow. For instance, in the glycerol droplet, the forcing terms act on only one side of the droplet and push the flow towards the rear of the droplet. In the case of water, the momentum source terms are more symmetrical and push the fluid in the two opposite directions. Each individual momentum source term results in two eddies, forming the four-swirl pattern. This is particularly visible in figure 3.

Interestingly, these caustics can be easily constructed from geometrical acoustics. This means that the flow can be, at least qualitatively, predicted from simple geometrical arguments. This assertion must be mitigated by the important role played by the viscosity and the attenuation of sound in the system.

5.2 The four-to-two eddies transition

Since the flow pattern heavily relies on the formation of caustics, it is strongly influenced by the ratio between the droplet diameter and the acoustic wave attenuation length in the fluid:

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}=\frac{D}{L_{a}}=\frac{D\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D708}b}{c_{0}^{3}}.\end{eqnarray}$$

For water sessile droplets with a 3.7 mm diameter, $\unicode[STIX]{x1D6EC}\simeq 0.07$ whereas for 90 % glycerol droplets $\unicode[STIX]{x1D6EC}\simeq 2.80$ . According to the numerical and experimental results, the progressive transition from two to four eddies is located within $0.29<\unicode[STIX]{x1D6EC}<1.0$ . In this regard, we understand why a transition of flow pattern happens for this range of viscosity. In glycerol droplets, the sound wave experiences very few reflections before fading whereas in water it should bounce at least sixteen times at the droplet interface. Hence, three regimes naturally emerge: $\unicode[STIX]{x1D6EC}\ll 1$ , $\unicode[STIX]{x1D6EC}\simeq 1$ and $\unicode[STIX]{x1D6EC}\gg 1$ . In the first one, the stable nature of the surface caustics favours the four eddies flow pattern as observed in water droplets. The intermediate regime happens for large glycerol concentrations. In this case, both surface and volume caustics are present as attested by figure 7(b). Thus, the flow in very viscous droplets is dominated by both the volume and surface caustics. Due to the large attenuation level, the surface caustic does not complete an orbit and the forcing is limited to only one side of the droplet, resulting in the two eddies flow pattern. Decreasing further $\unicode[STIX]{x1D6EC}$ was achieved previously by Beyssen et al. (Reference Beyssen, Le Brizoual, Elmazria, Alnot, Perry and Maillet2006) who worked at 40 MHz frequency with 90 wt.% glycerol mixtures. At such extreme regime, even the volume caustic vanishes and the flow pattern turns into a single vortex with a horizontal vorticity axis.