Oscillations of small bubbles and medium yielding in elastoviscoplastic fluids

We investigate the radial oscillations of small gas bubbles trapped in yield-stress fluid and driven by an acoustic pressure field. We model the rheological behavior of the yield-stress fluid using the recently developed elasto-visco-plastic (EVP) constitutive equation that takes into account the elastic and visco-plastic deformations of the material [P. Saramito, J. NonNewton. Fluid Mech. 158 (1-3) (2009) pp.154-161]. Assuming that the bubble remains spherical during the pressure driving, we reduce the problem to a set of ODEs and an integrodifferential equation, which we solve numerically for the case of two yield-stress fluids, a sot Carbopol gel and a stiffer Kaolin suspension. We find that, depending on the amplitude and frequency of the pressure field, the radial oscillations of the bubble produce elastic stresses that may or may not suffice to yield the surrounding material. We evaluate the critical amplitude of the acoustic pressure required to achieve yielding and we find a good agreement between numerical simulations and an analytical formula derived under the assumption of linear deformations. Finally, we examine the bubble oscillation amplitude for a very wide range of applied pressures both below and above the critical value to assess the impact of yielding on the bubble dynamics. This analysis could be used to identify a signature of yielding in experiments where the radial dynamics of a bubble is measured. More generally, these results can be used to rationalize the optimal conditions for pressure-induced bubble release from yield-stress fluids, which is relevant to various biomedical and industrial applications, including oil industry and food processing.

Significant bubble entrapment is desirable in food engineering to improve texture and slow down melting of ice cream [12] and improve the perceived crunchiness of snacks [13]. A small, controlled amount of bubbles is also allowed in concrete to improve workability and freeze-thaw resistance, at the expense of reduced compressive strength and concrete blisters [14]. Bubble entrapment in yield-stress fluids can be an undesirable byproduct of fluid processing: the presence of air pockets in polydimethylsiloxane (PDMS) used for microfluidics applications can severely reduce its quality and transparency [15]; bubbles trapped in fluids used in the oil industry can result in undesired permeable slurries [16] or lead to explosions in drilling mud, which may delay production and potentially inflict a huge burden on the ecosystem [17,18]. Air bubbles induce bacterial contamination and poor final appearance in personal care products [19].
Complete bubble removal can be achieved by centrifuging the sample or by using a vacuum pump to inflate bubbles. Recent experiments have confirmed that mechanical agitation is an effective method of removing gas voids from granular materials for which the yield stress is a consequence of the frictional contact network formed by the microscopic phase [20]. Stein and Buggisch [21] and Karapetsas et al. [22] showed that driving a bubble into volumetric oscillations using a dynamic external pressure can generate sufficiently large deformation and mechanical stresses to locally yield the material, thus promoting bubble release from yield-stress fluids. The same mechanism has been exploited by Iwata et al. [23,24] to enhance bubble removal from highly viscous shear-thinning and viscoelastic fluids.
In contrast, the literature on bubble dynamics in yield-stress fluids is very limited and focuses almost exclusively on the problem of bubble rise due to buoyancy. In such fluids, rising motion results from the interplay between the gravity, bubble shape, position of the yield surface, and rheology of the material both below and above the yield stress. Experiments require great care to suppress internal stresses and achieve repeatability [11,51]. Numerical investigation of the problem proves to be equally challenging. The existing works model the material rheology using either Bingham or Herschel-Bulkley constitutive equations [6,7,9,10,21,52]. Both of these models predict a discontinuity of the viscosity at the yield surface, the location of which is unknown in flows that are two dimensional (2D), three-dimensional (3D), or time dependent. To avoid this problem and make numerical solutions feasible, these equations are either regularized (e.g., [53]), which unfortunately reduces the solid region of the material to a liquid with very large viscosity, or solved via the augmented Lagrangian method [9], which converges extremely slowly, but recently has been improved [54]. More importantly, the Bingham and Herschel-Bulkley models implicitly assume that the unyielded material cannot deform, even if the yield surface may adjust to flow, especially in time-dependent problems, and leave undetermined the stress field there. While this problem is not critical for the case of bubble rise, bubble oscillation prescribes a nonzero strain field everywhere in the fluid. This leads to two unphysical results: first, it implies that any finite oscillation amplitude causes yielding of all the material; second, it means that the stress applied by the bubble is everywhere above the yield stress, and is therefore infinite. The former issue contrasts with the experimental findings of Stein and Buggisch [21], who report that a finite oscillation amplitude is required to achieve yielding. This issue is not discussed in Ref. [52] and is circumvented in Ref. [21] by prescribing the deformation field and thus the dynamics of the yield surface. Recently, Karapetsas et al. [22] performed a detailed theoretical and numerical analysis of a bubble rising in a Bingham fluid and driven by an acoustic field into volumetric oscillations. The authors developed a simplified 1D spherosymmetric model and also performed detailed numerical simulations that take into account axisymmetric deformations of the gas-liquid interface. Their results confirm that 073301-2 a Bingham material is yielded everywhere during the oscillations of the bubble, that a yield surface cannot be defined for this type of flow, and that an acoustic field promotes the release of bubbles that would be trapped by the yield stress otherwise.
In this paper, we investigate the radial oscillations of a microbubble trapped in an elastic yield-stress fluid and driven by an external pressure field. We focus on the case of ultrasonic fields in the frequency range ∼10−100 kHz, which is relevant to industrial equipment. Small bubbles with radii ∼30−300 μm are resonant in this frequency range, that is, they are efficiently excited into volumetric oscillations by the ultrasound field. We perform numerical simulations employing a generalized Rayleigh-Plesset equation and a recently developed elastoviscoplastic (EVP) constitutive equation that takes into account elastic and viscoplastic deformations of the material [55]. Employing this model, we resolve the conceptual difficulties introduced in Refs. [21,52] by the choice of the Bingham model. Using numerical simulations and an approximate linear theory, we evaluate the critical acoustic pressure required to yield the material and compute the dynamics of the yield surface. Finally, we explore the impact of yielding on the radial oscillations of the bubble. Our results represent a step towards the investigation of pressure-induced bubble release from yield-stress fluids and could potentially be used to identify the signature of yielding in experiments. The theoretical and numerical framework developed in this paper can be applied to validate constitutive equations for yield-stress fluids by comparing with experiments performed under controlled extensional deformation imparted via bubble oscillation.

