2.1. First-order scattered field

Before we proceed to calculations, a comment should be made. In the case of a single bubble undergoing shape modes, a scattered wave generated by a certain mode has the same angular dependence as the generating mode. For example, a shape mode depending on ${P_n}(\cos \theta )$, where ${P_n}$ is the Legendre polynomial of degree n, generates a scattered wave depending on ${P_n}(\cos \theta )$. In the case of two bubbles, the structure of the total scattered wave is different. Due to multiple rescattering between the bubbles, every oscillation mode generates a scattered wave with ${P_m}(\cos \theta )$ of all degrees, i.e. ${P_m}(\cos \theta )$. This effect can be seen in solutions obtained by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2015). Although Doinikov & Bouakaz (Reference Doinikov and Bouakaz2015) assume that two bubbles only undergo pulsation and translation, the scattered wave includes terms with ${P_m}(\cos \theta )$ of all degrees and the amplitudes of these terms are related to one another. This fact is taken into account in the derivation below.

We assume that each bubble undergoes shape modes, which can have different orders and different frequencies because some of them are excited parametrically. The modes can be divided into groups in accordance with their frequency. These groups can be considered separately because the interaction of modes having the same frequency can only contribute to acoustic microstreaming (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a,Reference Doinikov, Cleve, Regnault, Mauger and Inserrab; Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a,Reference Inserra, Regnault, Cleve, Mauger and Doinikovb). We consider modes with a frequency $\omega $, which can be equal to the driving frequency or any other sub- or ultra-harmonic of the driving frequency if the modes are excited parametrically.

We assume that the jth bubble undergoes ${N_j}$ modes with a frequency $\omega $. The modes have numbers $M_1^{(j)},M_2^{(j)},M_{{N_j}}^{(j)}$. In this case, the $\omega $-dependent deformation of the surface of the jth bubble can be represented as

(2.1)\begin{equation}r_s^{(j)} = {\textrm{e}^{ - \textrm{i}\omega t}}\mathop \sum \limits_{n = M_1^{(j)}}^{M_{{N_j}}^{(j)}} s_n^{(j)}{P_n}({\mu _j}),\end{equation}

where ${\mu _j} = \cos {\theta _j}$ and $s_n^{(j)}$ is the complex amplitude of the nth mode of the jth bubble, which is assumed to be known (for example, obtained from the experimental bubble dynamics).

The total first-order liquid velocity generated by both bubbles is given by

(2.2)\begin{equation}\boldsymbol{v} = {\boldsymbol{v}^{(1)}} + {\boldsymbol{v}^{(2)}},\end{equation}

where ${\boldsymbol{v}^{(j)}}$ is the first-order liquid velocity generated by the jth bubble. The value of ${\boldsymbol{v}^{(j)}}$ can be represented as

(2.3)\begin{equation}{\boldsymbol{v}^{(j)}} = \boldsymbol{\nabla }{\varphi ^{(j)}} + \boldsymbol{\nabla } \times {\boldsymbol{\psi }^{(j)}},\end{equation}

where ${\varphi ^{(j)}}$ and ${\boldsymbol{\psi }^{(j)}}$ are the scalar and the vector velocity potentials, respectively.

In view of axial symmetry, ${\varphi ^{(j)}}$ and ${\boldsymbol{\psi }^{(j)}}$ are written as (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a)

(2.4)\begin{gather}{\varphi ^{(j)}} = {\textrm{e}^{ - \textrm{i}\omega t}}\mathop \sum \limits_{n = 0}^\infty a_n^{(j)}{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)^{n + 1}}{P_n}({\mu _j}),\end{gather}

(2.5)\begin{gather}{\boldsymbol{\psi }^{(j)}} = {\psi ^{(j)}}({r_j},{\theta _j},t){\boldsymbol{e}_{\boldsymbol{\varepsilon j}}} = {\textrm{e}^{ - \textrm{i}\omega t}}{\boldsymbol{e}_{\boldsymbol{\varepsilon j}}}\mathop \sum \limits_{n = 1}^\infty b_n^{(j)}h_n^{(1)}({k_v}{r_j})P_n^1({\mu _j}),\end{gather}

