4.1. Rheology of the material

As a base material, we chose the 1.5 % polyox solution used in LLJ, which was characterized sufficiently in the work of Joseph et al. (Reference Joseph, Liu, Poletto and Feng1994) . The zero-shear rate viscosity, ${\tilde{\eta }_0}$, is equal to the viscosity of the material in the plateau shown for small shear rates. We determine the relaxation time, $\tilde{\lambda }$, from the intersection of the storage modulus, $G^{\prime}$, and loss modulus, $G^{\prime\prime}$, in the frequency sweep experiments. Finally, we find the ${\epsilon _{PTT}}$ value following a nonlinear regression on the flow curve.

We cannot estimate simultaneously the values of $\tilde{\lambda }$ and ${\varepsilon _{PTT}}$ from the flow curve alone because it is their product that appears in the constitutive equation, prohibiting their independent calculation only from one type of experiment. Data from extensional experiments are not reported, which limits determining fully the extensional behaviour of the rheological models. Also, we note that to achieve a perfect fit of the model, we should use a ${\varepsilon _{PTT}}$ value larger than 1, which is very uncommon. Hence, we select eventually a smaller one equal to 0.4. We will ascertain the effect of the value of ${\epsilon _{PTT}}$ in our problem in the following sections. In figure 4(a), we plot the predictions of the model superimposed on the experimental values and we summarize the parameter values in table 1. The value of the interfacial tension, $\tilde{\sigma }$, of LLJ and in our study is equal to $0.063\ \textrm{N}\ {\textrm{m}^{ - 1}}$.

Figure 4. (a) Steady-state flow curve for the 1.5 % polyox solution of Liu et al. (Reference Liu, Liao and Joseph1995). The symbols represent the experimental data and the continuous line represents the e-PTT model predictions. (b) Prediction of the steady shear response for various viscoelastic models. The black, continuous line corresponds to the e-PTT model; the blue, dashed-dotted line corresponds to the l-PTT model; and the red, dashed line corresponds to the FEG model.

Table 1. Rheological parameters of the e-PTT model found via a nonlinear regression for 1.5 % polyox solution.

It is known that the extensional rheology of the material plays a crucial role in bubble dynamics, especially downstream. Even though we choose the e-PTT model based on its previously successful numerical implementations, we consider it appropriate to examine its rheological behaviour under uniaxial extension and compare it with two other constitutive models that present a different behaviour. We select the linear PTT model, l-PTT, because it is another popular viscoelastic model, and the finite-extensibility Giesekus (FEG) model (Beris & Edwards Reference Beris and Edwards1994), because it presents a steady state extensional behaviour like the e-PTT but a different transient one. We fit their parameters matching the simple shear behaviour of e-PTT. The constitutive equations along with the values of the parameters are given in Appendix A. Although all three models show the same shear behaviour (figure 4b), their extensional one differs. We are interested not only in their steady-state response but also in their transient one because Varchanis et al. (Reference Varchanis, Pettas, Dimakopoulos and Tsamopoulos2021b) concluded that the transient response of the material is responsible for the appearance of the sharkskin instability in the extrusion process. The flow of polymers around the bubble is always transient in a Lagrangian sense due to the upper convective derivative, even if the flow field is steady in a Eulerian setting. So, the transient behaviour of the used models is relevant even in our case, where we solve the steady-state equations. In figure 5(a), we plot the steady-state uniaxial extensional viscosity and in figure 5(b), we plot the startup extensional viscosity for the three models. The I-PTT predicts extension rate hardening and the extensional viscosity reaches a plateau. Instead, the e-PTT and the FEG initially present extension rate hardening, followed by strong extension rate thinning. Regarding the transient response of these models, l-PTT and FEG predict a monotonic increase in the extensional viscosity until it reaches a steady-state value. In contrast, the e-PTT predicts an overshoot in the startup extensional viscosity prior to steady state. These observations will help elucidate the emergence of the knife-edge shape.

Figure 5. (a) Predictions of the steady-state uniaxial elongational viscosity, ${\tilde{\eta }_e}$, versus the extension rate. (b) Time evolution of the normalized startup extensional viscosity at $\widetilde {\dot{\varepsilon }} = 10\ {\textrm{s}^{ - 1}}$. The black, continuous line corresponds to the e-PTT model; the blue, dashed-dotted line corresponds to the l-PTT model; and the red, dashed line corresponds to the FEG model.

4.2. Comparison with bubble shapes reported by Liu et al.

In what follows, the blockage ratio, ${B_R}$, is 10, which is a representative value based on the experiments of LLJ. Only in § 4.7.3, the size of the cylindrical container is varied to show that it affects only the onset of the asymmetric solution, not the flow field. We commence the presentation of our results with a bubble of volume ${\tilde{V}_b} = 0.1\ \textrm{c}{\textrm{m}^3}$ or ${\tilde{R}_b} = 2.9\ \textrm{mm}$ in a 1.5 % polyox solution. The corresponding dimensionless numbers are given in table 2, as case 1.

Table 2. Bubble radius examined and corresponding dimensionless numbers.

Here, $Ar$ is nearly zero, which hints at the minimal inertial effects. Also, $Bo$, which quantifies surface tension effects, is approximately one, and we expect the deformation of the bubble not to be significant. In the bubble problem, the elastically induced deformations are hindered due to surface tension, which tends to preserve the spherical shape of the bubble. Figure 6 shows the experimental shape and on top of it, with a red line, our predicted shape, matching each other nicely. A slight difference is observed regarding the shape of the trailing edge, which is sharper and more pointed in the experiment, while simulations underestimate slightly the elastic response of the material. We attribute this deviation to the lack of rheological data in elongational flow, and the model parameters (table 1) underestimate slightly the elastic response of the material. We use figure 8 of LLJ and the effective radius of the bubble to determine the bubble dimensions. Only one side of the bubble is presented because the shape is axisymmetric. The bubble attains a prolate shape, which is a consequence of the material elasticity. The normal stresses developed at the trailing edge pull the interface downstream. For slightly larger $Eg$, it becomes an inverted teardrop. Note that the value of $Eg$ is smaller than unity but still high enough for elastic effects to manifest themselves. The prolate and the inverted teardrop shapes are reported in numerous experimental and numerical works dealing not only with viscoelastic fluids but also with elastoviscoplastic ones (Moschopoulos et al. Reference Moschopoulos, Spyridakis, Varchanis, Dimakopoulos and Tsamopoulos2021).

Figure 6. For a bubble of ${\tilde{R}_b} = 2.7\ \textrm{mm}$ in the 1.5 % polyox solution, comparison between the steady-state experimental shape by LLJ with our prediction indicated by the continuous red line.

