4.1.1. Bubble response neglecting mass transfer

Since in this subsection, mass transfer is not accounted for, the last two dimensionless numbers in table 5 are not relevant. With no mass transfer between the bubble and the surrounding material, the bubble is a closed system in a thermodynamic sense. Hence, the energy content of the bubble remains proportional to the product of its initial pressure and its initial volume.

Figure 2. Time evolution of (a) the imposed ambient pressure, (b) the corresponding radius of the bubble and (c) the displacement of the centre-of-volume of the bubble without accounting for mass transfer.

Figure 2(a) presents the imposed ambient pressure as a function of time. Each pressure plateau is maintained for 1 h, consistent with the experimental protocol, except for the first one, which lasts only 2 min. This reduced duration is due to the system reaching a steady state rapidly, with no further observable changes. Given that the rate of mass transfer before the first pressure-decrease is very slow, even when mass transfer is accounted for, the duration of this step remains equally small. At each stepwise pressure reduction, the bubble responds immediately, resulting in stepwise increases in its radius and location of its centre of volume, as shown in figures 2(b) and 2(c). The sharp radius increases verifies that the reductions in the ambient pressure are immediately felt in the bubble, so that the product of bubble pressure and volume remains constant.

The bubble volume increase takes place anisotropically. The upper part (front) of the bubble moves upwards, because of the combination of bubble expansion and lower hydrostatic pressure compared with its value in the lower part. The lower part moves slightly downwards during the first two pressure reductions and then it translates upwards during the pressure change, although ever so slightly. This variable direction in smaller translation is caused by the fact that here, the same two effects act in opposite directions. Thus, the shape changes from spherical to prolate. Consequently, the centre-of-volume translates upward, as shown in figure 2(c). The horizontal plateaus observed in the centre-of-volume plot indicate that the bubble remains completely entrapped between pressure steps. This constitutes the first deviation from the reported experimental findings of Daneshi & Frigaard (Reference Daneshi and Frigaard2023) (their figure 3), who reported small displacement of the bubble centre-of-volume even while pressure is maintained constant during each step. We will address this discrepancy in our forthcoming analysis.

Figure 3. Snapshots of bubble shapes and locations at the end of each pressure decrease along with the yielded/unyielded regions to the left and normal axial stresses to the right around the bubble at (a) $\tilde{t}=2\,{\rm min}$ , (b) $\tilde{t}=62.5\,{\rm min}$ , (c) $\tilde{t}=123\,{\rm min}$ , (d) $\tilde{t}=183.5\,{\rm min}$ , (e) $\tilde{t}=244\,{\rm min}$ and (f) $\tilde{t}=304.5\,{\rm min}$ in the absence of mass transfer. The corrugated stress isolines result from the gradually coarser mesh used as we move away from the bubble, but they do not affect in any way the accuracy of our predictions.

Figure 4. Comparison of the experimental values and numerical results for (a) Y calculated with the bubble radius at that time, (b) radius of the bubble versus the ambient pressure with no mass transfer.

Figure 3 shows snapshots of the bubble shapes at the final instant of each pressure plateau. The complete transient evolution of the bubble without mass transfer is presented in supplementary movie 1 available at https://doi.org/10.1017/jfm.2026.11298. Initially in figure 3(a), only minimal yielding occurs very close to the two poles of the bubble, where compression (negative axial stress) dominates in the front and elongation (positive axial stress) in the rear of the bubble. After the first pressure decrease in figure 3(b), the bubble expands outwards maintaining an almost spherical shape. The uniform outward motion of the interface compresses the surrounding material everywhere around it, denoted by the negative axial stress both at the front and the rear of the bubble. The magnitude of the stresses is large enough to yield entirely the material in a spherical-like shell, which is shifted more towards the front pole. The bubble is located within this fully yielded envelope of stagnant fluid and remains static. We will elaborate further on why there is no motion albeit the fully yielded region at the end of § 4.1.2. In figure 3(c), the same pattern as previously is observed, with an even larger, fully yielded region surrounding the slightly prolate bubble. In figure 3(d), the pattern changes. Now, unyielded material appears at the south pole of the visibly prolate bubble, while the front yielded region has grown to a larger extent. The volume has increased enough, enhancing buoyancy and thus the tendency of the bubble to rise. The previously compressed rear region of the bubble is now decompressing due to this tendency. Hence, the stress drops below the yield limit giving rise to unyielded material. Further volume increase generates slight extension at the rear pole in figure 3(e) and more pronounced extension in figure 3(f). Still, the intensity of the extension is not adequate to yield the material there and it is not expected to yield thereafter without an additional pressure intervention. It is impossible for rising to occur while the far-field unyielded material is in contact with the surface of the bubble. This condition has been shown to be sufficient for achieving entrapment of a bubble in a viscoplastic fluid (Tsamopoulos et al. Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008). In the elasticity-including models of yield stress fluids, the unyielded state is expected to transiently deform and eventually yield after a sufficient deformation (Kordalis et al. Reference Kordalis, Pema, Androulakis, Dimakopoulos and Tsamopoulos2023).

The magnitude of the dimensionless stress in the vicinity of the bubble is close to unity; hence, the dimensional stress scales with the buoyancy term, i.e. $| \tilde{\boldsymbol{\tau }}| \sim O(\tilde{\rho }\tilde{g}\tilde{R}_{b})$ , as originally claimed. All other stress components also follow this scaling. Since buoyancy is the driving force for the upward motion, the material generates stresses that counteract it exactly, maintaining the immobility in agreement with Newton’s third law.