where ${R_{j0}}$ is the equilibrium radius of the jth bubble, $P_n^1$ is the associated Legendre polynomial of the first order and degree n, $h_n^{(1)}$ is the spherical Hankel function of the first kind, ${k_v} = (1 + i)/\delta $ is the viscous wavenumber, $\delta = \sqrt {2\nu /\omega } $ is the viscous penetration depth, $\nu = \eta /\rho $ is the kinematic liquid viscosity, $\eta $ is the dynamic liquid viscosity, $\rho $ is the equilibrium liquid density and ${\boldsymbol{e}_{\boldsymbol{\varepsilon j}}}$ is the unit azimuth vector of the jth bubble. Note that axial symmetry allows us to set ${\varepsilon _1} = {\varepsilon _2}$ and ${\boldsymbol{e}_{\boldsymbol{\varepsilon }1}} = {\boldsymbol{e}_{\boldsymbol{\varepsilon }2}}$.

Substitution of (2.4) and (2.5) into (2.3) yields

(2.6)\begin{gather}v_r^{(j)} ={-} \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_j}}}\sum\limits_{n = 0}^\infty {(n + 1)} \left[ {a_n^{(j)}{{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)}^{n + 1}} + nb_n^{(j)}h_n^{(1)}({k_v}{r_j})} \right]{P_n}({\mu _j}),\end{gather}

(2.7)\begin{gather}v_\theta ^{(j)} = \frac{{{\textrm{e}^{ - \textrm{i}\omega t}}}}{{{r_j}}}\sum\limits_{n = 1}^\infty {\left\{ {a_n^{(j)}{{\left( {\frac{{{R_{j0}}}}{{{r_j}}}} \right)}^{n + 1}} - b_n^{(j)}[h_n^{(1)}({k_v}{r_j}) + {k_v}{r_j}h_n^{(1)/}({k_v}{r_j})]} \right\}P_n^1({\mu _j})} ,\end{gather}

where the prime denotes the derivative with respect to the argument in brackets. Note that when calculating (2.6) and (2.7), we have used mathematical properties of ${P_n}(\mu )$ and $P_n^1(\mu )$ (Abramowitz & Stegun Reference Abramowitz and Stegun1972).

To find $a_n^{(j)}$ and $b_n^{(j)}$, called linear scattering coefficients, boundary conditions at the bubble surfaces are applied. Since this calculation is cumbersome, it is performed in Appendix A.

As an example, figure 2 shows the magnitudes of $a_n^{(j)}$ and $b_n^{(j)}$ calculated for bubbles with ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ at the inter-bubble distance $d = 300\;\mathrm{\mu }\textrm{m}$ and the oscillation frequency $f = 30\;\textrm{kHz}$. The surrounding liquid is water. The bubbles are assumed to undergo the radial oscillation with an amplitude of 5 μm. As one can see, the scattering coefficients with n > 0 are not zero despite the fact that both bubbles only pulsate, but the magnitude of the scattering coefficients quickly decreases with increasing n. It is worth noticing that the coefficients $a_n^{(1)}$ are six orders of magnitudes larger than the coefficients $b_n^{(1)}$. This difference is explained by the weight of the radial functions ${({R_{j0}}/{r_j})^{n + 1}}$ and Hankel function $h_n^{(1)}({k_v}{r_j})$ that are respectively associated with the coefficients $a_n^{(1)}$ and $b_n^{(1)}$ in (2.6) and (2.7).

Figure 2.

Magnitude of the linear scattering coefficients for bubble 1. Due to the equivalence of the calculation parameters, the linear scattering coefficients for bubble 2 have the same magnitude.

2.2. Acoustic streaming

As shown by Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a), the Eulerian velocity of acoustic streaming can be written as

(2.8)\begin{equation}{\boldsymbol{v}_E} = \boldsymbol{\nabla } \times \boldsymbol{\varPsi },\end{equation}

where $\boldsymbol{\varPsi }$ obeys the vorticity equation derived by Westervelt (Reference Westervelt1953) for the case of solenoidal first-order motion,

(2.9)\begin{equation}{\varDelta ^2}\boldsymbol{\varPsi } ={-} \frac{1}{\nu }\boldsymbol{\nabla } \times \langle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\rangle ,\end{equation}

