2.1 Linear oscillatory field at first order

The fluid velocity and pressure $\boldsymbol{u}(\boldsymbol{x},t)$ and pressure $p(\boldsymbol{x},t)$ solve the Navier–Stokes equation

(2.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\\ \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$ being the kinematic viscosity. With the perturbation expansion (2.1), and using a time-periodic expression for the first-order field

(2.4a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{(1)}(\boldsymbol{x},t)=\boldsymbol{u}_{1}(\boldsymbol{x})\text{e}^{\text{i}\unicode[STIX]{x1D714}t},\quad p^{(1)}(\boldsymbol{x},t)=p_{1}(\boldsymbol{x})\text{e}^{\text{i}\unicode[STIX]{x1D714}t}, & & \displaystyle\end{eqnarray}$$

one obtains for the momentum equation at $O(\unicode[STIX]{x1D700})$

(2.5) $$\begin{eqnarray}\displaystyle \text{i}\unicode[STIX]{x1D714}\boldsymbol{u}_{1}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p_{1}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{1}. & & \displaystyle\end{eqnarray}$$

Equation (2.5) is solved using a Helmholtz decomposition

(2.6a-e ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{1}=\boldsymbol{u}_{\unicode[STIX]{x1D711}}+\boldsymbol{u}_{A},\quad \boldsymbol{u}_{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D711},\quad \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}=0,\quad \boldsymbol{u}_{A}=\unicode[STIX]{x1D735}\times A,\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }A=0. & & \displaystyle\end{eqnarray}$$

Note that the generality of the Nyborg formulation (Nyborg Reference Nyborg1958) is premised on finding the velocity in terms of values of the potential component $\boldsymbol{u}_{\unicode[STIX]{x1D711}}$ and its derivative at the boundary. The vortical part $\boldsymbol{u}_{A}$ satisfies

(2.7a-c ) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D6FB}^{2}+H^{2})\boldsymbol{u}_{A},\quad H^{2}=-\text{i}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D708},\quad H=(1-\text{i})\sqrt{\frac{\unicode[STIX]{x1D714}}{2\unicode[STIX]{x1D708}}}=(1-\text{i})\unicode[STIX]{x1D6FD}, & & \displaystyle\end{eqnarray}$$

where the sign of $H$ was chosen for decaying solution of $\exp (-\text{i}Hz)$ . The solution has the typical structure of Stokes boundary layer for an oscillatory outer driving flow $\boldsymbol{u}_{\unicode[STIX]{x1D711}}$ (of order ${\sim}U$ and varying in a large scale ${\sim}L$ ) near a wall with boundary layer thickness $\unicode[STIX]{x1D6FF}=1/\unicode[STIX]{x1D6FD}\ll L$ . We seek solution $(u,w)$ in an axisymmetric geometry. Accordingly, the solution of the radial component $u_{A}$ is straightforward and chosen to ensure a zero tangential velocity countering $u_{\unicode[STIX]{x1D711}}$ as

(2.8) $$\begin{eqnarray}\displaystyle u_{A}=-u_{\unicode[STIX]{x1D711}}\text{e}^{-\text{i}Hz}. & & \displaystyle\end{eqnarray}$$

It satisfies (2.7). The axial component $w_{A}$ is chosen as

(2.9a,b ) $$\begin{eqnarray}\displaystyle w_{A}=-\frac{\unicode[STIX]{x1D6E4}}{\text{i}h}\text{e}^{-\text{i}Hz},\quad \unicode[STIX]{x1D6E4}=-\frac{\unicode[STIX]{x2202}w_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}z}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(ru_{\unicode[STIX]{x1D711}}). & & \displaystyle\end{eqnarray}$$

Due to (2.6), $w_{\unicode[STIX]{x1D711}}$ and $\unicode[STIX]{x1D6E4}$ are harmonic, and therefore (2.9) satisfies (2.7). We note that

(2.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\unicode[STIX]{x1D711}}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(ru_{\unicode[STIX]{x1D711}})+\frac{\unicode[STIX]{x2202}w_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}z}=0, & & \displaystyle\end{eqnarray}$$

to obtain

(2.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{A}=-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(ru_{\unicode[STIX]{x1D711}})\text{e}^{-\text{i}Hz}-\frac{\unicode[STIX]{x2202}w_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}z}\text{e}^{-\text{i}Hz}+\frac{1}{\text{i}H}\frac{\unicode[STIX]{x2202}^{2}w_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}z^{2}}\text{e}^{-\text{i}Hz}\approx 0, & & \displaystyle\end{eqnarray}$$

the last term being higher order in the small quantity $\unicode[STIX]{x1D6FF}/L$ compared to the first two terms and therefore was neglected here as well as below (effectively $\unicode[STIX]{x1D6E4}$ being treated as approximately a constant). The velocity $\boldsymbol{u}_{\unicode[STIX]{x1D711}}+\boldsymbol{u}_{A}$ however does not satisfy the zero normal velocity condition at $z=0$ due to $w_{A}$ . Correcting for that, a modified total first-order velocity is found as

(2.12a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=\boldsymbol{u}_{\unicode[STIX]{x1D711}}+\boldsymbol{u}_{A}+w_{c}\hat{\boldsymbol{e}}_{z},\quad w_{c}=\unicode[STIX]{x1D6E4}/\text{i}H, & & \displaystyle\end{eqnarray}$$

keeping in mind that the non-decaying $w_{c}$ is only meaningful while considering the velocity field in the small boundary layer region. Therefore, one obtains the time-periodic first-order velocity field (superscript $^{(1)}$ indicates the corresponding term with periodic time dependence included):

(2.13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}\displaystyle u^{(1)}=u_{\unicode[STIX]{x1D711}}^{(1)}+u_{A}^{(1)}\; & =\; & u_{\unicode[STIX]{x1D711}}[\cos \unicode[STIX]{x1D714}t-\text{e}^{-\unicode[STIX]{x1D6FD}z}\cos (\unicode[STIX]{x1D714}t-\unicode[STIX]{x1D6FD}z)],\\ \displaystyle w^{(1)}\; & =\; & w_{\unicode[STIX]{x1D711}}^{(1)}+w_{c}^{(1)}+w_{A}^{(1)}=w_{\unicode[STIX]{x1D711}}\cos \unicode[STIX]{x1D714}t\\ \; & \; & \displaystyle +\,\frac{\unicode[STIX]{x1D6E4}}{\sqrt{2}\unicode[STIX]{x1D6FD}}[\cos (\unicode[STIX]{x1D714}t-\unicode[STIX]{x03C0}/4)-\text{e}^{-\unicode[STIX]{x1D6FD}z}\cos (\unicode[STIX]{x1D714}t-\unicode[STIX]{x1D6FD}z-\unicode[STIX]{x03C0}/4)].\end{array}\right\} & & \displaystyle\end{eqnarray}$$

This expression is consistent with Nyborg (Reference Nyborg1958).

2.2 Streaming at second order

At the second order, one observes the forcing term $\boldsymbol{F}$ on the right-hand side of the average Stokes equation (2.2). As was noted by Nyborg, quadratic product of the irrotational part $\boldsymbol{u}_{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$ does not contribute to streaming and can balance the pressure gradient term as in Bernoulli term

(2.14) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}p_{2}=\boldsymbol{u}_{\unicode[STIX]{x1D711}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{\unicode[STIX]{x1D711}}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D735}(|\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}|^{2}), & & \displaystyle\end{eqnarray}$$

reducing (2.2) into

(2.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\langle \boldsymbol{u}^{(2)}\rangle _{t}=\boldsymbol{F}=\unicode[STIX]{x1D70C}\langle (\boldsymbol{u}_{A}^{(1)}+w_{c}^{(1)}\hat{e}_{z})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{\unicode[STIX]{x1D711}}^{(1)}+\boldsymbol{u}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\boldsymbol{u}_{A}^{(1)}+w_{c}^{(1)}\hat{e}_{z})\rangle _{t}. & & \displaystyle\end{eqnarray}$$

Therefore, the forcing term in the $r$ -direction becomes