Not capturing bubble mobilisation is perfectly reasonable since the simulation does not achieve the necessary bubble volume, hence sufficient buoyancy, for the onset of motion. In contrast, the experiment shows bubble rising at the last step of pressure decrease. This constitutes the second deviation from the reported experimental findings of Daneshi & Frigaard (Reference Daneshi and Frigaard2023) and, as we will show later, it will be eliminated when mass transfer effects are included leading to increased bubble volume that initiates its rise.

In figure 4(a), we report a comparison between our numerical predictions and the experimental findings concerning the auxiliary transient yield number $ Y_{t}$ , i.e. the yield number calculated with the equivalent radius of the bubble at each time t. In figure 4(b), we compare our numerical predictions concerning the equivalent radius of the bubble with the experimental result at each ambient pressure level. The experimental work provides the values of $Y_{t}$ , which we convert to radii based on the definition of $Y$ . The computed $Y_{t}$ is systematically larger than the experimental value. This is translated into a systematic underestimation of the effective radius of the bubble. The small offset of experimental and numerical $Y_{t}$ at the first pressure level is attributed to the different yield stress values between our study and the experiments. This results from using the SHB to obtain $\tilde{\tau }_{y}=8.7\,{\rm Pa}$ , whereas the experimental study using the HB model results in $\tilde{\tau }_{y}=8.5\,{\rm Pa}$ . Although the initial difference is minor, the deviation increases progressively as the ambient pressure decreases. Figure 4(b) shows that the calculated and experimental radii are equal at the first pressure value but diverge thereafter. This constitutes the third deviation from the reported experimental findings of Daneshi & Frigaard (Reference Daneshi and Frigaard2023).

Examining the energy content of the bubble in the experimental work is essential. In figure 5, we show the $PV$ product of the bubble, calculated from the experimental data. The pressure of the bubble is calculated as the sum of the ambient pressure and the hydrostatic contribution. Although the bubble depth is not explicitly stated in the experimental work, a brief remark suggests that the hydrostatic pressure is approximately 2 kPa, corresponding to a depth of 2 $0\,{\rm cm}$ . The ambient pressure varies in the range $\tilde{P}_{a}\sim O(10^{4}{-}10^{5}\,{\rm Pa})$ and the hydrostatic part is $O(10^{3}\,{\rm Pa})$ . We neglect both capillarity and the fluid stress contribution because the former is $O(10{-}100\,{\rm Pa})$ and the latter is $O(\tilde{\rho }\tilde{g}\tilde{R}_{b})\sim O(10\,{\rm Pa})$ , as shown in figure 3. Finally, we neglect the minute upward displacement of the bubble. The experimental values of $PV$ progressively deviate from their initial value, i.e. the numerically imposed constant $PV$ in the simulation. The experimental values exceed their initial value by a factor larger than two. The $PV$ discrepancy is quite significant to be caused by experimental error during volume measurements. This suggests that the assumption of a constant $PV$ , as adopted by the authors, may not hold in their experiment. This constitutes the fourth deviation from the reported experimental findings of Daneshi & Frigaard (Reference Daneshi and Frigaard2023).

Figure 5. Calculated energy content of the bubble at each ambient pressure value in the experiment and the simulation with no mass transfer.

As a whole, the four deviations identified here indicate that the constant PV assumption is not valid under the experimental conditions. According to the ideal gas law, two primary explanations arise: either the temperature varies during the process or the bubble does not behave as a closed system and undergoes mass exchange with its environment. The former cause seems unlikely, because depressurisation occurs slowly and the heat capacity of the system is too large to allow any temperature change. At any rate, if we assume that the process is isentropic instead of isothermal, then (4.1) holds with a polytropic exponent equal to $1.4$ for air,

(4.1) \begin{equation}\tilde{P}_{b}\tilde{V}_{b}^{\gamma }=\tilde{P}_{b,o}\tilde{V}_{b,o}^{\gamma }.\end{equation}

Figure 6 presents the bubble volume data from the experiment, along with values derived under the assumptions of isothermal and isentropic variation, and the simulation values. The isentropic volumes are consistently smaller than the isothermal predictions, due to cooling of the gas during expansion. This discrepancy definitely rules out an isentropic process and, consequently, temperature variation as the primary cause of the observed $PV$ deviation. Thus, the only remaining explanation for the four deviations between the reported experimental findings and the simulation results is the mass exchange between the bubble and the surrounding material.

Figure 6. Bubble volumes from the experiment in comparison to those derived assuming isothermal or isentropic conditions in the absence of mass transfer.

4.1.2. Bubble response including mass transfer

The initial bubble radius is again $\tilde{R}_{b,o}=1.73\,{\rm mm}$ and the dimensionless numbers are listed in table 5. The value of ${\textit{Pe}}$ appears quite large, despite the slow motion of the bubble. Note that this nominal ${\textit{Pe}}$ is based on the characteristic velocity derived from the viscoplasticity-buoyancy scaling. While this scaling is the most appropriate choice, as it predicts the smallest possible velocity, it still overestimates its actual magnitude. Using the calculated velocity instead, which is of order $O(10^{-7} {\rm m}\,{\rm s}^{-1})$ yields the actual ${\textit{Pe}}$ of order $O(10^{-1})$ , which is physically reasonable. Despite the fact that convective terms could be neglected based on the low effective ${\textit{Pe}}$ value, we intentionally retain them for generality. Figure 7 presents simulations approximately up to $\tilde{t}=250\,{\rm min}$ . The onset of bubble motion occurs at the start of the sixth pressure level, marked by the sudden cessation of the blue dotted line in figure 7(a). The translation of the bubble is evident in figure 7(c), where the displacement exhibits a sharp increase. Note that the full displacement curve is not shown in figure 7(c), because the large values of $\Delta \tilde{z}_{b}$ after the bubble mobilisation distort the scale and obscure the earlier stasis phases. The omitted portion is provided separately in figure 12.