The most striking observation in the work of LLJ and Hassager (Reference Hassager1979) is the appearance of the knife-edge bubble shape. For a bubble with a volume ${\tilde{V}_b} = 0.8\ \textrm{c}{\textrm{m}^3}$ or ${\tilde{R}_b} = 5.7\ \textrm{mm}$ in the 1.5 % polyox solution, the values of the dimensionless numbers corresponding to this experiment are given in table 2, case 2. Note that in our simulations, we use a value of 4.5 for the Bond number for the reasons stated in § 3.2. Inertial effects are still unimportant $(Ar \ll 1)$. Also, one can argue that even the $Eg$ number is not particularly large. For example, $Eg$ was as high as 3 for the onset of the velocity jump discontinuity in the work of Fraggedakis et al. (Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016). Figure 7(a) depicts the predicted shapes with red lines superimposed on the experimental photo at the plane $x = 0$. From this side, the bubble exhibits an inverted teardrop shape with a tip, like those found in earlier works (Astarita & Apuzzo Reference Astarita and Apuzzo1965; Pilz & Brenn Reference Pilz and Brenn2007; Fraggedakis et al. Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016). However, when we turn our attention to the perpendicular side, namely at the plane $y = 0$ (figure 6b), the shape changes and exhibits the peculiar knife edge. Simulations and experiments match each other very well. Based on the terminology of LLJ, the predicted asymmetric shape resembles an axe. To our knowledge, the present simulations are the first ones to predict the breakup of axial symmetry in the rising bubble problem in viscoelastic materials and capture the knife-edge shape accurately. At a first glance, it might seem that we use smaller bubble volumes in simulations than experiments, but this is not true. The bubble shape is not axisymmetric anymore. We compare the numerical prediction against experiments in two perpendicular planes. We have no information of how the shape varies when we move from one to the other, where we expect that the simulation attains a bulkier shape than experiments. Still, the matching between simulations and experiments in the two perpendicular planes is satisfactory. The bottom of the shape presents only a very slight deviation from the horizontal direction. Experiments show a small change of curvature as we approach the knife edge (indicated with black arrows in figure 7b), which is captured by simulation to a smaller extent. The larger value of the experimental Bond number than that used in simulations permits significantly larger deformations of the free surface that are not achievable without increasing prohibitively the computational cost. Nevertheless, our choice of the maximum value of the Bond number does not affect our results so much. Let us comment also on the tip of the trailing edge. Even though it has a cusp-like appearance, as shown in figure 7(a), it is not a true cusp. When we zoom into it (figure 7c), we observe that it is rounded with a very small radius. Our numerical observations agree qualitatively with the theoretical analysis of Jeong & Moffatt (Reference Jeong and Moffatt1992) for cusping-free surfaces. They demonstrate that the radius of the tip becomes extremely small when the importance of surface tension fades away. In our simulations, the radius of the tip is approximately $5 \times {10^{ - 3}}$. Their conclusions illustrate the very small length scales we must capture accurately when the Bond number increases, leading to a tremendous computational cost. We investigate the effect of surface deformation on the knife-edge dynamics in a subsequent section.

Figure 7. For a bubble of ${\tilde{R}_b} = 5.7\ \textrm{mm}$ in the 1.5 % polyox solution, comparison between the steady-state experimental shape by LLJ with our prediction indicated by the continuous red line. (a) Cusp-like tip plane at $x = 0$. (b) Knife-edge plane at the plane $y = 0$. The black arrows indicate the location where the curvature of the bubble interface changes and the bubble becomes concave. (c) Closeup near the cusp-like tip, we identify the viscoelastic fluid with grey and the red line is the interface.

4.3. Flow kinematics

We investigate the spatial variation of all three velocity components, one at a time. Note that we transform the components of the velocity vector from cartesian coordinates to cylindrical ones in what follows. We select a bubble with ${\tilde{R}_b} = 4.9\ \textrm{mm}$, which corresponds to the dimensionless numbers given in table 3, as case 3. Figure 8(a) shows the bubble in an oblique view and the blue circle indicates the region on which we focus. We plot the contours of the radial velocity component, ${u_r}$, at the plane $y = 0$ (figure 8b). As we move from the equator of the bubble towards the knife edge, the radial velocity is negative, which means that material flows towards the $z$-axis. As we approach the point that differentiates between the flat edge and the rest of the bubble, at $z \approx{-} 1.7$, the magnitude of the radial component decreases, remaining negative. However, this observation only holds for a small area near the bubble interface where the material slips along the bubble interface. When we reach the vicinity of the flat edge, the flow kinematics changes rather dramatically. In the region just below the knife edge, the radial velocity becomes positive. A positive radial velocity denotes that material flows away from the $z$-axis, a striking difference from the flow of the material above the trailing edge. Also, the magnitude of the radial velocity which points away from the z-axis is much larger than the one that the material attains before it reaches the knife edge. The maximum positive radial velocity is located below the knife edge and slightly further to the right and left of it. We turn our attention to the ${u_Z}$ velocity component, shown in figure 8(c). Note that we changed the reference frame to the laboratory one. We visualize the velocity field in this reference frame, which means that the bubble moves upwards and the fluid is stationary far from it. The axial velocity is positive near the bubble surface, meaning that the material moves upwards. Along the flat part of the bubble interface, the axial velocity component is equal to the bubble rising velocity. As we move downstream, the axial velocity changes sign and becomes negative. The direction of the flow changes and the material moves away from the bubble. This strange flow configuration has been termed negative wake by Hassager (Reference Hassager1979). Harlen (Reference Harlen2002) determined that the negative wake results from the competition between shear and extensional stresses. To visualize better the spatial variation of ${u_z}$ inside the negative wake structure, we zoom in around the trailing edge and near the z-axis (figure 8d). As explained before, the axial velocity is positive on the bubble interface. Moving downwards, it decreases, becomes zero at the stagnation point of the negative wake structure and then attains a negative sign which identifies the negative wake. It is much more extended in the $x$-direction compared with the findings of axisymmetric solutions (Fraggedakis et al Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016; Niethammer et al. Reference Niethammer, Brenn, Marschall and Bothe2019), in which the negative wake is confined around the $z$-axis. Also, the maximum of the magnitude of the axial velocity does not lie on the z-axis, but is found marginally further to the right and left from it. We do not plot the azimuthal velocity, ${u_\theta }$, because it is zero on the $y = 0$ plane, by virtue of the symmetry condition.

Figure 8. (a) Oblique view of the bubble. The blue circle denotes the area of the knife-edge shape. (b) Contours of ${u_r}$ on the plane $y = 0$. (c) Contours of ${u_z}$ on the plane $y = 0$. (d) Contours of ${u_z}$ on the plane $y = 0$ and zoomed in on the trailing edge and around the z-axis.

Table 3. Critical $Eg$ based on ${A_p}$ for the onset of the instability for ${\varepsilon _{PTT}} = 0.4$ and $Ar = 0$.

Now, we turn our viewpoint to the perpendicular plane, namely $x = 0$. In figure 9(a), we present an oblique view of the bubble just like in figure 8(a), but now the location of the blue circle changes, and we investigate the kinematics around the cusp-like tip. We plot the radial component of the velocity, ${u_r}$, in figure 9(b). Here, we recover a flow field that resembles that found in axisymmetric solutions. The material flows toward the $z$-axis (${u_r}$ is negative), and its maximum value is located at the bubble interface, right before the tip. Note that small positive radial velocities are encountered in a tightly confined area below the trailing edge at the locality of the negative wake. Also, the variation of the axial velocity, ${u_z}$, around the tip in the $x = 0$ plane (figure 9c) resembles that found in axisymmetric simulations. Based on these observations, we conclude that the flow kinematics resembles the already known one of axisymmetric simulations when viewed from the tip side.