(2.16) $$\begin{eqnarray}\displaystyle F_{r} & = & \displaystyle \unicode[STIX]{x1D70C}\langle \boldsymbol{u}_{A}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{\unicode[STIX]{x1D711}}^{(1)}+w_{c}^{(1)}\hat{e}_{z}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{\unicode[STIX]{x1D711}}^{(1)}+\boldsymbol{u}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{A}^{(1)}\rangle _{t}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D70C}\langle \!u_{A}^{(1)}\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}^{(1)}+w_{A}^{(1)}\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}+w_{c}^{(1)}\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}\nonumber\\ \displaystyle & & \displaystyle +\,u_{\unicode[STIX]{x1D711}}^{(1)}\unicode[STIX]{x2202}_{r}u_{A}^{(1)}+u_{A}^{(1)}\unicode[STIX]{x2202}_{r}u_{A}^{(1)}+w_{\unicode[STIX]{x1D711}}^{(1)}\unicode[STIX]{x2202}_{z}u_{A}^{(1)}+w_{c}^{(1)}\unicode[STIX]{x2202}_{z}u_{A}^{(1)}+w_{A}^{(1)}\unicode[STIX]{x2202}_{z}u_{A}^{(1)}\rangle _{t}.\end{eqnarray}$$

Noting the order of terms

(2.17) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}u_{\unicode[STIX]{x1D711}}\sim u_{A}\sim U,\\ w_{A}\sim w_{c}\sim \unicode[STIX]{x1D6FF}U/L,\\ \unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}\sim \unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}\sim \unicode[STIX]{x2202}_{r}u_{A}\sim U/L,\\ \unicode[STIX]{x2202}_{z}u_{A}\sim U/\unicode[STIX]{x1D6FF},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

and neglecting higher order in $\unicode[STIX]{x1D6FF}/L$ , we obtain

(2.18) $$\begin{eqnarray}\displaystyle F_{r}=\unicode[STIX]{x1D70C}\langle u_{A}^{(1)}\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}^{(1)}+u_{\unicode[STIX]{x1D711}}^{(1)}\unicode[STIX]{x2202}_{r}u_{A}^{(1)}+u_{A}^{(1)}\unicode[STIX]{x2202}_{r}u_{A}^{(1)}+w^{(1)}\unicode[STIX]{x2202}_{z}u_{A}^{(1)}\rangle _{t}. & & \displaystyle\end{eqnarray}$$

Each term on the right-hand side is $O(U^{2}/L)$ . In an effort to express each term in terms of $u_{\unicode[STIX]{x1D711}}$ and $\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}^{(1)}$ , following Nyborg, we express the odd term $w_{\unicode[STIX]{x1D711}}$ in (2.13) as

(2.19) $$\begin{eqnarray}\displaystyle w_{\unicode[STIX]{x1D711}}\approx z\left(\frac{\unicode[STIX]{x2202}w_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}z}\right)_{z=0}=-\frac{z}{r}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(ru_{\unicode[STIX]{x1D711}})\right]_{z=0}. & & \displaystyle\end{eqnarray}$$

Note that with (2.19) $w_{\unicode[STIX]{x1D711}}\sim \unicode[STIX]{x1D6FF}U/L$ and therefore $w_{\unicode[STIX]{x1D711}}^{(1)}\unicode[STIX]{x2202}_{z}u_{A}^{(1)}$ in (2.18) is also $O(U^{2}/L)$ consistent with the other terms there. After substituting (2.13) in (2.18) averaging over time, one obtains

(2.20) $$\begin{eqnarray}\displaystyle F_{r} & = & \displaystyle \frac{\unicode[STIX]{x1D70C}}{2} (\!u_{\unicode[STIX]{x1D711}}\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}(\text{e}^{-2\unicode[STIX]{x1D6FD}z}-2\text{e}^{-\unicode[STIX]{x1D6FD}z}\cos \unicode[STIX]{x1D6FD}z)\nonumber\\ \displaystyle & & \displaystyle -\,u_{\unicode[STIX]{x1D711}}(1/r)\unicode[STIX]{x2202}_{r}(ru_{\unicode[STIX]{x1D711}})\text{e}^{-\unicode[STIX]{x1D6FD}z}\{\unicode[STIX]{x1D6FD}z(\cos \unicode[STIX]{x1D6FD}z+\sin \unicode[STIX]{x1D6FD}z)-\sin \unicode[STIX]{x1D6FD}z\}\!).\end{eqnarray}$$

In the governing equation (2.15) in the second order, we note that the vertical ( $z$ ) derivative is larger than the transverse ( $r$ ) derivative by $O(L/\unicode[STIX]{x1D6FF})$ to obtain

(2.21) $$\begin{eqnarray}\displaystyle \left\langle \frac{\unicode[STIX]{x2202}^{2}u^{(2)}}{\unicode[STIX]{x2202}^{2}z}\right\rangle _{t}=-\frac{F_{r}}{\unicode[STIX]{x1D707}}. & & \displaystyle\end{eqnarray}$$

Integrating and noting (2.9), we get

(2.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \langle u^{(2)}\rangle _{t}=\frac{1}{\unicode[STIX]{x1D714}}\left(\frac{1}{2}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}(u_{\unicode[STIX]{x1D6FC}}-u_{\unicode[STIX]{x1D6FD}})-\frac{u_{\unicode[STIX]{x1D711}}^{2}}{r}u_{\unicode[STIX]{x1D6FD}}\right),\quad \text{where}\\ \displaystyle u_{\unicode[STIX]{x1D6FC}}=\frac{1}{4}\text{e}^{-2\unicode[STIX]{x1D6FD}z}+\text{e}^{-\unicode[STIX]{x1D6FD}z}\sin \unicode[STIX]{x1D6FD}z-\frac{1}{4},\\ \displaystyle u_{\unicode[STIX]{x1D6FD}}=\frac{1}{2}\unicode[STIX]{x1D6FD}z\text{e}^{-\unicode[STIX]{x1D6FD}z}(\cos \unicode[STIX]{x1D6FD}z-\sin \unicode[STIX]{x1D6FD}z)-\text{e}^{-\unicode[STIX]{x1D6FD}z}\left(\sin \unicode[STIX]{x1D6FD}z+\frac{1}{2}\cos \unicode[STIX]{x1D6FD}z\right)+\frac{1}{2},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

as was also found by Nyborg (Reference Nyborg1958).

The vertical component of streaming velocity $\langle w^{(2)}\rangle _{t}$ – not found in the literature but critical for plotting streamlines – is obtained by using the equation of mass conservation and taking into account the no-slip condition on the rigid wall:

(2.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle w^{(2)}\rangle _{t}=\int _{0}^{z}-\frac{1}{r}\frac{\unicode[STIX]{x2202}(r\langle u^{(2)}\rangle )}{\unicode[STIX]{x2202}r}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad \times \left(-\frac{1}{8\unicode[STIX]{x1D6FD}}\text{e}^{-2\unicode[STIX]{x1D6FD}z}-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-\unicode[STIX]{x1D6FD}z}(6\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+8\sin \unicode[STIX]{x1D6FD}z+14\cos \unicode[STIX]{x1D6FD}z)-\frac{7}{4}z+\frac{29}{8\unicode[STIX]{x1D6FD}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad \times \left(-\frac{1}{8\unicode[STIX]{x1D6FD}}\text{e}^{-2\unicode[STIX]{x1D6FD}z}-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-\unicode[STIX]{x1D6FD}z}(2\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+4\sin \unicode[STIX]{x1D6FD}z+6\cos \unicode[STIX]{x1D6FD}z)-\frac{3}{4}z+\frac{13}{8\unicode[STIX]{x1D6FD}}\right).\qquad \quad\end{eqnarray}$$

Here the $z$ -dependence of $\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}/\unicode[STIX]{x2202}r$ and $\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}/\unicode[STIX]{x2202}r^{2}$ is neglected in the boundary layer similar to what was assumed in (2.20). The acoustic streaming velocity field is therefore known in terms of the outer irrotational velocity field. One can compute the shear stress on the wall (Nyborg Reference Nyborg1958) as

(2.24) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{wall}=\left.\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\right|_{z=0}=\left.\frac{\unicode[STIX]{x1D70C}_{0}}{4\unicode[STIX]{x1D6FD}}u_{\unicode[STIX]{x1D711}}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}r}\right|_{z=0}. & & \displaystyle\end{eqnarray}$$

2.3 Lagrangian streaming velocity

For comparison with the average path of a tracer particle, it was noted before that one would also need to find the Lagrangian or particle-averaged streaming velocity (Nyborg Reference Nyborg1958) which contains an additional Stokes drift term