Figure 7. Time evolution of (a) the ambient imposed pressure, (b) equivalent radius of the bubble and (c) displacement of the centre-of-volume of the bubble, with comparison between the no mass transfer case and the mass transfer case. The solid red lines are identical to those in figure 2.

In figure 7(b), the bubble radius evolves differently from the case without mass transfer. Here, the radius gradually increases during each constant pressure interval, unlike the horizontal plateau observed when no mass transfer takes place. The growth due to mass transfer becomes evident by the end of the second pressure period at $\tilde{t}=62.5\,{\rm min}$ , where the bubble radius exceeds that of the no mass transfer case in the subsequent step, $\tilde{t}\gt 62.5\,{\rm min}$ . At each successive pressure level, the difference between the two curves steadily increases, exceeding $1\,{\rm mm}$ at $\tilde{t}\approx 244\,{\rm min}$ , while stasis is still maintained in both cases. This continuous, slow growth causes a gradual upward shift in the centre-of-volume during each constant pressure period, as observed in figure 7(c). The first deviation identified earlier has been removed.

Figure 8 shows the bubble growth, along with the yield surface evolution at each pressure level and the reduced concentration profile. The full transient response of the bubble with mass transfer is illustrated in supplementary movie 2. Initially in figure 8(a), the reduced concentration is positive around the bubble, indicating that the surrounding medium is slightly undersaturated. This is because the dimensional far-field concentration, $\tilde{k}_{H}\tilde{P}_{a,o}$ , is lower than the dimensional concentration around the bubble, $\tilde{k}_{H}\tilde{P}_{b}\cong \tilde{k}_{H}(\tilde{P}_{a,o}+\tilde{\rho }\tilde{g}\tilde{h}_{o})$ . As a result, gas is very slowly transferred from the bubble to Carbopol, which is consistent with the small positive value of ${\textit{Sa}}_{t}=3\times 10^{-4}$ . This trend is swiftly reversed after the first pressure reduction, shown in figure 8(b). The decreased gas concentration in the bubble makes the surrounding material oversaturated and the reduced concentration becomes negative. Although the dimensional concentration of gas in Carbopol far away from the bubble remains $\tilde{k}_{H}\tilde{P}_{a,o}$ , the new local dimensional concentration around the bubble drops to $\tilde{k}_{H}\tilde{P}_{b}\cong \tilde{k}_{H}(0.33\tilde{P}_{a,o}+\tilde{\rho }\tilde{g}\tilde{h}_{o})$ . Consequently, air is transferred from the medium to the bubble, with the current transient saturation number taking the value ${\textit{Sa}}_{t}=-0.036$ . As we explained in § 2.2, the negative sign denotes mass diffusing towards the bubble. The observed growth following this pressure drop in figure 8(b) suggests that the critical absolute value of the transient saturation number has been exceeded, i.e. $| Sa_{t,cr}| \lt 0.036$ . Further pressure reductions in figure 8(c–e) lower the reduced concentration even more. In figure 8(a–c), the negligible bubble movement creates a diffusion-dominated concentration profile with nearly concentric contours around the bubble. In figure 8(d), the motion is still small, but we observe the reduced concentration levels being slightly distorted, extended in the rear and squeezed in the front. This distortion is more pronounced in figure 8(e) since the bubble is very close to mobilising and exhibits a very slow upward motion.

Figure 8. Snapshots of bubble shapes and locations along with the yielded/unyielded regions (left half of each panel) and the reduced concentration profile (right half of each panel) around the bubble at (a) $\tilde{t}=2\,{\rm min}$ , (b) $\tilde{t}=62.5\,{\rm min}$ , (c) $\tilde{t}=123\,{\rm min}$ , (d) $\tilde{t}=183.5\,{\rm min}$ , and (e) $\tilde{t}=244\,{\rm min}$ , for when mass transfer is accounted.

Figure 9 presents a direct comparison between bubble shapes and the yielded areas when mass transfer is not accounted for (left half) and when it is (right half). Initially, in figure 9(a), both bubbles are spherical and exhibit similar yielded regions confined around their two poles. In figure 9(b), the difference in size becomes evident. The mass-transfer bubble expands more than the no-mass transfer bubble, leading to visibly larger fluidised areas. In both cases, however, the surrounding material is fully fluidised while stasis prevails. The additional mass elongates the bubble with mass transfer further, but generates unyielded material at its lower pole in figure 9(c), whereas the lack of mass transfer allows for full fluidisation of the material around the smaller bubble on the left. In figure 9(d), the mass transfer bubble has created significant extension at its south pole, sufficient to yield the material locally. Nevertheless, it is still in contact with unyielded material and static. Finally, in figure 9(e), the unyielded material has barely lost contact with the surface of the mass transfer bubble. Also, its shape is at the most elongated state during this expansion. Prolate entrapped bubbles are more agile than the spherical ones and require a greater plastic resistance to remain stationary (Pourzahedi et al. Reference Pourzahedi, Chaparian, Roustaei and Frigaard2022). According to these remarks, the shape of the bubble with mass transfer and its yielded areas in figure 9(e) suggest that it is ready to start rising, which is initiated at the next reduced level of pressure. The second deviation identified earlier has been removed.

