3.1 Limiting cases

Before proceeding any further, it is useful to examine first two important limiting cases.

3.1.1 Rising bubble in the absence of an acoustic field

We first consider the case which corresponds to a bubble that rises under the effect of buoyancy having a constant volume in the presence of a constant pressure field. Thus setting ${\dot{R}}=\ddot{R}=0$ and $R=R_{o}$ while assuming ${\dot{z}}=\ddot{z}=0$ it is also possible, using (3.22), to estimate the critical Bingham number, $Bn_{c}$ , below which the bubble becomes entrapped in the viscoplastic material:

(3.25) $$\begin{eqnarray}\displaystyle Bn_{c}=4\unicode[STIX]{x1D712}/3\approx 0.175. & & \displaystyle\end{eqnarray}$$

It is important to note that this value is very close to the value that has been reported by Tsamopoulos et al. (Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008) (i.e. $Bn_{c}=0.143$ ) and Dimakopoulos et al. (Reference Dimakopoulos, Pavlidis and Tsamopoulos2013) (i.e. $Bn_{c}=0.129$ ) who performed axisymmetric finite element simulations fully taking into account the deformability of the bubble using the regularized Papanastasiou model and the augmented Lagrangian method, respectively. The remarkable agreement with such detailed calculations corroborates the validity of the present analysis.

3.1.2 Pulsating bubble without translation and negligible buoyancy

Next, we consider the case which corresponds to a bubble that pulsates in the presence of an acoustic field without any translation. Thus setting $U=\dot{U}=0$ and neglecting buoyancy we obtain the following dimensional equation:

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle P_{r}[1-(R_{o}/R)^{3k}+a\sin (ft)]+2\unicode[STIX]{x1D6FE}R^{-1}[1-(R_{o}/R)^{3k-1}]\qquad \qquad \qquad \qquad & \displaystyle \nonumber\\ \displaystyle & \displaystyle +\;4\unicode[STIX]{x1D707}{\dot{R}}R^{-1}+\unicode[STIX]{x1D70F}_{y}[-2/\sqrt{3}+\sqrt{3}\ln (R_{\infty }/R)]\text{sign}({\dot{R}})+\unicode[STIX]{x1D70C}(R\ddot{R}+{\textstyle \frac{3}{2}}{\dot{R}}^{2})=0. & \displaystyle\end{eqnarray}$$

We note that this equation is almost identical to the equation that has been derived by Chan & Yang (Reference Chan and Yang1969); the only difference is a factor of 2 in front of the logarithm term.

Assuming small amplitude of pressure and radial oscillations (i.e. assuming $a\ll 1$ and $R=R_{o}+\unicode[STIX]{x1D700}R^{(1)}(t)$ with $\unicode[STIX]{x1D700}R^{(1)}(t)\ll R_{o}$ ), the above equation can be linearized, and the following expression can be derived:

(3.27) $$\begin{eqnarray}\displaystyle & & \displaystyle \ddot{R}^{(1)}(t)+\frac{4\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}R_{o}^{2}}{\dot{R}}^{(1)}(t)+\frac{[3kP_{r}+2(3k-1)\unicode[STIX]{x1D6FE}/R_{o}-2\sqrt{3}\unicode[STIX]{x1D70F}_{y}\text{sign}({\dot{R}}^{(1)})]}{\unicode[STIX]{x1D70C}R_{o}^{2}}R^{(1)}(t)\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{P_{r}}{\unicode[STIX]{x1D70C}R_{o}}\sin (ft).\end{eqnarray}$$

Thus, we have obtained an equation of the generic form of a forced harmonic oscillator with damping factor $4\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}R_{o}^{2}$ and natural frequency

(3.28) $$\begin{eqnarray}\displaystyle f_{n}=\sqrt{\frac{3kP_{r}+2(3k-1)\unicode[STIX]{x1D6FE}/R_{o}-2\sqrt{3}\unicode[STIX]{x1D70F}_{y}\text{sign}({\dot{R}})}{\unicode[STIX]{x1D70C}R_{o}^{2}}-\frac{8\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D70C}^{2}R_{o}^{4}}}. & & \displaystyle\end{eqnarray}$$

We should note that the presence of $\text{sign}({\dot{R}}^{(1)})$ in the numerator of the first term implies that the spring constant of the oscillator varies depending on the phase of the oscillation. The average natural frequency, though, matches that of a Newtonian liquid ( $\unicode[STIX]{x1D70F}_{y}=0$ ) for which (3.26) reduces to the well-known Rayleigh–Plesset equation. Thus, the natural frequency for our system is given by (see Brennen Reference Brennen2014)