Figure 9. (a) Oblique view of the bubble. The blue, dotted outline denotes the area of the cusp-like tip. (b) Contours of ${u_r}$ on the plane $x = 0$. (c) Contours of ${u_z}$ on the plane $x = 0$.

How are the different flow kinematics possible in the two perpendicular planes? The missing part of the puzzle is uncovered when we explore the azimuthal component of the velocity. In figure 10(a), we plot the trailing section of the bubble with a semi-transparent, opaque pink colour and an isosurface of the magnitude of ${u_\theta }(|{u_\theta }|= 0.43)$ on top of it, in blue. Four lobes of finite ${u_\theta }$ are found, signalling that the material rotates below the bubble, breaking the axial symmetry. Note that the transparency of the bubble shape allows us to discern the created lobes behind the bubble and present the three-dimensional (3-D) perspective of the flow field. To investigate further the rotational character of the flow, we plot the contours of the azimuthal velocity (figure 10b) in the plane with $y = 0.02$. The radial and the axial velocity components (not shown here) follow the same pattern found in figure 7. Note that the right-hand side of the figure depicts the simulation results and the left-hand side shows the mirrored ones. Of course, when we mirror our results, we interchange the sign of ${u_x}$ which leads to a change of sign of ${u_\theta }$. The large magnitude of ${u_\theta }$, which is comparable with the magnitude of the other velocity components, does not surprise us. The material rotates near the knife edge and between the two perpendicular planes, namely $x = 0$ and $y = 0$, and the observed change of the flow kinematics becomes possible.

Figure 10. (a) Oblique view of the trailing edge of the bubble, shown in opaque pink, and isosurface of the magnitude of ${u_\theta }$ shown in blue. The isosurface corresponds to $|{u_\theta }|= 0.43$. (b) Contours of ${u_\theta }$ at the plane $y = 0.02$. (c) Contours of ${u_\theta }$ at the plane $z ={-} 1.57$. Black arrows indicate the direction of the flow.

The rotational motion of the fluid is even better visualized if we plot ${u_\theta }$ at the plane $z={-} 1.69$ (figure 10c), located right above the trailing edge. We mention that we visualize the $yx$ plane from above and the azimuthal coordinate is positive in the counterclockwise direction. The superimposed black arrows to the contours of ${u_\theta }$ indicate the direction of the flow. We conclude from figures 10(a) and 10(b) that ${u_\theta }$ is extremely localized near the knife edge. The material approaches the interface of the bubble perpendicularly to the long (knife) side and is deflected by it. This initiates rotation of the liquid, which now moves parallel to the knife side in opposite directions stretching the interface. This rotational motion resembles that found in the four-mill arrangement (Feng & Leal Reference Feng and Leal1997). In another study, Joseph et al. (Reference Joseph, Nelson, Renardy and Renardy1991) investigated experimentally the cusping of a free surface by two counter-rotating cylinders half-immersed in the liquid. The cylinders drag material from the reservoir at their exterior sides and then plunge it in the area between them. For viscoelastic materials and a critical threshold of the rotational velocity of the two cylinders, they reported that the rounded interface, formed at the location where the material re-enters the reservoir, transitions to a cusped one. The created flow field approximates nicely that found in the bubble rise problem, concurring with our argument that this rotational fluid motion generates the knife edge.

Having investigated the different velocity components, we can summarize our findings regarding flow kinematics. In figure 11(a), we present schematically the flow encountered in smaller bubbles translating through viscoelastic materials and just below the trailing edge. The fluid converges towards the $z$-axis from both the x and y directions (blue vectors in figure 11a), squeezing the bubble interface in both directions, and a stagnation point lies at the $z$-axis. At the same time, fluid flows away from the stagnation point in the z direction (yellow arrows in figure 11a). The fluid is compressed in the first and second directions ($x$ and y, respectively) and is elongated in the third one ($z$ direction). So, the flow establishes a uniaxial elongation stretching.

Figure 11. (a) Schematic of the flow field behind the trailing edge corresponding in axisymmetric solution for a smaller bubble. (b) Schematic of the flow field behind the trailing edge corresponding in asymmetric solution for a larger bubble.

In figure 11(b), we present schematically the flow field obtained from our novel numerical results for larger bubbles and in the case of asymmetric solutions. The material flows towards the stagnation point along the y direction, but it flows away from it in the perpendicular one, namely the x direction. This major flow field alteration differentiates the asymmetric from the axisymmetric cases. It is accompanied by the emergence of a non-zero azimuthal velocity component, which is depicted with the rotating, green and brown arrows. Unquestionably, an elastic instability takes place and induces this transition of the flow field. Regarding the flow in the z direction, fluid moves away from the stagnation point, similar to the axisymmetric configuration.

Consequently, the fluid is compressed in one direction $(y)$, while it is elongated in the remaining two ($x$ and $z$). The flow becomes asymmetric and transitions from uniaxial to biaxial stretching after a critical $Eg$. The material continues to squeeze the bubble interface from one side, but it diverts and drags simultaneously the interface in the perpendicular direction. As a result, the bubble forms a tip at one side, which points in the direction of the positive radial velocity component, and the knife-edge shape develops on the perpendicular side, where the radial velocity component is negative. Note that the direction of the tip side is not unique but depends on the type and size of the cross-section of the container, as analysed by LLJ, but this is out of the scope of the present study.

4.4. Flow dynamics

It is imperative to investigate the flow dynamics, specifically how the spatial variations of the stresses evolve at the trailing edge when the bifurcating solution is just formed, which will help us construct the mechanism for the onset of the instability. In figure 12, we plot the contours of the stress components, namely ${\tau _{zz}}$, ${\tau _{rz}}$, ${\tau _{\theta \theta }}$ and ${\tau _{rr}}$ on the plane $y = 0$ for two slightly different bubble radii. The left-hand side of all panels concerns a bubble with ${\tilde{R}_b} = 3.54\ \textrm{mm}$, with dimensionless numbers given in table 2 as case 4. The right-hand side concerns a bubble with ${\tilde{R}_b} = 3.57\ \textrm{mm}$ with dimensionless numbers given in table 2 as case 5, just after the onset of the instability. The former bubble size induces an axisymmetric flow field whereas the latter results in an asymmetric one. Between the two, the latter bubble elongates slightly more, and the extent of its asymmetric deformation is too small to measure. As shown in figure 12(a), the maximum value of ${\tau _{zz}}$ increases from 6 to 9 approximately, which constitutes a 50 % increase. However, it is not localized along the $z$-axis but extends along the $x$-axis also. At the same time, the shear stress, ${\tau _{rz}}$, increases significantly (figure 12b) when $Eg$ is increased due to the deformation of the bubble interface. Its absolute value steps up from 0.6 to 2, a 233 % increase.