II. EQUATIONS GOVERNING THE BUBBLE DYNAMICS
We consider a bubble of equilibrium radius R 0 suspended in an incompressible yield-stress fluid. The bubble is driven by a time-dependent pressure p ∞ (t ) imposed far from the bubble. We assume that the Cauchy stress tensor T is given by where p denotes the pressure, v is the velocity field, η s is the shear viscosity of the solvent, and τ is a nontraceless and non-Newtonian contribution to the total stress. In the present work, we use the elastoviscoplastic model developed by Saramito [55] that gives excellent predictions when compared to experimental results in shear [56] and for the case of a sedimenting sphere [57]. This model predicts a neo-Hookean elastic behavior before yielding and a viscoelastic behavior afterwards, with an instantaneous transition from solid-to liquidlike behavior when the second invariant of the deviatoric part of τ is larger than the yield stress. The evolution of τ is governed by where G is the elastic modulus of the material and the symbol above a tensor X denotes its upper convected derivative, defined as with D/Dt the material derivative. In Eq. (2), |τ d | denotes the square root of the second invariant of the deviatoric part of the stress tensor τ, and K is the consistency parameter of the yielded phase with n its power-law index. For stresses |τ d | < τ y , the material is unyielded and experiences no viscoplastic deformation; for |τ d | > τ y , the material behaves as a viscoelastic fluid, thus undergoing both elastic and viscoplastic deformation. In the limit G → ∞, the constitutive model given by Eq.  [58,59]. However, recent experiments considering the extensional deformation of yield-stress fluids [59,60] suggest that the third invariant of τ d may also be required in the yielding criterion. Given the small size of the examined bubbles (∼100 μm), the relevant Bond number is less than 10 −3 ; hence we assume that the center of volume of the bubble is stationary. We define a spherical coordinate system with its origin at the center of the bubble, with r, θ , and φ the radial, azimuthal, and polar coordinates, respectively. Under the assumption that the bubble undergoes spherical oscillations only, the problem is spherically symmetric, implying v = v r (r)r, (4a) The incompressibility of the material implies that the radial velocity v r is related to the timedependent radius of the bubble, R, through In the case ofṘ > 0, the bubble expands and the material undergoes a nonuniform spherosymmetric compression. A nonuniform spherosymmetric extension is applied to the medium during bubble compressionṘ < 0. In this deformation field, phenomena characteristic of shear deformation such as wall slip and shear banding are not expected. We assume that the pressure at infinity changes due to the acoustic driving as p ∞ (t ). The time evolution of the bubble radius is governed by the generalized Rayleigh-Plesset equation that is obtained by integrating the radial component of the momentum balance from R to infinity [61]: where ρ is the density of the medium and we have assumed that the stress tensor τ vanishes at infinity because the rate of strain goes to zero far from the bubble. In Eq. (6), p(R) and τ rr (R) are the pressure and the radial component of τ evaluated at the surface of the bubble, respectively. These quantities are related through the normal stress balance, where γ is the surface tension of the interface between the gas and the yield-stress fluid. In principle, the surface tension between the gas and the yield-stress fluid could depend on the rheological state of the material. In this paper, we assume that the surface tension is the same, regardless of the stresses inside the yield-stress fluid. We denote the pressure inside the bubble with p gas and we neglect the viscous stresses in the gas phase. We assume that the bubble undergoes isothermal compression, p gas = (p 0 + 2γ /R 0 )(R 0 /R) 3 , and that the driving pressure is periodic, p ∞ (t ) = p 0 + p sin(ω t ), with the equilibrium pressure p 0 , the driving pressure amplitude p, and the angular frequency ω.
The dimensionless forms of Eqs. (6) and (2) are obtained introducing the following dimensionless quantities, noted with stars * :  Substituting the pressure p(R) in Eq. (6) and dropping the stars, we obtain De De In Eqs. (9) and (10), we have introduced the relevant dimensionless numbers. The ratio of the solvent to the viscosity of the yield-stress material is given by α = η s ω 1−n /K. The Bingham number, Bn = τ y /Kω n , denotes the ratio of the yield stress to the material viscous stress. The Reynolds number, Re = ρR 2 0 ω 2−n /K, expresses the relative importance of inertial to viscous stresses. The dimensionless number, P stat = R 0 p 0 /2γ , gives the ratio of the static pressure to surface tension. The Deborah number, De = Kω n /G, expresses the ratio of the viscoelastic relaxation time to the characteristic flow timescale. The dimensionless ambient pressure is given by P 0 = p 0 /Kω n and the dimensionless pressure amplitude, P = p/Kω n , compares the amplitude of acoustic pressure to the viscous stresses.
The values of the constitutive parameters for yield-stress fluids can vary over a wide range. In this work, we consider two yield-stress fluids: a soft Carbopol gel and a stiff kaolin suspension. The physical properties of the two materials are taken from the measurements of Lacaze et al. [62] and Luu and Forterre [63] and are summarized in Table I. Lacaze et al. used a Carbopol 940 gel at a concentration of 0.1% in weight. The kaolin used by Luu and Forterre is a colloidal suspension of clay in water at 55% in weight, which was supplied by Imerys Ceramics France. Both yieldstress fluids display a pronounced shear-thinning response. The large difference in the elastic moduli between the two fluids (2.0×10 5 Pa for kaolin versus 81.5 Pa for Carbopol) allows us to explore the effect of the ratio τ y /G for two materials with yield stress that are not very small (91 Pa for kaolin versus 8.6 Pa for Carbopol). Since the solvent used in the Carbopol gel and in the kaolin suspension is water, we assume that the Newtonian viscosity in Eq. (1) is given by η s = 0.001 Pa s. This makes α small; typically less than 0.01. In the Appendix, we report the rheological response predicted by the EVP constitutive equation for the two yield-stress fluids. Finally, throughout this paper, we fix the ambient pressure p 0 = 1.13×10 5 Pa. We neglect residual stresses that are potentially present in the material after its preparation and that decay over long timescales [51,64]. Thus, we assume equilibrium initial conditions: R(0) = R 0 ,Ṙ(0) = 0, and τ rr (r, 0) = τ θθ (r, 0) = 0.