with $\left\langle {} \right\rangle$ denoting the time average.

Equation (2.9) is shown by Doinikov et al. (Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a) to be transformed to

(2.10)\begin{equation}{\varDelta ^2}\boldsymbol{\varPsi } = \frac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [(\varDelta {\boldsymbol{\psi }^ \ast }) \times \boldsymbol{v}]\} ,\end{equation}

where Re means “the real part of” and the asterisk denotes the complex conjugate. Substituting $\boldsymbol{v} = {\boldsymbol{v}^{(1)}} + {\boldsymbol{v}^{(2)}}$ and $\boldsymbol{\psi } = {\boldsymbol{\psi }^{(1)}} + {\boldsymbol{\psi }^{(2)}}$, one obtains

(2.11)\begin{equation}\begin{array}{ccccc} {\varDelta ^2}\boldsymbol{\varPsi } & = \dfrac{1}{{2\nu }}Re \{ \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(1) \ast }}) \times {\boldsymbol{v}^{(1)}}] + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(2) \ast }}) \times {\boldsymbol{v}^{(2)}}]\\ & \quad + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(1) \ast }}) \times {\boldsymbol{v}^{(2)}}] + \boldsymbol{\nabla } \times [\varDelta ({\boldsymbol{\psi }^{(2) \ast }}) \times {\boldsymbol{v}^{(1)}}]\} . \end{array}\end{equation}

Equation (2.11) can be represented in the coordinates $({r_1},{\theta _1})$ or $({r_2},{\theta _2})$. It is shown in Appendix B that $\boldsymbol{\varPsi }$ in the coordinates $({r_1},{\theta _1})$ is written as

(2.12)\begin{equation}\boldsymbol{\varPsi }({r_1},{\theta _1}) = \varPsi ({r_1},{\theta _1}){\boldsymbol{e}_{\varepsilon 1}},\end{equation}

where $\varPsi ({r_1},{\theta _1})$ obeys the following equation:

(2.13) \begin{align} & {\left( {\Delta_{r\theta }^{(1)} - \dfrac{1}{{r_1^2{{\sin }^2}{\theta_1}}}} \right)^2}\varPsi ({r_1},{\theta _1}) = \dfrac{1}{{2\nu R_{10}^4}}\left\{\displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty {{P_n}({\mu_1})P_m^1({\mu_1})[K_{1nm}^{(1)}({r_1}) + L_{1nm}^{(1)}({r_1}) + M_{1nm}^{(1)}({r_1})]}\right.\nonumber\\ & \quad \left. + \sqrt {1 - \mu_1^2} \displaystyle\sum_{n,m = 1}^\infty {P_n^1({\mu_1})P_m^{1/}({\mu_1})[K_{2nm}^{(1)}({r_1}) + K_{2mn}^{(1)}({r_1}) + L_{2nm}^{(1)}({r_1}) + L_{2mn}^{(1)}({r_1}) + M_{2nm}^{(1)}({r_1}) + M_{2mn}^{(1)}({r_1})]}\right\}. \end{align}

Here, the operator $\varDelta _{r\theta }^{(1)}$ is defined by (B5) and the ${r_1}$-dependent functions are calculated by (B9), (B10), (B13), (B14), (B20) and (B21).

Equation (2.13) is valid for ${r_1} < d$ since this restriction is imposed on coordinate transformations (Appendix A) that were used to derive (2.13). Therefore, the solution of (2.13) describes acoustic streaming in the domain ${D_1}$ shown in figure 3.

Figure 3.

Regions where (2.13) and (2.15) are valid.

As shown in Appendix B, $\boldsymbol{\varPsi }$ in the coordinates $({r_2},{\theta _2})$ is written by

(2.14)\begin{equation}\boldsymbol{\varPsi }({r_2},{\theta _2}) = \varPsi ({r_2},{\theta _2}){\boldsymbol{e}_{\varepsilon 2}},\end{equation}

where $\varPsi ({r_2},{\theta _2})$ is calculated by