Figure 12. Contours of stress components: (a) ${\tau _{zz}}$; (b) ${\tau _{rz}}$; (c) ${\tau _{\theta \theta }}$; (d) ${\tau _{rr}}$. In all panels, the left-hand side shows stresses for ${\tilde{R}_b} = 3.54\ \textrm{mm}$, case 4 in table 2, and the right-hand side shows stresses for ${\tilde{R}_b} = 3.57\ \textrm{mm}$, case 5 in table 2.

They both relax further downstream. In both bubbles, the negative wake structure emerges which is a result of the competition between ${\tau _{zz}}$ and ${\tau _{rz}}$. Still, the most crucial difference between the axisymmetric and the asymmetric cases is found in the hoop stresses, ${\tau _{\theta \theta }}$, in figure 12(c) downstream. Positive values of ${\tau _{\theta \theta }}$ denote that the fluid resists its extension in the azimuthal direction and negative ones that it opposes its compression, just like positive values of ${\tau _{zz}}$ indicates extension of the polymeric chains in the axial direction. Here, $\boldsymbol{\nabla }{u_{\theta \theta }} = (1/r)(\partial {u_\theta }/\partial \theta ) + {u_r}/r$ is responsible for the development of azimuthal stresses. In axisymmetric solution, the first term is almost zero, but hoop stresses are generated due to the radial term. To understand the creation of hoop stresses due to the radial motion, we can image a ribbon consisting of small fluid parcels that form a circle due to axial symmetry. For example, when two parcels have a positive radial velocity, the angle that is formed between them does not change, but their distance with respect to the axis of symmetry increases. As a result, the arclength that connects them, or alternatively their azimuthal distance, increases. So, the material is stretched in the azimuthal direction and positive hoop stresses develop. In contrast, when ${u_r} < 0$, their distance with respect to the axis of symmetry decreases and so does their azimuthal distance; thus, the material is compressed.

Note that a small expansion of the material also develops for the axisymmetric case right after the negative wake, but the positive values of ${\tau _{\theta \theta }}$ are very small and confined near the axis of symmetry (left-hand side in figure 12c), whereas they increase three-fold in the asymmetric case (right-hand side). When the solution becomes asymmetric, we have (a) an azimuthal velocity component, so $\partial {u_\theta }/\partial \theta $ is finite and not zero; and (b) ${u_r}/r$ is larger compared with the axisymmetric case. Both lead to larger $\boldsymbol{\nabla }{u_{\theta \theta }}$ and consequently to larger hoop stresses. Moreover, the maximum ${\tau _{\theta \theta }}$ does not lie along the z-axis but arises away from it. Finally, radial, normal stresses, ${\tau _{rr}}$, increase and become as shown in the right-hand side of figure 12(d), indicating a radial expansion.

Next, we examine the stress distribution downstream of the semi-transparent bubble when the knife edge is readily visible for ${\tilde{R}_b} = 5.7\ \textrm{mm}$, case 2 in table 2. We plot the isosurfaces of ${\tau _{zz}}$ (figures 13a and 13b), ${\tau _{rr}}$ (figure 13c) and ${\tau _{\theta \theta }}$ (figure 13d) on top of the bubble shape. Notice the striking asymmetric stress distribution. In figure 13(a), the isosurface of ${\tau _{zz}}$ (with value ${\tau _{zz}} = 2.2$) attains a butterfly-like shape. Along the z-axis, a birefringent stand is formed, albeit very small, due to the strain-rate thinning nature of the material. Two wings appear, confined along the $x$-axis and parallel to the knife side. These stresses pull the interface downwards, flattening the lower part of the broad side of the knife edge. Here, ${\tau _{zz}}$ is better visualized if we plot its isosurface from an angle that shows the bottom part of the bubble, as shown in figure 13(b), where the nearly planar shape of the stress isosurface is readily viewed. Note that the flow direction points into the plot. Regarding ${\tau _{rr}}$ (figure 13c), we plot an isosurface for ${\tau _{rr}}( = 0.12)$, and we find two larger lobes attached to the bubble surface, at the tip side, and two smaller islands below the bubble. The two small islands correspond to the expansion of the material at the negative wake, which is also observed in the axisymmetric solution (figure 12d). The two larger lobes touching the bubble interface are formed due to radial expansion of the material after its rotation downstream, which extends the bubble along the $x$-axis. Finally, an isosurface of positive hoop stresses $({\tau _{\theta \theta }} = 0.18)$ is plotted in figure 13(d) where four lobes are created. The off-axis hoop stresses shown in figure 12(d) expand even more, strengthen and move further away from the z-axis. Their shape demonstrates that the material expands asymmetrically along the azimuthal direction.

Figure 13. 3-D view of the trailing edge of the bubble, shown in opaque pink and isosurfaces of stress components in blue and from different views: (a) oblique view of ${\tau _{zz}}$ (value $2.2$); (b) bottom view of ${\tau _{zz}}$ (value $2.2$); (c) oblique view of ${\tau _{rr}}$ (value $0.78$) and (d) oblique view of ${\tau _{\theta \theta }}$ (value $0.41$). Case 2 in table 2 with $Ar = 6.2 \times {10^{ - 3}}$, $Eg = 1.37$, $Bo = 5.14$.

4.5. Determination of the onset and the magnitude of the bifurcating solution

Clearly, the instability that causes the flow to transition from axisymmetric to asymmetric is controlled by elastic effects, characterized by the value of $Eg$. We propose two different quantities to indicate whether the instability is triggered or not. The first one is the average of the absolute value of the azimuthal velocity, which we call the azimuthal parameter, ${A_p}$, defined as

(4.1)\begin{equation}{A_p} = \int_{{V_{flow}}} {|{u_\theta }|\,\textrm{d}V} ,\end{equation}

where the integral is calculated in the whole flow domain, ${V_{flow}}$. Here, ${u_\theta }$ is non-zero only just below the trailing edge and it is zero elsewhere.

The second is the absolute value of the difference between the arclength of the entire bubble interface measured at $x = 0$ with the arclength measured at $y = 0$, denoted as ${A_L}$, and called arclength parameter, which reads as

(4.2)\begin{equation}{A_L} = \left|{\int_{\textrm{d}l{|_{x = 0}}} {\textrm{d}l} - \int_{\textrm{d}l{|_{y = 0}}} {\textrm{d}l} } \right|.\end{equation}

In figure 14, we plot ${A_p}$ and ${A_L}$ as functions of $Eg$. When the flow is axisymmetric, both quantities are zero. When the elastic instability is triggered, or equivalently at $Eg = E{g_{cr}}$, they start increasing and a supercritical pitchfork bifurcation is encountered. We encounter a pitchfork bifurcation because a stable solution, namely the axisymmetric solution for $Eg < E{g_{crit}}$, forks into two stable ones, i.e. the asymmetric solutions, and the axisymmetric one becomes unstable. Note that we calculate the azimuthal parameter using the absolute value of the azimuthal velocity so only one bifurcated branch is depicted in figure 14 and the similar ones that follow. The second branch is the mirror image of the first one, i.e. ${A_p} < 0$. In one branch, the knife edge is in the x-axis and in the other, it is in the y-axis, which is the only difference between the two branches. Following the experiments of LLJ, we categorize the asymmetric solution as the stable one for $Eg > E{g_{crit}}$. Also, we categorize the pitchfork bifurcation as supercritical because the bifurcated solution evolves towards larger values of $Eg$ and the asymmetric solutions are stable. Still, a formal stability analysis is needed which is beyond the scope of the present study.