III. NUMERICAL APPROACH
The solution of the partial differential equations (9) and (10) in the present form is complicated by the motion of the bubble surface that makes the domain time dependent. To avoid this problem and immobilize the boundary, we follow previous works [38,65] and transform the radial coordinate 073301-5 into the Lagrangian coordinate, This coordinate transformation reduces Eqs. (9) and (10) to a system of first-order ordinary differential and integrodifferential equations defined on y, where the surface of the bubble is given by the point y = 0: Re De We discretize the coordinate y into a set of N points, y i = y 1 , . . . , y N , with the last point describing the conditions far from the bubble. Following Kafiabad and Sadeghy [36], we discretize the spatial integral in Eq. (13) using a Gauss-Laguerre quadrature method, which is suitable for integrals of the type To apply the Gauss-Laguerre method, we first rewrite the spatial integral in Eq. (13) , and then we approximate I with the sum where w i are the weight factors defined as is the Laguerre polynomial of order N + 1 evaluated at the grid point y i . The positions of the grid points are given by the zeros of the Laguerre polynomial of the order of N and are found by solving the implicit equation L N (y i ) = 0 for y i . This procedure automatically divides the domain in elements of different size, with smaller elements near the surface of the bubble where larger gradients are expected. The resulting system of ordinary differential equations is solved using a fourth-order Runge-Kutta implicit scheme with a variable time-step size giving a maximum relative error of 10 −9 . We found that the choice N = 150 yields numerically convergent results and guarantees that the far-field conditions are met at the last point y N . Numerical simulations performed with more refined grids, larger domains, and smaller error threshold gave indistinguishable results.