(2.15) \begin{align}& {\left( {\Delta_{r\theta }^{(2)} - \dfrac{1}{{r_2^2{{\sin }^2}{\theta_2}}}} \right)^2}\varPsi ({r_2},{\theta _2}) = \dfrac{1}{{2\nu R_{20}^4}}\left\{ \displaystyle\sum_{\substack{n = 0\\ m = 1}}^\infty {{P_n}({\mu_2})P_m^1({\mu_2})[K_{1nm}^{(2)}({r_2}) + L_{1nm}^{(2)}({r_2}) + M_{1nm}^{(2)}({r_2})]}\right.\nonumber\\ & \quad \left. { + \sqrt {1 - \mu_2^2} \displaystyle\sum_{n,m = 1}^\infty {P_n^1({\mu_2})P_m^{1/}({\mu_2})[K_{2nm}^{(2)}({r_2}) + K_{2mn}^{(2)}({r_2}) + L_{2nm}^{(2)}({r_2}) + L_{2mn}^{(2)}({r_2}) + M_{2nm}^{(2)}({r_2}) + M_{2mn}^{(2)}({r_2})]} } \right\}. \end{align}

Here, the operator $\varDelta _{r\theta }^{(2)}$ is defined by (B5) and the ${r_2}$-dependent functions are calculated by (B26), (B27), (B30), (B31), (B34) and (B35).

Equation (2.15) is valid for ${r_2} < d$ and hence its solution describes acoustic streaming in the domain ${D_2}$ shown in figure 3. Equations (2.13) and (2.15) must give the same results in the space between the bubbles so the domains ${D_1}$ and ${D_2}$ can be combined to get acoustic streaming around both bubbles, as shown in figure 3.

It should be emphasised that the above restrictions on the solutions of (2.13) and (2.15) are not related to the ratio of the bubble radii to the viscous penetration depth or to the inter-bubble distance. The restrictions shown in figure 3 result from geometrical transformations of coordinates that are used in our derivation and hence are of purely geometrical nature. They are necessary to provide the convergence of infinite sums used in the solutions.

Equations (2.13) and (2.15) are solved in Appendix C. The result is given by

(2.16)\begin{equation}\varPsi ({r_j},{\theta _j}) = \sum\limits_{l = 1}^\infty {\varPsi _l^{(j)}({r_j})P_l^1({\mu _j})} ,\quad j = 1,2,\end{equation}

where $\varPsi _l^{(j)}({r_j})$ is calculated by (C8). Substituting (2.16) into (2.8), one obtains the components of the Eulerian streaming velocity,

(2.17)\begin{gather}{v_{Er}}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {l(l + 1)\varPsi _l^{(j)}({r_j}){P_l}({\mu _j})} ,\end{gather}

(2.18)\begin{gather}{v_{E\theta }}({r_j},{\theta _j}) ={-} \frac{1}{{{r_j}}}\sum\limits_{l = 1}^\infty {[\varPsi _l^{(j)}({r_j}) + {r_j}\varPsi _l^{(j)/}({r_j})]P_l^1({\mu _j})} ,\end{gather}

where $\varPsi _l^{(j)/}({r_j})$ is calculated by (C17).

The components of the Lagrangian streaming velocity are defined by

(2.19)\begin{gather}{v_{Lr}}({r_j},{\theta _j}) = {v_{Er}}({r_j},{\theta _j}) + {v_{Sr}}({r_j},{\theta _j}),\end{gather}

(2.20)\begin{gather}{v_{L\theta }}({r_j},{\theta _j}) = {v_{E\theta }}({r_j},{\theta _j}) + {v_{S\theta }}({r_j},{\theta _j}),\end{gather}

where ${v_{Sr}}({r_j},{\theta _j})$ and ${v_{S\theta }}({r_j},{\theta _j})$ denote the components of the Stokes drift velocity (Longuet-Higgins Reference Longuet-Higgins1998), which are calculated in Appendix E and given by (E15) and (E16).

Combining the solutions in the coordinates $({r_1},{\theta _1})$ and $({r_2},{\theta _2})$ provides acoustic microstreaming in the area shown in figure 3. As a consequence, when the distance between the bubbles is reduced, the boundary of the spatial domain where our second-order solutions are valid shifts to the bubble interfaces. It is also worth mentioning that the domain $D = \{ {r_1} < d \cup {r_2} < d\} $ is a theoretical upper limit for the spatial domain where (2.19) and (2.20) are valid. Since in the process of numerical calculations one has to truncate the infinite sums in (2.17) and (2.18), one can only approach this upper limit by increasing the number of terms in the truncated sums, which naturally leads to more time-consuming computations.