Figure 9. Comparison of the bubble shapes and locations along with the yielded areas between the no-mass-transfer bubble (left half of panels) and the mass-transfer bubble (right half of panels), while the bubble remains static at (a) $\tilde{t}=2\,{\rm min}$ , (b) $\tilde{t}=62.5\,{\rm min}$ , (c) $\tilde{t}=123\,{\rm min}$ , (d) $\tilde{t}=183.5\,{\rm min}$ and (e) $\tilde{t}=244\,{\rm min}$ .

Figure 10 presents the evolution of $Y_{t}$ and of the bubble radius at each pressure level, now that mass transfer effects are incorporated. In both figures 10(a) and 10(b), the mass transfer predictions align closely with the experimental results. In figure 10(a), the blue curve consists of segments with positive slopes followed by vertical ones. The former segments correspond to periods of pressure decrease, while the latter to periods of constant pressure. As shown in figure 10(b), the vertical segments increase the radius due to mass transfer while the pressure remains constant. These vertical parts, which last 60 minutes, adjust the trajectory of the curve towards the experimental values. During the remaining segments, the bubble radius satisfies the constant $PV$ equation. Physically, this occurs because the pressure reduction occurs much faster and abruptly, not allowing enough time for mass transfer to become noticeable. We will provide a detailed investigation and supporting evidence for this mechanism in § 4.4.2. Furthermore, the simulation manages to capture the onset of the bubble mobilisation. We find numerically that the critical yield number for mobilisation is $Y_{t,cr}=0.2$ , which is close to the experimental value, $Y_{cr,exp}=0.18$ . Overall, mass transfer is identified as the key missing mechanism to capture the faster increase of the bubble volume. The third deviation identified earlier has been removed.

Figure 10. Comparison of the experimental values and the numerical predictions with and without mass transfer for (a) $Y$ calculated via the instantaneous bubble radius, (b) the bubble radius versus the ambient pressure for protocol I.

Figure 11 presents the mass-transfer adjusted energy content ( $PV$ ). The predicted results closely align with the experimental ones. The final energy content exceeds its initial value by more than a factor of two. Although Daneshi & Frigaard (Reference Daneshi and Frigaard2023) assumed no mass transfer effects in their analysis, they observed indirectly their presence. In the validation of their experimental method, they noted that smaller $\Delta \tilde{t}_{\textit{step}}$ cases require an additional pressure reduction step to induce bubble rise. This is readily explained by noting that a smaller $\Delta \tilde{t}_{\textit{step}}$ allows less time for inward mass diffusion, resulting in a smaller bubble volume at a given pressure compared with a larger $\Delta \tilde{t}_{\textit{step}}$ . Therefore, the bubble cannot rise without a further decrease in the ambient pressure. Hence, the fourth deviation identified earlier has been removed.

Figure 11. Estimated energy content of the bubble at each ambient pressure value.

Figure 12. Time evolution of (a) the bubble velocity and (b) the bubble displacement during the rising phase of the bubble for when mass transfer is accounted. The experimental results correspond to figure 6 of Daneshi & Frigaard (Reference Daneshi and Frigaard2023). The dimensionless displacement from the experimental work is multiplied by the plateau value (5.9 mm) of the simulation in panel (b), while the time instant $\tilde{t}=0$ of the experimental work is set approximately at 237 min of the simulation.

Finally, we present the eventual rising of the previously immobile bubble when the sixth step of pressure is applied at $\tilde{t}\sim 244\,{\rm min}$ . This induces the last velocity peak shown in figure 12(a) followed by a mild velocity decrease to non-zero values between $0.05\,{\rm mm}\,{\rm s}^{-1}$ and $0.1\,{\rm mm}\,{\rm s}^{-1}$ . The bubble clearly translates upwards. In figure 12(b), we report the displacement of the bubble during the same period compared with the respective experimental data. The comparison is remarkable with the two data sets practically superimposing, an observation that reassures us on the validity of the current implementation. At the end of the simulation, the bubble has travelled more than $25\,{\rm mm}$ , a distance which exceeds its initial radius by 15 times. The transition to the rising phase is presented in figure 13.

Figure 13. Snapshots during the initiation of motion of the bubble presenting bubble shapes and locations along with the yielded/unyielded regions to the left and normal axial stress to the right around the bubble at (a) $\tilde{t}=244\,{\rm min}$ and (b) $\tilde{t}=246\,{\rm min}$ . In panel (b), the path lines are also reported.

In figure 13(a), we show the last instant during the penultimate pressure value. The far-field unyielded material has barely lost contact with the bubble surface and no motion has been detected yet. Then, the final pressure value is imposed and, shortly after the resulting bubble enlargement, the bubble is mobilised in figure 13(b). The material’s resistance to motion has been overcome, with the far-field unyielded material having a visible distance from the bubble surface. The entrapment criterion of Tsamopoulos et al. (Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008) is no longer satisfied, leading to rising.