Validation of the code
We validate the model and the numerical implementation by studying the dynamics of a bubble in the linear regime. Under the assumption of linearity, the material behaves as a Kelvin-Voigt viscoelastic solid. For small driving amplitudes, the radial oscillations of the bubble can be expressed as , where x(t ) follows the same dynamics as the forcing, x(t ) = ( R/R 0 )sin(ω t + φ), with a phase shift φ. In this regime, Eqs. (9) and (10) reduce to a damped harmonic oscillator driven by the external pressure [66], The damping coefficient β and the resonance frequency ω 0 are given by [47,48] In the linear regime, the medium behaves as a Kelvin-Voigt material, and thus the yield-stress does not enter in the expression for the damping coefficient and the resonance frequency. The solution of Eq. (16) gives the amplitude of oscillation, R, and the phase, φ, as a function of the frequency, We test the numerical solution of Eqs. (12)- (14) against the predictions of Eq. (18). To avoid transient effects, we run simulations for 2000 periods and we compute R as the maximum radial excursion over the last period. In Fig. 1, we compare the amplitude of the radial oscillations predicted by the linear theory given by Eq. (18) with that obtained from the numerical solution. We considered the case of a bubble with equilibrium radius R 0 = 100 μm, which is of the size used in the experiments performed by Jamburidze et al. [48], driven at p = 100 Pa in the Carbopol gel [ Fig. 1(b)] and at p = 1 Pa in the kaolin suspension [ Fig. 1(a)]. Such a small acoustic pressure is required to ensure that the dynamics remains in the linear regime. One might question the need for such a smaller p in kaolin as opposed to Carbopol resulting in R/R less than 10 −3 in order to remain in the linear regime. The reason is that nonlinearity in this material is not induced by an increased amplitude of the radial oscillations, but by its yielding at smaller radial deformations caused by its larger elastic modulus; see related discussion in Sec. IV. In Fig. 1, the frequency at which the bubble experiences the largest radial excursion is approximately given by the resonance frequency ω 0 because the damping coefficient is small. The numerical and analytical solutions show excellent agreement for all the frequencies investigated, thus showing that the numerical implementation of the EVP model correctly reduces to the a Kelvin-Voigt model for small deformations.

A. Analysis of bubble dynamics at resonance
As a result of the harmonic change in pressure due to the acoustic driving, the bubble undergoes periodic compression and expansion. The dynamics of the bubble radius is strongly dependent on 073301-7 the amplitude of the pressure applied by the acoustic field and, for sufficiently small p, it is given by the same harmonic function of the forcing.
To explore the transition to the nonlinear regime, we report in Fig. 2(a) the dynamics of a 100 μm bubble driven at its resonance frequency, ω = ω 0 , in the Carbopol gel. We discard the initial transient response of the bubble and we focus on its dynamics sufficiently far from the initial condition. As the acoustic pressure is increased, the amplitude of the radial excursion becomes larger and the dynamics deviate from the single harmonic response expected in the linear regime. At the largest driving pressure, the bubble spends more time in its expanded than its contracted state because (i) in the former state, the increased liquid inertia decreases its acceleration [67], and (ii) the bubble pressure varies as P g ∼ R −3 , which changes much faster when the bubble radius is minimized. Both observations lead to an effect equivalent to the added mass effect in a translating bubble. To obtain a quantitative insight on the transition from the linear to the nonlinear regime, we plot in Fig. 2(b) the amplitude of the Fourier modes |R n | of R(t )/R 0 . Figure 2(b) shows that the first harmonic is the mode with the largest amplitude for all driving investigated. In the case of p = 0.1 k Pa, the amplitude of the modes higher than 1 decays very fast with the mode number n. Conversely, in the case of p = 2.5 k Pa, the amplitude of the high-order modes decays slower with n and multiple harmonics play a role in the radial response. The significant coupling between different modes |R n | is a signature of the nonlinear dynamics of the bubble at large driving pressures.
The radial dynamics of a 100 μm bubble driven at its resonance frequency in the kaolin suspension is reported in Fig. 3(a). In contrast to the case of a bubble driven in the Carbopol gel, Fig. 3(a) suggests that the radial dynamics is given by a single harmonic for all the driving amplitudes. This is confirmed by Fig. 3(b), which shows that the amplitude of the modes n > 1 is significantly smaller than the mode n = 1 for all the p investigated. Nevertheless, the bubble dynamics are nonlinear, as can be seen in the change of mode amplitude as a function of the driving amplitude. In the linear regime, one expects the first mode to be linearly proportional to p. Figure 3(b) shows that increasing p by ten times, from 1 k Pa to 10 k Pa, results in a fivefold increase of |R 1 |, which is a signature of the nonlinear dynamics of the bubble, although this signature is weaker than in Carbopol. In other words, in kaolin, even when much higher pressure amplitudes are used than in Carbopol, nonlinearity cannot be detected either by the amplitude of the radial oscillations or by the (quite small) Fourier modes that are hardly present.