(2.25) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{u}_{T}\rangle _{t}=\left\langle \left(\int \boldsymbol{u}^{(1)}\,\text{d}t\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{(1)}\right\rangle _{t}, & & \displaystyle\end{eqnarray}$$

expressed in terms of the first-order velocity given in (2.13). For the radial component, one can obtain

(2.26) $$\begin{eqnarray}\displaystyle \langle u_{T}^{(1)}\rangle _{t} & = & \displaystyle \left\langle \left(\int u^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{r}u^{(1)}\right\rangle _{t}+\left\langle \left(\int w_{\unicode[STIX]{x1D711}}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}\right\rangle _{t}+\left\langle \left(\int w_{A}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}\right\rangle _{t}\nonumber\\ \displaystyle & & \displaystyle +\,\left\langle \left(\int w_{c}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}\right\rangle _{t}+\left\langle \left(\int w_{\unicode[STIX]{x1D711}}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{A}^{(1)}\right\rangle _{t}+\left\langle \left(\int w_{A}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{A}^{(1)}\right\rangle _{t}\nonumber\\ \displaystyle & & \displaystyle +\,\left\langle \left(\int w_{c}^{(1)}\,\text{d}t\right)\unicode[STIX]{x2202}_{z}u_{A}^{(1)}\right\rangle _{t}.\end{eqnarray}$$

Using (2.19) and (2.17), we realize that the second, third and fourth terms involving $\unicode[STIX]{x2202}_{z}u_{\unicode[STIX]{x1D711}}^{(1)}$ are higher order in the small quantity $\unicode[STIX]{x1D6FF}/L$ than the others. The first term is identically zero. We obtain

(2.27) $$\begin{eqnarray}\displaystyle \langle u_{T}\rangle _{t} & = & \displaystyle \frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{1}{2}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}+\frac{u_{\unicode[STIX]{x1D711}}^{2}}{r}\right)u_{\unicode[STIX]{x1D700}},\nonumber\\ \displaystyle & & \displaystyle \quad \text{where }u_{\unicode[STIX]{x1D700}}=\unicode[STIX]{x1D6FD}z\text{e}^{-\unicode[STIX]{x1D6FD}z}(\cos \unicode[STIX]{x1D6FD}z-\sin \unicode[STIX]{x1D6FD}z)-\text{e}^{-\unicode[STIX]{x1D6FD}z}\cos \unicode[STIX]{x1D6FD}z+\text{e}^{-2\unicode[STIX]{x1D6FD}z}.\end{eqnarray}$$

Nyborg (Reference Nyborg1958) stressed that the vertical component is much smaller. However, it is critical for computing the streamlines and can be found by using continuity:

(2.28) $$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!\langle w_{T}\rangle _{t}=\int _{0}^{z}-\frac{1}{r}\frac{\unicode[STIX]{x2202}(r\langle u_{T}\rangle )}{\unicode[STIX]{x2202}r}\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{4\unicode[STIX]{x1D714}\unicode[STIX]{x1D6FD}}\left(\frac{1}{2}\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}+\frac{3}{2r}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}\right)(-2\unicode[STIX]{x1D6FD}z\text{e}^{-\unicode[STIX]{x1D6FD}z}\sin \unicode[STIX]{x1D6FD}z-2\text{e}^{-\unicode[STIX]{x1D6FD}z}\cos \unicode[STIX]{x1D6FD}z+\text{e}^{-2\unicode[STIX]{x1D6FD}z}+1).\qquad \qquad\end{eqnarray}$$

Note that Raney, Corelli & Westervelt (Reference Raney, Corelli and Westervelt1954) found that adding this correction improved the comparison of theory with experimental observation of streaming near a cylinder.

2.4 Potential velocity $u_{\unicode[STIX]{x1D711}}$ due to oscillating bubble above a rigid surface

The problem of a bubble of initial radius $R_{0}$ near a plane rigid wall at a distance $h$ from the wall is equivalent to two bubbles of the same initial radius separated by a distance of $2h$ ; it satisfies the impermeability condition in the plane of the rigid wall.

The velocity potential $\unicode[STIX]{x1D719}$ of the fluid around the microbubble is

(2.29) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}=-\left(\frac{1}{S_{1}}+\frac{1}{S_{2}}\right){\dot{R}}R^{2}, & & \displaystyle\end{eqnarray}$$

where $S_{1}$ and $S_{2}$ are the distances from the centre of the real and image microbubbles to the desired location in the fluid and ${\dot{R}}$ and $R$ are the velocity and instantaneous radius of the microbubble. The radial and vertical components of the irrotational velocity $\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$ are

(2.30) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u^{(\unicode[STIX]{x1D711})}=u_{\unicode[STIX]{x1D711}}\cos \unicode[STIX]{x1D714}t={\dot{R}}R^{2}\left(\frac{r}{(r^{2}+(z-h)^{2})^{3/2}}+\frac{r}{(r^{2}+(z+h)^{2})^{3/2}}\right),\\ \displaystyle w^{(\unicode[STIX]{x1D711})}=w_{\unicode[STIX]{x1D711}}\cos \unicode[STIX]{x1D714}t={\dot{R}}R^{2}\left(\frac{z-h}{(r^{2}+(z-h)^{2})^{3/2}}+\frac{z+h}{(r^{2}+(z+h)^{2})^{3/2}}\right).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Here $u_{\unicode[STIX]{x1D711}}$ and $w_{\unicode[STIX]{x1D711}}$ are time-independent parts of the potential velocity components in radial and vertical directions and $r$ and $z$ are the radial and vertical coordinates. The instantaneous bubble radius $R$ and velocity ${\dot{R}}$ are described by the Rayleigh–Plesset-type equation

(2.31) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}\left(R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}\right)=P_{b}-P_{0}+P_{ex}\sin \unicode[STIX]{x1D714}t,\\ \displaystyle P_{b}=P_{g_{0}}\left(\frac{R_{0}}{R}\right)^{3\unicode[STIX]{x1D705}}-\frac{4\unicode[STIX]{x1D707}{\dot{R}}}{R}-\frac{2\unicode[STIX]{x1D6FE}}{R}-P_{sc}(h,t),\\ \displaystyle P_{sc}(h,t)=\unicode[STIX]{x1D70C}\frac{R}{2h}(R\ddot{R}+2{\dot{R}}^{2}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Here $P_{b}$ is the fluid pressure adjacent to the microbubble, $P_{0}$ is the ambient pressure (100 kPa) and $P_{ex}$ is the ultrasound excitation amplitude. Parameter $\unicode[STIX]{x1D6FE}$ is the gas–fluid surface tension, $P_{g_{0}}$ is the initial gas pressure inside the microbubble and $\unicode[STIX]{x1D705}$ is the polytropic constant. Note that the impermeability at the wall is accounted for by the effect of the wall being considered as a pressure $P_{sc}(h,t)$ scattered from the image bubble located at a distance $2h$ from the real bubble.

2.5 Linearized bubble dynamics and streaming fields

Assuming that the microbubble is pulsating with a small amplitude $R=R_{0}(1+x)$ , equation (2.31) is linearized in $x$ . Non-dimensionalizing $t$ and $P_{ex}$ as $t^{\ast }=t\unicode[STIX]{x1D714}$ and $P_{ex}^{\ast }=P_{ex}/P_{0}$ results in an equation of damped harmonic oscillator

(2.32) $$\begin{eqnarray}\displaystyle \ddot{x}+\left(\frac{4}{Re_{b}(1+R_{0}/2h)}\right){\dot{x}}+\left(\frac{3\unicode[STIX]{x1D705}Eu_{b}}{1+R_{0}/2h}+\frac{6\unicode[STIX]{x1D705}-2}{We_{b}(1+R_{0}/2h)}\right)x=\frac{P_{ex}^{\ast }Eu_{b}\sin (\unicode[STIX]{x1D714}t)}{1+R_{0}/2h}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where

(2.33a-c ) $$\begin{eqnarray}\displaystyle Re_{b}=\frac{\unicode[STIX]{x1D70C}R_{0}^{2}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D707}},\quad Eu_{b}=\frac{P_{0}}{\unicode[STIX]{x1D70C}R_{0}^{2}\unicode[STIX]{x1D714}^{2}},\quad We_{b}=\frac{\unicode[STIX]{x1D70C}R_{0}^{3}\unicode[STIX]{x1D714}^{2}}{\unicode[STIX]{x1D6FE}}, & & \displaystyle\end{eqnarray}$$

are the characteristic Reynolds (we have taken $R_{0}\unicode[STIX]{x1D714}$ as the velocity scale), Euler and Weber numbers. And