3.1. Convergence analysis

In the present modelling, the interaction between bubbles is considered both when the linear velocity field is calculated and when the second-order velocity field of acoustic microstreaming is calculated. Physically, the interaction is realised by re-scattering of sound between the bubbles. One can see that in the course of the calculation of the linear velocity, multiple scattering effects are included into the derivation of the linear scattering coefficients $a_n^{(j)}$ and $b_n^{(j)}$, which are given by the system of (A41) and (A42) for $j = 1$ and $n \ge 1$. The infinite summation appearing in (A41) and (A42) must be truncated at some finite value of a given number n of the linear scattering coefficients, while the coefficients appearing in this system of equations, which are given by (A43)–(A48) also involve infinite sums, which determine the number of re-scattering events that are allowed for. These sums must also be truncated to some finite value ${m_{lim}}$. The numerical case under consideration is the interaction between two bubbles of identical size ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ separated by the inter-bubble distance $d = 150\;\mathrm{\mu }\textrm{m}$ that oscillate radially in phase, with the amplitude 5 μm at the frequency 30 kHz. In order to estimate the number of the linear scattering coefficients required to accurately describe multiple scattering effects, figure 4(a) shows the evolution of the total energy $\sum\nolimits_{i = 0}^n {||a_i^{(1)}|{|^2}}$ carried by the linear scattering coefficients $a_i^{(1)}$ up to the value n. The energy reaches a plateau at $n = 5$, providing the limiting number of first-order interactions that should be taken into account. The linear scattering coefficient of the highest amplitude $a_i^{(1)}$, in the case $n \ne 0$, is plotted in figure 4(b) as a function of the limit value ${m_{lim}}$. Setting ${m_{lim}} = 5$ in our case is enough to ensure the convergence of the linear scattering coefficients.

Figure 4.

Convergence analysis of the numerical modelling for the case of two spherically oscillating bubbles of identical radii $({R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m)}$ driven by a 30 kHz ultrasound wave in water. (a) Convergence of the energy spectrum of the linear scattering coefficients as a function of the number n of scattered modes. (b) Convergence of the highest scattering coefficient $a_1^{(1)}$ a function of the finite upper limit ${m_{lim}}$ of the sums appearing in (A41)–(A48). (c,d) Convergence of the radial ${V_r}$ and tangential ${V_\theta }$ components of the Lagrangian streaming velocity as a function of the finite upper limit ${l_{lim}}$ of the sums appearing in (2.17) and (2.18).

In the second-order streaming velocities, the velocity field is also defined by an infinite sum of the angular functions ${P_n}(\mu )$ and $P_n^1(\mu )$, such as in the system of (2.17) and (2.18) for the components of the Eulerian streaming velocity. The numerical truncation of these infinite sums to a finite value ${l_{lim}}$ is quantified on the radial and tangential components of the Lagrangian streaming velocity in figure 4(c,d). The velocity components are computed at the location $({r_1} = {R_{01}} + 5\;\mathrm{\mu }\textrm{m,}{\theta _1} = {\rm \pi}/2)$ in the reference frame of bubble 1. The velocity components have converged for ${l_{lim}} \ge 6$. In the following, numerical results are only considered for which the convergence has been achieved.

3.2. Effect of the inter-bubble distance

A single bubble that pulsates radially in an unbounded fluid does not produce any microstreaming (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a). Due to the presence of a second bubble, multiple sound scattering leads to a non-zero microstreaming flow. Figures 5(a)–5(c) shows streamline patterns, for different inter-bubble distances, around two bubbles with identical equilibrium radii ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$, oscillating radially in phase with the amplitude $5\;\mathrm{\mu }\textrm{m}$ at the frequency 30 kHz. The bubbles are indicated by full circles and the flow pattern is only displayed within the region of validity $D = \{{{r_1} < d \cup {r_2} < d} \}$.

Figure 5.