(3.29) $$\begin{eqnarray}\displaystyle f_{n}=\sqrt{\frac{3kP_{r}+2(3k-1)\unicode[STIX]{x1D6FE}/R_{o}}{\unicode[STIX]{x1D70C}R_{o}^{2}}-\frac{8\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D70C}^{2}R_{o}^{4}}}, & & \displaystyle\end{eqnarray}$$

or in dimensionless form

(3.30) $$\begin{eqnarray}\displaystyle Sr_{n}=\frac{Ar}{\sqrt{Ar(3kp_{r}+2(3k-1)Bo^{-1})-8}}. & & \displaystyle\end{eqnarray}$$

As will be shown below, the mobilization of the bubble becomes maximized when the imposed frequency of the acoustic field, or Sr in dimensionless terms, is close to the natural frequency that is calculated by (3.29) or its equivalent (3.30).

When it comes to the flow of a viscoplastic material, it is always important to investigate the existence of an unyielded region and the position/shape of the possible yield surface. This limiting case can be used to draw some interesting conclusions regarding the flow of the viscoplastic material that surrounds the bubble. In the absence of any translation, the velocity potential is simply given by $\unicode[STIX]{x1D711}={\dot{R}}R^{2}/r$ , which can be introduced into (3.3) to evaluate the stress tensor, and subsequently its dimensional second invariant is given by

(3.31) $$\begin{eqnarray}\displaystyle \Vert \unicode[STIX]{x1D64F}\Vert =2\sqrt{3}\unicode[STIX]{x1D707}\left|\frac{R^{2}{\dot{R}}}{r^{3}}\right|+\unicode[STIX]{x1D70F}_{y}. & & \displaystyle\end{eqnarray}$$

According to the Von Mises criterion, an unyielded region exists in areas where the following inequality is satisfied: $\Vert \unicode[STIX]{x1D64F}\Vert \leqslant \unicode[STIX]{x1D70F}_{y}$ . When the bubble expands or contracts, i.e. ${\dot{R}}\neq 0$ , $\Vert \unicode[STIX]{x1D64F}\Vert >\unicode[STIX]{x1D70F}_{y}$ , and therefore it becomes evident that in the case of sphero-symmetric bubble oscillation a true unyielded region cannot exist except for the time instants that ${\dot{R}}$ is equal to zero, i.e. the time of maximum expansion or contraction of the bubble. However, since the yielding of the material is typically related to the buildup or breakdown of a specific structure which requires a finite amount of time, we conclude that practically unyielded regions do not exist in the case of a sphero-symmetric bubble oscillation. The non-existence of an unyielded region in this type of flow can be attributed to the incompressibility of the material. As the bubble expands sphero-symmetrically, there will always be flow in the radial direction as the bubble displaces the material while it expands or due to the fact that the liquid follows the contraction of the bubble to occupy the empty space.

5.1 Sphero-symmetric bubble oscillations: one-dimensional reduced-order model

We begin our study with a discussion about the predictions of the sphero-symmetric model that we developed in § 3. Note that for all simulations presented below, we have considered $R_{\infty }=10^{4}$ ; considering values one order of magnitude lower or higher than this value led to negligible differences. First, we examine in figure 2 the effect of the yield stress by varying the value of Bn. In this figure we plot the evolution of the dimensionless velocity of the bubble, $\text{d}z_{cm}/\text{d}t$ (see figure 2 a), along with the position of its centre of mass, $z_{cm}$ (see figure 2 b). It is already known for the case of steady bubble rise under a static pressure field that motion takes place below a critical value of Bn and the bubble rise velocity decreases with increasing viscoplasticity of the material because of the increasing effective viscosity. As shown in figure 2, a similar behaviour is encountered when the bubble is subjected to an acoustic field resulting in deceleration of the bubble with increasing yield stress. For $Bn=0.1$ , the bubble covers a distance equal to 0.3 times the bubble radius after approximately 10 cycles of pressure oscillation (at $t=61.53$ ), whereas in the case of a Newtonian liquid the same distance is covered at $t=26.75$ ; assuming that the imposed frequency is 200 Hz the corresponding time in dimensional terms would be 0.308 and 0.134 s, respectively.

Figure 2.

Evolution of (a) bubble rise velocity and (b) position of bubble centre of mass for different values of Bn. The other parameter values are given in table 1.