(2.34) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x1D714}_{0}^{2}}{\unicode[STIX]{x1D714}^{2}}=\frac{3\unicode[STIX]{x1D705}Eu_{b}+(6\unicode[STIX]{x1D705}-2)/We_{b}}{(1+R_{0}/2h)},\\ \displaystyle \unicode[STIX]{x1D6FF}_{t}=\frac{4\unicode[STIX]{x1D714}}{Re_{b}(1+R_{0}/2h)\unicode[STIX]{x1D714}_{0}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{0}$ and $\unicode[STIX]{x1D6FF}_{t}$ are the circular natural frequency and the damping term of the microbubble. Note that the natural frequency is modified due to the presence of the wall by the factor $(1+R_{0}/2h)^{-1}$ . The analytical solution of equation (2.32) in the steady region is

(2.35) $$\begin{eqnarray}\displaystyle x=\frac{Eu_{b}P_{ex}^{\ast }\sin (\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D719})}{Z_{m}(1+R_{0}/2h)}=\unicode[STIX]{x1D709}_{m}\sin (\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D719}), & & \displaystyle\end{eqnarray}$$

where $Z_{m}=\sqrt{(\unicode[STIX]{x1D6FF}_{t}\unicode[STIX]{x1D714}_{0}/\unicode[STIX]{x1D714})^{2}+(1/\unicode[STIX]{x1D714}^{4})(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D714}^{2})^{2}}$ is the non-dimensional absolute value of impedance and $\unicode[STIX]{x1D719}$ is the oscillation phase that can be absorbed by redefining time as $t^{\prime }=t^{\ast }+\unicode[STIX]{x1D719}/\unicode[STIX]{x1D714}$ . As a result the irrotational velocity $u_{\unicode[STIX]{x1D711}}$ from equation (2.30), when non-dimensionalized by $R_{0}\unicode[STIX]{x1D714}$ , becomes

(2.36) $$\begin{eqnarray}\displaystyle u_{\unicode[STIX]{x1D711}}^{\ast } & = & \displaystyle \frac{\unicode[STIX]{x1D709}_{m}}{r^{\ast 2}}\left.\left(\frac{1}{\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}+\frac{1}{\left(1+\left(\displaystyle \frac{z^{\ast }+h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}\right)\right|_{z\ast =0}\nonumber\\ \displaystyle & = & \displaystyle \frac{2\unicode[STIX]{x1D709}_{m}}{r^{\ast 2}\left(1+\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}},\end{eqnarray}$$

(2.37) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{\ast }}{\unicode[STIX]{x2202}r^{\ast }} & = & \displaystyle \frac{\unicode[STIX]{x1D709}_{m}}{r^{\ast 3}\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}\left.\left(\frac{3}{\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)}-2\right)\right|_{z=0}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D709}_{m}}{r^{\ast 3}\left(1+\left(\displaystyle \frac{z^{\ast }+h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}\left.\left(\frac{3}{\left(\displaystyle \frac{z^{\ast }+h^{\ast }}{r^{\ast }}\right)^{2}\left(1+\left(\displaystyle \frac{z^{\ast }+h^{\ast }}{r^{\ast }}\right)^{2}\right)}-2\right)\right|_{z=0}\nonumber\\ \displaystyle & = & \displaystyle \frac{2\unicode[STIX]{x1D709}_{m}}{r^{\ast 3}\left(1+\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}\left(\frac{3}{\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\left(1+\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\right)}-2\right),\end{eqnarray}$$

with a similar expression for $w_{\unicode[STIX]{x1D711}}^{\ast }$ . Note that in conformity with the approximation that $u_{\unicode[STIX]{x1D711}}^{\ast }$ varies in a larger scale ( $h$ and $R_{0}$ ) compared to the boundary layer thickness $\unicode[STIX]{x1D6FF}$ , the expressions were evaluated at the wall $z=0$ . Correspondingly, shear stress (2.24) non-dimensionalized with respect to $P_{0}$ becomes

(2.38) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}^{\ast }=\frac{2\unicode[STIX]{x1D709}_{m}^{2}r^{\ast }}{Eu_{b}(2Re_{b})^{1/2}h^{\ast 5}}\frac{(1-2r^{\ast 2}/h^{\ast 2})}{(1+r^{\ast 2}/h^{\ast 2})^{4}}. & & \displaystyle\end{eqnarray}$$

2.6 Effects of translation of microbubble

An oscillating microbubble translates towards the wall due to Bjerknes force as has been carefully analysed by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2014). The velocity potential is modified by an additional dipole term

(2.39) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}=-\left(\frac{1}{S_{1}}+\frac{1}{S_{2}}\right){\dot{R}}R^{2}-b\left(\frac{h-z}{S_{1}^{3}}+\frac{h+z}{S_{2}^{3}}\right). & & \displaystyle\end{eqnarray}$$

Corresponding linearized non-dimensional potential velocity in the radial direction is also modified as

(2.40) $$\begin{eqnarray}\displaystyle u_{\unicode[STIX]{x1D711}}^{\ast } & = & \displaystyle \frac{\unicode[STIX]{x1D709}_{m}}{r^{\ast 2}}\left.\left(\displaystyle \frac{1}{\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}+\frac{1}{\left(1+\left(\displaystyle \frac{z^{\ast }+h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}\right)\right|_{z=0}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{b^{l}}{r^{\ast 4}}\left.\left(\frac{3(h^{\ast }-z^{\ast })}{\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)^{5/2}}+\frac{3(h^{\ast }+z^{\ast })}{\left(1+\left(\displaystyle \frac{z^{\ast }-h^{\ast }}{r^{\ast }}\right)^{2}\right)^{5/2}}\right)\right|_{z=0}\nonumber\\ \displaystyle & = & \displaystyle \frac{2\unicode[STIX]{x1D709}_{m}}{r^{\ast 2}\left(1+\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\right)^{3/2}}+\frac{6h^{\ast }b^{l}}{r^{\ast 4}\left(1+\left(\displaystyle \frac{h^{\ast }}{r^{\ast }}\right)^{2}\right)^{5/2}},\end{eqnarray}$$

where $b^{l}$ is given by (Doinikov & Bouakaz Reference Doinikov and Bouakaz2014)

(2.41) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle b^{l}=\left|\frac{\unicode[STIX]{x1D709}_{m}}{4h^{\ast 2}}\left(\frac{\unicode[STIX]{x1D6FC}^{3}+3\text{i}\unicode[STIX]{x1D6FC}^{2}-6\unicode[STIX]{x1D6FC}-6\text{i}}{-\unicode[STIX]{x1D6FC}^{3}-3\text{i}\unicode[STIX]{x1D6FC}^{2}+18\unicode[STIX]{x1D6FC}+18\text{i}}\right)\right|,\\ \displaystyle \unicode[STIX]{x1D6FC}=(1+\text{i})R_{0}\unicode[STIX]{x1D6FD}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

We note that the faster decay with $r^{\ast }$ characteristic of a dipole makes the translational part much smaller at larger radial distance. Also at smaller $r^{\ast }\ll 1$ , the translational part is smaller than the radial part by a factor of $3/4h^{\ast 3}$ . For separation distances of the microbubble from the wall $h\geqslant 2R_{0}$ , this factor is quite small, and therefore the effect of translation here has not been considered, except briefly. However, note that it can easily be accounted for using the expression (2.41).

2.7 Effects of the shell on a contrast microbubble

Microbubbles (size ${\sim}1{-}10~\unicode[STIX]{x03BC}\text{m}$ ) used for contrast-enhanced ultrasound imaging (Goldberg, Raichlen & Forsberg Reference Goldberg, Raichlen and Forsberg2001; Paul et al. Reference Paul, Nahire, Mallik and Sarkar2014) are coated by a monolayer of lipids, proteins or other surface-active molecules to stabilize them against premature dissolution due to gas diffusion (Katiyar, Sarkar & Jain Reference Katiyar, Sarkar and Jain2009; Sarkar, Katiyar & Jain Reference Sarkar, Katiyar and Jain2009). There have been numerous models of contrast microbubble coating (de Jong et al. Reference de Jong, Hoff, Skotland and Bom1992; Church Reference Church1995; Hoff, Sontum & Hovem Reference Hoff, Sontum and Hovem2000; Chatterjee & Sarkar Reference Chatterjee and Sarkar2003; Marmottant et al. Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005; Sarkar et al. Reference Sarkar, Shi, Chatterjee and Forsberg2005; Tsiglifis & Pelekasis Reference Tsiglifis and Pelekasis2008; Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010). We have recently shown that most of these models can be expressed as an interfacial rheological model with an effective surface tension $\unicode[STIX]{x1D6FE}(R)$ and an effective dilatational viscosity $\unicode[STIX]{x1D705}^{s}(R)$  (Katiyar & Sarkar Reference Katiyar and Sarkar2011). Here we have used a strain-softening model called the exponential elasticity model (Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010):

(2.42a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}(R)=\unicode[STIX]{x1D6FE}_{0}+\unicode[STIX]{x1D6FD}^{s}E^{s},\quad E^{s}=E_{0}^{s}\exp (-\unicode[STIX]{x1D6FC}^{s}\unicode[STIX]{x1D6FD}^{s}),\quad \unicode[STIX]{x1D6FD}^{s}=\left(\frac{R}{R_{E}}\right)^{2}-1,\quad \unicode[STIX]{x1D705}^{s}\text{ constant}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $E^{s}$ is the shell dilatational elasticity, $\unicode[STIX]{x1D6FD}^{s}$ is the change in area fraction from the equilibrium radius $R_{E}=R_{0}[1+(1-\sqrt{1+4\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D6FC}^{s}/E_{0}^{s}})/2\unicode[STIX]{x1D6FC}^{s}]^{-1/2}$ , and $\unicode[STIX]{x1D6FE}_{0}$ , $\unicode[STIX]{x1D6FC}^{s}$ and $E_{0}^{s}$ are material properties of the coating. We have demonstrated how these characteristic properties for a microbubble contrast agent can be measured through acoustic experiments (Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010, Reference Paul, Russakow, Rodgers, Sarkar, Cochran and Wheatley2013; Kumar & Sarkar Reference Kumar and Sarkar2015; Xia, Porter & Sarkar Reference Xia, Porter and Sarkar2015; Kumar & Sarkar Reference Kumar and Sarkar2016). For example, for contrast agent Sonazoid, they are $\unicode[STIX]{x1D6FE}_{0}=0.019~\text{N}~\text{m}^{-1}$ , $E_{0}^{s}=0.55~\text{N}~\text{m}^{-1}$ , $\unicode[STIX]{x1D6FC}^{s}=1.5$ and $\unicode[STIX]{x1D705}^{s}=1.2\times 10^{-8}~\text{N}~\text{s}~\text{m}^{-1}$ . As a result of the coating stresses, equations (2.34) become

(2.43) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \!\!\frac{\unicode[STIX]{x1D714}_{0}^{2}}{\unicode[STIX]{x1D714}^{2}}=\frac{1}{(1+R_{0}/2h)}\left(3\unicode[STIX]{x1D705}Eu_{b}+\frac{2(\sqrt{1+4\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D6FC}^{s}/E_{0}^{s}}/\unicode[STIX]{x1D6FC}^{s})(1+2\unicode[STIX]{x1D6FC}^{s}-\sqrt{1+4\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D6FC}^{s}/E_{0}^{s}})}{We_{b}^{s}}\right),\\ \displaystyle \unicode[STIX]{x1D6FF}_{t}=\frac{4\unicode[STIX]{x1D714}}{(1+R_{0}/2h)\unicode[STIX]{x1D714}_{0}}\left(\frac{1}{Re_{b}}+\frac{1}{Re_{b}^{s}}\right),\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $We_{b}^{s}=\unicode[STIX]{x1D70C}R_{0}^{3}\unicode[STIX]{x1D714}^{2}/E_{0}^{s}$ , $Re_{b}^{s}=\unicode[STIX]{x1D70C}R_{0}^{3}\unicode[STIX]{x1D714}/\unicode[STIX]{x1D705}^{s}$ , $\unicode[STIX]{x1D6FE}_{0}/E_{0}^{s}$ and $\unicode[STIX]{x1D6FC}^{s}$ are the non-dimensional numbers related to the coating parameters. We have shown that different models of coating often give rise to qualitatively similar behaviours (Kumar & Sarkar Reference Kumar and Sarkar2015; Xia et al. Reference Xia, Porter and Sarkar2015; Kumar & Sarkar Reference Kumar and Sarkar2016) even for nonlinear response. We expect that they would vary very little in their predictions of linear dynamics.