The most important observation in figure 13(b) is related to the formation of the streamlines. It is a well-known fact that the buoyancy-driven motion of a body, either rigid or deformable, in any incompressible fluid, Newtonian or complex, generates a recirculation at the side of the object in an inertial reference frame, similar to the one of figure 13(b) (Arigo & McKinley Reference Arigo and McKinley1998; Brücker Reference Brücker1999; Leal, Skoog & Acrivos Reference Leal, Skoog and Acrivos1971; ten Cate et al. Reference ten Cate, Nieuwstad, Derksen and Van den Akker2002; Ninomiya & Yasuda Reference Ninomiya and Yasuda2006; Ortiz et al. Reference Ortiz, Lee, Figueroa-Espinoza and Mena2016). As the body moves due to buoyancy, it displaces the fluid at its front to occupy its new position, leaving a void at its prior position. Then, the fluid recirculates to fill this void and maintain continuity. The size of the recirculation depends on the viscosity of the fluid, among other properties. This fluid rotation generates velocity gradients, which in turn affect the stress field. In the specific case of yield stress materials, the recirculation is present (Mougin et al. Reference Mougin, Magnin and Piau2012; Holenberg et al. Reference Holenberg, Lavrenteva, Liberzon, Shavit and Nir2013; Putz et al. Reference Putz A.M., Burghelea, Frigaard and Martinez2008) and its size depends also on plasticity. We attribute the characteristic shape of an ‘umbrella’ in the front yielded area (Moschopoulos et al. Reference Moschopoulos, Spyridakis, Varchanis, Dimakopoulos and Tsamopoulos2021; Esposito, Dimakopoulos & Tsamopoulos Reference Esposito, Dimakopoulos and Tsamopoulos2024; Yazdi and D’Avino Reference Yazdi and D’Avino2023) to the presence of the recirculation, the size of which is such as to allow the turning in the trajectory of the fluid parcels. However, if plasticity dominates and the yielded space is not adequate to encompass the mandatory-for-rising recirculation, then the object will remain trapped in place. Since the surrounding material cannot move to fill the void from the movement of the body, then the only way for continuity to be sustained is the fluid to remain static. This is the reason why a bubble inside a spherical-like shell of fully yielded material, as shown in figure 9, can remain static, irrespective of mass transfer being accounted for or not. Lastly, we notice the reversal of the flow below the bubble, the elastic phenomenon named ‘negative wake’. In Appendix C, we discuss the effect of the initial concentration of dissolved air on the bubble dynamics regarding the present problem.

4.3.1. Mass transfer properties

We simulate the material used under protocol II with three values of the diffusion coefficient and three values of Henry’s constant, each one including the respective base value. As we concluded in the previous section, the absence of mass transfer eradicates hysteresis.

Figure 18(a) shows the effect of the diffusion coefficient on the variation of the bubble radius with the imposed ambient pressure. The values we use match the range of existing values in real gas–water systems. The intermediate value is the one we have already examined thoroughly. There are two key points to consider during the evaluation of protocol II, which are the maximum value of the radius at the lowest pressure during the ramp down and the last radius at the final pressure during the ramp up, quantifying the magnitude of the hysteresis. We observe that the whole ramp-down phase, and most importantly the maximum radius value, is better approximated by the largest diffusion coefficient. However, this diffusion coefficient results in the largest deviation from the experimental values during the ramp up phase. Conversely, the lowest diffusion coefficient accurately predicts the final value of the radius in the ramp-up phase; it almost superimposes to the experimental one. Nevertheless, both the ramp-down and the ramp-up phase significantly differ from the experiment. Therefore, the base value strikes a balance between the ramp-down and ramp-up phases, providing a relatively accurate approximation of both the peak value and the final bubble radius. Overall, increasing the diffusion coefficient increases the difference between the two branches of the hysteresis, primarily affecting the ramp-up phase, although these increases are rather mild.

Figure 18. Radius versus imposed ambient pressure for the parametric analysis while varying (a) the diffusion coefficient and (b) the Henry constant. The orange curves correspond to the base values. In each case, the rest of the properties are those of the base case.

Figure 18(b) shows the effect of increasing the Henry constant. The values we used are at the lower end of gas-fluid pair solubility. Starting with the base value, which is the lowest one for the air–water pair and doubling it, we reach a medium–low value. The Henry constant at the upper end of the range can reach one to two orders of magnitude greater than our largest value. Yet, the trend from increasing this property is satisfactorily captured. The only value that can capture the experimental findings is the base one. Already using the second value of $\tilde{k}_{H}$ , leads to predictions that far exceed the measurements of the ramp-up phase, resulting in a much larger hysteresis at the final pressure. Naturally, increasing $\tilde{k}_{H}$ further leads to unacceptable predictions. Clearly, $\tilde{k}_{H}$ has a more pronounced effect on the predictions than $\tilde{D}$ .

4.3.2. Rheological properties

We examine the effect on the predictions of protocol II keeping all properties the same as the base case, except for the shear modulus and the yield stress, which we change separately. In figure 19, we present these results. First, we double the shear modulus and then we increase the yield stress to the largest value we could simulate. Most interestingly, all curves superimpose with the base case and no observable deviation can be seen.

Figure 19. Radius versus imposed ambient pressure for the material properties shown in the figure inset.

The material properties directly affect the stresses generated around the bubble. To understand why this superposition occurs, we need to trace the connection between stress and volume, which regulates the radius of the bubble,

(4.2) \begin{align} \tilde{P}_{b}&=\tilde{P}_{a}+\tilde{\rho }\tilde{g}\tilde{h}+\frac{2}{\tilde{R}_{b}}\tilde{\sigma }-\tilde{\tau }_{\textit{nn}}, \end{align}

(4.3) \begin{align} \tilde{\tau }_{\textit{nn}}&=\frac{\int _{{\tilde{S}_{b}}}\tilde{\tau }_{\textit{nn}}\,\mathrm{d}\tilde{S}}{\int _{{\tilde{S}_{b}}}\,\mathrm{d}\tilde{S}}. \end{align}