The presence of the acoustic field, though, causes volume oscillations of the bubble (see figure 3 a for the evolution of the bubble radius), due to its compressibility. The volume oscillation modifies the stress in the liquid surrounding the bubble and thus allows overcoming the yield stress of the material. Therefore, even though the present simplified model predicts that in the presence of a static pressure field the bubble will become entrapped for $Bn=0.175$ (see (3.26)), we find that in the presence of an acoustic field bubble motion is still possible beyond this critical value, i.e. for $Bn=0.18$ . For this value of Bn, the bubble rises during the expanding phase of the bubble whereas remains stagnant when the bubble contracts. Eventually, for even higher values of Bn, e.g. for $Bn=0.2$ , the model predicts that the rise velocity becomes equal to zero and thus the bubble is entrapped.

Figure 3.

Evolution of (a) bubble radius and (b) position of bubble centre of mass for different values of Bo. The other parameter values are given in table 1.

It is noteworthy that the effect of Bn on the amplitude of volume oscillations is negligible (not shown here for conciseness). However, as shown in figure 3(a), the amplitude of radial oscillations depends strongly on the surface tension of the liquid–air interface, and thus on the value of the Bond number, which acts as a restoring force to maintain the bubble shape. Higher values of Bo result in bubbles which are susceptible to larger volume oscillations leading to a slight increase of their mobility; the latter is not significant for the specific parameter values that we have considered (see figure 3 b).

The amplitude of volume oscillations also depends significantly on the amplitude of the imposed pressure oscillations, $a$ , far from the bubble (see figure 4 a). For small values of $a$ , the difference in the bubble volume is not significant between the two phases of the flow, i.e. contraction and expansion. On increasing the amplitude, $a$ , this symmetry breaks and the bubble undergoes substantial growth which is followed by limited contraction. For amplitude $a=0.7$ , i.e. 70 % of the reference pressure, during expansion the bubble radius increases by approximately 33 % which corresponds to an increase in bubble volume of ${\sim}235\,\%$ . On the other hand, during contraction the minimum bubble radius reaches ${\sim}88\,\%$ of its initial value, i.e. the volume decreases at most by ${\sim}32\%$ . As expected, increasing the amplitude of bubble oscillations leads to an increase of the mobility of the bubble (see figure 4 b).

Figure 4.

Evolution of (a) bubble radius and (b) position of bubble centre of mass for different amplitudes of pressure oscillation. The other parameter values are given in table 1.

In figure 5 we examine the effect of Ar, which is a measure of the importance of inertia in our system, on the bubble rise velocity. At this point, it is important to mention that the particular selection of our scaling poses an inherent difficulty in the interpretation of our results when varying Ar while keeping fixed the value of Sr. These two numbers are related to each other according to the following expression: $Ar=Sr^{2}f^{2}R_{o}/g$ . Therefore, assuming that the characteristic time is fixed, i.e. a constant imposed frequency, would imply that bubbles of different radii are examined with fixed Sr and increasing values of Ar. Since lengths are scaled with $R_{o}$ it would be difficult to evaluate the dependence of the actual bubble rise velocity on Ar. To overcome this problem, we introduce an alternative scaling for the length, using the viscous length scale $l_{v}=(\unicode[STIX]{x1D707}^{2}/\unicode[STIX]{x1D70C}^{2}g)^{1/3}$ , and for time, using the viscous time scale $t_{v}=(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}g^{2})^{1/3}$ , which are solely based on liquid properties. Assuming a liquid with density $2400~\text{kg}~\text{m}^{-3}$ and viscosity 1 Pa s, we get that $t_{v}=0.0163~\text{s}$ , $l_{v}=2.605~\text{mm}$ and $l_{v}/t_{v}=0.16~\text{m}~\text{s}^{-1}$ . We denote the new dimensionless position of the centre of mass and dimensionless time as $z_{cm}^{\prime }$ and $t^{\prime }$ , respectively; accordingly, the bubble rise velocity, $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ , has been scaled with $l_{v}/t_{v}=(\unicode[STIX]{x1D707}g/\unicode[STIX]{x1D70C})^{1/3}$ . The new scaling allows a direct comparison of the bubble rise velocity for different values of Ar. As depicted in figure 5, increasing the value of Ar results in a considerable speeding up of the bubble rising motion due to the increased effect of buoyancy. Note that in figure 5(a) $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ is plotted against $t$ to clearly depict the dependence of the rise velocity for all values of Ar for the same number of pressure cycles, whereas in figure 5(b) we prefer to plot $z_{cm}^{\prime }$ versus $t^{\prime }$ since from the slope of these curves it is possible to infer the rise velocity of the bubble, $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ . Moreover, as shown in figure 5(a), nonlinear effects become prominent for the highest value of Ar that we have considered, i.e. for $Ar=10$ .