B. Conditions for oscillation-induced yielding
As the driving pressure is increased, the large-amplitude oscillations of the bubble generate considerable elastic stresses, potentially yielding the surrounding material. The periodic expansion and compression of the bubble generates extensional and compressive strains in the yield-stress fluid. If the amplitude of the acoustic pressure, p, is sufficiently large, the periodic elongational stresses due to the radial oscillations of the bubble can be larger than the yield stress. For a fixed set of constitutive parameters and a given bubble equilibrium size, there exists a frequency-dependent critical driving pressure, p crit , above which the material around the bubble yields during an oscillation cycle. The maximum normal stress difference, τ rr − τ θθ , occurs at the bubble surface and decays to zero at infinity. Thus, p crit is defined as the minimum pressure amplitude, p, for which the von Mises criterion, is satisfied at the bubble surface at least at one instant during a cycle. To find p crit , one has to solve the system of Eqs. (12)- (14) numerically for different acoustic pressure amplitudes and frequencies and find the minimum p for which Eq. (19) is satisfied. An estimate of p crit can be obtained by assuming that the dynamics of the bubble and that of the yield-stress fluid are linear until yielding occurs. The validity of this assumption is verified later through numerical simulations. In the linear regime, it is τ θθ = −τ rr /2 [38,48] and the elastic stress is linearly related to the strain, We evaluate Eq. (20) at r = R, with the assumption of small radial oscillations R Since in the linear regime x(t ) = R/R 0 sin(ωt + φ), the maximum amplitude during each cycle is given by R/R 0 . To evaluate R/R 0 , we use Eq. (18), resulting in the following maximum value of the radial stress at the bubble surface during each cycle: with β and ω 0 defined in Eq. (17). By inserting the maximum radial stresses given by Eq. (22) in the yielding criterion given by Eq. (19), and considering that τ θθ = −τ rr /2, we obtain an equation for the critical driving pressure for yielding: Equation (23)  is desired, e.g., to promote release of bubbles of a known size from a yield-stress fluid, Eq. (23) can be used as a rule of thumb to select the acoustic pressure and its frequency.
The dynamics of the bubble can be nonlinear even before yielding, due to the nonlinear inertial and elastic terms in Eqs. (13) and (14). It follows that the assumptions used to derive Eq. (23) might break down. To verify its relevance, we compare the linear yielding criterion p crit obtained from Eq. (23) with the numerical results obtained from the solution of the nonlinear Eqs. (12)- (14). We run the simulations for 2000 driving periods, always starting from the rest state, and discard the initial transient response. We consider the material yielded if at any instant during the last ten cycles Eq. (19) is satisfied. For a fixed frequency, we run simulations at increasing p until yielding is detected at a single instant during a cycle. This value of p is considered the p crit for that particular frequency. By repeating this process for different frequencies, we construct the curve p crit (ω).
In Fig. 4, we plot the critical pressure amplitude as a function of the angular frequency as predicted by Eq. (23) and by the numerical simulations. We considered the case of a bubble with an equilibrium radius R 0 = 100 μm suspended in the two yield-stress fluids considered, which have very different elastic moduli. In Fig. 4(a) we show the results obtained with the stiff kaolin suspension, and in Fig. 4(b) we consider the rheological parameters for the soft Carbopol gel. Figure 4 shows that p crit depends strongly on the driving frequency and shows a pronounced minimum at ω = ω 0 that is the resonance frequency of the bubble. At resonance, the driving pressure that yields the medium can be orders of magnitude smaller than that obtained off resonance. At frequencies that are much smaller than the resonance frequency of the bubble, p crit approaches a constant value given by p crit ≈ τ y (4G + 3p 0 + 4γ /R 0 )/(2 √ 3G). This expression and Eq. (23) indicate that the yielding pressure amplitude is inversely proportional to the elastic modulus. At frequencies larger than the resonance frequency, p crit is a linearly increasing function of ω. It is noteworthy that a stiffer EVP material generally requires smaller p to yield, i.e., it yields easier. In the case of a stiff yield-stress fluid, Fig. 4(a) shows that the linear theory gives a very good prediction of the p crit for all the angular frequencies explored and that very low pressure amplitudes are required to yield the material. The good agreement between the linear theory and the simulations is a consequence of the small strains, R/R 0 ≈ τ y /G, at which a very stiff material yields. Since for the stiff kaolin suspension it is τ y /G = 0.00045, prior to yielding the oscillations of the bubble are very small and the dynamics is in the linear regime.
In the case of the soft Carbopol gel, shown in Fig. 4(b), the critical pressure computed with the numerical simulations is in good agreement with the linear approximation for almost all the frequencies, except for ω ≈ ω 0 /2. The discrepancy observed at ω ≈ ω 0 /2 is a signature of the 073301-10 Amplitude of the Fourier modes |R n | of the radial dynamics. Linear theory predicts a single harmonics response, n = 1. Numerical simulations reveal that the dominant mode is |R 2 |, which corresponds to oscillations at twice the driving frequency. The total radial excursion is larger than that predicted by the linear theory.
weakly nonlinear response of the bubble that is neglected in the linear theory. The characteristic yield strain τ y /G = 0.1 of the Carbopol gel is much larger than that of the kaolin suspension and, at this strain, the dynamics of the bubble can be weakly nonlinear. In Fig. 5, we show that for ω = ω 0 /2, the radial oscillations depart from a single harmonic response predicted in the linear regime and the bubble experiences multiple harmonics, with the dominant one being twice the angular frequency of the acoustic driving ω. This is confirmed by investigating the amplitude of the Fourier modes reported in Fig. 5(b), which show that the largest mode |R n | is given by n = 2. The additional harmonics shown in Fig. 5(b) induce larger radial excursions compared to those predicted by the linear theory, hence resulting in larger strains and a smaller p crit . In summary, the results obtained with Carbopol and kaolin suggest that Eq. (23) is a very good estimate of the critical pressure for materials with characteristic yield strain, i.e., τ y /G smaller than one.