Evolution of the streaming flow as a function of the inter-bubble distance d in the case of two spherically oscillating bubbles (with the same parameters as in figure 4); (a) $d = 150\;\mathrm{\mu }\textrm{m}$, (b) $d = 300\;\mathrm{\mu }\textrm{m}$, (c) $d = 750\;\mathrm{\mu }\textrm{m}$. (d) Evolution of the magnitude of the Lagrangian streaming velocity as a function of the inter-bubble distance. The velocity is measured at the locations $({r_1} = 1.3{R_{01}},{\theta _1} = {\rm \pi}/2)$ and $({r_2} = 1.5{R_{01}},{\theta _2} = {\rm \pi})$.

The streamline patterns consist of large lobes located above and below each bubble. The form and the size of the lobes are changed with increasing inter-bubble distance, namely the centre of rotation is shifted away from the bubbles and the extent of the vortices increases. Along the horizontal direction ${\theta _1} = 0$, the streaming flow is directed towards the bubbles in the inter-bubble space. Two pairs of recirculation vortices are observed in the distance from the bubbles along the horizontal axis. The size of these vortices decreases when the inter-bubble distance increases. Figure 5(d) depicts the evolution of the magnitude of the Lagrangian streaming velocity as a function of the inter-bubble distance at the locations $({r_1} = 1.3{R_{01}},{\theta _1} = {\rm \pi}/2)$ and $({r_2} = 1.5{R_{01}},{\theta _2} = {\rm \pi})$. The velocity magnitude rapidly decreases to values of the order of 10 μm s⁻¹. Such small velocities cannot be observed experimentally. It is worth noticing that these results are in good qualitative agreement with the experimental observations of Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016). Although their study is devoted to the analysis of streaming flows induced by pancake-like bubbles, the same kind of long-range streaming pattern is observed around the interacting bubble pair.

3.3. Contribution of the present model

In order to illustrate the contribution of the present model, we perform a comparison between the microstreaming pattern that is obtained by the superimposition of two streaming flows calculated by the single-bubble theory and the pattern predicated by the present model. The microstreaming pattern induced by each bubble in the non-interacting case is calculated by the theoretical approach of Inserra et al. (Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020b). Numerical results are shown in figure 6 for two identical bubbles with ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ at the inter-bubble distance $d = 200\;\mathrm{\mu }\textrm{m}$, exhibiting radial $(s_0^{(1,2)} = 5\;\mathrm{\mu }\textrm{m)}$ and mode 2 $(s_2^{(1,2)} = 10\;\mathrm{\mu }\textrm{m)}$ oscillations at the driving frequency 30 kHz in water. As one can see, the microstreaming patterns in figures 6(a) and 6(b) look markedly different. Note, for example, that the direction of streamlines in the areas above and below the bubbles is different in figures 6(a) and 6(b).

Figure 6.

Comparison of (a) the microstreaming pattern obtained by the superimposition of two flows given by the single-bubble theory and (b) the pattern predicted by the present model.

In order to reinforce the validity of these results, a computational fluid dynamics simulation has been performed using StarCCM+ software. The computational fluid dynamics (CFD) simulation is based on the interface motion of a deformable spheroid that sets in motion the fluid in its surroundings. The interface motion is described by (2.1) for each bubble. The interface deformation is achieved by the morphing module implemented in StarCCM+. This method was first validated successfully when investigating the microstreaming flow induced by a single axisymmetric bubble, whose theoretical modelling is known (Doinikov et al. Reference Doinikov, Cleve, Regnault, Mauger and Inserra2019a,Reference Doinikov, Cleve, Regnault, Mauger and Inserrab; Inserra et al. Reference Inserra, Regnault, Cleve, Mauger and Doinikov2020a,Reference Inserra, Regnault, Cleve, Mauger and Doinikovb). In the present study, the numerical modelling has been extended to two bubbles undergoing axisymmetric shape oscillations. The computational domain consists of a quarter of a millimetre sphere (liquid medium) containing two quarters of sub-millimetre spheroids (i.e. two bubbles), see Appendix F. The axes of each quarter of the sphere are aligned together. The Navier–Stokes and continuity equations are resolved through a laminar, incompressible and unsteady model using a pressure boundary applied at the far field and a slip condition at the interfaces. The smallest mesh size at the interface is 1 μm. To compare the fluid field given by the CFD simulation with that given by our analytical model, the velocity field is averaged over one period. Figure 7 displays the microstreaming pattern obtained by the CFD simulation. The microstreaming pattern is very similar to the analytic one shown in figure 6(b). The location and the size of the lobes are well captured and the flow direction is identical to that obtained analytically. Figure 7(b) shows a partial three-dimensional view of the streamlines induced by the bubble oscillations. This figure highlights the generation of vortex rings around the bubbles, which correspond to the lobes shown in the different cross-view.