Figure 5.

Evolution of (a) bubble rise velocity and (b) position of bubble centre of mass for different values of Ar; $t$ and $t^{\prime }$ denote the time scaled with $f^{-1}$ and $t_{v}$ , respectively. The other parameter values are given in table 1.

The effect of Sr, over a wide range of values, is examined in figure 6. For very small values of Sr, i.e. $Sr=0.0005$ , which correspond to high-frequency pressure oscillations, the bubble does not have enough time to adjust to the changing pressure field and therefore the amplitude of volume oscillations is reduced (see figure 6 a) while the rise velocity (see figure 6 b) soon approaches the steady rise velocity for a static pressure field $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }=(1/9)(1-(3/4)\unicode[STIX]{x1D712}Bn)Ar^{2/3}\approx 0.0102$ ; the latter expression has been derived using (3.24) assuming ${\dot{R}}=0$ , $\ddot{z}=0$ and ${\dot{z}}>0$ . The amplitude of radial oscillations increases with increasing Sr indicating that the resonance frequency of the system should be in the regime $0.001<Sr<0.01$ . For even higher values of Sr, i.e. $Sr\geqslant 0.1$ , an asymptotic regime is reached where the behaviour of the system is determined only by the amplitude of the imposed pressure oscillations.

Figure 6.

Evolution of (a) bubble radius and (b) bubble rise velocity for different values of Sr. The other parameter values are given in table 1.

Figure 7.

Evolution of (a) bubble rise velocity and (b) bubble radius for values of Sr which correspond to frequencies in the region of the resonance frequency. The other parameter values are given in table 1.

As discussed in § 3.1.2, Sr which corresponds to the natural frequency of a pulsating bubble in the absence of buoyancy effects is given by (3.30). For the parameter values that we consider as our base case, this equation predicts that resonance should take place for $Sr\approx 0.002185$ . Our numerical simulations indicate that this is a good approximation of the natural frequency even in our case where buoyancy is important. This is clearly demonstrated in figure 7 where we plot the evolution of the rise velocity, $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ , and bubble radius, $R$ , with time for three different values of Sr. Near resonance the volume oscillations become maximized affecting the mobility of the bubble in the viscoplastic material through two different mechanisms. The first mechanism is related to the increased effect of buoyancy due to the maximized bubble volume during the expansion phase that can be attained near resonance. The second mechanism is related to the rapid growth of the bubble size inducing high shear rates in the liquid surrounding the bubble which results in decreased effective viscosity of the material.

These two mechanisms are always present when a bubble is subjected to an acoustic field. We have already seen in figure 2 that even for conditions far from resonance the bubble is able to rise for values of Bn slightly higher than the critical value. However, the effect of these mechanisms on the stress field is significantly more near resonance, and, as shown in figure 8, our simplified model predicts that for $Sr=0.002185$ the range of Bn for which bubble motion is possible widens considerably, i.e. even for $Bn=0.25$ the bubble rises slowly whereas bubble entrapment is predicted for $Bn=0.175$ in the case of static pressure field.

Figure 8.

Evolution of the position of bubble centre of mass, $z_{cm}$ , for different values of Bn. The value of $Sr=0.002185$ corresponds to an imposed frequency close to the natural frequency of the system. The other parameter values are given in table 1.

5.2 Axisymmetric deformable bubble: two-dimensional model

So far, we have examined in detail the predictions of a simplified sphero-symmetric model. In this part of our work, we perform simulations using a more rigorous model that fully takes into account bubble deformability and the flow field surrounding the bubble, assuming axial symmetry only. The details of the model are presented in § 4. Note that for all the simulations that will be presented below, we have considered $R_{\infty }=50$ ; considering values higher than this value (e.g. $R_{\infty }=100$ ) led to negligible differences (see appendix B). Moreover, the exponential factor, $N$ , which arises in (4.3) was considered to be equal to $10^{4}$ for all the simulations that will be presented below, in line with our previous work (Tsamopoulos et al. Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008); considering values higher than this value (e.g. $N=5\times 10^{4}$ ) led to negligible differences in the dynamics of bubble motion (see appendix B).