The average stress on the surface of the bubble, defined by (4.3), contributes to the bubble pressure as shown in (4.2). So, the bubble volume multiplied by this stress-including pressure equals the right-hand side in the ideal gas law. It is therefore necessary to assess the magnitude of this additional stress. Figure 20 shows the normal axial stress for the cases with $\tilde{\tau }_{y}=12 \,{\rm Pa}$ and $\tilde{G}=200 \,{\rm Pa}$ . The dimensionless stress magnitude is $O(1)$ , so the dimensional stress is $O(\tilde{\rho }\tilde{g}\tilde{R}_{b})$ . This agrees with our earlier observation in § 4.1.1. As a result, the stress contribution is significantly smaller than the dominant ambient and hydrostatic pressure, and has a negligible effect on the bubble volume. This explains why the curves for different material properties in figure 19 coincide. In essence, the rheological properties of the material do not influence the magnitude of the stress field around the bubble when the bubble is immobilised, but rather the stresses arise primarily to counteract the bubble’s buoyancy to maintain it nearly stationary. The material’s rheological properties only contribute to suspending the bubble.

Figure 20. Normal axial stress field around the bubble for $\tilde{\tau }_{y}=12\,{\rm Pa}$ to the left half of the panel and $\tilde{G}=200\,{\rm Pa}$ to the right half of the panel at $\tilde{t}=40.5\,{\rm min}$ at the minimum pressure. The solid black line is the yield surface. In each case, the rest of the properties are those of the base case.

4.4.1. Diffusion around the bubble under constant pressure

Our previous detailed simulations indicate that convective effects may be neglected, especially when the bubble is entrapped. We take advantage of this observation to develop a semi-analytical model to approximate the mass transfer induced variation of the bubble radius during a constant pressure level following an abrupt pressure change. First, we assume that the bubble retains its initial spherical shape, allowing adoption of a spherical coordinate system and reducing the problem to a one-dimensional one, and then, we account for the deviation from the spherical shape. The starting point is the species interfacial balance, (2.8) in dimensional form. Upon neglecting the two convective terms and solving for the rate of change of the radius, we obtain

(4.4) \begin{equation}\frac{\mathrm{d}\tilde{R}^{*}}{\mathrm{d}\tilde{t}}=\frac{\tilde{D}}{\left(\tilde{c}_{b}-\left.\tilde{c}\right| _{{\tilde{R}^{*}}}\right)}\left.\frac{\partial \tilde{c}}{\partial \tilde{r}}\right| _{{\tilde{R}^{*}}}.\end{equation}

The assumed spherical symmetry dictates that the bubble radius is a function of time only $\tilde{R}^{*}(\tilde{t})$ , while species concentration depends on time and radial coordinate. This simpler dependence of the bubble radius is indicated with an asterisk.

The solution to this diffusion problem in an unbounded domain around a bubble of instantaneous radius $\tilde{R}^{*}$ is given by (Carslaw & Jaeger Reference Carslaw and Jaeger1959)

(4.5) \begin{align}\tilde{c}&=\tilde{c}_{\infty }+\left(\left.\tilde{c}\right| _{{\tilde{R}^{*}}}-\tilde{c}_{\infty }\right)\frac{\tilde{R}^{*}}{\tilde{r}}erfc\left(\frac{\tilde{r}-\tilde{R}^{*}}{2\sqrt{\tilde{D}\tilde{t}} }\right)\!, \end{align}

(4.6) \begin{align}\left.\frac{\partial \tilde{c}}{\partial \tilde{r}}\right| _{{\tilde{R}^{*}}}&=-\left(\left.\tilde{c}\right| _{{\tilde{R}^{*}}}-\tilde{c}_{\infty }\right)\left(\frac{1}{\tilde{R}^{*}}+\frac{1}{\sqrt{\pi \tilde{D}\tilde{t}}}\right)\!, \end{align}

where $\tilde{c}_{\infty }$ is the species concentration in the far-field. Introducing this expression for the concentration in (4.6) yields

(4.7) \begin{equation}\frac{\mathrm{d}\tilde{R}^{*}}{\mathrm{d}\tilde{t}}=-\tilde{D}\frac{\left.\tilde{c}\right| _{{\tilde{R}^{*}}}-\tilde{c}_{\infty }}{\tilde{c}_{b}-\left.\tilde{c}\right| _{{\tilde{R}^{*}}}}\left(\frac{1}{\tilde{R}^{*}}+\frac{1}{\sqrt{\pi \tilde{D}\tilde{t}}}\right)\!.\end{equation}

Equation (4.7) is a generalisation of the model reported by Epstein & Plesset (Reference Epstein and Plesset1950) (henceforth, indicated as EP), because of the appearance of the term $\left.\tilde{c}\right| _{{\tilde{R}^{*}}}$ in the denominator. Krieger, Mulholland & Dickey (Reference Krieger, Mulholland and Dickey1967) used the EP model for an immobilised spherical bubble to determine the diffusion coefficient of common gases in Newtonian fluids. Others have used the EP equation with and without the square-root term in the parenthesis and observed significant deviations (Cable Reference Cable1967) or solved it numerically to determine the evolution of the bubble radius (Cable & Evans Reference Cable and Evans1967).

We use Henry’s law and the ideal gas law to eliminate concentrations in favour of pressures in (4.7) and perform a variable change, $\tilde{t}-\tilde{t}_{i}$ , to allow the starting point for mass transfer to be any time $\tilde{t}_{i}$ :

(4.8) \begin{equation}\frac{\mathrm{d}\tilde{R}^{*}}{\mathrm{d}\tilde{t}}=-\tilde{D}\frac{\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}}{1-\left(\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}\right)}\left(1-\frac{\tilde{P}_{a,o}}{\tilde{P}_{b,i}}\right)\left[\frac{1}{\tilde{R}^{*}}+\frac{1}{\sqrt{\pi \tilde{D}\left(\tilde{t}-\tilde{t}_{i}\right)}}\right]\!.\end{equation}