Figure 7.

Microstreaming flow obtained by CFD simulation. (a) Cross-section view of the flow. (b) Three-dimensional view of the streamlines illustrating the vortex ring generation around the bubbles.

It is worth noting that the computational time of the CFD simulation is of the order of 24 hours, while the calculation of streamlines by the analytical model takes approximately 2 hours.

We have also made a comparison with the theory developed by Doinikov & Bouakaz (Reference Doinikov and Bouakaz2016). Recall that this theory, first, assumes the distance between the bubbles to be large compared with the bubble radii and, second, ignores the interaction between the bubbles in second-order solutions. Figure 8(a) demonstrates a streamline pattern predicted by this theory, while figure 8(b) shows a pattern given by the present model. The simulations were made for bubbles with ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ at $d = 300\;\mathrm{\mu }\textrm{m}$, assuming that both bubbles undergo radial oscillations with the amplitude 5 μm at the frequency 30 kHz in water. As one can see, figure 8(a) predicts an opposite direction of the streaming flow, which is likely to be a consequence of the fact that the theory of Doinikov & Bouakaz (Reference Doinikov and Bouakaz2016) only allows for the interaction of the bubbles in the first-order solutions.

Figure 8.

Comparison of (a) the streamline pattern given by the theory of Doinikov & Bouakaz (Reference Doinikov and Bouakaz2016) and (b) the pattern predicted by the present model.

3.4. Effect of viscosity

Figures 9–12 illustrate the effect of the liquid viscosity on the behaviour of the microstreaming. The simulations were made for bubbles with ${R_{10}} = {R_{20}} = 50\;\mathrm{\mu }\textrm{m}$ at the frequency 30 kHz, three values of the inter-bubble distance d and three values of the viscosity $\eta = 0.001$, 0.01 and 0.1 Pa s.

Figure 9.

Streamline patterns produced by two radially oscillating bubbles of equal radii at different inter-bubble distances and different values of the liquid viscosity: (a,d,g) $\eta = 0.001\;\textrm{Pa}\;\textrm{s}$; (b,e,h) $\eta = 0.01\;\textrm{Pa}\;\textrm{s}$; (c,f,i) $\eta = 0.1\;\textrm{Pa}\;\textrm{s}$.

Figure 10.

The magnitude of the Lagrangian streaming velocity $|{\boldsymbol{v}_L}|$, which is generated by two radially oscillating identical bubbles, as a function of the distance from the bubble surface along two directions. Direction 1 is counted from the surface of bubble 2 at ${\theta _2} = 0$, while direction 2 is counted from the surface of bubble 1 at ${\theta _1} ={-} {\rm \pi}/2$. The other parameters are as in figure 9.

Figure 11.

Streamline patterns produced by bubbles of equal radii, undergoing radial and mode 2 oscillations, at different inter-bubble distances and different values of the liquid viscosity: (a,d,g) $\eta = 0.001\;\textrm{Pa}\;\textrm{s}$; (b,e,h) $\eta = 0.01\;\textrm{Pa}\;\textrm{s}$; (c,f,i) $\eta = 0.1\;\textrm{Pa}\;\textrm{s}$.

Figure 12.

The magnitude of the Lagrangian streaming velocity $|{\boldsymbol{v}_L}|$, which is generated by two identical bubbles undergoing radial and mode 2 oscillations, as a function of the distance from the bubble surface along two directions. Direction 1 is counted from the surface of bubble 2 at ${\theta _2} = 0$, while direction 2 is counted from the surface of bubble 1 at ${\theta _1} ={-} {\rm \pi}/2$. The other parameters are as in figure 11.