Below, we determine the yield surface as the location where $\Vert \unicode[STIX]{x1D64F}\Vert =Bn=0.1$ , employing the Von Mises criterion. It is important to keep in mind, though, that according to the Papanastasiou model, even for $\Vert \unicode[STIX]{x1D64F}\Vert <Bn$ , finite motion of the material is allowed, albeit with very large viscosity, and therefore using this model it is not possible to predict a true zone of unyielded material. In fact, in our case, given the radial character of the flow and the fact that the liquid is considered to be incompressible, we expect that the radial velocity will acquire small but finite values even very far from the bubble and therefore we expect a material following the Papanastasiou model, and thus exhibiting no elasticity, should actually yield. On the other hand, in the case of elastoviscoplastic materials, the material could exhibit true unyielded zones, since in this case the material could also experience an elastic solid deformation. Therefore, in the discussion that follows, we will still refer to a yield surface but with the understanding that this does not correspond to a true boundary between truly yielded and unyielded domains. It will rather be used as a convenient way to study the evolution and visualize regions with relatively low viscosity (which hereafter we will refer to as yielded regions) and regions with much higher viscosity (which hereafter we will refer to as unyielded regions). It should also be noted that the position of this pseudo-yield surface is quite sensitive to the value of $N$ but its evolution at different time instants is qualitatively similar.

Figure 9 shows the evolution of the bubble shape and the yielded (red) and unyielded (grey) regions for the ‘base’ case while we depict four different isolines of the second invariant of the stress tensor with values 0.1, 0.2, 0.3 and 1. We observe that the bubble at all times acquires a nearly spherical shape while the stresses monotonically decrease away from the bubble resulting in the presence of unyielded material there; according to the Von Mises criterion the yield surface arises where $\Vert \unicode[STIX]{x1D64F}\Vert =Bn=0.1$ . When the bubble rate of expansion (at $t=31.42$ and $t=37.70$ ) or contraction (at $t=34.56$ ) is at its peak, the stress field in the material surrounding the bubble is maximized leading thus to the maximum size of the yielded regions. At these time instants the yield surface appears to be nearly symmetric. For time instants where the bubble volume has reached a maximum (at $t=32.99$ ) or minimum (at $t=36.13$ ), the radial velocity decreases (see figure 10) and thus the stress field reduces significantly resulting in a decrease of the size of the yielded regions. Another characteristic of this system is that we do not observe at any time instant the formation of an unyielded region near the equatorial plane. Such formations have been reported by Tsamopoulos et al. (Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008) in the case of a static pressure field, due to the locally zero shear and normal stresses. However, in our case, the strong hoop stresses that the liquid experiences adjacent to the liquid–air interface due to the radial oscillations of the bubble do not allow such formations. Note that the destruction of unyielded zones around the bubble equator has also been reported in a freely oscillating bubble by Tripathi et al. (Reference Tripathi, Sahu, Karapetsas and Matar2015).

Figure 9.

Yielded (red) and unyielded (grey) domains for the ‘base’ case at different time instants: (a) $t=31.42$ , (b) $t=32.99$ , (c) $t=34.56$ , (d) $t=36.13$ , (e) $t=37.70$ . Four isolines of the second invariant of the stress tensor with values 0.1, 0.2, 0.3 and 1 are shown. The yield surface is evaluated as the location where $\Vert \unicode[STIX]{x1D70F}\Vert =Bn$ with $N=10^{4}$ . The other parameter values are given in table 1.

Figure 10.

Contour plots of the velocity field around the bubble: (a) $t=31.42$ , (b) $t=32.99$ , (c) $t=34.56$ , (d) $t=36.13$ . We present the radial component in spherical coordinates of the velocity; the origin of the coordinate system is taken at the bubble centre of mass. The other parameter values are given in table 1.

The velocity field around the bubble is depicted in figure 10. In this figure, we assume a spherical coordinate system, the origin of which is placed at the bubble centre of mass, and plot the radial velocity component for the same time instants as in figure 9; the velocity profile for $t=37.70$ is not shown since it looks very similar to the one for $t=31.42$ due to the periodic character of the flow. Clearly, the velocity field around the bubble is non-symmetric which can be attributed to its rising motion; a symmetric profile would be expected in the absence of buoyancy. During the phase of expansion, the bubble increases its pace, due to the continuously increasing effect of buoyancy. As it rises, it encounters unyielded material at the top, whereas the material below has been yielded at earlier times. As a result, the bubble, as it expands, encounters an area with high effective viscosity at the top and low effective viscosity at the bottom which leads in turn to preferential expansion towards the bottom part of the bubble. On the other hand, at $t=34.56$ (for which the bubble reaches its maximum rate of contraction), the bubble at a previous time instant creates an unyielded region around it. Due to its contraction though at a next time instant the bubble acquires smaller volume and thus decelerates covering less distance leaving more yielded material at the top and resulting again in an asymmetry of the effective viscosity of the material surrounding bubble. Thus, the fore–aft symmetry, which would be expected in the case of a spherical bubble rising steadily in the limit of low values of Ar, breaks and the flow field becomes non-symmetric for the case of a pulsating bubble.