2.8 Circulation of the ring vortex in microstreaming flow

Microstreaming gives rise to an axisymmetric vortex ring near the wall (see figure 5). We compute the circulation of the vortex ring, a single quantitative scalar measure of the vortex, in the $r$ – $z$ plane using the Stokes theorem. However, it requires us to objectively define the spatial extent of the vortex. The problem of a quantitative and objective identification of vortices has been a major issue in experimentally measured velocity data in turbulent flows. In fact although vortices are ubiquitous, there is no universally accepted definition of a ‘vortex’ (Chakraborty, Balachandar & Adrian Reference Chakraborty, Balachandar and Adrian2005). The problem can be appreciated by noting that simple shear has vorticity, yet does not contain a vortex. A local method of vortex identification obtains a point function that classifies the point to be outside or inside a vortex. The criterion being completely defined by the local velocity gradient is Galilean invariant. There have been many local methods proposed in the literature and their relative merits, specifically their utility in identifying vortical structures in complex turbulent flow measurements, have been investigated in detail (Chakraborty et al. Reference Chakraborty, Balachandar and Adrian2005). Here, we apply the d2 method (discriminant of non-real eigenvalues of velocity gradient tensor) (Vollmers Reference Vollmers2001) principally for its ease of use in two-dimensional data. It is a Galilean-invariant method that distinguishes vortical structures from boundary layers and shear layers (Najjari & Plesniak Reference Najjari and Plesniak2016). A negative value of $d_{2}$ corresponds to the non-real eigenvalue of velocity gradient and indicates a vortex region:

(2.44) $$\begin{eqnarray}\displaystyle d_{2}=\left(\frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}+\frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\right)^{2}-4\left(\frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}-\frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}\right), & & \displaystyle\end{eqnarray}$$

where

(2.45) $$\begin{eqnarray}\displaystyle \left.\!\!\!\begin{array}{@{}rcl@{}}\displaystyle \frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}\; & =\; & \displaystyle \frac{1}{\unicode[STIX]{x1D714}}\left(\frac{1}{2}\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}(u_{\unicode[STIX]{x1D6FC}}-u_{\unicode[STIX]{x1D6FD}})-\left(\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{r\unicode[STIX]{x2202}r}-\frac{u_{\unicode[STIX]{x1D711}}^{2}}{r^{2}}\right)u_{\unicode[STIX]{x1D6FD}}\right),\\ \displaystyle \frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\; & =\; & \displaystyle \frac{1}{\unicode[STIX]{x1D714}}\left(\frac{1}{2}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}\unicode[STIX]{x1D6FD}\text{e}^{-\unicode[STIX]{x1D6FD}z}(\cos \unicode[STIX]{x1D6FD}z-2\sin \unicode[STIX]{x1D6FD}z+\unicode[STIX]{x1D6FD}z\cos \unicode[STIX]{x1D6FD}z-0.5\text{e}^{-\unicode[STIX]{x1D6FD}z})\right.\\ \; & \; & -\displaystyle \left.\frac{u_{\unicode[STIX]{x1D711}}^{2}}{r}(\unicode[STIX]{x1D6FD}\text{e}^{-\unicode[STIX]{x1D6FD}z}(\sin \unicode[STIX]{x1D6FD}z-\unicode[STIX]{x1D6FD}z\cos \unicode[STIX]{x1D6FD}z))\right),\\ \displaystyle \frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}\; & =\; & \displaystyle -\frac{1}{2\unicode[STIX]{x1D714}}\left(-\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}+\frac{1}{r}\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right)\\ \; & \; & \displaystyle \times \left(-\frac{1}{8\unicode[STIX]{x1D6FD}}\text{e}^{-2\unicode[STIX]{x1D6FD}z}-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-\unicode[STIX]{x1D6FD}z}\left(6\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+8\sin \unicode[STIX]{x1D6FD}z+14\cos \unicode[STIX]{x1D6FD}z\right)-\frac{7}{4}z+\frac{29}{8\unicode[STIX]{x1D6FD}}\right)\\ \displaystyle \; & \; & -\,\displaystyle \frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{\unicode[STIX]{x2202}^{3}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{3}}\right)\\ \; & \; & \displaystyle \times \left(-\frac{1}{8\unicode[STIX]{x1D6FD}}\text{e}^{-2\unicode[STIX]{x1D6FD}z}-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-\unicode[STIX]{x1D6FD}z}\left(2\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+4\sin \unicode[STIX]{x1D6FD}z+6\cos \unicode[STIX]{x1D6FD}z\right)-\frac{3}{4}z+\frac{13}{8\unicode[STIX]{x1D6FD}}\right),\\ \displaystyle \frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\; & =\; & \displaystyle -\frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}\right)\\ \; & \; & \displaystyle \times \left(\frac{1}{4}\text{e}^{-\unicode[STIX]{x1D6FD}z}(\text{e}^{-\unicode[STIX]{x1D6FD}z}+6\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+16\sin \unicode[STIX]{x1D6FD}z+6\cos \unicode[STIX]{x1D6FD}z-6\unicode[STIX]{x1D6FD}z\cos \unicode[STIX]{x1D6FD}z)-\frac{7}{4}\right)\\ \displaystyle \; & \; & -\,\displaystyle \frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right)\\ \; & \; & \displaystyle \times \left(\frac{1}{4}\text{e}^{-\unicode[STIX]{x1D6FD}z}(\text{e}^{-\unicode[STIX]{x1D6FD}z}+2\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+8\sin \unicode[STIX]{x1D6FD}z+2\cos \unicode[STIX]{x1D6FD}z-2\unicode[STIX]{x1D6FD}z\cos \unicode[STIX]{x1D6FD}z)-\frac{3}{4}\right).\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We compute the total circulation by detecting the vortex region using the d2 method, and then numerically integrating:

(2.46) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\mathop{\sum }_{i}\unicode[STIX]{x1D6EC}_{i}=\mathop{\sum }_{i}\unicode[STIX]{x1D6FA}_{i}\unicode[STIX]{x0394}r\unicode[STIX]{x0394}z=\mathop{\sum }_{i}\left(\frac{\unicode[STIX]{x2202}\langle w^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x2202}\langle u^{(2)}\rangle _{t}}{\unicode[STIX]{x2202}z}\right)\unicode[STIX]{x0394}r\unicode[STIX]{x0394}z, & & \displaystyle\end{eqnarray}$$

which can be non-dimensionalized by  $R_{0}^{2}\unicode[STIX]{x1D714}$ to obtain the non-dimensional circulation $\unicode[STIX]{x1D6EC}^{\ast }$ .

Figure 2. The radial pulsation of the free microbubble when $Re_{b}=48$ , $Eu_{b}=0.11,We_{b}=20,h/R_{0}=2$ .

3.1 Microstreaming due to free bubble oscillation

Although the analysis delineated here is valid for different bubble radius and excitation parameters, we are primarily interested in microbubble-based contrast-enhanced ultrasound imaging and drug delivery applications with $R_{0}$ of the order of micrometres and $f_{ex}$ in the megahertz range. Specifically, we consider the case of $R_{0}=1.6~\unicode[STIX]{x03BC}\text{m}$ , the average radius of the contrast agent Sonazoid (Sontum et al. Reference Sontum, Ostensen, Dyrstad and Hoff1999; Sontum Reference Sontum2008), at an excitation frequency of 3 MHz. The driving amplitude was chosen to be small enough, $P_{ex}=70~\text{kPa}$ ( $P_{ex}^{\ast }=0.7$ ), to ensure linear oscillation, giving rise to $Re_{b}=48$ , $Eu_{b}=0.11$ and $We_{b}=20$ . Figure 2 shows the radial pulsation of the microbubble with non-dimensional time located at $h=2R_{0}$ obtained solving the full Rayleigh–Plesset equation (2.31) using a standard Matlab solver, ode15s. The eventual periodic oscillation with a non-dimensional amplitude $\unicode[STIX]{x1D700}=0,11$ seen here is also predicted by the linearized form (2.35) demonstrating that the excitation is small enough for linear dynamics, and the linear analysis is valid. Note that the bubble streaming Reynolds number $Re_{bs}=\unicode[STIX]{x1D700}^{2}(2Re_{b})^{1/2}=0.12$ . We also check the approximation made in neglecting the $z$ -variation of $u_{\unicode[STIX]{x1D711}}^{\ast }$ (varying in a scale $R_{0}$ ) inside the boundary layer.

Figure 3. The non-dimensional irrotational velocity at $z=0$ and $z=\unicode[STIX]{x1D6FF}$ for $h=1.5R_{0},h=2R_{0},h=3R_{0}$ when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $P_{ex}^{\ast }=0.7$ .

In figure 3, we plot $u_{\unicode[STIX]{x1D711}}^{\ast }$ at $z=0$ and at $z=\unicode[STIX]{x1D6FF}$ for three different separations $h=1.5R_{0}$ ,  $h=2R_{0}$ and $h=3R_{0}$ . We note that the variation is small for all cases, being largest for the smallest separation $h=1.5R_{0}$ . Therefore the approximation adopted here – $u_{\unicode[STIX]{x1D711}}^{\ast }(z=0)$ and $\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{\ast }(z=0)$ in (2.36) and (2.37) – following Nyborg (Reference Nyborg1958) is justified. The average streaming velocity (2.22) and (2.23) can be written in terms of the irrotational component of the first-order oscillating velocity $u_{\unicode[STIX]{x1D711}}^{\ast }$ . Figure 4 shows the streamlines of $u_{\unicode[STIX]{x1D711}}^{\ast }$ surrounding the microbubble at two different non-dimensional times corresponding to microbubble expansion and compression.

Figure 4. The irrotational velocity around the free microbubble when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20,h/R_{0}=2$ , $P_{ex}^{\ast }=0.7$ : (a) expansion; (b) compression.

We plot the streamlines of microstreaming velocities $\langle u^{(2)}\rangle _{t}$ and $\langle w^{(2)}\rangle _{t}$ according to equations (2.22) and (2.23) in figures 5(a) and 5(b) for two locations of the microbubble at $h=2R_{0}$ and $h=3R_{0}$ , respectively. This case is the same one discussed before for $\text{Re}_{b}=48,Eu_{b}=0.11,We_{b}=20$ . They show the cross-sectional plot of the axisymmetric ring vortex over the wall. Inside the vortex, the flow near the wall is directed radially outward, while it is directed inward beyond the vortex length. The vortex ring extends up to $r^{\ast }=1.41$ and $r^{\ast }=2.12$ defining $l_{vortex}$ along the wall for the cases $h=2R_{0}$ and $h=3R_{0}$ respectively. Length $l_{vortex}$ increases with $h^{\ast }=h/R_{0}$ .

To understand the generation of the vortex we examine (2.2) showing $\unicode[STIX]{x1D70C}\langle \boldsymbol{u}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{(1)}\rangle _{t}$ driving the Stokes equation for the streaming velocity. In particular, one sees that $F_{r}$ in (2.20) is primarily dominated by the first term proportional to $u_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}^{\ast }$ (the second term proportional to $\unicode[STIX]{x1D6FD}z$ is small and zero at the wall). We plot $u_{\unicode[STIX]{x1D711}}^{\ast }$ and $\unicode[STIX]{x2202}_{r}u_{\unicode[STIX]{x1D711}}^{\ast }$ in figures 5(c) and 5(d) for the same cases as in figures 5(a) and 5(b). We see that the former achieves a peak and the latter changes sign at $r^{\ast }=1.41$ and $r^{\ast }=2.12$ , resulting in $F_{r}$ changing sign at those radial distances. This in turn changes the direction of $\langle u^{(2)}\rangle _{t}$ along the wall at these locations signalling the radial extent of the vortex $l_{vortex}$ . In fact, one finds from (2.36)

(3.1) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{\ast }}{\unicode[STIX]{x2202}r^{\ast }}\right|_{z^{\ast }=0}=0\rightarrow \frac{l_{vortex}}{R_{0}}=\frac{h}{\sqrt{2}R_{0}}. & & \displaystyle\end{eqnarray}$$

Note that this expression is independent of all other parameters and primarily results from the functional dependence in (2.36) arising due to the reflection from the wall. It matches with the values found above: $l_{vortex}=1.41R_{0}$ and $2.12R_{0}$ for $h=2R_{0}$ and $h=3R_{0}$ respectively. In contrast to the vortex length, its vertical extent $d_{vortex}$ remains unchanged with varying $h$ , the microbubble separation from the wall ( $0.32R_{0}$ for both $h^{\ast }=h/R_{0}=2$ and $h^{\ast }=3$ ). It is of the order ${\sim}\unicode[STIX]{x1D6FF}$ , the boundary layer thickness. One can compute it by considering the vertical component of streaming velocity close to the axis of symmetry being zero $\langle w^{(2)}\rangle _{t}(r^{\ast }\rightarrow 0,d_{vortex}/R_{0})=0$ . Near the axis of symmetry where $r=0$ , one can approximate