Figure 14. Azimuthal parameter, ${A_p}$, as a function of $Eg$ to be read on the left vertical axis as the black arrow indicates and the arclength parameter, ${A_L}$, to be read on the right vertical axis as the blue arrow indicates.

Let us comment further on the transition from the axisymmetric to the asymmetric solution. If we zoom in the bifurcation point, we notice that the intersection of the upper branch is not perpendicular to the x-axis. This observation is not obvious in figure 14, but becomes evident in the following ones (e.g. figure 20b). This slight deviation in the slope is attributed to the numerical error that the solution method introduces inevitably. Such observation is also reported in the work of Poole, Alves & Oliveira (Reference Poole, Alves and Oliveira2007) (see their figure 3a, for example), who examine the pitchfork bifurcation encountered in the viscoelastic flows through the cross-slot geometry. Furthermore, one should consider that the moving boundary nature of our problem leads to azimuthal disturbances which grow gradually after some value of $Eg$ which is a numerical artifact. This behaviour decreases as we refine the mesh. Regarding the base case, namely the axisymmetric solution, we capture it even for $Eg > E{g_{crit}}$, but not much larger. To allow the iterations to converge to the base solution, we do not perturb it and we do not use arc-length continuation. However, the mesh deformation introduces, in effect, the needed perturbation in the third direction and forces the Newton method to converge to the asymmetric solution. Nevertheless, the axisymmetric solution for $Eg > E{g_{crit}}$ can be obtained enforcing axial symmetry.

The rate of increase of the two parameters determining the magnitude of the asymmetry is different. The flow field becomes asymmetric, but the rotational motion is not strong enough to deform greatly the bubble interface asymmetrically. When the azimuthal velocity increases due to larger $Eg$, the parameter that measures the asymmetric deformation of the bubble interface grows and the knife-edge shape is formed. The asymmetry of the flow field is detected more easily when $Eg$ is slightly larger than $E{g_{cr}}$. However, the use of the arclength parameter is not always practical or even possible in experiments. For instance, the bubble deformation might be too weakly asymmetric to make the required optical measurements. The difference between the two arc lengths should become large enough to identify the asymmetric deformation which occurs in $Eg > E{g_{cr}}$. Based on these arguments, if we try to calculate experimentally the value of $E{g_{cr}}$ for the onset of the instability using only the difference in the two arc lengths or possibly some other geometrical criterion, we may overestimate its value.

4.6. Dependence of the bubble rise velocity on its effective radius

Before we move on to a parametric study, we deem it worthwhile to examine the dependence of the speed of the bubble on its volume. Following the work of LLJ, we define the capillary number, which depends on the velocity of the bubble, and it reads as

(4.3)\begin{equation}Ca = Bo\ast {U_b}.\end{equation}

We visualize its dependence on the volume of the bubble in figure 15(a). The experimental data denoted with symbols show a clear jump for a volume $0.1\ \textrm{c}{\textrm{m}^3}$ accompanied by the appearance of the knife-edge shape. However, the simulations do not capture it. The velocity predicted from simulations lies between the upper and the lower branches of the experimental results. Still, the difference between the two is not that large and simulations follow experiments, excluding the presence of the jump discontinuity. After a volume of $0.4\ \textrm{c}{\textrm{m}^3}$, the slope of the predicted velocity increases faster and starts to catch up with the values of the experiments. Around the volume $0.9\ \textrm{c}{\textrm{m}^3}$, a turning point is approached, but our simulations cannot continue further due to convergence issues. We assume that this turning point is the beginning of the hysteresis loop that precedes the velocity jump discontinuity as identified by Fraggedakis et al. (Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016). Moreover, experiments determine that the knife edge arises approximately for a volume $0.1\ \textrm{c}{\textrm{m}^3}$ or ${\tilde{R}_b} = 0.29\ \textrm{cm}$ (star symbol in figure 15a), whereas simulations indicate a critical volume of $0.18\ \textrm{c}{\textrm{m}^3}$ or ${\tilde{R}_b} = 0.35\ \textrm{cm}$ (triangular symbol in figure 15a). We attribute this deviation to two reasons: first, the extensional nature of the material is not correctly characterized because extensional, rheological experiments were not performed. Thus, we may underestimate its elastic response as we argue in relation to figure 6(b). Second, the distance between the bubble and the sidewall of the container becomes smaller when the volume of the bubble increases in experiments, whereas we keep it constant in simulations. As we explain in § 5.7.3, the proximity of the container to the bubble affects only the onset of the instability.

Figure 15. (a) Capillary number as a function of the bubble volume and (b) drag coefficient as a function of the bubble volume. Comparison between our theoretical results, denoted by the line, and the experimental results, denoted by the square symbols, for the 1.5 % polyox solution of LLJ. Critical volumes for the occurrence of the knife edge in experiments and simulations are indicated by the star symbol (experiments) and by the triangular symbol (simulations).

The drag coefficient, ${C_d}$, for a steadily rising bubble is given by

(4.4)\begin{equation}{C_d} = \frac{{{{\tilde{F}}_d}}}{{\left( {\dfrac{{{{\tilde{\eta }}_0}{{\tilde{U}}_b}}}{{{{\tilde{R}}_b}}}} \right)({\rm \pi}\tilde{R}_b^2)}} = \frac{4}{{3{U_b}}},\end{equation}

where ${\tilde{F}_d}$ is the drag force on the bubble, which is balanced by buoyancy in equilibrium. In figure 15(b), we compare the calculated ${C_d}$ against the experimental one as a function of the bubble volume. The same mismatch found in figure 15(a) is observed also in the case of ${C_d}$. Experiments predict that the bubble obtains the knife-edge shape and simultaneously its velocity jumps one order of magnitude, which translates in the sudden decrease of ${C_d}$. We predict a milder decrease of ${C_d}$ until the volume $0.9\ \textrm{c}{\textrm{m}^3}$, where an abrupt decrease is indicated by simulations. Thus, the rapid decrease of ${C_d}$ happens after the onset of the instability, a result which holds throughout our analysis.

4.7. Parametric analysis

Up to this point, we examined different bubble radii, compared our novel numerical results with the experiments of LLJ, and analysed the flow kinematics and dynamics captured by numerical simulations. These cases resulted in axisymmetric and asymmetric bubble shapes, where the teardrop bubble shape arose in smaller bubbles, but the peculiar knife-edge bubble shape was obtained in larger bubbles.