C. Dynamics of the yield surface
In the case of a driving pressure larger than p crit , part of the material surrounding the bubble is yielded and behaves as a liquid and the remaining part behaves as a solid. As a result, the bubble oscillates in a cavity with a time-dependent radius, filled by a viscoelastic liquid and surrounded by an elastic solid. This situation has been studied by Vincent et al. [68,69] in the context of cavitation in trees. The liquid and the solid regions are separated by the yield surface whose instantaneous position, r y (t ), is defined as the radial coordinate at which the von Mises criterion is satisfied: |τ rr (r y (t ), t ) − τ θθ (r y (t ), t )| = √ 3τ y . The periodic compression and expansion of the bubble generate cyclic elongational stresses that result in a time-dependent yield surface r y (t ). Since a bubble trapped in a yield-stress fluid can only rise when the surrounding material is yielded, it is interesting to investigate the evolution of the yield surface as its dynamics could have a strong impact on the rising velocity of the bubble.
In Fig. 6, we show snapshots of the bubble radius and of the position of the yield surface at three different instants during one cycle: minimum radius, equilibrium radius, and maximum expansion. We consider the case of a bubble of equilibrium radius R 0 = 100 μm trapped in a Carbopol gel and driven at ω = 0.9 ω 0 and p = 10 k Pa. To avoid transient effects, the period shown in Fig. 6 is chosen sufficiently far away from t = 0. In Fig. 6(a), the bubble is compressed, R < R 0 , and the Carbopol gel surrounding the bubble is yielded. As the bubble radius increases to reach its equilibrium value, the strain and the elastic stresses decrease; thus the yield surface moves towards the surface of the bubble. Eventually, the Carbopol unyields for R = R 0 ; see Fig. 6(b). During the subsequent expansion of the bubble, the strain and the elastic stresses increase again, the Carbopol gel yields, and the yielded region grows; see Fig. 6(c). To get a more detailed insight of the dynamics of the yield surface, we plot r y (t ) and the bubble radius R(t ) in Fig. 7 during one period. When the 073301-11 material is unyielded, the yield surface is not defined, which explains why r y (t ) in Fig. 7 is clipped for certain time intervals. The Carbopol gel unyields and then yields twice per cycle during the compression and the expansion phases as R goes through R 0 . This is a consequence of the change of sign of the normal stress difference between the compression and expansion phases, which implies that the normal stress difference must go through zero. It follows that for R ≈ R 0 , the deviatoric part of the stress tensor is smaller than the yield stress and the material unyields everywhere.
The dynamics of the yield surface and of the bubble radius in a kaolin suspension are shown in Fig. 8. We consider the case of a bubble of equilibrium radius R 0 = 100 μm driven at ω = 0.9 ω 0 and p = 0.16 k Pa, which corresponds to 5.5 times the critical pressure. Due to the large elastic modulus of the kaolin suspension, the elastic stresses are sufficiently large to yield part of the material, despite the very small oscillations of the bubble. In contrast to the case of a bubble oscillating in the Carbopol gel, the yielded region in the kaolin suspension behaves essentially as a viscous fluid because De is small. Despite this difference, Fig. 8 reveals that r y (t ) obtained in the case of a kaolin suspension is qualitatively similar to that obtained for a Carbopol gel. FIG. 7. Numerical simulations of the dynamics of the yield surface r y (t ) and of the bubble radius R(t ) during one cycle for the same parameters used in Fig. 4. The material unyields and yields twice per cycle when the bubble radius is close to its equilibrium radius and the elastic strains are small. We derive an approximate expression for the dynamics of the yield surface under the assumptions that the dynamics of the bubble is linear and that the material behaves as a Kelvin-Voigt solid everywhere. These assumptions are reasonable if the characteristic yield strain is small, τ y /G 1, if the material behaves mostly elastically in the yielded region, De 1, and if the driving pressure is close to the critical pressure, p ≈ p crit . The stresses in the material are then given by The position of the yield surface is given by the radial coordinate at which the von Mises yielding criterion is satisfied, giving the implicit equation subjected to the constraint that r y (t ) > R 0 . If at any instant it is r y (t ) R 0 , the material is unyielded. Substitution of the stresses given by Eq. (24) into Eq. (25) gives an equation for r y (t ): with R/R 0 and φ given by Eq. (18). Since Eq. (26) is derived under the assumption of linear bubble dynamics, we expect it to break down for driving pressures much larger than the critical pressure p crit . We explore the range of validity of Eq. (26) by finding the position of the yield surface through numerical simulations at different p/ p crit , with p crit computed from Eq. (23). In Fig. 9, we report the evolution of the yield surface predicted by the numerical simulations and by Eq. (26) for a bubble with equilibrium radius R 0 = 100 μm driven by an acoustic field at ω = 0.9ω 0 in the Carbopol gel. The Deborah number corresponding to this case is De = 7.36. In Fig. 9(a), the bubble is driven at a pressure close to the critical pressure, p = 1.62 p crit , and the dynamics of the yield surface given by the linear approximation given by Eq. (26) is very close to that obtained in the numerical simulation. As expected, Fig. 9(b) shows that by increasing the acoustic driving to p = 5.6 p crit , the dynamics becomes pronouncedly nonlinear and the linear theory fails to predict the evolution of the yield surface quantitatively. The linear theory systematically overpredicts r y (t ) in the first half of the cycle 073301-13 and underpredicts it in the second half. This is a consequence of the linearization, which neglects the advection of the yield surface due to the displacement of the bubble surface.