Figure 11.

Evolution of (a) bubble rise velocity and (b) position of bubble centre of mass for different values of Bn. The other parameter values are given in table 1.

Next, we examine the effect of Bn on the rise velocity of the bubble. In figure 11 we plot the evolution of the dimensionless velocity of the bubble, $\text{d}z_{cm}/\text{d}t$ (figure 11 a), along with the position of its centre of mass, $z_{cm}$ (figure 11 b). For the case of a Newtonian liquid we find a behaviour similar to that described by our sphero-symmetric model, i.e. a sinusoidal dependence of the rise velocity that follows the pressure oscillations. We should note, though, that although the sphero-symmetric model is able to qualitatively describe the flow, it provides an underestimation of bubble rise velocity. In the case of a viscoplastic material, the picture changes radically since the sphero-symmetric model fails to capture certain qualitative characteristics of the flow. As shown in figure 11(a), for each cycle of the acoustic pressure field, the bubble velocity exhibits a double minimum at two different time instants which coincide with the times where the bubble reaches its maximum and minimum volume. As explained above, at these time instants the stresses decrease significantly resulting in the decrease of the size of the yielded region around the bubble causing a severe deceleration. Possibly, the failure of the sphero-symmetric model to capture this behaviour can be attributed to our assumption of a simplified flow field, as given by (3.5). The latter also leads to an underestimation of the second invariant of the rate of strain tensor which in turn results in an underestimation of the stress in the surrounding material. Therefore, the simplified sphero-symmetric model also provides an underestimation of the critical value of Bn for bubble entrapment since we clearly see in figure 11 that the bubble rises for values of Bn that can be as high as 0.5; note that according to Tsamopoulos et al. (Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008) the critical value for Bn for steady bubble rise has been evaluated to be approximately equal to 0.143.

Figure 12.

Evolution of (a) bubble rise velocity and (b) bubble radius for different values of Bo. The other parameter values are given in table 1.

In figure 12, we show how the bubble rise velocity and the bubble radius are affected by the value of the Bond number. As already noted in the discussion of figure 3, increasing the value of Bo results in volume oscillations with larger amplitude; similar evolution of the bubble radius is also found here. We observe negligible differences for the two highest values that are considered in this figure, i.e. for $Bo=0.01$ and $Bo=0.1$ , while a small decrease in the rise velocity is observed for $Bo=0.001$ . For the lowest value of Bo (0.0001), though, we notice that besides the considerable decrease of the overall rise velocity, its average value is significantly higher during the expansion phase of the bubble and lower during the compression phase. As noted above, even for the highest values of Bo that we have examined in this study the bubble retains its spherical shape. To evaluate whether any deviations from its spherical shape are important we may consider a spherical coordinate system with origin the centre of mass of the bubble and decompose its shape in spherical harmonics. Thus, we may evaluate the Legendre coefficients using the following expression:

(5.1) $$\begin{eqnarray}c_{i}=\frac{2k+1}{2}\int _{0}^{\unicode[STIX]{x03C0}}f(\unicode[STIX]{x1D703})P_{k}(\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ denotes the azimuthal angle and $f(\unicode[STIX]{x1D703})$ is the radial distance to the bubble surface from the centre of the local spherical coordinate system. The evolution of the Legendre decomposition of the shape of the bubble is shown in figure 13. As expected, we find that the bubble shape is dominated by volume oscillations and $c_{2}$ , $c_{3}$ and $c_{4}$ remain small for all times while higher modes are negligible.

Figure 13.

Variation with time of selected Legendre coefficients of the shape of the bubble for the ‘base’ case. The decomposition is done in a spherical coordinate system located on the instantaneous centre of mass of the bubble. The other parameter values are given in table 1.

The effects of the amplitude of the imposed pressure oscillations, $a$ , on the bubble rise velocity and position of the centre of mass are examined in figure 14. We observe in figure 14(a) that for large values of $a$ the local minimum in the rise velocity, which corresponds to the time where the bubble volume reaches its maximum point, tends to disappear. This is due to the fact that even though the rate of change of the bubble radius is very small, the bubble has a large enough size so that only the effect of buoyancy is sufficient to overcome the yield stress of the material and thus maintain its motion. On the other hand, the local minimum which corresponds to the point of minimum bubble volume also increases due to the higher level of stresses that the fluid experiences with increasing amplitude of the pressure oscillations also maintaining the overall stress in the material above its yield stress. Moreover, as discussed in § 5.1, increasing the amplitude, $a$ , leads to a break of symmetry with the bubble undergoing substantial growth during the expansion phase which is followed by limited contraction; the evolution of the bubble radius looks very similar to the one shown in figure 4(a) and it is not shown here for conciseness. This is also clearly reflected in the evolution of the bubble rise velocity which is shown in figure 14(a). We have also examined the effect of different values of the reference pressure while keeping $a$ constant. Such a situation would arise, for example, if the bubble was located at larger depths; however, we found that varying $P_{r}$ between 2000 and 10 000 led to negligible differences in the dynamics of the bubble (not shown here for brevity).

Figure 14.

Evolution of (a) bubble rise velocity and (b) position of bubble centre of mass for different amplitudes of pressure oscillation. The other parameter values are given in table 1.

To examine the effect of Ar and Sr, we again employ the alternative scaling for the length and time using the viscous scales, as done in § 5.1, and derive the new dimensionless position of the centre of mass, $z_{cm}^{\prime }$ , time, $t^{\prime }$ , and bubble rise velocity, $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ . As shown in figure 15, on increasing the value of Ar the rise velocity increases due to the increased effect of buoyancy. For low values of Ar, the average rise velocities of the compression and expansion phases are similar since the bubble has adequate time to adjust to the changing pressure field. With increasing values of Ar, the effect of buoyancy becomes more important and the flow becomes dominated by the rising motion of the bubble causing a substantial area around the bubble to remain yielded even at times where the rate of expansion or compression is very small. As noted above, increasing the value of Ar, while keeping fixed the value of Sr also implies that we consider bubbles with increasing radius which also explains why the amplitude of rise velocity oscillations increases, i.e. due to the increased size of the amplitude of volume oscillations. The asymmetry between the compression and expansion phase is also apparent in this figure, especially for the case of moderate values of Ar, i.e. for $Ar=1$ , while for higher value of the local minimum corresponding to the point of maximum volume it disappears.

Figure 15.

Evolution of the bubble rise velocity for different values of Ar. The other parameter values are given in table 1.

Finally, we examine the effect of excitation frequency in figure 16. As discussed in § 5.1, in order to maximize the effect of the acoustic field on the bubble motion and to minimize the energy damping, it is reasonable to choose a frequency in the region of the resonance frequency of the system. Our axisymmetric simulations indicate that, in accordance with the prediction of the sphero-symmetric model and (3.29), resonance takes place in the region of $Sr=0.0022$ (see figure 16 c). To examine how effective acoustic excitation can be in mobilizing a bubble in viscoplastic liquids, we have performed simulations at resonance frequency for a wide range of Bn as shown in figure 16(d). Interestingly we find that motion of the bubble takes place for the highest value of Bn that we have examined, i.e. for $Bn=20$ , a value which is more than one order of magnitude higher than the critical Bn for bubble entrapment in a static viscoplastic liquid. It is therefore evident that acoustic excitation can be an extremely efficient way to mobilize bubbles inside a viscoplastic material. For oscillations of very high and very low frequency (see figure 16 a,b) the steady rise velocity approaches that for a static pressure field. In the former case, the reason is because the viscoplastic liquid does not have enough time to adjust to the rapidly changing pressure field, whereas in the latter case the change in the pressure field is so slow that the liquid experiences a reduced level of stress due to the radial oscillation and the stress field is dominated by the rising motion due to buoyancy.

Figure 16.

Evolution of bubble rise velocity for (a) low values of Sr, (b) high values of Sr and (c) near resonance. (d) Evolution of the bubble centre of mass for different values of Bn near resonance. The other parameter values are given in table 1.

So far, we have performed a thorough parametric study and have investigated in detail the effect of all dimensionless groups of our system. However, in experiments it is often difficult to control each parameter separately while keeping the rest fixed. An important parameter, though, which can be easily fixed is the size of the bubble at equilibrium. Therefore, it is also of value to examine the effect of the bubble radius on the particular system, while keeping the liquid properties fixed. To this end, we define the following dimensionless groups: the Morton number, $Mo=g\unicode[STIX]{x1D707}^{4}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}^{3}$ , the dimensionless frequency, $\unicode[STIX]{x1D6FA}=f(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}g^{2})^{1/3}$ , the yield stress parameter, $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D70F}_{y}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}$ , and dimensionless reference pressure, $\unicode[STIX]{x1D6F1}=P_{r}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}$ . These dimensionless groups depend solely on the liquid properties and for our ‘base’ case acquire the following values: $Mo=32.7$ , $\unicode[STIX]{x1D6FA}=3.26$ , $\unicode[STIX]{x1D6F6}=0.16$ and $\unicode[STIX]{x1D6F1}=2000$ . Keeping the values of these groups constant and varying only the bubble radius, we end up with the values of the other dimensionless groups that are shown in table 2.