Still, some of the questions posed in the introduction remain unanswered. The following parametric study assists in elucidating them, after which we will propose a physical mechanism for the observed instability, which leads to the knife-edge bubble shape. Also, we modify the solution procedure. We do not vary ${R_b}$ anymore because it affects all dimensionless numbers and does not allow us to isolate the effect of the different parameters on the flow. Instead, we fix the values of $Ar$ and $Bo$, and perform an arclength continuation on the $Eg$ number. In all forthcoming cases, the inertial effects are minimal, because $Ar$ is smaller than ${10^{ - 3}}$. Thus, we set $Ar$ equal to zero for all, unless we state otherwise. Also, ${B_R} = 10$ for the majority of the following results, except for § 4.7.3, where we investigate the effect of the confinement.

4.7.1. Effect of the Bond number, formation of a bubble with cross-edge instead of a knife edge

We investigate the effect of surface tension by varying the Bond number. In figure 16(a), we plot the azimuthal parameter versus $Eg$ for five different Bond values. The azimuthal parameter starts from 0, indicative of an axisymmetric flow field, and it increases suddenly at the onset of the flow instability, which occurs at a critical $Eg$ value. Larger $Eg$ is needed to trigger the instability when $Bo$ is smaller, because surface tension opposes the flow-induced deformation. When $Eg$ increases, the stresses exerted on the trailing edge perturb the spherical shape. The bubble obtains an inverted teardrop shape and finally the flow field becomes asymmetric. It seems that only after the inverted teardrop shape is formed does the instability emerge. We summarize the critical $Eg$ values for the onset of the instability in table 3.

Figure 16. (a) Azimuthal parameter, ${A_p}$ as a function of $Eg$ for different Bond values, (b) the arclength parameter, ${A_L}$, for different Bond values.

In figure 16(b), we plot the arc length difference defined in (4.2) for different Bond numbers. As stated, the arclength parameter increases less and attains smaller values than the azimuthal parameter. When $Bo < 1.5$, ${A_L}$ attains very small values even when $Eg > 1.3$. Experiments do not easily capture these small values because the determination of the bubble interface is prone to image processing errors. For $Bo \ge 2$, the arclength parameter increases and should be measurable experimentally. However, it shows a different dependence on $Eg$ for $Bo = 1.5$. It increases slightly (not observable in figure 16(b) due to the ordinate scale) and remains almost constant for the range of $Eg$ we investigated. To elucidate this behaviour, we examine the flow kinematics.

In figure 17(a), we plot the contours of the radial velocity where the left-hand side corresponds to a slice at $x = 0$ and the right-hand side corresponds to a slice at $y = 0$. Positive radial velocities are found in both slices, and the material stretches the interface along both the x and y planes simultaneously, which indicates a change of the flow field compared with the case shown in figures 8 and 9. There we found that the material pulled the interface only along only one plane ($x$ or $y$) when the bubble attains the knife-edge shape. In figure 17(b), we plot side by side the plane at $y = 0$(right-hand side) and a slice that corresponds to an angle $\theta = 45^\circ $ (left-hand side); positive radial velocities are found in the former, but negative in the latter. So, the dihedral angle between the two planes in which the flow kinematics changes is $45^\circ $, corresponding to a reduction of the angle by 45° compared with that previously shown in figure 9. To illustrate this change better, we also plot the azimuthal velocity, ${u_\theta }$, at a plane z = −1.428, shown in figure 17(c). Here, the plane of $\theta = 45^\circ $ serves as an internal plane of symmetry and the fluid converges towards the bubble surface with zero azimuthal velocity. Hence, the material rotates simultaneously in the x and $y$ directions, which results in minimal differences in the two-interface arc-lengths calculated at the planes with $x = 0$ and $y = 0$. The odd stretching of the bubble interface prevents the easy distinction of the asymmetric shape. Still, small differences arise if we compare the arc-length of the bubble interface at plane $x = 0$ and the plane that corresponds to $\theta = 45^\circ $. However, the asymmetric shape is better discerned when we view it from a z plane, like in figure 17(c). As a result, the knife edge is not formed and the shape resembles a deformed cross. Still, the deformation is not as pronounced as in the other cases and the asymmetry in the bubble shape might be overlooked in experiments, as opposed to the flow around a sphere in a similar fluid (Varchanis, Haward & Shen Reference Varchanis, Haward and Shen2023).

Figure 17. (a) Contours of the radial velocity, ${u_r}$. The right-hand side of the panel corresponds to the plane $y = 0$ and the left-hand side corresponds to the plane $x = 0$. (b) Contours of the radial velocity, ${u_r}$. The right-hand side corresponds to the plane $y = 0$ and the left-hand side corresponds to the plane $\theta = 45^\circ $ or $x = y$. (c) Contours of the azimuthal velocity, ${u_\theta }$, at the plane $z ={-} 1.428$. Black dotted arrows indicate the direction of the flow and the dotted lines denote the two planes in which the fluid moves towards the $z$-axis with zero azimuthal velocity. $Bo = 1.5$, $Eg = 1.3$.

Following these observations, we illustrate ${\tau _{zz}}$ and ${\tau _{\theta \theta }}$ using isosurfaces on top of the bubble shape in figure 18. The spatial variation of stresses is different compared with that in figure 13. The ${\tau _{zz}}$ transitions from a butterfly-like to a cross-like configuration, stretching the interface downwards along two perpendicular planes. These planes form an angle of 45° with the $x$- and $y$-axis and their common contact point lies along the $z$-axis. Also, a very small birefringent strand exists and extends along the z-axis (figure 18a). Regarding ${\tau _{\theta \theta }}$, four more lobes are created perpendicular to the already existing ones (figure 13c).

Figure 18. 3-D view of the trailing edge of the bubble, shown in opaque red, and isosurfaces of stress components in blue and from different views: (a) oblique view of ${\tau _{zz}}$ (value $0.52$); (b) bottom view of ${\tau _{zz}}$ (value $0.52$); (c) oblique view of ${\tau _{\theta \theta }}$ (value $0.18$). $Bo = 1.5$, $Eg = 1.3$.

4.7.2. Effect of extensional viscosity

As discussed previously, when the knife-edge shape is formed, the flow type changes from uniaxial elongational to biaxial stretching behind the trailing edge. Even though the onset of the elastic instability renders the flow field three-dimensional, the flow character remains predominantly extensional. Also, the experiments of LLJ demonstrate that bubbles obtain the knife-edge shape not only in shear-thinning materials, like polyox, but also in Boger-type fluids, which present a constant shear viscosity. So, shear thinning does not affect the creation of the knife edge and we focus safely on the extensional viscosity. To examine the effect of the extensional viscosity on the dynamics, we use the different viscoelastic models mentioned in § 4.1 and in more detail in Appendix A.