D. Impact of yielding on the radial dynamics
In this section, we explore the impact of yielding on the radial dynamics of a bubble. To highlight the effects of viscoplastic deformations, we compare the dynamics of a 100 μm bubble in an EVP fluid and in a neo-Hookean solid with the same elastic modulus. If any difference between the two behaviors is observed, it must be due to the yielding of the medium. In Fig. 10, we plot the dynamics of a bubble driven at p = 2 k Pa and at ω = ω 0 in the Carbopol gel and in the kaolin suspension, compared to its dynamics in a neo-Hookean solid. In the case of the Carbopol gel, Fig. 10(a) shows that the dynamics of the bubble is indistinguishable from that predicted in a neo-Hookean solid. Due to the small elastic modulus of the Carbopol gel compared to the kaolin suspension, the relaxation time of the liquid in the fluidized region is much larger than the driving frequency and De = 7.73. It follows that the yield-stress material behaves as an elastic solid both in the yielded and in the unyielded region, thus making the dynamics of the bubble identical to that predicted by the neo-Hookean model. Conversely, the dynamics of a bubble oscillating in the kaolin suspension is markedly different from that predicted by a neo-Hookean model, with the oscillations being significantly damped. In this case, it is De = 0.018 and the yielded region behaves as a viscous fluid. It follows that yielding of the material manifests itself as a larger damping compared to that expected for a neo-Hookean solid. These findings suggest that in the case De < 1, it should be possible to experimentally verify if the material has yielded by observing the dynamics of the bubble. Ideally, if, by increasing the acoustic pressure above p crit , a qualitative change in the dynamics of the bubble is observed due to yielding of the material, one might be able to measure the yield stress or at least identify yielding at high frequencies.
We test this hypothesis by computing the maximum radial excursion of a bubble in a kaolin suspension, for the same parameters used in Fig. 10 but varying the driving amplitude. The comparison with the neo-Hookean model reported in Fig. 11 shows that the amplitude of oscillations is significantly lower than that predicted in an elastic solid for driving pressures larger than p crit . This is a consequence of the additional dissipative processes taking place in the yielded region. The onset of the additional damping is not sharp at p = p crit because, for pressure slightly larger than the critical pressure, the size of the yielded region is small and the viscous stresses do not significantly impact the dynamics of the bubble. Finally, Fig. 11 shows that for pressures p > p crit , the oscillations of the bubble in the EVP fluid grow less than linearly with the driving amplitude. The sublinear increase of the oscillation amplitude is a consequence of the increase of viscous dissipation as the yielded region grows. The significant signature that yielding can have on the dynamics of the bubble suggests that a potential protocol for investigating yielding in experiments using acoustically driven microbubbles is to increase the driving power at a fixed frequency progressively. These findings have implications for bubble removal: as the yielded region grows, most of the power input by the pressure waves is lost to viscous dissipation. As a consequence, there might be an optimal choice for the power that maximizes the efficiency of the bubble release process.