Figure 17.

Evolution of the position of bubble centre of mass for bubbles of different size. We consider a liquid with $Mo=g\unicode[STIX]{x1D707}^{4}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}^{3}=32.7$ and $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D70F}_{y}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}=0.16$ . The dimensionless reference pressure and angular frequency are $\unicode[STIX]{x1D6F1}=P_{r}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}=2000$ and $\unicode[STIX]{x1D6FA}=f(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}g^{2})^{1/3}=3.26$ , respectively.

Table 2.

Typical values of the dimensionless groups for bubbles of different size for a liquid with $Mo=g\unicode[STIX]{x1D707}^{4}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}^{3}=32.7$ and $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D70F}_{y}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}=0.16$ . The dimensionless reference pressure and angular frequency are $\unicode[STIX]{x1D6F1}=P_{r}/(\unicode[STIX]{x1D70C}g^{2}\unicode[STIX]{x1D707}^{2})^{1/3}=2000$ and $\unicode[STIX]{x1D6FA}=f(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}g^{2})^{1/3}=3.26$ , respectively, and $l_{v}=(\unicode[STIX]{x1D707}^{2}/\unicode[STIX]{x1D70C}^{2}g)^{1/3}$ denotes the viscous length scale.

In figure 17 we plot the evolution of the bubble rise velocity, $\text{d}z_{cm}^{\prime }/\text{d}t^{\prime }$ , for six different bubble radii and it can be clearly seen that the viscoplastic character of the fluid is most prominent for bubbles of smaller radii where we observe the two local minima of the rise velocity that arise inside each pressure cycle. For bubble of large radius, the flow is dominated by buoyancy effects and the system tends to a more Newtonian behaviour.

Figure 18.

Yield regions and bubble shapes: (a) $t=62.83$ , (b) $t=63.62$ , (c) $t=64.40$ , (d) $t=65.19$ , (e) $t=65.97$ , (f) $t=66.76$ , (g) $t=67.54$ , (h) $t=68.33$ . The yield surface is evaluated as the location where $\Vert \unicode[STIX]{x1D70F}\Vert =Bn$ with $N=10^{4}$ . The dimensionless parameters are the same as for figure 17 and we consider a bubble with nominal radius at equilibrium equal to 7.5 mm, i.e. $R_{o}/l_{v}=2.878$ .

In the latter case, the bubble also is more deformable and acquires the shape of a spherical cap as shown in figure 18 where we also plot the yielded (red) and unyielded (grey) regions of the viscoplastic liquid. During a period of pressure oscillation the bubble grows and contracts retaining its deformation which also affects the shape of the yield surface, although the size of the indentation decreases as the bubble reduces in size tending to recover its oblate shape. It is also noted that an island of unyielded material arises either above the bubble, during the expansion phase, or below the bubble, during the contraction phase. The size of this island grows at times where the rate of expansion or contraction decreases, i.e. at times where the bubble volume reaches its minimum or maximum point, and eventually merges with the outer unyielded domain which grows in size at the same time.

We have also examined the case of a liquid with the same properties and bubble of the same size as in figure 18 but for a frequency that corresponds to $Sr=0.2$ which is close to the natural frequency of the system ( $Sr_{n}\approx 0.09$ ). As expected, the deformation of the bubble is much more marked in this case and even leads to rupture at the end of the eleventh cycle (see figure 19). Moreover, the large deformations of the bubble and high amplitude of the volume oscillations that the bubble experiences cause a significant increase in the size of the yielded region that surrounds the bubble; the yield surface, though, is within the region of our computational domain at all times.

Figure 19.

Yield regions and bubble shapes for $\unicode[STIX]{x1D6FA}=0.738$ : (a) $t=31.42$ , (b) $t=32.20$ , (c) $t=32.99$ , (d) $t=33.77$ , (e) $t=34.56$ , (f) $t=35.34$ , (g) $t=35.80$ . The yield surface is evaluated as the location where $\Vert \unicode[STIX]{x1D70F}\Vert =Bn$ with $N=10^{4}$ . The other dimensionless parameters are the same as for figure 18.