(3.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\simeq \text{Lim}_{r\rightarrow 0}\frac{1}{r}\left(\left.\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right|_{r}-\left.\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right|_{0}\right)=\text{Lim}_{r\rightarrow 0}\left.\frac{1}{r}\frac{\unicode[STIX]{x2202}^{2}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r^{2}}\right|_{r}. & & \displaystyle\end{eqnarray}$$

Therefore, close to the axis, one can write (2.23) as

(3.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle w^{(2)}\rangle _{t}\simeq -\frac{1}{2\unicode[STIX]{x1D714}}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}\right)\nonumber\\ \displaystyle & & \displaystyle \quad \times \left(-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-2\unicode[STIX]{x1D6FD}z}-\frac{1}{4\unicode[STIX]{x1D6FD}}\text{e}^{-\unicode[STIX]{x1D6FD}z}(8\unicode[STIX]{x1D6FD}z\sin \unicode[STIX]{x1D6FD}z+12\sin \unicode[STIX]{x1D6FD}z+20\cos \unicode[STIX]{x1D6FD}z)-\frac{10}{4}z+\frac{42}{8\unicode[STIX]{x1D6FD}}\right).\qquad\end{eqnarray}$$

One then obtains

(3.4) $$\begin{eqnarray}\displaystyle \langle w^{(2)}\rangle _{t}=0\rightarrow \frac{d_{vortex}}{R_{0}}=\frac{1.6}{R_{0}\unicode[STIX]{x1D6FD}}=1.6\frac{\unicode[STIX]{x1D6FF}}{R_{0}}=\frac{1.6}{R_{0}}\sqrt{\frac{2\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D714}}}=\frac{1.6\sqrt{2}}{\sqrt{Re_{b}}}. & & \displaystyle\end{eqnarray}$$

Here, it results in $d_{vortex}=0.32R_{0}$ as seen in figure 5(a,b). The width therefore is typical of boundary layer scaling and depends only on the excitation frequency and kinematic viscosity of the fluid.

Figure 5. Microstreaming streamlines near the plane rigid wall when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $P_{ex}^{\ast }=0.7$ and the microbubble is located at (a)  $h/R_{0}=2$ and (b)  $h/R_{0}=3$ . (c,d) Corresponding non-dimensional amplitude of radial potential velocity (black solid line) on the rigid wall and its derivative (red dashed line) with respect to radial distance along the wall. (e,f) Corresponding non-dimensional shear stress on the wall using limiting streaming velocity and microstreaming velocity.

Figures 5(e) and 5(f) show the shear stress on the wall induced by the microstreaming flow shown in figures 5(a) and 5(b), respectively. It changes sign at $r=l_{vortex}$ where the radial first-order irrotational velocity achieves its maximum. The maximum shear stress on the wall, as expected, decreases when the microbubble is excited further away from the wall. The value of the maximum shear stress appears at a distance $r_{max\,stress}^{\ast }$ . One can take the expression (2.38) and find its maximum at $r_{max\,stress}^{\ast }=0.2865h^{\ast }$ as was also found by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2010a ,Reference Doinikov and Bouakaz b ). For the conditions investigated here ( $R_{0}=1.6~\unicode[STIX]{x03BC}\text{m}$ , $f_{ex}=3~\text{MHz}$ , $P_{ex}=70~\text{kPa}$ ), the maximum shear stress is computed to be 12.8 and 2.2 Pa (note that $P_{0}=10^{5}~\text{Pa}$ ) from figure 5(e,f). We noted before that an average streaming shear stress expression $\unicode[STIX]{x1D70F}_{L}=\unicode[STIX]{x1D707}u_{L}/\unicode[STIX]{x1D6FF}$ has been widely used in the literature to estimate the maximum shear stress on the wall induced by microstreaming. In figure 5(e,f) the non-dimensional average streaming shear stress $\unicode[STIX]{x1D70F}_{L}^{\ast }$ (non-dimensionalized with respect to the ambient pressure) on the wall has been calculated using the limiting streaming velocity $u_{L}$ ( $z\rightarrow \infty$ in (2.38))

(3.5) $$\begin{eqnarray}\displaystyle u_{L}=-\frac{1}{\unicode[STIX]{x1D714}}\left(\frac{3}{8}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D711}}^{2}}{\unicode[STIX]{x2202}r}+\frac{u_{\unicode[STIX]{x1D711}}^{2}}{2r}\right). & & \displaystyle\end{eqnarray}$$

Note also that the direction of the radial limiting velocity is opposite to that of the radial microstreaming velocity close to the wall as can also be observed from the microstreaming streamlines in the regions close to and far from the wall. Therefore, opposite sign has been chosen for computing $\unicode[STIX]{x1D70F}_{L}^{\ast }$ and compared to the wall stress. In any event, $\unicode[STIX]{x1D70F}_{L}^{\ast }$ which is often used in the literature is computed to be much higher (almost three times in figure 5 e,f) than the wall shear stress.

Observing the motion of the tracer particles to delineate microstreaming, one obtains the Lagrangian streamlines. Figure 6 compares Lagrangian with Eulerian streamline patterns demonstrating that for the same flow field they are significantly different. Specifically, the vertical width of the vortical structure is narrower in the Lagrangian field causing the Lagrangian velocity to be of opposite sign at the same position. In figure 7, we plot the radial and vertical velocity components in both Eulerian and Lagrangian descriptions at three different vertical locations $z=0.03R_{0}$ , $z=0.16R_{0}$ and $z=0.32R_{0}$ from the wall. We find that the velocities in the two different descriptions are similar and in the same direction closest to the wall, but at $z=0.125R_{0}$ , while Eulerian radial velocity is positive, i.e. outward close to the axis in the bottom region of the vortex, its Lagrangian counterpart is mostly negative in the top part of the Lagrangian vortex (figure 7 a). The vertical velocity in figure 7(b) shows similar differences in the two descriptions. Note that the highest magnitude of the radial Eulerian velocity at $z=0.32R_{0}$ (top of the vortex) is ${\sim}3~\text{mm}~\text{s}^{-1}$ .

Figure 6. Comparing microstreaming streamlines near the plane rigid wall due to a free microbubble using Lagrangian and Eulerian streaming when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $h/R_{0}=2$ , $P_{ex}^{\ast }=0.7$ .

Figure 7. Comparing (a) radial component and (b) vertical component of Lagrangian and Eulerian microstreaming velocity at three different vertical distances $z=0.03R_{0}$ , $z=0.16R_{0}$ and $z=0.32R_{0}$ from the wall when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $h/R_{0}=2$ , $P_{ex}^{\ast }=0.7$ .

We have examined the effects of bubble translational motion on the theoretical expression of microstreaming in § 2.6 noting that the effects are small. In figure 8, we validate this conclusion by noting that accounting for translation leads to very little change in the streamline pattern.

Figure 8. Streamlines due to microstreaming near the plane rigid wall due to a free microbubble with and without accounting for translational motion when $Re_{b}=48,Eu_{b}=0.11$ ,  $We_{b}=20,h/R_{0}=2$ , $P_{ex}^{\ast }=0.7$ .

Figure 9. Vorticity contours due to the ring vortex induced by free microbubble when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ at $h/R_{0}=2$ , $P_{ex}^{\ast }=0.7$ .

Figure 9 plots the vorticity. The region of high negative vorticity (on the right-hand side of the figure; the left-hand side is antisymmetric) in the area adjoining the wall is generated by the strong shear due to the no-slip condition, coinciding with the flatter streamlines in the lower part of the vortex (figure 8). The vorticity increases with vertical distance from the wall changing sign to become positive and achieving the highest positive value near the edge of the vortex flow due to the strong shear layer there. The vorticity field, although generated by the vortex flow, does not provide an easy identification of the vortex, embodying the problem of vortex identification described in § 2.8. We identify the vortex region using the d2 method – inside the region enclosed by the $d_{2}=0$ curve in figure 9, where $d_{2}$ according to (2.44) is negative – and compute circulation of the vortex. We note that the region identified by the d2 method is smaller than what can be estimated from the streamline plot in figure 8. The closely spaced straight streamlines near the boundaries of the vortex in figure 8 were identified by the method as due to the boundary shear rather than due to the vortex. As described before, such difficulties of vortex identification methods have been discussed before in the literature (Chakraborty et al. Reference Chakraborty, Balachandar and Adrian2005). However, the vortex identified allows us to compute the circulation. Figure 10(a) shows the non-dimensional circulation $\unicode[STIX]{x1D6EC}^{\ast }$ of the axisymmetric vortex ring (non-dimensionalized by $R_{0}^{2}\unicode[STIX]{x1D714}$ ) as a function of separation distance $h$ of the microbubble from the wall. We also plot the maximum shear stress $\unicode[STIX]{x1D70F}_{max}^{\ast }$ in the same figure. They both decrease with increasing separation of the microbubble from the wall with a power-law variation with $h/R_{0}$ .