V. CONCLUSIONS
We have investigated the dynamics of a bubble driven by an oscillating pressure field in an incompressible and elastic yield-stress fluid using numerical simulations and an approximate linear theory. We modeled the rheological behavior of the fluid using a recently developed constitutive model [55] that takes into account both elastic and viscoplastic deformations. By assuming that the bubble remains spherical during the pressure driving, we reduced the problem to a set of integrodifferential equations that we solve numerically using a Gauss-Laguerre method for the spatial integral and a fourth-order implicit Runge-Kutta time integration method. To explore the 073301-15 effects of different rheological parameters, we considered the case of a bubble driven by an acoustic field in a soft Carbopol gel and in a stiff kaolin suspension.
For a given bubble, there exists a frequency-dependent critical pressure at which the oscillations of the bubble yield the material. The critical pressure varies significantly with the frequency and it shows a pronounced minimum at the resonance frequency of the bubble. The critical pressure is very well approximated by an analytical formula derived under the assumption of linear bubble dynamics. In the case of an acoustic pressure larger than the critical pressure, a dynamic yield surface is developed inside the yield-stress fluid in the immediate environment of the bubble. We found that the position of the yield surface evolves significantly during one period, both in the Carbopol gel and kaolin suspension. The material unyields and then subsequently yields twice per period as the bubble goes through its equilibrium configuration. This is a consequence of the small elastic stress imparted to the yield-stress fluid by a bubble that is close to its equilibrium configuration. We developed an equation for the dynamics of the yield surface based on a linear approximation of the bubble oscillations. The linear theory is in good agreement with the fully nonlinear numerical simulations for p ≈ p crit , but deviates for larger driving amplitudes for which the assumption of linearity breaks down.
Finally, we explored the impact of yielding of the medium on the radial oscillations of the bubble. In the case of soft yield-stress fluids with elastic modulus of the order of G ≈ 100 Pa, yielding of the medium has negligible effects on the dynamics of the bubble. These materials have relaxation times that are much larger than typical inverse ultrasonic frequencies; thus, the yielded region behaves as an elastic solid. Conversely, we found that yielding has a significant impact on a bubble oscillating in the stiff kaolin suspension. In this case, the yielded region behaves as a viscous fluid, which is responsible for an extra oscillation damping. It sets in driving pressure larger than the critical pressure and it induces a sublinear dependence of the oscillation amplitude with the driving pressure.
Our results show that considering the elastic behavior of the yield-stress fluid is crucial to predict yielding of the material and bubble oscillations due to a finite pressure driving. Numerical simulations of bubble dynamics in stiff yield-stress fluids suggest that the onset of an additional damping at a critical pressure amplitude, p ≈ p crit , could be exploited to identify the signature of yielding in experiments, which would be cumbersome to assess otherwise. Finally, the theoretical and numerical framework presented in this paper can support experimental investigation of yielding under extensional deformation, which is relatively unexplored compared to yielding under shear deformation. Experiments are currently underway and will be the subject of a forthcoming paper.

APPENDIX: RHEOLOGICAL PREDICTIONS OF THE EVP MODEL
In this Appendix, we report the shear and extensional rheology of the Carbopol gel and of the kaolin suspension predicted by the EVP model. The constitutive parameters are given in Table I. We perform shear rheology simulations by fixing the shear rateγ and computing the shear stress τ xy and the first normal stress difference N 1 = 0.5(τ xx − τ yy ). The second normal stress difference is zero for the EVP constitutive model considered in the present work [56]. Figures 12 and 13 show the steady-state shear stress and first normal stress difference as a function of the shear rate for the case of the Carbopol gel and the case of the kaolin suspension, respectively. Both yield-stress fluids considered are shear thinning. The Carbopol gel shows much larger normal stresses than the kaolin suspension. This behavior is a consequence of its longer relaxation time, which results in larger elastic stresses. The difference in relaxation times between the two yield-stress fluids is better visualized in the transient shear stress response reported in Fig. 14. In contrast to the case of the Carbopol gel, the shear stress in the kaolin suspension reaches its steady-state value over ≈10 −3 s. The extensional rheology simulations are performed by applying a uniaxial extension rate˙ and computing the steady-state extensional viscosity η e = (τ xx − τ yy )/3˙ . Figure 15(a) shows that the Carbopol gel is extensional thinning, in agreement with the measurements performed by Louvet et al. [70]. Figure 15(b) shows that the EVP model predicts extensional thinning also for the case of the kaolin suspension. Measurements of extensional viscosity of kaolin suspensions show extensional thinning behavior at low extension rates and extensional thickening at large extension rates [71]. Since the EVP model predicts extensional thinning [see Fig. 15(b)], it is apparent that it is unable to correctly predict the rheological behavior of kaolin suspensions at large extension rates. Nevertheless, the choice of the EVP constitutive equation to model the kaolin suspension is justified because we explore extension rates for which experiments report extensional thinning.