The subscript $i$ indicates that the constant pressure $\tilde{P}_{b,i}$ was first set at time $\tilde{t}_{i}$ . It is noteworthy that the transient saturation number, ${\textit{Sa}}_{t,i}$ defined as $\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}(1-({\tilde{P}_{a,o}}/{\tilde{P}_{b,i}}))$ , naturally arises in this expression. Its sign signals either growth or shrinkage of the bubble depending on the size of $({\tilde{P}_{a,o}}/{\tilde{P}_{b,i}})$ . Another interesting observation is the dimensionless Henry constant appearing in the denominator. The low solubility of air in water leads to $\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}=0.02$ , which is small compared with unity. In contrast, our tests have shown that its presence becomes progressively more important when its value exceeds 0.05, resulting from a Henry constant $\tilde{k}_{H}\gt 2\times 10^{-5} \mathrm{mol}\,\mathrm{m}^{-3}\,\mathrm{Pa}^{-1}$ at the given temperature. A fairly large value arises with $\mathrm{CO}_{2}$ in an aqueous solution, where $\tilde{k}_{H}=3.3\times 10^{-4} \mathrm{mol}\,\mathrm{m}^{-3}\,\mathrm{Pa}^{-1}$ resulting in $\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}\approx 0.82$ . This makes plain the importance of the generalisation of (4.8) over the EP model. In all cases, we keep this term in our analysis.

Equation (4.8) may be modified to apply even for ellipsoidal bubbles. To this end, the spherical radius $\tilde{R}^{*}$ is substituted by the equivalent spherical radius $\tilde{R}_{b}$ divided by a correction factor $A$ , which accounts for the increased surface of an ellipsoidal bubble over a spherical one with the same volume. More details are given in Appendix E. For an ellipsoidal bubble, the rate of change of the bubble radius is given by

(4.9) \begin{equation}\frac{\mathrm{d}\tilde{R}_{b}}{\mathrm{d}\tilde{t}}=-A\tilde{D}\frac{Sa_{t,i}}{1-\left(\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}\right)}\left[A\frac{1}{\tilde{R}_{b}}+\frac{1}{\sqrt{\pi \tilde{D}\left(\tilde{t}-\tilde{t}_{i}\right)} }\right]\!.\end{equation}

Application of (4.9) requires the working conditions, i.e. the pressure used to saturate the medium before the experiment, $\tilde{P}_{a,o}$ , the induced pressure, $\tilde{P}_{b,i}$ , at time $\tilde{t}_{i}$ , and the temperature. It provides a correlation between the experimentally measured evolution of the radius and the properties of mass transfer, $\tilde{D}$ and $\tilde{k}_{H}$ .

In the absence of independent data for yield-stress fluids, we cannot validate it. Instead, we will use it to extract these properties by fitting the radius versus time results acquired from our earlier detailed simulations, which accurately replicated the experimental values. After reducing the ambient pressure to bring the bubble pressure at $\tilde{P}_{b,i}$ , the corresponding initial radius $\tilde{R}_{i}=\tilde{R}(\tilde{t}=\tilde{t}_{i})$ is used as an initial condition. We perform this procedure using the second constant pressure value of protocol I to obtain the fitted values of the properties and compare them with the actual values we imposed in the detailed model.

In figure 21, we observe that the simulation curve and the one resulting from (4.9) are in very good agreement. The latter returns values of the predicted properties with less than 5 % relative error, as shown in table 7. This finding seems very promising for capturing both the Henry constant and the diffusion coefficient with a single experiment. Fitting the radius in $\mathrm{mm}$ returns a diffusion coefficient in $\mathrm{mm}^{2}\,\mathrm{s}^{-1}$ , as shown in table 7.

Table 7. Mass transfer properties either used in the detailed simulation or resulting from fitting (4.9).

Figure 21. Time evolution of the radius at the second constant pressure level of protocol I either from the detailed simulation or from solving (4.9) with the two fitted properties.

4.4.2. Simulation of the protocols

As we established in § 4.3.2, the material’s rheology has a negligible effect on the bubble volume. The volume evolution is determined by two effects: (a) the abrupt pressure change during the transition to the new pressure level and (b) the mass transfer during the constant pressure period. These two effects correspond exactly to the two terms in (2.13) and become alternately dominant.

Figure 22 supports this claim. It presents the two contributions along with their sum, which is the time derivative of the volume during protocol II. Figure 22(a) shows the three quantities during the whole protocol. The variation of the ambient pressure causes the spikes in the time-derivative of the volume. We mainly focus on the period between the last pressure decrease and the first pressure increase, as indicated by the grey dotted perimeter and magnified in figure 22(b). During the constant pressure period, within $3300\leq t\leq 3700$ in figure 22(b), only the mass variation modifies the bubble volume. The orange short-dashed perimeter denotes the last pressure decrease and is magnified in figure 22(c), while the brown long-dashed perimeter denotes the first pressure increase and is magnified in figure 22(d). During these two spikes, the volume variation follows the abrupt change in the pressure contribution, which exceeds by two to three orders of magnitude the mass contribution.

Figure 22. Time evolution of the bubble volume along with its two generating factors, the pressure contribution and the mass contribution during (a) the entire duration of protocol II, (b) the magnified period between the last pressure decrease and the first pressure increase, (c) the magnified last pressure decrease spike, and (d) the magnified first pressure increase spike.