For all the investigated cases in this section, we use $Bo = 2$. In figure 19(a), we plot the azimuthal parameter as a function of $Eg$. We observe that the cases of e-PTT and FEG demonstrate the same behaviour qualitatively, albeit the difference in the values of the critical $Eg$ for the onset of the instability. The e-PTT predicts $E{g_{crit}} = 0.78$, whereas FEG predicts $E{g_{crit}} = 1.11$. Their steady-state, extensional viscosity shows similar dependence on the extension rate, but the strain-rate thinning response of FEG arises for a larger extension rate. Furthermore, their transient extensional viscosity presents different behaviour. The e-PTT model predicts an overshoot in the startup extensional viscosity before a steady-state is attained, in contrast to the FEG model, which predicts a monotonic increase of the startup extensional viscosity. As seen in the work of Varchanis et al. (2021), this overshoot serves as a destabilizing agent. In our case, this overshoot, coupled with the earlier strain-rate thinning response of the e-PTT material, seems to promote the instability to occur at smaller $Eg$ values. To our surprise, the l-PTT fluid does not introduce any three-dimensional character to the flow field; the asymmetry parameter remains equal to zero for the whole range of the investigated $Eg$ numbers and the solution is axisymmetric. All models predict slight strain-hardening effects from $\widetilde {\dot{\varepsilon }} = {10^{ - 2}}$ to $1\ {\textrm{s}^{ - 1}}$, but only the l-PTT predicts the extensional viscosity to be constant and independent of the extension rate for larger values $(\widetilde {\dot{\varepsilon }} \ge 10\ {\textrm{s}^{ - 1}})$, see figure 5. Furthermore, we perform simulations using two more models: (a) the FENE-CR model, which exhibits constant shear viscosity, extension-rate hardening and bounded extensional viscosity; and (b) the Oldroyd-B model, which predicts an infinite increase of the extensional viscosity when the strain rate becomes large enough. The FENE-CR model is presented in Appendix A. Nonetheless, axisymmetric solutions (not shown here) are obtained with both models. Thus, extension-rate thinning is necessary for the onset of the instability.

Figure 19. (a) Dependence of the azimuthal parameter, ${A_p}$, on $Eg$ for different constitutive models. (b) Dependence of the azimuthal parameter, ${A_p}$, on $Eg$ for different values of ${\varepsilon _{PTT}}$. $Bo = 2$.

We vary the ${\varepsilon _{PTT}}$ parameter to assess the effect of the strain-hardening and we find that the critical $Eg$ number depends strongly on it (figure 19b). As we decrease ${\varepsilon _{PTT}}$, the rate-of-strain thinning starts at larger extension rates. So, the elastic stresses (or equivalently, the extensional viscosity) become greater and deform the bubble interface more. This increased deformation triggers the instability sooner, just like in the case of larger $Bo$ (figure 16a). Also, ${\varepsilon _{PTT}}$ affects the growth of the instability. For large values, the slope of ${A_p}$ is steeper and rises almost vertically, which demonstrate that extension-rate thinning amplifies the asymmetry of the flow. The changes that ${\varepsilon _{PTT}}$ induces in the values of $E{g_{crit}}$ are summarized in table 4.

Table 4. Critical $Eg$ based on ${A_p}$ for the onset of the instability for $Bo = 2$ and $Ar = 0$.

We note that You et al. (Reference You, Borhan and Haj-Hariri2009) investigated the bubble rise in a viscoelastic material by performing a stability analysis of axisymmetric simulations. They determined that the flow becomes unstable and a three-dimensional configuration is stable when they increased $Eg$ (termed $Wi$ in their paper). Still, they use a FENE-CR material, which contradicts our observations. However, we believe this is not a true outcome of the governing equations, but results from the deterioration of the accuracy of their simulations. The extreme mesh deformation at the trailing pole of the bubble necessitates a very fine discretization at the bubble tip, like that shown in figure 3. If the mesh is not adequate, then the side of the last element on the bubble interface is not perpendicular to the axis of symmetry but forms a small angle. Thus, the assumption of axial symmetry is not imposed accurately and wrong results may be obtained.

4.7.3. Effect of confinement

Studies of viscoelastic flows demonstrate that the behaviour of the flow dynamics is sensitive to relevant dimensionless geometric parameters (Varchanis et al. Reference Varchanis, Hopkins, Shen, Tsamopoulos and Haward2020b). In our case, the relevant scaling measure is the bubble-cylinder ratio or the blockage ratio, ${B_R}$, defined in (2.1). To investigate the scaling of the critical $Eg$ number with the flow geometry, we perform additional simulations with a series of containers of smaller radii. Following the previous observations, we set $Bo = 2.5$ to promote the onset of the instability. The azimuthal parameter demonstrates a non-intuitive dependence on $Eg$. Previous works report that the increase of the blockage ratio promotes the destabilization of the flow and instabilities should occur at smaller $Eg$ numbers, which is opposite to what we observe in figure 20(a). We note that the confinement changes only the onset of the instability and the flow field remains qualitatively the same. This stems from the circular cross-section of the container. As LLJ state, if the cross-section of the container is rectangular then the orientation of the knife-edge shape depends on the proximity of the bubble to the walls. Nonetheless, we do not pursue further the flow alteration that a rectangular container might induce.

Figure. 20. (a) Azimuthal parameter, ${A_p}$, as a function of $Eg$ for different blockage ratios. (b) Azimuthal parameter, ${A_p}$, as a function of $Wi$ for different blockage ratios. $Bo = 2.5$.

We realize that the bubble velocity or alternatively the upstream velocity of the fluid depends on the blockage ratio in our simulations because we do not impose a prescribed value a priori, but we determine the bubble speed along with the solution of the problem. End-wall effects affect the motion of bubbles and spheres (Sigli & Coutanceau Reference Sigli and Coutanceau1977). As the boundaries are placed closer to the bubble, its velocity decreases. So, $Eg$ is not an appropriate measure of elasticity in the present situation because it does not account for the flow intensity. To this end, we define the Weissenberg, $Wi$, which depends on the bubble velocity and relates to $Eg$ as

(4.5)\begin{equation}Wi = Eg\ast {U_B} \equiv \tilde{\lambda }\frac{{{{\tilde{U}}_b}}}{{\tilde{R}}}.\end{equation}

When we examine the dependence of the azimuthal parameter on $Wi$ (figure 20b), our results concur with those made in earlier investigations of flows exhibiting elastic instabilities (Varchanis et al. Reference Varchanis, Hopkins, Shen, Tsamopoulos and Haward2020b). For larger blockage ratios, the onset of the elastic instability occurs at smaller $Wi$. Nonetheless, the blockage ratio does not induce significant differences in the critical value of $Wi$. This result is also supported by LLJ, who show that the velocity and the cusp-like shape of the bubble do not depend much on the size of the container. Indeed, the container shape and size affect primarily the orientation of the knife edge, as observed experimentally, but we will not investigate this in the present study.

The results presented here correspond to $Ar = 0$ (creeping flow). Hence, the instability we investigate is categorized as purely elastic. In their pioneering works, McKinley, Pakdel & Öztekin (Reference Mckinley, Pakdel and Öztekin1996) and Pakdel & McKinley (Reference Pakdel and Mckinley1996) show that both the geometric (i.e. blockage ratio) and rheological (i.e. $Wi$) scaling can be unified into a single criterion to determine the onset of purely elastic instabilities. The proposed scaling incorporates the non-zero tensile stresses and the magnitude of the streamline curvature in the flow, and can be expressed as

(4.6) \begin{equation}\left[ \frac{{\tilde{\lambda }{{\tilde{U}}_b}}}{{\tilde{R}}}\frac{{{{\tilde{\tau }}_{ss}}}}{{\tilde{\eta }}_0\dot{\tilde{\gamma}}} \right] \ge M_{crit}^2,\end{equation}