We investigate these variations in detail. In order to accommodate widely different magnitudes of these quantities in the same figure, we plot in figure 10(b) these quantities – circulation, maximum shear stress – by normalizing them by their highest values occurring at $h/R_{0}=1$ , i.e. $\unicode[STIX]{x1D6EC}^{\ast }(h/R_{0}=1)$ and $\unicode[STIX]{x1D70F}_{max}^{\ast }(h/R_{0}=1)$ . We also plot maximum vorticity $\unicode[STIX]{x1D6FA}_{max}$ normalizing it by $\unicode[STIX]{x1D6FA}_{max}(h/R_{0}=1)$ . We note similar scaling in shear stress and vorticity $\unicode[STIX]{x1D70F}_{max}^{\ast },\unicode[STIX]{x1D6FA}_{max}\sim (h/R_{0})^{-4.3}$ . This can be explained by noting that the term $\unicode[STIX]{x2202}\langle u_{2}\rangle /\unicode[STIX]{x2202}y$ responsible for $\unicode[STIX]{x1D70F}$ is the dominant part in $\unicode[STIX]{x1D6FA}$ . We find $\unicode[STIX]{x1D6EC}^{\ast }\sim (h/R_{0})^{-3.4}$ . Approximate relation $\unicode[STIX]{x1D6EC}^{\ast }\sim \unicode[STIX]{x1D6FA}_{max}(h/R_{0})$ can be explained by noting the Stokes theorem $\unicode[STIX]{x1D6EC}=\int \unicode[STIX]{x1D6FA}\,\text{d}r\,\text{d}z\sim \unicode[STIX]{x1D6FA}_{max}l_{vortex}d_{vortex}$ and the relations found before (3.1) and (3.4). These relations are also seen for other microbubble radii $R_{0}=3~\unicode[STIX]{x03BC}\text{m}$ and $R_{0}=5~\unicode[STIX]{x03BC}\text{m}$ (correspondingly different $Re_{b}$ , $Eu_{b}$ and $We_{b}$ ) in figures 10(c) and  10(d). The exact power-law indices differ but $\unicode[STIX]{x1D70F}_{max}^{\ast }$ and $\unicode[STIX]{x1D6FA}_{max}$ have the same index and $\unicode[STIX]{x1D6EC}^{\ast }\sim \unicode[STIX]{x1D6FA}_{max}(h/R_{0})$ approximately holds in all cases.

Figure 10. (a) Non-dimensional circulation of the ring vortex and non-dimensional maximum shear stress as a function of the separation distance of the microbubble from the wall when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $P_{ex}^{\ast }=0.7$ . (b) Circulation, maximum vorticity and maximum wall shear stress normalized by their values at $h/R_{0}=1$ as a function of $h/R_{0}$ for the same condition as in (a). (c) The same quantities as plotted in (b) for $Re_{b}=170,Eu_{b}=0.03,We_{b}=133$ . (d) The same quantities as plotted in (b) for $Re_{b}=471$ , $We_{b}=617$ , $Eu_{b}=0.011$ .

3.2 Effects of microbubble coating on microstreaming

Here we consider the effects of the stabilizing coating on a microbubble. As noted before, we consider the contrast agent Sonazoid and use its property values determined using the exponential elasticity model (Paul et al. Reference Paul, Katiyar, Sarkar, Chatterjee, Shi and Forsberg2010): $\unicode[STIX]{x1D6FE}_{0}=0.019~\text{N}~\text{m}^{-1}$ , $E_{0}^{s}=0.55~\text{N}~\text{m}^{-1}$ , $\unicode[STIX]{x1D6FC}^{s}=1.5$ and $\unicode[STIX]{x1D705}^{s}=1.2\times 10^{-8}~\text{N}~\text{s}~\text{m}^{-1}$ . Note that we found qualitatively similar results with other models (not shown here) such as that of Marmottant et al. (Reference Marmottant, van der Meer, Emmer, Versluis, de Jong, Hilgenfeldt and Lohse2005). We choose bubble radius of $1.6~\unicode[STIX]{x03BC}\text{m}$ , average size of Sonazoid $h=2R_{0}$ and identical excitation parameters (70 kPa, 3 MHz) as before giving rise to $Re_{b}=48,Eu_{b}=0.11$ , $We_{b}=20$ , $Re_{b}^{s}=6.4,We_{b}^{s}=2.6$ . Two material non-dimensional parameters are $\unicode[STIX]{x1D6FE}_{0}/E_{0}^{s}=0.0345$ and $\unicode[STIX]{x1D6FC}^{s}=1.5$ . Resulting non-dimensional periodic oscillation amplitude is $\unicode[STIX]{x1D700}=0.08$ , smaller than the free bubble. The bubble streaming Reynolds number $Re_{bs}=0.06$ is also correspondingly smaller. Figure 11(a) shows the microstreaming streamlines due to a pulsating Sonazoid microbubble for this condition. The streamline pattern is similar to that due to a free microbubble without any coating shown in figure 4(a). The dimensions of the vortex $l_{vortex}=1.41R_{0}$ and $d_{vortex}=0.32R_{0}$ are the same as those in free vortex (as determined, for example, by equation (3.1)). In figure 11(b), we show the wall shear stress for the same condition. Shear stress profiles are similar to those for a free bubble, reaching zero value at $r^{\ast }=l_{vortex}/R_{0}$ . The value of the maximum shear stress for this coated bubble is 6 Pa, smaller than the case of a free bubble for the same conditions (size, excitation parameters and separation) with identical position and the shear stress expectedly decreases and with increasing bubble separation; at the same time, $l_{vortex}$ increases as per equation (3.1).

Figure 11. (a) Streaming velocity field near the plane rigid wall due to coated microbubble when $Re_{b}=48$ , $Eu_{b}=0.11$ , $We_{b}=20$ , $Re_{b}^{s}=6.4$ , $We_{b}^{s}=2.6$ , $h/R_{0}=2$ , $\unicode[STIX]{x1D6FE}_{0}/E_{0}^{s}=0.0345$ and $\unicode[STIX]{x1D6FC}^{s}=1.5$ . (b) Wall shear stress for the same conditions as well as for two other values of $h/R_{0}$ .

Figure 12. Non-dimensional circulation of half the ring vortex and maximum wall shear stress due to coated microbubble (a) with respect to dilatation viscosity and (b) with respect to coating elasticity when $Re_{b}=48,Eu_{b}=0.11,We_{b}=20$ , $h/R_{0}=2$ and $\unicode[STIX]{x1D6FC}^{s}=1.5$ .

To examine the effects of the shell parameters, in figures 12(a) and 12(b) we plot the non-dimensional circulation of the vortex $\unicode[STIX]{x1D6EC}^{\ast }$ and non-dimensional maximum shear stress on the wall with respect to shell viscosity ratio $\unicode[STIX]{x1D705}^{s}/\unicode[STIX]{x1D705}_{(sonazoid)}^{s}$ and shell elasticity ratio $E_{0}^{s}/E_{0(sonazoid)}^{s}$ , respectively, with other parameters held the same as before. The values for Sonazoid coating agent are $E_{0(sonazoid)}^{s}=0.55~\text{N}~\text{m}^{-1}$ and $\unicode[STIX]{x1D705}_{sonazoid}^{s}=1.2\times 10^{-8}~\text{N}~\text{s}~\text{m}^{-1}$ . Here $\unicode[STIX]{x1D6EC}^{\ast }$ and $\unicode[STIX]{x1D70F}_{max}^{\ast }$ show identical variations with coating parameters. Coating parameters affect the bubble pulsation amplitude and that in turn affects $\unicode[STIX]{x1D6EC}^{\ast }$ and $\unicode[STIX]{x1D70F}_{max}^{\ast }$ . The variations in $\unicode[STIX]{x1D6EC}^{\ast }$ and $\unicode[STIX]{x1D70F}_{max}^{\ast }$ are then identical for the same reasons as for the free bubble. Figure 12(a) shows that circulation and maximum shear stress decrease with increasing dilatational viscosity, as can be expected from the fact that at higher dilatational viscosities, the coated microbubble has smaller pulsation amplitudes and thereby lower streaming velocity and lower shear stress as well as circulation. The change of circulation and shear stress with shell elasticity (figure 12 b) shows a maximum corresponding to $E_{0}^{s}=0.27~\text{N}~\text{m}^{-1}$ . The maximum signals a resonance of the coated microbubble at the excitation frequency; changing the shell elasticity changes the spring constant and correspondingly the resonance frequency of the coated bubble according to equation (2.43), i.e. the oscillation is maximum at this value of  $E_{0}^{s}$ .