Therefore, the radius evolution of a spherical bubble during both the constant and varying pressure segments may be given by a two-branch function. In the first branch, the radius evolution is caused by mass transfer and it is determined by (4.9) with $A=1$ , while in the second branch, it is caused by the pressure jump, which is related to the bubble volume via the ideal gas law for fixed mass in the bubble,

(4.10) \begin{align}\frac{\mathrm{d}\tilde{R}_{b}}{\mathrm{d}\tilde{t}}=\left\{\begin{array}{rl} -\tilde{D}\dfrac{Sa_{t,i}^{\textit{eff}}}{1-\left(\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}\right)}\left[\dfrac{1}{\tilde{R}_{b}}+\dfrac{1}{\sqrt{\pi \tilde{D}\left(\tilde{t}-\tilde{t}_{i}\right)}}\right], & \tilde{P}_{b}=\tilde{P}_{b,i}=\textit{const}.,\\[12pt] -\dfrac{\tilde{R}_{b}}{3\tilde{P}_{b}}\dfrac{d\tilde{P}_{b}}{d\tilde{t}}, & \tilde{P}_{b}=f\left(\tilde{t}\right)\!. \end{array}\right.\end{align}

We observe that the first branch of (4.10) contains ${\textit{Sa}}_{t,i}^{\textit{eff}}$ . The solution of the transient mass transfer around a spherical object (see (4.5)), used to derive the first branch, by definition assumes that the surrounding medium has constant concentration $\tilde{c}_{\infty }$ (no spatial variation) when $\tilde{t}=\tilde{t}_{i}$ . So, each time the first branch is reactivated, a new concentration field is used with $\tilde{c}|_{\tilde{t}={\tilde{t}_{i}}}=\tilde{c}_{\infty }$ everywhere in the fluid. Each cycle of this simplified model is connected to its previous cycle only through the value of the radius. Therefore, the concentration gradient used is sharper than its actual value, overestimating mass transfer.

(4.11) \begin{align}&\qquad\qquad\qquad\qquad\qquad {\textit{Sa}}_{t,i}^{\textit{eff}}=\tilde{k}_{H}\tilde{R}_{\textit{gas}}\tilde{T}_{o}\left(1-\frac{\tilde{P}_{a,o}^{\textit{eff}}}{\tilde{P}_{b,i}}\right)\!, \end{align}

(4.12) \begin{align}& \tilde{P}_{a,o}^{\textit{eff}}=\frac{\tilde{c}\big(\tilde{r}=\beta \tilde{R}_{b}\big)}{\tilde{k}_{H}}=\tilde{P}_{a,o}+\beta \left(\tilde{P}_{b,i}-\tilde{P}_{a,o}\right)\textit{erfc}\left(\frac{\left(\beta -1\right)\tilde{R}_{b}}{2\sqrt{\tilde{D}\tilde{t}} }\right)\!,\quad \beta \geq 1.4. \end{align}

The definition of ${\textit{Sa}}_{t,i}^{\textit{eff}}$ in (4.11) contains the newly defined effective saturation pressure in (4.12). The latter is based on the concentration, (4.5) at a distance $\beta \tilde{R}_{b}$ from the bubble. Note that the argument of the error function does not contain the quantity $\tilde{t}-\tilde{t}_{i}$ , but $\tilde{t}$ . As time passes and the bubble pressure is reduced, the effective saturation pressure decreases. This macroscopic reduction to the effective driving force imitates the diffusion around the bubble that dampens the gradient on the bubble interface during each constant pressure value. Hence, using $\tilde{t}$ in this expression provides the evolution of the concentration profile during the process. Through trial and error, we determined that the effective distance $\beta$ should exceed 1.4. Giving it a much larger value nullifies the complementary error function and the actual saturation pressure is retrieved; hence, ${\textit{Sa}}_{t,i}^{\textit{eff}}=Sa_{t,i}$ . However, setting $\beta$ to unity leads to ${\textit{Sa}}_{t,i}^{\textit{eff}}=0$ , ‘switching off’ the mass transfer. We concluded that the effective distance coefficient $\beta$ lies in the range of 1.4–1.7. We give more details of the model in Appendix F.

A value of $\beta$ close to unity produces a small driving force, practically switching off the mass transfer. In the case of protocol II, it tends to collapse the ramp-down and ramp-up branches, hindering the hysteresis. In contrast, a value above 2.5 produces a driving force equivalent to $\tilde{P}_{a,o}$ . Figure 23 shows that the predictions of our new model agree quite nicely with the experimental results by fine tuning $\beta$ . In figure 23(c), the simplified model approaches more accurately the experimental data than our two-dimensional (2-D) simulations. Since the β parameter regulates the apparent concentration gradient and, in this case, a rather small value of $\beta$ has been used, we tend to believe that the initial solute concentration used in the 2-D simulation has probably been overestimated, as stated earlier. It is noteworthy that our simplified model could be used to simulate the work of Miele et al. (Reference Miele, Abate, Taki and Di Maio2024) during the so-called pressure quenching protocol for foam formation. This protocol consists of a bubble-generating phase through an initial pressure decrease in a polymeric melt, followed by a pressure increase to regulate the rate of growth of the bubble. There are strong similarities between this experimental work and the examined protocol II, especially during the inversion from the ramp-down to the ramp-up phase.

Figure 23. Comparison of the predictions using the simplified model with those using our detailed simulations and the experimental data for (a) protocol I with $\beta =1.7$ , (b) protocol II for $\Delta \tilde{t}_{\textit{step}}=5\,{\rm min}$ with $\beta =1.4$ and (c) protocol II for $\Delta \tilde{t}_{\textit{step}}=1\,\textrm{h}$ with $\beta =1.5$ .