Figure 9 shows streamline patterns in the case that both bubbles undergo radial oscillations with the amplitude 5 μm. As one can see, if d is not large, with increasing $\eta $, the form of vortices is changed noticeably. The vortices become increasingly detached from the bubble surfaces and their extent increases. Figure 10 shows how the magnitude of the Lagrangian streaming velocity $|{\boldsymbol{v}_L}|$ depends on the liquid viscosity. $|{\boldsymbol{v}_L}|$ depicted as a function of the distance from the bubble surface at two values of the inter-bubble distance along two directions. Direction 1 is counted from the surface of bubble 2 at ${\theta _2} = 0$, while direction 2 is counted from the surface of bubble 1 at ${\theta _1} ={-} {\rm \pi}/2$. As one can see, increasing $\eta $ leads to increasing $|{\boldsymbol{v}_L}|$ and the peak value of $|{\boldsymbol{v}_L}|$ is shifted away from the bubble surface. It is interesting to note that, along direction 2, the increase of $|{\boldsymbol{v}_L}|$ at $\eta = 0.01\;\textrm{Pa}\;\textrm{s}$ is higher than at $\eta = 0.1$ but the region where $|{\boldsymbol{v}_L}|$ is relatively high remains wider at $\eta = 0.1$.

Figure 11 shows streamline patterns in the case that both bubbles undergo radial and mode 2 oscillations. The amplitude of both modes is 5 μm and the phase shift between modes 0 and 2 is zero. As one can see, with increasing $\eta $, the vortices become increasingly detached from the bubble surfaces. At $d = 600\;\mathrm{\mu }\textrm{m}$ and $\eta = 0.1\;\textrm{Pa}\;\textrm{s}$, a complicated structure of streamlines is observed in the inter-bubble space. The dependence of the magnitude of the Lagrangian streaming velocity on the liquid viscosity is illustrated by figure 12, which is analogous to figure 10. As one can see, along direction 1, $|{\boldsymbol{v}_L}|$ at $\eta = 0.01\;\textrm{Pa}\;\textrm{s}$ is higher than at $\eta = 0.001\;\textrm{Pa}\;\textrm{s}$ but further increasing the viscosity to 0.1 Pa s results in a $|{\boldsymbol{v}_L}|$ that is lower than at $\eta = 0.001\;\textrm{Pa}\;\textrm{s}$. Other effects are similar to those in figure 10.

3.5. Comparison with existing experimental observations

Streaming flows induced in the vicinity of an interacting bubble pair have been reported by Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016) for flattened bubbles undergoing non-spherical oscillations in a microchannel. They consider the case of two bubbles of radii ${R_{10}} = {R_{20}} = 38\;\mathrm{\mu }\textrm{m}$, separated by $d = 200\;\mathrm{\mu }\textrm{m}$, undergoing spherical, translational and mode 2 oscillations. The associated mode dynamics, meaning the values of the complex amplitudes $s_{0,1,2}^{(1,2)}$, is provided. The resulting streaming pattern, which is shown in figure 13(a), consists of two anti-fountain vortices on both sides of the bubble pair and two pairs of smaller vortices that are located between the bubbles. Along the line joining the bubble centres, the streaming flow is directed towards the bubbles. Experimental streaming velocities are found to be of the order of millimetres per second. The measured input data for $s_{0,1,2}^{(1,2)}$ provided by Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016) allow one to perform a qualitative comparison with our model. Although our theory implies that bubbles are initially of spherical shape and located in an unbounded fluid, one can expect that their dynamics should exhibit similar qualitative features as in the case investigated by Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016). The theoretical streaming pattern is shown in figure 13(b) and is in good agreement with the experimental one: both the large lobes and the small vortices between the bubbles are observed. Also, our theory provides velocity magnitudes of the same order as in the experimental results.

Figure 13.

Comparison between (a) experimental observations of the steaming flow induced by a bubble pair undergoing radial, translational and mode 2 oscillations (reproduced from Mekki-Berrada et al. (Reference Mekki-Berrada, Combriat, Thibault and Marmottant2016), see the copyright form) and (b) the theoretical results given by the present model.