where ${\tilde{\tau }_{ss}}$ is the tensile stress along a streamline and $\tilde{R}$ is the characteristic radius of curvature of the streamline. If the left-hand side of (4.6) exceeds the $M_{crit}^2$ at some area of the flow field, the flow is prone to an instability that originates from that area.

In figure 21, we plot the inverse of the critical $Wi$ versus the blockage ratio, which shows clearly the linear dependence of $1/W{i_{cr}}$ on ${B_R}$, as argued by McKinley et al. (Reference Mckinley, Pakdel and Öztekin1996) and Pakdel & McKinley (Reference Pakdel and Mckinley1996). Thus, any initial perturbation in the flow field behind the trailing edge of the bubble is amplified by the presence of strong elastic normal stresses along curved streamlines in the vicinity.

Figure 21. Inverse of the critical $Wi$ for the onset of the instability as a function of the blockage ratio.

Following the earlier works, we relate the radius of curvature through a reciprocal linear combination of the radius of the bubble, ${\tilde{R}_b}$, and the radius of the cylindrical container, ${\tilde{R}_c}$, which reads as follows:

(4.7)\begin{equation}\frac{1}{{\tilde{R}}} = \frac{\alpha }{{{{\tilde{R}}_b}}} + \frac{\beta }{{{{\tilde{R}}_c}}}.\end{equation}

As argued by McKinley et al. (Reference Mckinley, Pakdel and Öztekin1996), the radius of curvature of the streamline is a more appropriate length scale of the flow field, rather than the radius of the bubble or of the cylindrical container. However, streamline curvature should scale with a combination of the two, and the simplest relationship is that employed in (4.7). Suppose now that in the vicinity of the trailing edge, $\dot{\tilde{\gamma}} $ can be approximated as $\dot{\tilde{\gamma}} \approx {\tilde{U}_b}/\tilde{R}$ and that ${\tilde{\tau }_{ss}} = 2\tilde{\lambda }{\tilde{\eta }_0}{\dot{\tilde{\gamma}}^2}$, which is the normal stress for an Olroyd-B fluid in shear flow. We combine (4.7) and the former arguments in the dimensionless instability criterion (4.6). After certain rearrangements, the following expression is obtained:

(4.8)\begin{equation}\frac{1}{{W{i_{cr}}}} = \frac{{\sqrt 2 \alpha }}{{{M_{crit}}}} + 2{B_R}\frac{\beta }{{{M_{crit}}}} \equiv \frac{{\bar{\alpha }}}{{{M_{crit}}}} + {B_R}\frac{{\bar{\beta }}}{{{M_{crit}}}}.\end{equation}

We wish to determine the value of ${M_{crit}}$ along with the numerical factors $\bar{\alpha }$ and $\bar{\beta }$. To do so, we use a linear polynomial to fit the simulation data depicted in figure 20. Note that we can obtain only two constants from the regression procedure for the slope and the intercept. Without loss of generality, we set $\bar{\beta } = 1$ and in the limit of ${B_R} \to 0$, the characteristic radius of curvature reduces to $\tilde{R} = {\tilde{R}_c}$ and the shear flow reads as $\dot{\tilde{\gamma}} \approx {U_b}/{\tilde{R}_c}$. The fitted values of ${M_{crit}}$, $\bar{\alpha }$ and $\bar{\beta }$ are summarized in table 5.

Table 5. Fitted values of the parameters in (4.8)

4.8. Physical mechanism leading to the asymmetric shapes

In figure 22, we give a visual representation of the physical mechanism that leads to asymmetric flow fields downstream of the bubble and the appearance of the knife-edge shape. We choose to split it into three consecutive stages. During stage one, the increase of the bubble volume induces a stronger elastic response of the fluid due to larger induced velocities and strains, resulting in stronger extensional forces in the trailing edge of the bubble. The elastic stresses pull downwards the rear pole deforming the bubble, and the inverted teardrop shape is formed.

Figure 22. Schematic representation of the sequence of events that lead to the creation of the knife-edge shape. The first two stages are relevant also in axisymmetric shapes while the third stage corresponds to fully 3-D simulations.

In the second stage, the inverted teardrop shape permits the development of large shear stresses close to the bubble interface, even though the shear-free condition is applied along it. A delicate balance between the tensile and shear viscoelastic stresses may develop in the region right after the bubble, which gives rise to the negative wake structure. Harlen (Reference Harlen2002) argues that the tensile stresses should not be much larger than the shear ones. Otherwise, their balance is not established and the negative wake does not appear. As a consequence of the emergence of the negative wake, small and positive radial velocity arises, which is confined in the region of the negative wake (figure 9c).

Stages one and two are encountered in both axisymmetric and 3-D shapes. Due to the shape of the bubble, the streamlines are curved at the trailing edge and large tensile stresses are exerted along them, a combination that is known to destabilize flows of viscoelastic materials (Pakdel & McKinley Reference Pakdel and Mckinley1996). Typically, the destabilization begins when a small radial perturbation thrusts the polymer chain not to lie along a single streamline and the shearing motion between the streamlines stretches the chain nonuniformly amplifying the non-Newtonian ‘hoop stresses’, namely the $\theta \theta $ component of the stress tensor. This necessarily positive radial thrust is present when the negative wake is formed behind the bubble. Thus, all the prerequisites for triggering an elastic instability exist, making stage three possible.

However, not all bubbles rising in viscoelastic materials obtain asymmetric shapes, despite the presence of a negative wake. So, all these ingredients are not sufficient to form the knife-edge shape. Based on the discussion presented in § 4.6.2, the fluid should also exhibit extension-rate hardening for intermediate strain rates and thinning for larger ones for the onset of the elastic instability, which allows the knife edge to form. Let us expand on the necessity of both strain hardening and thinning effects in the constitutive model. The strength of the tensile stresses should increase initially in response to the induced strain and deform the bubble. The prolate or inverted teardrop shape of the bubble prompts the development of shear stresses at the trailing edge. As the bubble size increases (or alternatively the strain rate increases), tensile stresses should decrease compared with the shear stresses and allow the destabilizing effect of the latter to amplify the azimuthal perturbations, as indicated by Pakdel & McKinley (Reference Pakdel and Mckinley1996).

We believe that this very subtle interplay between tensile and shear stresses is the reason behind the limited experimental evidence of the knife-edge bubble shape and the non-existence of numerical simulations that demonstrate it. This mechanism explains clearly the effect of $Bo$ and ${\varepsilon _{PTT}}$ on the instability. The more the tensile stresses deform the bubble interface, the more the shear stresses increase near the trailing edge, which are convected by the flow and trigger the instability sooner. In general, the deformation of the bubble is a result of either weaker surface tension effects (or equivalently larger Bond number) or stronger tensile stresses (smaller ${\varepsilon _{PTT}}$) that pull the interface downward. Note that after the emergence of small azimuthal perturbations, ${\varepsilon _{PTT}}$ amplifies their growth and larger values facilitate the steeper increase of ${A_p}$ (figure 18b).