1. Introduction
The rise velocity of a bubble in a surfactant solution can be reduced to half its velocity in pure water (Magnaudet & Eames Reference Magnaudet and Eames2000; Takagi & Matsumoto Reference Takagi and Matsumoto2011), which is the same as that of a solid sphere. This behaviour can be traced back to mechanisms first elucidated by Frumkin & Levich (Reference Frumkin and Levich1947) and Levich (Reference Levich1962). As bubbles rise at a high Péclet number
$Pe$
(
$={\rm d}U/D$
, where
$d$
is the bubble diameter,
$U$
is the bubble velocity and
$D$
is the diffusion coefficient), surfactants at the bubble’s surface move from the front to the rear owing to surface advection. Consequently, a higher surfactant concentration occurs at the rear of the bubble, creating a gradient of surface tension along the bubble’s surface. Unlike the typical gas–liquid boundary conditions that have zero shear stress, here, the stress occurs owing to the surface tension gradient, known as the Marangoni effect. This modifies the flow around the bubble, increasing its drag. In many scenarios, even a minimal surfactant contamination of the order of a few ppm can cause the drag force on the bubble to be comparable to that of a rigid sphere, as reported by Harper (Reference Harper1972) and Clift, Grace & Weber (Reference Clift, Grace and Weber1978).
The variables influencing the velocity of a bubble were studied by Bachhuber & Sanford (Reference Bachhuber and Sanford1974), who observed that a bubble for Reynolds number
$\textit{Re}$
(
$=\rho {\rm d}U/\mu$
, where
$\rho$
is the liquid density,
$\mu$
is the dynamic velocity)
$\lt$
40 in water begin its ascent with the velocity of a sphere under a free-slip boundary condition, but as the bubble moved upwards, it rises with the velocity of a solid sphere. This indicates that the surfactant concentration on the bubble’s surface increases with distance travelled, resulting in an increase in drag. Notably, Zhang & Finch (Reference Zhang and Finch2001) showed that, at low surfactant concentrations, the distance to reach the terminal velocity of a solid sphere can span several metres because of the slow pace of adsorption and desorption of the surfactant.
The aforementioned mechanism is modelled in the so-called stagnant cap model (Savic Reference Savic1953; Davis & Acrivos Reference Davis and Acrivos1966; Harper Reference Harper1982; Sadhal & Johnson Reference Sadhal and Johnson1983). In this model, two boundary conditions are considered, separated by the cap angle
$\theta _C$
(
$\theta _C$
corresponds to the leading edge of the stagnant region and is defined here from the front stagnant point): zero shear stress in the bubble front up to
$\theta _C$
and a no-slip condition afterwards. Sadhal & Johnson (Reference Sadhal and Johnson1983) formulated a drag equation as a function of cap angle in creeping flow, demonstrating significant drag force variations within a specific angle range (
$60^\circ \lt \theta _{C} \lt 150^\circ$
). The flow structure in the vicinity of a bubble has been investigated primarily using numerical analysis (e.g. Cuenot, Magnaudet & Spennato Reference Cuenot, Magnaudet and Spennato1997; Takemura & Yabe Reference Takemura and Yabe1999; Ponoth & McLaughlin Reference Ponoth and McLaughlin2000; Dani et al. Reference Dani, Cockx and Guiraud2006, Reference Dani, Cockx, Legendre and Guiraud2022; Pesci et al. Reference Pesci, Weiner, Marschall and Bothe2018; Kentheswaran et al. Reference Kentheswaran, Dietrich, Tanguy and Lalanne2022, Reference Kentheswaran, Antier, Dietrich and Lalanne2023). Intriguingly, many studies reported that this drag increase as a function of the cap angle remains consistent even at intermediate Reynolds numbers (
${\textit{Re}}$
), except for the specific angle range where the slope is slightly steeper than that predicted by the model (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997).
The change in the boundary conditions owing to the surfactant has other important effects in addition to its effect on the velocity of the bubble. Tagawa, Takagi & Matsumoto (Reference Tagawa, Takagi and Matsumoto2014) experimentally investigated the effects of surfactants on the path instability of a single bubble. They observed that when
${\textit{Re}}$
is fixed at 400, the bubble motion transitions from rectilinear to a helical spiral and further to a zigzag motion with an increase in the drag coefficient, corresponding to increased surfactant concentration. This mechanism is explained by changes in vorticity generation owing to altered boundary conditions and vorticity stretching/tilting. Therefore, the extent and location of boundary condition changes, and the effect of these changes on the flow field in the vicinity of an object, should be clarified.
However, experimental observations on stress fields near a bubble/sphere moving in a dilute surfactant solution are quite limited. The measurement of a flow field in the vicinity of a bubble remains a significant challenge owing to the thin boundary layer on the bubble’s surface. The classical visualisation technique, i.e. particle image velocimetry (PIV), is particularly difficult because it uses micron-sized fluid tracers. Hosokawa et al. (Reference Hosokawa, Masukura, Hayashi and Tomiyama2017) successfully extracted the Marangoni effect of a settling droplet using spatio-temporal filter velocimetry with boundary-fitted measurement. Their method connected the internal and external flows of the droplet and enabled detailed analysis of interfacial stress transfer. However, the direct visualisation of the Marangoni effect around a gas bubble remains challenging, especially because internal flow tracking is not feasible. More importantly, to locate the location of the cap angle, we would prefer to visualise the stress field. Although fluid stresses can be estimated from PIV data, the obtained stress data are often noisy because estimating stresses from PIV requires spatial derivatives of the velocity field.
In this study, we experimentally examined the stress field around a rising bubble in a dilute surfactant solution using photoelastic measurement, a recently developed flow visualisation method employing nano-sized rods (Nakamine et al. Reference Nakamine, Yokoyama, Worby, Muto and Tagawa2024). To validate this technique, we first analysed the stress field around a solid sphere sinking in a quiescent liquid with terminal velocity, comparing our results with well-established numerical simulation. This comparison also served to characterise the stress field simulating a fully contaminated bubble. Subsequently, we investigated a bubble rising rectilinearly in a quiescent surfactant solution at intermediate
${\textit{Re}}$
(
$20 \lt {\textit{Re}} \lt 220$
) and high Péclet numbers
$Pe = O(10^3{-}10^5)$
. By observing bubbles at various distances travelled, we characterised the stress field near the bubble’s surface, where the boundary condition transitioned from a clean to a contaminated state. The cap angle was then estimated from the experimental result. Finally, we discuss the relevance of the cap angle in the stagnant cap model and previously reported and previous numerical study by Cuenot et al. (Reference Cuenot, Magnaudet and Spennato1997).
2. Experimental and numerical methods
2.1. Experimental methods
2.1.1. Experimental set-up
We observed the flow structure around a bubble and a solid particle in a cellulose nanocrystal (CNC) suspension that exhibits photoelasticity (Lane et al. Reference Lane, Rode and Rösgen2022; Nakamine et al. Reference Nakamine, Yokoyama, Worby, Muto and Tagawa2024; Worby et al. Reference Worby, Nakamine, Yokoyama, Muto and Tagawa2024). Flow birefringence occurs when non-spherical nanoparticles align in a certain orientation under the effects of shear (Calabrese, Haward & Shen Reference Calabrese, Haward and Shen2021). By controlling the polarisation state before flow birefringence and measuring the polarisation state after, we can determine the phase retardation
$\varDelta$
and azimuth
$\phi$
. According to the stress-optic law,
$\varDelta$
corresponds to the magnitude of the principal stress difference, and
$\phi$
to its orientation.
Figure 1 shows a schematic of the experimental set-up. Back light generated by a light source (SOLIS-525C, Thorlabs) with a wavelength
$\lambda$
of 520 nm was circularly polarised through a linear polariser and a quarter-wave plate. The circularly polarised light passed through the flow around the object, resulting in the emission of elliptically polarised light with
$\varDelta$
and
$\phi$
. A high-speed polarisation camera (CRYSTA PI-1P, Photron Ltd.) was used to measure
$\varDelta$
and
$\phi$
. The camera had linear polarisers whose angles differed by
$45^\circ$
at neighbouring pixels (Onuma & Otani Reference Onuma and Otani2014). Using the intensity distribution of the four pixel set (superpixels), we derived the phase retardation and azimuth. Applying the four-step phase-shifting method, we described
$\varDelta$
and
$\phi$
using the intensities measured through each linear polariser, denoted as
$I_0, I_{45}, I_{90}$
and
$I_{135}$
(Lane et al. Reference Lane, Rode and Rösgen2022; Yokoyama et al. Reference Yokoyama, Mitchell, Nassiri, Kinsey, Korkolis and Tagawa2023; Nakamine et al. Reference Nakamine, Yokoyama, Worby, Muto and Tagawa2024)
\begin{align}& \varDelta = \frac{\lambda }{2\pi }\sin ^{-1}{\frac{\sqrt {(I_{90}-I_0)^2 + (I_{45}-I_{135})^2}}{I/2}}, \end{align}
\begin{align}&\qquad\qquad \phi = \frac{1}{2}\tan ^{-1}{\frac{I_{90}-I_0}{I_{45}-I_{135}}}, \end{align}
The schematic of the experimental set-up. (a) Overview. (b) The schematic diagram of the measurement principle.
where
$I$
is the sum of the four intensities, i.e.
$I_0, I_{45}, I_{90}$
and
$I_{135}$
. The image resolution for the phase retardation and azimuth images was 512
$\times$
512 pixels. The camera captured images at 2000 f.p.s. and 10–20
$\mathrm{\unicode{x03BC}m\,pix^{-1}}$
. A time-averaged values over 10 frames was applied to all experimental data.
We used a CNC (Cellulose Lab Ltd.) suspension. A 0.5 wt
$\%$
of CNC and 0.4 ppm of Triton X-100 were mixed with ultrapure water using a magnetic stirrer and then sonicated for 15 min using a homogeniser (UX-300, Mitsui Electric Co., Ltd). This treatment increases the transparency of a suspension, and the suspension behaves as a Newtonian fluid. This rheological uniformity enables direct comparison with Newtonian Direct Numerical Simulation results. The viscosity of the CNC suspension
$\mu$
was
$1.7+0.1\;\mathrm{mPa}\boldsymbol{\cdot }\mathrm{s}$
at
$10\, ^\circ \mathrm{C}$
measured using a rheometer (MCR302, Anton Paar Co. Ltd.). The 0.5 wt % CNC suspension is essentially Newtonian under the present conditions: at 10
$^\circ$
C the viscosity is
$\mu = 1.7 + 0.1$
mPa s over the shear–rate range
$20 \leqslant \dot {\gamma } \leqslant 1000\,\mathrm{s}^{-1}$
. The highest shear rate
$\dot {\gamma }\sim 6U/D$
estimated from the bubble diameter and velocity observed in the experiment is
$O(10^2)$
, which is within the range of viscosity measurement, further supporting the validity of the Newtonian approximation adopted in this study. The suspension exhibits moderate light attenuation, with a noticeable decrease in optical transmittance for a single-pass path length exceeding approximately 0.1 m. Therefore, the present implementation is best suited for vessels with optical path lengths of
$O(0.1)$
m or less. The surface tension of the CNC and Triton X-100 mixed suspension was measured using the pendant drop method (Berry et al. Reference Berry, Neeson, Dagastine, Chan and Tabor2015). The surface tension 5 s after droplet formation was found to be
$74\pm 4$
mN m–1, which agrees, within experimental uncertainty, with the commonly reported value for pure water at
$10^\circ \mathrm{C}$
(74.23
$\pm$
0.37 mN m−1; Vargaftik, Volkov & Voljak Reference Vargaftik, Volkov and Voljak1983).
The CNC suspension was filled into a vessel (90 mm
$\times$
90 mm
$\times$
450 mm, glass), in which the effects of the wall on the flow structure around an object rising/falling at the centre of the vessel could be neglected (Clift et al. Reference Clift, Grace and Weber1978). A single bubble was generated from the bottom of the vessel using a sound speaker (FF125K, Fostex) (Kusuno, Yamamoto & Sanada Reference Kusuno, Yamamoto and Sanada2019). The sound speaker was activated using an input signal via a function generator (WF1974, NF Corp.) and an audio amplifier (AP20d, Fostex). The static pressure at the outlet of the bubble generator was controlled using a pressure controller (PACE 5000, Baker Hughes). The measurement height varied between 5 and 400 mm from the outlet of the bubble generator. We observed a bubble rising in a rectilinear path within the measurement range, indicating that the flow around the bubble was axisymmetric. The Reynolds number was within the range of
$20 \lt {\textit{Re}} \lt 220$
.
To verify the performance of the polarisation measurement, we used a solid sphere (
$1770\mathrm{\:kg\, m^{-3}, 6 \:mm \:in\:diameter, Plastic}$
). The sphere was quietly dropped from the surface of the suspension. We confirmed that no gas was adsorbed at the solid–liquid interface. The camera was positioned in an area with a sphere falling distance of 120 mm, which was sufficiently long to assume a steady-state condition to be reached at
${\textit{Re}} = 120$
.
2.1.2. Integrated photoelasticity
Considering a two-dimensional stress field of infinitesimal thickness
${\rm d}z$
around an object, the phase retardation
$\varDelta$
, azimuth
$\phi$
and stress
$\sigma$
are related by
\begin{equation} \varDelta \cos 2\phi = C(\sigma _{\textit{xx}}-\sigma _{\textit{yy}}){\rm d}z,\\[-13pt] \nonumber\end{equation}
\begin{equation}\varDelta \sin 2\phi = 2C\sigma _{\textit{xy}}{\rm d}z, \end{equation}
where
$\sigma _{\textit{xx}}$
,
$\sigma _{\textit{yy}}$
and
$\sigma _{\textit{xy}}$
are the stress components in the Cartesian coordinate system shown in figure 1(b), and
$C$
is the material-specific constant. In Mohr’s circle terms,
$\varDelta$
is the circle diameter, while
$\phi$
gives the angular position of that point. Because the isotropic pressure contribution cancels,
$\varDelta$
and
$\phi$
isolate the viscous stresses. Here,
$C$
was determined by fitting the numerical birefringence fields around a sphere to the experimentally measured retardation profiles (see in § 3.1.1). Based on the proportionality between the phase retardation and the principal stress difference, a value of
$C=1.65\times 10^{-5}\mathrm{\:Pa^{-1}}$
provided the best agreement. This value is comparable to previously reported values, such as
$C=1.5\times 10^{-5}\mathrm{\:Pa^{-1}}$
in Nakamine et al. (Reference Nakamine, Yokoyama, Worby, Muto and Tagawa2024). Given that the rheological and optical properties of CNC suspensions vary depending on preparation conditions, we consider this agreement satisfactory and within an acceptable range for quantitative birefringence analysis. Because we consider a two-dimensional stress field, the stress component in the optical axis direction, i.e.
$z$
-direction, does not affect the photoelastic measurement results. However, our experimental results for the phase retardation and azimuth were three-dimensional stress fields in which the fluid stress field was integrated along the optical path. Here, we adopted the method of integrated photoelasticity (Ramesh Reference Ramesh2000), where the stress field is sliced into a collection of virtual thin plates, which are sufficiently thin to be assumed as a two-dimensional stress fields (figure 1
b). In this case, the polarisation state through the three-dimensional stress field can be calculated by multiplying the Stokes parameter
$\boldsymbol{S}$
by the Mueller matrix
$\boldsymbol{X^{(i)}}$
corresponding to each plate
\begin{equation} \boldsymbol{S'}=\boldsymbol{X^{(A_\theta )}X^{(N)} \ldots X^{(i)} \ldots . X^{(2)}X^{(1)}X^{(Q_{45})}X^{(P_0)}S}, \end{equation}
where
$\boldsymbol{X^{(i)}}$
is the number of the plate from which the fluid stress field is sliced (including the linear polariser
$\boldsymbol{X^{(P_0)}}$
, quarter-wave plate
$\boldsymbol{X^{(Q_{45})}}$
and oriented analyser
$\boldsymbol{X^{(A_\theta )}}$
) and is expressed in terms of phase retardation and azimuth (Nakamine et al. Reference Nakamine, Yokoyama, Worby, Muto and Tagawa2024)
\begin{equation} \boldsymbol{X^{(i)}}= \left [ \begin{array}{ccccc} 1 & 0 & 0 & 0 \\ 0 & 1-(1-\cos \varDelta ^{(i)})\sin ^2 2\phi ^{(i)} & (1-\cos \varDelta ^{(i)})\sin 2\phi ^{(i)} \cos 2\phi ^{(i)} & -\sin \varDelta ^{(i)}\sin 2\phi ^{(i)} \\0 & (1-\cos \varDelta ^{(i)})\sin 2\phi ^{(i)} \cos 2\phi ^{(i)} & 1-(1-\cos \varDelta ^{(i)})\cos ^2 2\phi ^{(i)} & \sin \varDelta ^{(i)}\cos 2\phi ^{(i)} \\0 & \sin \varDelta ^{(i)} \sin 2\phi ^{(i)} & -\sin \varDelta ^{(i)}\cos 2\phi ^{(i)} & \cos \varDelta ^{(i)} \\ \end{array} \right ] . \end{equation}
The Mueller matrices
$\boldsymbol{X^{(i)}}\ldots \: \boldsymbol{X^{(N)}}$
are calculated using (2.3), (2.4) and (2.6), based on stress field obtained from the numerical results explained in § 2.2. Because the numerical results are velocity–pressure fields, viscous stresses are obtained from the velocity fields and then transformed into integrated phase retardation and azimuth using the Stokes parameters and Mueller matrices. Finally, using (2.5), we can obtain the integrated phase retardation and azimuth, which can be directly compared with experimental results (for more details, see Yokoyama et al. Reference Yokoyama, Mitchell, Nassiri, Kinsey, Korkolis and Tagawa2023).
2.1.3. Surfactant effects
In aqueous surfactant solutions, the drag coefficient of the bubble changes from that with a free-slip condition to that of a no-slip condition (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997). At high Péclet numbers and concentrations below the critical micelle concentration (CMC), surfactants on the surface are transported towards the rear of the rising bubble, resulting in a concentration gradient along the surface. Consequently, the surface tension gradient induces the shear stress.
In this study, Triton X-100 was used, and we set the bulk concentration
$C$
to 0.4 ppm (7
$\times 10^{-4}$
mol
$\, \rm{ m^{-3}}$
). Because the concentration in this study was much smaller than the CMC of Triton X-100 (0.23 mol
$\mathrm{m^{-3}}$
), micelles did not form. Furthermore, when the Péclet number is of the order of
$10^3$
or higher, advection occurs on a much shorter time scale than diffusion. Under such conditions, the sudden change in surfactant concentration on the surface, or stress jump, is determined by the advection of surfactant and the kinetic balance of adsorption and desorption.
Adsorption/desorption models of interfaces have been proposed previously (Levich Reference Levich1962; Lin, McKeigue & Maldarelli Reference Lin, McKeigue and Maldarelli1990). Following Langmuir adsorption kinetics, the surfactant exchange
$j_{s}$
can be expressed by the following equation:
\begin{equation} j_{s}=k_{a}C_{\!s}(\varGamma _\infty -\varGamma )-k_{d}\varGamma , \end{equation}
where
$C_{\!s}$
is the sublayer concentration,
$\varGamma _\infty$
is the maximum adsorbed surface concentration,
$\varGamma$
is the surface concentration,
$k_{a}$
is the adsorption constant and
$k_{d}$
is the desorption constant. For Triton X-100,
$\varGamma _\infty$
=
$2.9\times 10^{-6}$
mol
$\mathrm{m^{-2}}$
,
$k_{a}$
=
$50\,\mathrm{m^3\, (mol\,s})^{-1}$
and
$k_{d}$
=
$3.3\times 10^{-2}\,\mathrm{s}^{-1}$
(Lin et al. Reference Lin, McKeigue and Maldarelli1990). The equilibrium interfacial concentration
$\varGamma _{\textit{max}}$
at
$j_{s} = 0$
is
\begin{equation} \varGamma _{\textit{max}} = \frac{k_{a}C_{\!s}}{k_{a}C_{\!s}+k_{d}}\varGamma _\infty . \end{equation}
We estimated that
$\varGamma _{\textit{max}}/\varGamma _\infty$
is 0.52, it indicates that the interface is not saturated, and Marangoni effects are expected to remain active under the current conditions. For initially clean interfaces
$(\varGamma = 0)$
, the equilibrium surface concentration can be scaled using the following equation:
\begin{equation} \varGamma _{\textit{max}} = \int _0^{t_{eq}}j_{s}(t)dt\sim k_{a}C_{\!s}\varGamma _\infty t_{eq}, \end{equation}
where
$t_{eq}$
is the time required to reach the equilibrium concentration. From (2.8) and (2.9), the order of
$t_{eq}$
is
$(k_{a}C_{\!s}+k_{d})^{-1}$
. For a bubble rise velocity of 260
$\mathrm{mm\, s^{-1}}$
with a radius of 0.73
$\mathrm{mm}$
, the distance required to reach equilibrium is
$z_{eq} \sim t_{eq}U \sim 1\mathrm{\:m}$
. Therefore, we can infer that the velocity decreases gradually in the measurement region used in this study. Note that, in this study, the time-averaged value over 10 frames
$(5\mathrm{\:ms})$
was used. Because the distance travelled by the bubble during this period was
$O(10^{-4})\mathrm{\:m}$
, which was sufficiently short compared with the characteristic distance
$O(1)\mathrm{\:m}$
over which interfacial adsorption varies, the effect of averaging was not considered to be significant.
Furthermore, we evaluated whether the present experimental conditions satisfy the criteria for the stagnant cap regime, based on the theoretical framework presented in Palaparthi, Papageorgiou & Maldarelli (Reference Palaparthi, Papageorgiou and Maldarelli2006). In this model, the validity of the stagnant cap regime can be assessed by comparing the time scales of adsorption and diffusion with that of surface convection using the dimensionless numbers
${Bi}(1 + k)$
, where
$Bi = k_{d}R/U$
is the Biot number and
$k=k_{a}C/k_{d}$
is the non-dimensional bulk concentration. Substituting the parameters used in this study, we obtain
$Bi(1 + k) \simeq 2 \times 10^{-4}$
. This value is much smaller than unity, confirming that the surfactant transport is dominated by surface convection rather than adsorption or bulk diffusion. Therefore, the present system is indeed within the stagnant cap regime.
2.2. Numerical methods
A numerical simulation of the flow around the solid sphere was performed using the open source code Basilisk (Popinet Reference Popinet2015) to validate the phase retardation and azimuth obtained from the experimental results, which corresponded with the fluid stress field. A flow around a fixed sphere with uniform inflow was computed. The equations of motion and divergence-free conditions are expressed as
\begin{align}& \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{uu}) =\frac{1}{\rho } \left (-\boldsymbol{\nabla } p +\boldsymbol{\nabla }\boldsymbol{\cdot }2\mu \boldsymbol{D} \right ), \end{align}
\begin{align}&\qquad\qquad\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0 \end{align}
where
$\boldsymbol{D}=(\boldsymbol{\nabla }\boldsymbol{u}+(\boldsymbol{\nabla }\boldsymbol{u})^{\mathrm{T}})/2$
is the strain-rate tensor. The sphere was described using an embedded boundary method. For the solid sphere, the boundary condition was no slip. The grid size near the interface was 1/128 of the diameter, which was sufficient to resolve the boundary layer, of the order of
$O({\textit{Re}}^{-0.5})$
in thickness (e.g. Lamb Reference Lamb1945). The computational domain was 40 times larger than the diameter. The effect of boundary conditions at the wall of the computational domain could be neglected. The length of the standing eddy at
${\textit{Re}}$
= 100 was determined to be 0.82
$d$
with a separation angle of 128
$^\circ$
, which agreed with previous studies (Pruppacher, Le Clair & Hamielec Reference Pruppacher, Le Clair and Hamielec1970; Clift et al. Reference Clift, Grace and Weber1978; Magnaudet, Rivero & Fabre Reference Magnaudet, Rivero and Fabre1995), thus validating the numerical results. In addition, at very low Reynolds number
${\textit{Re}} = 0.1$
, the computed drag coefficient is 249 (with pressure and viscous contributions of 81.0 and 168, respectively), which agrees well with theoretical and numerical results in the Stokes regime. This supports the validity of our numerical results also in the creeping-flow regime.
3. Results and discussion
3.1. Flow around a solid sphere
3.1.1. Measurement validation: stress field around a solid sphere at
${\textit{Re}}=120$
In this subsection, we validate the photoelastic measurement technique by comparing experimental results with established numerical simulations. This comparison was conducted both qualitatively and quantitatively to demonstrate the method’s accuracy and reliability in capturing stress fields around bluff bodies moving through a quiescent liquid. Specifically, we focused on the flow around a solid sphere at
${\textit{Re}} = 120$
, representing an intermediate
${\textit{Re}}$
regime similar to that of bubbles.
Figures 2(a) and 2(b) show the experimental results for phase retardation and azimuth fields around a solid sphere falling in a CNC suspension at
${\textit{Re}}$
= 120. The white region corresponds to the sphere itself. The phase retardation, which is proportional to the principal stress difference, and the azimuth, which reflects the stress component alignment, were observed during the motion of the sphere.
Flow surrounding a solid sphere at
${\textit{Re}}$
= 120. (a) Experimental phase retardation. (b) Experimental azimuth. (c) Numerical phase retardation. (d) Numerical azimuth. (e) Phase retardation near the surface. (f) Azimuth near the surface. The dashed line of (a) is the region of
$r = 1.05R$
to
$1.30R$
.
Figures 2(c) and 2(d) present the corresponding numerical simulations for the same system. A direct comparison revealed reasonable agreement between the experimental and numerical results, indicating the accuracy and reliability of the polarisation technique. This validation highlighted the potential of this method in capturing stress fields quantitatively.
Figures 2(e) and 2(f) depict detailed phase retardation and azimuth profiles in the near-surface region, specifically from
$r = 1.05R$
to
$1.30R$
. Remarkably, the experimental measurements exhibited excellent agreement with numerical results within this range for
$20^\circ \lt \theta \lt 160^\circ$
, confirming the quantitative accuracy of the stress-optic law as applied through polarisation techniques. This agreement also demonstrated the capability of this method to accurately resolve stress fields in the boundary layer, whose thickness was approximately
$\delta /R \sim \sqrt {1/{\textit{Re}}} \sim O(10^{-1})$
. The ability to resolve stress fields in this region supports analysis of flow behaviour around the sphere, including detailed discussions on the boundary layer dynamics and wake structures, as discussed in the following sections.
In areas close to the stagnation points (both front and rear, i.e.
$\theta \sim 0^\circ$
and
$180^\circ$
), experimental azimuth measurements exhibited larger errors. This discrepancy resulted primarily from two factors: (i) the azimuth angle shifted sharply from
$180^\circ$
to
$0^\circ$
across the symmetry axis, where shear stress changed sign, and (ii) the phase retardation in these regions was minimal, resulting in reduced measurement accuracy. In addition, data acquisition techniques, particularly at the interface (
$r \lt 1.05R$
), were critical for further reducing errors by noise and data resolution near the interface owing to pixel shadowing and fluctuations.
Despite the inherent challenges, the polarisation technique demonstrated its capability to capture key features of the stress field. These results established its utility for analysing particle-laden flows, including the flow around spheres and bubbles examined in this study. Furthermore, the findings suggest the potential to extend the method to more complex multiphase and non-Newtonian flow systems, subject to further investigation and refinement.
3.1.2. Physical insights into the stress field around a sphere
In this section, we explore the physical characteristics of the stress field around a solid sphere, based on the experimental validation from the previous section (§ 3.1.1). Figure 3 provides a schematic and three-dimensional illustration of the flow structure around the sphere at
${\textit{Re}} = 120$
, highlighting key features such as the boundary layer, wake and standing eddy formation. To deepen our understanding, figure 4 systematically presents stress fields at various
${\textit{Re}}$
, including two limiting cases: Stokes and potential flows. This comparison enables us to explore the evolution of the stress field across different regimes, from low to intermediate
${\textit{Re}}$
and potential flow limits.
(a) Schematic of the flow field around a solid sphere at
${\textit{Re}} = 120$
. Key features of the flow are labelled: (a) mushroom-shaped region of vanishing vorticity (
$\omega \sim 0$
), (b) region where the normal stress gradient is negligible (
$\partial u_y / \partial x \sim 0$
) and (c) regions in the wake where (c i) the tangential stress gradient at the surface vanishes (
$\partial u_t / \partial n \sim 0$
) and (c ii) the normal stress gradient is nearly zero (
$\partial u_y / \partial y \sim 0$
). (b) Three-dimensional visualisation of the flow structure around the sphere, showing the boundary layer, wake and standing eddy formation. Experimental results related to this flow field are shown in figure 2.
Numerical results of flow around a solid sphere for various Reynolds numbers (
${\textit{Re}} = 0.1, 1, 10, 100$
) and analytical results for two limiting cases (Stokes flow and potential flow). (a) Phase retardation
$\varDelta$
and (b) azimuth
$\phi$
distribution. These results illustrate the evolution of stress fields and their dependence on
${\textit{Re}}$
and flow regimes.
To enhance readability and understanding, we first examine the stress fields across different
${\textit{Re}}$
using figure 4. For the potential solution, the phase retardation decreases with increasing distance from the sphere’s surface and exhibits a slight dependence on the angular coordinate,
$\theta$
. In contrast, the azimuth profile strongly depends on
$\theta$
rather than on the radial distance,
$r$
, aligning with the potential flow velocity distributions,
$u_r = U(1 - R^3 / r^3)\cos \theta$
and
$u_\theta = -U(1 + R^3 / 2r^3)\sin \theta$
. In Stokes flow, where the velocities in the fluid are
$u_r = U(1 + R^3/2r^3- 3R / 2r)\cos \theta$
and
$u_\theta = -U(1 - R^3 / 4r^3 + 3R/4r)\sin \theta$
, the phase retardation also decreases with increasing distance, but the azimuth is influenced by both
$\theta$
and
$r$
, resulting in a profile that clearly deviates from the potential solution. The phase and azimuth profiles in Stokes flow align well with behaviour observed in the low
${\textit{Re}}$
regime (
${\textit{Re}} \leqslant 1$
). At higher
${\textit{Re}}$
(
$10 \leqslant {\textit{Re}}$
), a boundary layer develops around the sphere, and the azimuth profile in the outer flow begins to resemble that of the potential solution. At intermediate
${\textit{Re}}$
(10, 100), the flow structure transitions from a viscous-dominated regime to an inertial-dominated regime, with noticeable deviations from symmetry in the phase retardation fields. These trends provide a framework for understanding the specific characteristics of the flow at
${\textit{Re}} = 120$
, which are discussed in detail in the following paragraphs.
The flow field at
${\textit{Re}} = 120$
is depicted in figure 3, highlighting key regions of interest:
-
(a) the mushroom-shaped phase retardation at the front of the sphere, observed despite not being the area of maximum stress, with negligible vorticity (
$\omega \sim 0$
); -
(b) the boundary layer, extending until the normal stress gradient becomes minimal (
$\partial u_y / \partial x \sim 0$
); and (c i) and (c ii), corresponding to regions in the wake in which tangential and normal stress gradients vanish. In the following, we discuss the distinct characteristics of each region.
In the region shown in figure 3(a), which corresponds to the area in front of the sphere and outside the boundary layer, the flow exhibits characteristics of potential flow, with negligible vorticity (
$\omega \sim 0$
). The azimuth profile shown in figure 2(b–d) aligns with the theoretical potential flow solution, and the phase retardation decreases smoothly with distance from the sphere. However, this alone does not explain the mushroom-shaped phase retardation pattern observed at
${\textit{Re}} = 120$
, which requires consideration of the combined effects of potential flow and the boundary layer.
In the outer region (figure 3
a), the flow follows potential flow behaviour, characterised by negligible vorticity and a stress distribution primarily dependent on the angular coordinate,
$\theta$
, rather than the radial distance,
$r$
. This corresponds with the velocity distribution predicted for potential flow, where phase retardation decreases gradually with increasing distance from the sphere, setting a baseline for the retardation profile in the outer flow region.
Closer to the sphere, within the boundary layer, the no-slip condition at the solid surface causes the velocity to sharply drop to zero, creating steep velocity gradients and significant vorticity. This leads to a pronounced increase in phase retardation near the sphere’s surface. At the boundary layer’s edge, where the transition to potential flow occurs, the velocity gradient reverses, reducing viscous stress and causing the phase retardation to drop abruptly. This sharp transition at the boundary layer edge significantly contributes to the mushroom-shaped retardation pattern.
Thus, the combined effects of smooth potential flow behaviour in the outer region and the steep velocity gradients within the boundary layer result in the formation of the mushroom-shaped phase retardation at
${\textit{Re}} = 120$
. The distinct drop in retardation at the boundary layer edge reveals the interplay between these two flow regimes, highlighting the unique structure of the stress field around the sphere at intermediate
${\textit{Re}}$
.
In the region shown in figure 3(b), the boundary layer forms as a result of the no-slip condition at the sphere’s surface. Within this region, phase retardation increases significantly owing to steep velocity gradients near the surface, whereas the azimuth profile reflects the orientation of the stress components induced by these gradients. A large phase retardation, indicating high stress, is observed from the front stagnation point to the equator of the sphere. In this region, the azimuth remains relatively constant, with values ranging from
$135^\circ$
to
$180^\circ$
on the left side of the sphere and
$0^\circ$
to
$45^\circ$
on the right side. This distinction in the azimuth profile highlights the difference in flow structure between the near-solid region and the surrounding flow, marking the boundary layer as a region dominated by vorticity (or viscous stress).
At the boundary layer’s edge, phase retardation decreases sharply, eventually becoming nearly zero as the transition to potential flow occurs. This transition is accompanied by a reversal in the velocity gradient, which reduces the viscous stress. The azimuth profile further illustrates this transition, with azimuth values of approximately
$135^\circ$
near the solid surface and
$45^\circ$
immediately outside the boundary layer. This trend is particularly pronounced at higher
${\textit{Re}}$
, where the thinner boundary layer sharpens the velocity gradient and stress distribution, emphasising the difference between the near-solid and outer flow regions.
In figure 3(b,c), areas of significant phase retardation do not extend to the rear stagnation point but instead spread downstream, forming the wake. Similarly, regions with characteristic azimuth extend into the wake. Near the rear stagnation point, phase retardation is minimal, and the azimuth is
$180^\circ$
(or
$0^\circ$
). However, as the distance from the rear stagnation point increases (approximately
$0.25d$
from the solid surface), the azimuth transitions to
$90^\circ$
. The region where the azimuth remains
$180^\circ$
represents part of the standing eddy but does not encompass its entire structure.
Overall, the unique interplay between the steep velocity gradients in the boundary layer and the transition to potential flow at the boundary layer edge contributes to the observed stress field characteristics. The phase retardation and azimuth profiles, particularly their transitions, reflect the underlying flow physics, emphasising the role of vorticity and velocity gradient reversals in defining the boundary layer and wake structure.
In the region shown in figures 3(c i) and 3(c ii), the wake exhibits minimal phase retardation, with the azimuth transitioning from
$180^\circ$
near the rear stagnation point to
$90^\circ$
further downstream. These regions highlight the formation of a standing eddy, whose size and intensity are determined by the
${\textit{Re}}$
, as depicted in figure 4. The stress gradients in this wake region reflect the transition from the recirculating flow behind the sphere to the external flow. Note that determining the presence and size of the standing eddy purely from phase retardation measurements is challenging because the retardation values in this region are small, as shown in figures 2 and 4. Instead, the azimuth profile provides more reliable insights into the structure of the standing eddy.
First, consider the separation point on the solid surface (figure 3
c i), where tangential stress becomes approximately zero as the velocity gradient reverses before and after separation. In the presence of a standing eddy, the azimuth increases with
$\theta$
. Next, in the standing eddy region along the symmetry axis (figure 3
c ii), the azimuth transitions from
$180^\circ$
near the rear stagnation point to
$90^\circ$
further downstream. The distance from the solid surface to this azimuth transition point is approximately
$0.24d$
, which does not represent the full length of the standing eddy (approximately
$0.94d$
at
${\textit{Re}} = 120$
). This transition corresponds to the condition
$\partial u_y / \partial y = 0$
in the
$(x, y)$
coordinate system, indicating the location of maximum velocity within the standing eddy. If no standing eddy is present, then
$\partial u_y / \partial y \lt 0$
behind the object, suggesting that the azimuth profile can be used to infer the presence or absence of a standing eddy.
These results suggest that the narrow region between the separation point (c i) and the maximum velocity point in the standing eddy (c ii) is critical for understanding the stress profile in the wake. This region provides valuable insights into the interplay between recirculating flow and the transition to the external flow.
As demonstrated above, analysing the phase retardation and azimuth provides valuable insights into distinguishing critical flow regions, such as the boundary layer and wake, and understanding their underlying structures. These measurements enable us to link stress distributions with the flow dynamics, revealing the intricate interplay between viscous and inertial effects. While pressure visualisation remains a challenge for future studies, the square of the phase retardation correlates with viscous energy dissipation, offering a potential pathway to relate stress field measurements to pressure drop and overall flow behaviour.
Drag coefficient of the bubble for various rising heights. The colours show the bubble sizes in each experiment. Closed symbols represent experiments conducted in the presence of surfactants, while open symbols indicate those performed in surfactant-free solutions. (i)–(iv) Correspond to figure 6.
3.2. Flow around a contaminated bubble
3.2.1. Surfactant-induced stress variation around a contaminated bubble
We observed an evolution in the drag coefficient of bubbles rising in a surfactant solution (ultrapure water + CNC 0.5 wt
$\%$
+ Triton X-100 0.4 ppm). Figure 5 shows the drag coefficient ratio for various heights (
$h=$
5, 20, 50, 100, 200, 300 and 400 mm), with the bubble radius in the range
$0.35 \lt R \lt 0.8$
mm. The drag coefficient ratio
$C_{D}^*$
is expressed as
\begin{equation} C_{D}^* = \frac{C_{D}-C_{D_{b}}}{C_{D_{s}}-C_{D_{b}}}, \end{equation}
where
$C_{D_{b}} = 16/{\textit{Re}}(16+3.315{\textit{Re}}^{0.5}+3{\textit{Re}})/(16+3.315{\textit{Re}}^{0.5}+{\textit{Re}})$
is the drag coefficient of the spherical clean bubble (Mei, Klausner & Lawrence Reference Mei, Klausner and Lawrence1994), and
$C_{D_{s}} = 24/{\textit{Re}}(1+0.27{\textit{Re}})^{0.43} + 0.47[1-\exp {(-0.04{\textit{Re}}^{0.38})}]$
is the drag coefficient of a solid sphere (Cheng Reference Cheng2009). The reason for selecting spherical bubbles was that the bubble aspect ratio
$\chi$
is less than 1.25 and almost less than 1.05 (see Appendix B for shape of gradually contaminated bubble). In fact, the Eötvös number
$Eo=\rho g d^2/\sigma$
in the present experiments was less than 0.32, indicating that deformation was sufficiently small. This, combined with the fact that the Reynolds number range of this study was in the intermediate regime, allowed for an accurate evaluation of the drag coefficient. As various drag coefficient models are available for solid spheres (Goossens Reference Goossens2019), we had to choose our model carefully. In the low
${\textit{Re}}$
region,
${\textit{Re}} \sim 20$
for the smallest bubble in this study, a discrepancy of up to 6
$\%$
was observed between the models. In this study, the model proposed by Cheng (Reference Cheng2009) was used as the drag coefficient model closest to the previous numerical results (Magnaudet et al. Reference Magnaudet, Rivero and Fabre1995) for
${\textit{Re}} = 20$
. The drag coefficient in the experiment was estimated based on the force balance described by the equation of motion (Magnaudet & Eames Reference Magnaudet and Eames2000)
\begin{equation} \rho _{\textit{GV}}\frac{{\rm d}\boldsymbol{U}}{{\rm d}t} = \boldsymbol{F}_{\!b} + \boldsymbol{F}_{\!d} + \boldsymbol{F}_{\!a}, \end{equation}
where
$\boldsymbol{F}_{\!b} = -\rho V \boldsymbol{g}$
is the buoyancy force,
$\boldsymbol{F}_{\!d} = - ({1}/{2}) \pi \rho C_{D} R^2 \boldsymbol{U} |\boldsymbol{U}|$
is the drag force and
$\boldsymbol{F}_{\!a} = -C_M \rho V {\rm d}\boldsymbol{U}/{\rm d}t$
is the added mass force. Since the mass of the bubble is negligible compared with that of the surrounding liquid, the added mass term effectively represents the inertial response of the system. When accounting for bubble deformation, the added mass coefficient
$C_M$
was determined using a linear interpolation model
$C_M = 0.62 \chi - 0.12$
(Klaseboer et al. Reference Klaseboer, Chevaillier, Maté, Masbernat and Gourdon2001). Note that the size of the bubble was controlled in the experimental set-up, but the camera was at fixed position and thus all results were for different individual bubbles. In the early stage of bubble displacement (
$h=5$
mm), the drag coefficient of the bubbles in the CNC suspension was close to that of clean bubbles (
$C_{D}^*= 0$
), but it deviated slightly. For bubbles that travelled a longer distance, the drag force increased around a height of
$h\sim$
300 mm, and the drag coefficient for all bubbles was closer to that of a solid sphere (
$C_{D}^* = 1$
). We argue that the bubbles became progressively more contaminated and the transition of the boundary conditions occurred owing to the Marangoni effect. The drag coefficient of some contaminated bubbles was even higher than that of solids. This overshoot has been also observed in numerical results (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997; Kentheswaran et al. Reference Kentheswaran, Antier, Dietrich and Lalanne2023). We emphasise that the primary cause of the observed drag increase is the presence of Triton X-100. Although it has been reported that certain tracer materials can exhibit weak surface activity and reduce bubble velocity (e.g. Weiner et al. Reference Weiner, Timmermann, Pesci, Grewe, Hoffmann, Schlüter and Bothe2019), such effects are generally limited in strength and occur only under specific conditions. In contrast, Triton X-100 is a well-known surfactant with established surface activity and a slow desorption dynamics, which lead to progressive interfacial immobilisation as the bubble rises. In our experiments, we confirmed that Triton X-100 substantially reduces interfacial mobility after a rise distance of approximately 300 mm. By comparison, the effect of CNC particles on the interface is significantly weaker. While CNC can adsorb at the interface to a certain extent, its surface activity is lower and their influence on interfacial mobility develops more gradually. This limited interfacial activity is consistent with the surface tension measurements presented in Appendix C, which show negligible differences between Triton X-100 solutions with and without CNC. This distinction is quantitatively supported by our measurements: the dimensionless drag coefficient
$C_{D}^*$
reaches values of
$0.2{-}0.4$
at a rise height of
$h=50$
mm in Triton X-100 solutions, whereas similar values are observed only at
$h=300$
mm in CNC suspensions. This sixfold difference in rise distance indicates that the interfacial immobilisation caused by CNC is substantially less pronounced under the same conditions. We acknowledge that CNC particles may introduce minor quantitative changes in parameters such as the drag coefficient or surface stress. However, the essential features reported in this study – particularly the abrupt transition in interfacial stress – remain clearly observable across all tested conditions. Therefore, we conclude that the central findings of this work are robust and not significantly influenced by the limited interfacial effects associated with CNC particles.
Experimental conditions for the gradually contaminated bubble case. (i)–(iv) Correspond to figure 6.
Retardation (a) and azimuth (b) field around the bubble, whose radius was 0.59
$\pm$
0.01 mm, at (i)
$C_{D}$
= 0.39 (
$h=5$
mm), (ii)
$C_{D}$
= 0.69 (
$h=50$
mm), (iii)
$C_{D}$
= 1.18 (
$h=100$
mm), (iv)
$C_{D}$
= 1.22 (
$h = 400$
mm). The position at which retardation jump occurred is indicated by a dashed circle. The value of
$C_{D}$
and
$h$
are also indicated in figure 5.
Let us examine the flow structure around the bubble at every height (
$h=5,50,100$
and
$400$
mm). Figure 6(i–iv) show the polarisation measurement results around the bubble with a diameter of 0.59
$\pm$
0.01 mm at heights of
$h$
= 5, 50, 100 and 400 mm, respectively, where the results are time averaged (5 ms, 10 frames). The conditions corresponding to each panel (i)–(iv) are summarised in table 1. Note that the averaging time was much shorter than the transient accumulation of the surfactant and much shorter than the bubble velocity change, as discussed in § 2.1.3. The phase retardation and azimuth near the bubble (
$r$
= 1.05
$R$
) are shown in figures 7(a) and 7(b).
(a) Phase retardation near the bubble’s surface (
$r = 1.05R$
). (b) Azimuth. The results for
$C_{D}$
= 0.39, 0.69, 1.18 and 1.22 correspond to those for figure 6(i–iv), respectively. The dashed circles indicate retardation jumps.
A bubble at
$h$
= 5 mm (figure 6
i) behaved similarly to a clean bubble. First, for a nearly clean bubble (for more detail see figure 7
a), the largest phase retardation was observed at the front stagnant point in the entire field. This indicated a completely different flow structure from the solid case, which had the largest phase retardation around the equator. The phase retardation was very small at the side of the bubble. We observed that this was because the boundary condition was a free-slip condition. Slightly further from the rear stagnant point (
$\theta \sim$
150
$^\circ$
), a large phase retardation was observed. This region was certainly a wake region, but the phase retardation was larger than in the solid-sphere case. Figure 6(i) shows that little azimuth change occurred in the radial direction, which was not observed in the solid sphere in figure 2. In addition, the azimuth around the bubble exhibited near fore–aft symmetry, similar to a potential flow, indicating that the vorticity generated at the surface was small. This indicated that the bubble’s surface largely maintained a free-slip condition, which resulted in a small deviation from the potential flow owing to insufficient vorticity generation. However, the azimuth deviated from the potential in the vicinity of the rear part indicated by the blue dashed circle. Therefore, the reason for the deviation from the potential flow was considered to be a transition from free slip to no slip in the flow structure owing to contamination in the vicinity of the rear part. This could be observed more clearly as the rise height increased with accumulating surfactants on the bubble’s surface.
The bubble at
$h$
= 50 mm (figure 6
ii) was expected to be at intermediate contamination based on its drag coefficient. Because its rise velocity was smaller than that of the clean bubble, the magnitude of the phase retardation around the bubble at
$h$
= 50 mm was smaller than that of the clean bubble because the phase retardation was proportional to the velocity in a similar flow, as can be observed from (2.3) and (2.4). Nevertheless, in the front part of the bubble, both phase retardation and azimuth exhibited the same trend as that of a clean bubble. Note that, in the region
$\theta \sim 105^\circ$
, a jump occurred in the phase retardation and azimuth near the bubble’s surface, indicating that the boundary conditions have transitioned. The jump in
$\phi$
signifies that, in the Cartesian frame, the ratio between the in-plane normal stress (
$\sigma _{\textit{xx}} - \sigma _{\textit{yy}}$
) and the shear stress
$\sigma _{\textit{xy}}$
changes abruptly. Assuming that the flow upstream of the bubble resembles a clean interface, where the normal component dominates, the sudden appearance of the jump implies an increase in the shear component at this polar angle, i.e. the interfacial condition locally shifts toward a no-slip condition. The phase retardation and azimuth in the region below the jump were similar to those of the solid sphere, except the non-zero phase retardation at the rear stagnant point. At the rear of the bubble, we estimated that detachment occurred because
$\phi$
increased with the increase in
$\theta$
.
The bubble at
$h$
= 100 mm (figure 6
iii) had a drag coefficient similar to a solid sphere and could be assumed to be a ‘fully contaminated’ bubble from the drag coefficient perspective. The magnitude of the phase retardation was even smaller than for
$h$
= 50 mm owing to the velocity reduction. Now, in the region
$\theta \sim 70^\circ$
, the phase retardation and azimuth changed drastically, i.e. a jump occurred near the bubble’s surface. The boundary condition was the no-slip condition in a larger area than the bubble at
$h$
= 50 mm because of the accumulated contamination, although the bubble was not a fully contaminated one from the flow structure perspective. When the contamination advanced to this stage, the boundary condition for
$0^\circ \leqslant \theta \lesssim 70^\circ$
was similar to that of the clean bubble in the front, and that for
$70^\circ \lesssim \theta \leqslant 180^\circ$
was close to that of the solid sphere in the rear. Unlike the
$h$
= 50 mm bubble, the phase retardation at the rear stagnation point was zero. The fact that the drag coefficient agreed with that of the solid sphere although the entire flow field did not match that of the solid sphere can be attributed to the manner in which normal and shear stress contributions offset each other. Specifically, although the local normal stress around
$\theta \sim 70^\circ$
differed from that of a fully no-slip sphere, the front region (
$0^\circ \leqslant \theta \lesssim 70^\circ$
) retained a free-slip boundary condition, which balanced out this difference such that the overall normal force remained comparable to that of the solid sphere. For the shear stress, its contribution to drag was weighted by
$\sin \theta$
. Because
$\sin \theta$
was small in the front part, the difference between free-slip and no-slip conditions had only a minor effect on the total shear contribution. Consequently, although the local flow field deviated from the fully no-slip cases, the integrated normal and shear stresses produced nearly the same drag coefficient as that of a solid sphere, as indicated by previous numerical studies (e.g. Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997).
The bubble at
$h$
= 400 mm (figure 6
iv) was the most contaminated bubble in this study. Compared with the bubble at
$h$
= 100 mm, the boundary condition transition occurred more forward (
$\theta \sim 40^\circ$
). The evolution of the angle where the jump in the boundary condition occurred 30
$^\circ$
ahead of the bubble of
$h$
= 100 mm, although the rising trajectory was four times longer. This indicated that a very long time was required for the front part to become a no-slip condition.
The polarisation measurement experiments described in this section showed that, as the rising distance of the bubbles increased, the transition from the free-slip to no-slip condition (flow around a clean bubble to flow around a solid sphere) occurred in sequence from the rear of the bubbles. The point at which the boundary conditions transitioned appeared as a jump in the phase retardation and azimuth.
At the point where the boundary condition shifts from free slip to no slip, a very strong vorticity (shear stress) must be generated at the surface to satisfy the continuity equation (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997).
(a) Normalised viscous shear stress and (b) normalised normal stress near the bubble’s surface (
$r = 1.05\,R$
). Stresses are normalised by
$\mu U/2 R$
. For conditions see figure 6.
Here, we reconstructed the normalised viscous shear stress
$\sigma _{nt}$
and normal stress
$\sigma _{nn}$
near the surface from the retardation and azimuth results (see Appendix A for axisymmetric reconstruction method). Note that the stress components are modified to the (
$r$
,
$\theta$
) coordinates along the surface from (
$x$
,
$y$
) coordinate system. Figure 8(a) shows that a local large shear stress acted on the point where the stress jumped at each height. In the front,
$\theta \sim 0^\circ$
, artificial stress oscillation occurred owing to reconstruction using the fluctuating azimuth. In the front part, excepting the artificial fluctuation region, a similar flow field was established irrespective of the bubble height. Figure 8(b) shows that the normal stress component also exhibited the spikes. This was due to the no-slip behaviour that caused the velocity to suddenly drop to zero. In the free-slip condition, the shear stress was zero and the vorticity was
$O(1)$
at the surface; therefore, such a large shear stress near the surface did not occur in a clean bubble. Therefore, this revealed that the shear stress (or the vorticity proportional to it) was caused by the Marangoni stress (Atasi et al. Reference Atasi, Ravisankar, Legendre and Zenit2023). This local shear stress indicated the localised generation of a surface tension gradient (or surfactant concentration gradient), and we inferred that a small amount of surfactant was present in the forward part from the jump point. Nevertheless, a small shear stress was still acting on the front part of the bubble, and the possibility of the effect of surfactant adsorption during advection could not be ruled out.
3.2.2. Measured cap angle comparison with the stagnant cap model and numerical simulations
In this section, we compare the measured cap angle with both the stagnant cap model (Sadhal & Johnson Reference Sadhal and Johnson1983) and the numerical results from Cuenot et al. (Reference Cuenot, Magnaudet and Spennato1997). The stagnant cap model, which was often used for comparative verification of previous numerical analyses, has not been compared with experimental results. Using the stagnant cap model, we can describe
$C_{D}^*$
as a function of
$\theta _C$
(see (53) of Sadhal & Johnson Reference Sadhal and Johnson1983)
\begin{equation} C_{D}^*(\theta _C) = \frac{1}{2\pi }\left [ 2(\pi -\theta _C) + \sin \theta _C+\sin 2\theta _C-\frac{1}{3}\sin 3\theta _C \right ]. \end{equation}
Figure 9 shows the relationship between the normalised drag coefficient and cap angle. The experimental results for various bubble sizes exhibited reasonable agreement with the model overall. This result strongly confirmed that the phase retardation and azimuthal jumps in the experiment corresponded to the shift of the boundary condition from free slip to no slip. However, two discrepancies were observed. The first was that the experimental drag at
$140^\circ \lt \theta _{\mathrm{C}}$
was larger than the theory. This was because the bubble had a spheroidal shape and the drag was larger than that of spherical shape (Blanco & Magnaudet Reference Blanco and Magnaudet1995; Sanada et al. Reference Sanada, Sugihara, Shirota and Watanabe2008). The second was that the slope of the experimental drag was slightly steeper than the theory near the equator (
$\theta _C \sim 90^\circ$
). A similar trend was observed in the numerical analyses performed at a finite
${\textit{Re}}$
of 100 (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997; Dani et al. Reference Dani, Cockx, Legendre and Guiraud2022), indicating that the drag transition is more sensitive than for the creeping flow. To the best of the authors’ knowledge, this is the first experimental result that is consistent with previous numerical studies of the flow structure around bubbles in surfactant solutions with finite
${\textit{Re}}$
.
Comparison of drag coefficient ratios of experimental results and stagnant cap model: —
$\!$
—, stagnant cap model (Sadhal & Johnson Reference Sadhal and Johnson1983); •, experimental results (cap angle is obtained from phase retardation jump);
$*$
, numerical results (Cuenot et al. Reference Cuenot, Magnaudet and Spennato1997). (i)–(iv) Correspond to figure 6.
Note that, although the drag coefficient of an apparently contaminated bubble approaches that of a solid sphere, the flow field differs, particularly for
$\theta _C \lt 60^\circ$
. Drag contributions primarily result from shear stress in the range
$45^\circ \lt \theta _C \lt 135^\circ$
and normal stress elsewhere. As shown in figure 8, viscous normal stress may counterbalance shear stress and the pressure gradient, even when the bubble front remains clean. While this flow field difference likely has a slight impact on the drag coefficient, it can influence other phenomena, such as bubble path instability (Tagawa et al. Reference Tagawa, Takagi and Matsumoto2014; Pesci et al. Reference Pesci, Weiner, Marschall and Bothe2018; Farsoiya et al. Reference Farsoiya, Popinet, Stone and Deike2024), lift reversal (Fukuta, Takagi & Matsumoto Reference Fukuta, Takagi and Matsumoto2008; Hayashi & Tomiyama Reference Hayashi and Tomiyama2018) and stability of multi-bubbles (Harper Reference Harper2008; Atasi et al. Reference Atasi, Ravisankar, Legendre and Zenit2023).
3.2.3. Maximum shear stress in the vicinity of the bubble
In this section, we examine how the drag coefficient
$C_{D}$
correlates with the maximum vorticity near the bubble’s surface. As discussed in recent numerical investigations (Atasi et al. Reference Atasi, Ravisankar, Legendre and Zenit2023; Hayashi et al. Reference Hayashi, Motoki, Legendre and Tomiyama2025), the maximum vorticity generated at the interface of contaminated bubbles is an important factor that determines the bubble dynamics such as drag and lift. For a clean spherical bubble, Legendre (Reference Legendre2007) showed that the drag coefficient scales with the maximum vorticity
$\omega _{\textit{max}}^b$
as
$C_{D}{\textit{Re}}=16\omega _{\textit{max}}^b$
. This gives
$C_{D}{\textit{Re}}=16$
for Hadamard–Rybczynski solution
$(\omega _{\textit{max}}^b=1)$
and
$C_{D}{\textit{Re}} = 48$
for the potential solution
$(\omega _{\textit{max}}^b=3)$
. For a solid sphere, Hayashi et al. (Reference Hayashi, Motoki, Legendre and Tomiyama2025) reported a linear relationship
$C_{D}{\textit{Re}}=11\omega _{\textit{max}}^s+7.5$
, where
$\omega _{\textit{max}}^s=3/2(1+16/11{\textit{Re}}^{0.687})$
. In the case of contaminated bubble, such relationships hold when diffusion dominates. However, under inertia-dominant conditions, a consistent scaling of
$C_{D}{\textit{Re}}$
with maximum vorticity has not yet been established. These existing relations are based on numerical results; in this section, we re-examine them experimentally.
Can vorticity nevertheless be estimated from polarisation measurements? The shear stress in the vicinity of the spherical bubble is given by
$\mu (\partial u_t/\partial r-u_t/R)$
, whereas the vorticity is
$(\partial u_t/\partial r+u_t/R)$
. For a free-slip condition, the shear stress is zero and we obtain
$\omega =2u_t/R$
, but
$u_t$
is unknown in our experiments. In contrast, for a no-slip condition, the vorticity is
$\omega =\partial u_t/\partial r$
which is directly proportional fo the shear stress. Assuming that the dominant contribution to the local vorticity arises from the tangential velocity gradient near the no-slip region, the maximum normalised shear stress obtained via birefringence can be considered an approximate indicator of the local vorticity. However, since the measurement is taken in the vicinity of the bubble
$(r = 1.05R)$
, and not on the surface, the result is still an approximation and requires careful interpretation.
Maximum shear stress in the vicinity of the contaminated bubble (
$r = 1.05R$
) and
$C_{D}{\textit{Re}}$
;
$\bigcirc$
, experimental results; - -, maximum vorticity for the solid sphere proposed by Hayashi et al. (Reference Hayashi, Motoki, Legendre and Tomiyama2025).
The product
$C_{D} {\textit{Re}}$
was plotted against the maximum normalised shear stress
$-\sigma _{nt}^{{\textit{max}}}$
, as shown in figure 10. It shows that higher shear stress corresponds to the vorticity correlating with higher drag, consistent with numerical predictions. However, the experimental vorticity tends to lie above the linear predicted trend
$(C_{D} {\textit{Re}} - 7.5)/11$
for a solid sphere. This deviation supports the numerical findings (Hayashi et al. Reference Hayashi, Motoki, Legendre and Tomiyama2025), and no generalised relationship between maximum vorticity and drag coefficient is identified in the present experiments.
The discrepancy from previously established
$C_{D}{\textit{Re}}-\omega$
relationships for clean bubbles and a solid sphere is likely attributable to localised spikes in the shear stress profile of a partially contaminated bubble. Interestingly, as the contamination level increases and the boundary condition approaches that of a solid sphere
$(C_{D}^*\rightarrow 1)$
, the maximum shear stress tends to align with the solid-sphere prediction. In contrast, bubble with lower
$C_{D}^*$
often exhibit much higher maximum shear stress. Returning to figure 8, at
$h=400$
mm, the tangential shear stress profile becomes smooth, where
$\theta _C$
is approximately
$40^\circ$
, and no sharp spike is observed. The corresponding maximum shear stress is comparable to that of a solid sphere. On the other hand, for bubble with lower contamination levels (i.e. smaller
$C_{D}^*$
), sharp spikes appear. In some cases, the maximum shear stress exceeds that of a more contaminated bubble, even though the
$C_{D}^*$
is smaller. These observations suggest that lower contaminated bubble may exhibit stronger vorticity generation than would be expected based on a solid-sphere model. It may be possible to correct the overestimation of local vorticity using the vorticity distribution of solid spheres at the cap angle, but such a correction requires a dedicated flow model around a solid sphere.
4. Conclusion
We experimentally determined the stress field of the flow around a bubble rising in a dilute surfactant solution using a newly developed unsteady polarisation measurement method. We first validated our method by measuring the flow around a solid sphere falling at terminal velocity in a quiescent solution without a surfactant effect. The measurements and numerical results were in quantitative agreement, confirming that the photoelastic method was successful in measuring the stress around the sphere. The stress field around the solid sphere revealed distinct regions, including the outer potential flow region, boundary layer and wake. Interestingly, the mushroom-shaped phase retardation observed at the front of the sphere was attributed to the interplay between potential flow characteristics and the boundary layer dynamics. These features, combined with azimuth profiles reflecting stress orientations, provided a detailed understanding of the stress distribution around the sphere and the mechanisms underlying flow structures at intermediate Reynolds numbers.
We then measured the flow around a bubble rising in a quiescent surfactant solution with a high Péclet number and intermediate Reynolds number. The results showed that the flow structure near the front of the bubble resembled that of a clean bubble, whereas the rear flow was similar to that of a solid sphere. However, the stress field in front of the bubble differed even when the drag coefficient was nearly the same, particularly when the cap angle was less than 60
$^\circ$
. Additionally, when the drag coefficient deviated from that of a solid sphere, slight differences in the rear stress field were observed. Sudden jumps in phase retardation and azimuth, indicated by the cap angle, were detected between the front and rear of the bubble. These jumps evolved with increasing distance travelled, reflecting the time-dependent adsorption of surfactants. By reconstructing the stress field under axisymmetric assumptions, we experimentally demonstrated the presence of localised stresses acting in the phase retardation jump region. The measured cap angle as a function of the normalised drag coefficient had reasonable agreement with the predictions of the stagnant cap model proposed by Sadhal & Johnson (Reference Sadhal and Johnson1983) and quantitatively aligned with numerical results at intermediate Reynolds numbers by Cuenot et al. (Reference Cuenot, Magnaudet and Spennato1997), marking the first experimental confirmation of previous numerical studies on flow structures in the vicinity of a bubble in a dilute surfactant solution at finite Reynolds numbers.
As a future prospect, the method employed in this study may create new opportunities for revealing the stress field in a wide range of dispersion flows. The observed evolution of the cap angle along the rise distance in the present study may offer an opportunity to validate theoretical models of surfactant transport that incorporate adsorption kinetics as the rate-limiting step (He, Maldarelli & Dagan Reference He, Maldarelli and Dagan1991). Although only a simple system was used in this study, various systems exist in the bubble dynamics, such as multi-bubbles and bubble motion in non-Newtonian fluids. For low Péclet number conditions and mixed surfactants, the boundary conditions become half-slip, in which both the shear stress and the tangential velocity are non-zero. In addition, as the amount of vorticity generation increases under slip conditions, unstable behaviour is caused by vorticity evacuation and accumulation generated at the surface. Moreover, lift reversal (Fukuta et al. Reference Fukuta, Takagi and Matsumoto2008; Hayashi & Tomiyama Reference Hayashi and Tomiyama2018) and bubble train stabilisation (Liger-Belair et al. Reference Liger-Belair, Marchal, Robillard, Dambrouck, Maujean, Vignes-Adler and Jeandet2000; Harper Reference Harper2008; Atasi et al. Reference Atasi, Ravisankar, Legendre and Zenit2023) are caused by vorticity stretching/tilting, respectively; both are organised by the maximum vorticity (Magnaudet & Mougin Reference Magnaudet and Mougin2007; Yang & Prosperetti Reference Yang and Prosperetti2007; Cano-Lozano et al. Reference Cano-Lozano, Martínez-Bazán, Magnaudet and Tchoufag2016). Polymer stretching generates a negative wake behind a bubble (Zenit & Feng Reference Zenit and Feng2018). Furthermore, since the present measurement captures the stress field in the bulk fluid near the interface, and not the interfacial stress itself, the technique may remain applicable even in the presence of additional interfacial stress components such as interfacial viscosity (Singh & Narsimhan Reference Singh and Narsimhan2022; Herrada et al. Reference Herrada, Ponce-Torres, Rubio, Eggers and Montanero2022). Therefore, the extension of the results of this study to viscous stress, boundary layer and wake regions will contribute to the elucidation of the detailed physical mechanism of these phenomena, as well as in complex flows.
Acknowledgements
The authors would like to thank Professor R. Zenit for carefully reading the manuscript and for giving constructive comments.
Funding
This work was supported by JSPS KAKENHI Grant Numbers JP22K20403, JP23K13251 and JP24H00289. We also acknowledge support from The Ihara Science Nakano Memorial Foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Axisymmetric reconstruction
In this section, we introduce the axisymmetric reconstruction procedure from the integrated phase retardation and azimuth to normal/shear stress. Axisymmetric reconstruction has been conventionally used for transparent solids (Aben, Ainola & Errapart Reference Aben, Ainola and Errapart2010; Yokoyama et al. Reference Yokoyama, Mitchell, Nassiri, Kinsey, Korkolis and Tagawa2023). Inverse Abel transform and onion-peeling methods have been used for the reconstruction, and we adopted the onion-peeling method based on them. For axisymmetric reconstruction, we must first transform the stresses in the Cartesian coordinate system to the cylindrical coordinate system
$(r,y,\varTheta )$
, assuming a steady-state condition, as follows:
\begin{align}&\qquad\qquad\qquad \sigma _{\textit{xy}} = \sigma _{\textit{ry}}\cos \varTheta , \end{align}
\begin{align}& \sigma _{\textit{xx}}-\sigma _{\textit{yy}} = \sigma _{\textit{rr}}\cos ^2\varTheta + \sigma _{\varTheta \varTheta }\sin ^2\varTheta -\sigma _{\textit{yy}}. \end{align}
From (2.4), with the phase retardation and azimuth obtained from the polarisation measurement, the shear stress is independent of the other stress components, indicating that it can be estimated by a procedure similar to that previously used in the analysis of solids
\begin{equation} \varDelta \sin 2\phi =2C\int \sigma _{\textit{ry}}\cos \varTheta {\rm d}z. \end{equation}
In contrast, normal stress cannot be estimated in the same manner as for solids. In the steady state of a solid, the derivation is based on the equilibrium equation
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } = 0$
, i.e. the shear and normal stresses are balanced, whereas no such limitations exist in a fluid owing to advection
$(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } = \rho \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{(uu)})$
. Therefore, we use the relation between the continuity equation and axisymmetric stresses assuming a Newtonian fluid. We treat only the viscous stresses neglecting pressure because the pressure cancels out in (A2). Assuming a Newtonian fluid, the continuity equation yields
$\sigma _{\textit{rr}}^{vis} + \sigma _{\varTheta \varTheta }^{vis} + \sigma _{\textit{yy}}^{vis} = 2\mu (\partial u_r/\partial r+u_r/r + \partial u_y/\partial y) = 0$
, where the subscript denotes the viscous stress component. Additionally, the axisymmetric stress relation
$\sigma _{\textit{rr}}^{vis} = \partial (r\sigma _{\varTheta \varTheta }^{vis})/\partial r$
, from (2.3) and (A1), the relation between
$\varDelta$
,
$\phi$
and
$\sigma _{\varTheta \varTheta }^{vis}$
becomes
\begin{equation} \varDelta \cos 2\phi = C\int 3\sigma _{\varTheta \varTheta }^{vis} + r(1+\cos ^2\varTheta )\frac{\partial \sigma _{\varTheta \varTheta }^{vis}}{\partial r}{\rm d}z. \end{equation}
When the
$\varTheta$
component is known, the other normal stress components are also uniquely determined. Once the full stress tensor
$\boldsymbol{\sigma }$
is reconstructed, it can be projected onto the local interface direction to obtain interfacial components. Specifically, the normal and tangential surface stress are given by
$\sigma _{nn}=\boldsymbol{n}^{\mathrm{T}}\boldsymbol{\sigma }\boldsymbol{n}$
and
$\sigma _{nt}=\boldsymbol{n}^{\mathrm{T}}\boldsymbol{\sigma }\boldsymbol{t}$
, where
$\boldsymbol{n}$
and
$\boldsymbol{t}$
are the unit normal and tangential vectors to the bubble’s surface, respectively.
Stress near the solid sphere’s surface (
$r = 1.05\,R$
);—
$\!$
—, stress obtained from the velocity;
$\bigcirc$
, stress obtained from the phase retardation and azimuth. For conditions see figure 2.
We performed tests on numerical data (figure 2) to validate our reconstruction method. Figure 11 shows the stress fields obtained from the velocity and phase retardation and azimuth fields using the onion-peeling method. The stresses in the far from the symmetry axis were calculated as 0. The figure shows that the stress fields obtained from the velocity and phase and azimuth were consistent over a wide range. Therefore, we conclude that this reconstruction method is applicable. Note that the values near the axis of symmetry and near the interface cannot be guaranteed owing to the orientation fluctuation as discussed in § 3.1.1.
Appendix B. Aspect ratio of gradually contaminated bubble
The deformation of the bubbles was further analysed by comparing the experimentally measured aspect ratios with theoretical predictions. We employed the correlation proposed by (Kentheswaran et al. Reference Kentheswaran, Antier, Dietrich and Lalanne2023), which accounts for the effect of interfacial contamination via the drag coefficient ratio
$C_{D}^*({\textit{Re}},Eo)$
, in addition to
${\textit{Re}}$
and
$Eo$
\begin{align}&\qquad\,\, \chi ^{cont}=(1-\alpha )\chi ^{\textit{solid}}+\alpha \chi ^{\textit{clean}},\\[-8pt]\nonumber \end{align}
\begin{align}&\,\, \chi ^{\textit{clean}}({\textit{Re}},Eo)=\left(1+0.016Eo^{1.12}{\textit{Re}}\right)^{0.388},\\[-8pt]\nonumber \end{align}
\begin{align}& \chi ^{\textit{solid}}({\textit{Re}},Eo)=\left(1+0.024Eo^{1.17}{\textit{Re}}^{0.44}\right)^{0.57},\\[-8pt]\nonumber \end{align}
\begin{align}&\qquad\qquad \alpha = \left(1-C_{D}^*({\textit{Re}},\chi )^{2.5}\right)^5,\\[-8pt]\nonumber \end{align}
\begin{align}& C_{D}^*({\textit{Re}},\chi )=\frac{C_{D}-C_{D}^{\textit{clean}}({\textit{Re}},\chi )}{C_{D}^{\textit{solid}}({\textit{Re}},\chi )-C_{D}^{\textit{clean}}({\textit{Re}},\chi )},\\[-8pt]\nonumber \end{align}
\begin{align}& C_{D}^{\textit{clean}}({\textit{Re}},\chi )=\frac{16}{{\textit{Re}}}\big(1+0.25{\textit{Re}}^{0.32}\chi ^{1.9}\big),\\[-8pt]\nonumber \end{align}
\begin{align}& C_{D}^{\textit{solid}}({\textit{Re}},\chi )=\frac{24}{{\textit{Re}}}\big(1+0.17{\textit{Re}}^{0.687}\chi ^{0.48}\big).\\[-10pt]\nonumber \end{align}
Here,
$\chi ^{{\textit{clean}}}$
was proposed by Aoyama et al. (Reference Aoyama, Hayashi, Hosokawa and Tomiyama2016),
$\chi ^{{\textit{solid}}}$
was proposed by Aoyama et al. (Reference Aoyama, Hayashi, Hosokawa and Tomiyama2018), empirical correlation for
$C_{D}^{{\textit{clean}}}$
was developed by Chen et al. (Reference Chen, Hayashi, Hosokawa and Tomiyama2019) and that for
$C_{D}^{{\textit{solid}}}$
was given by Kentheswaran et al. (Reference Kentheswaran, Antier, Dietrich and Lalanne2023) as a modified version of the model proposed by Aoyama et al. (Reference Aoyama, Hayashi, Hosokawa and Tomiyama2018). Figure 12 shows the comparison between the predicted aspect ratio
$\chi ^{p}$
and the experimentally measured aspect ratio
$\chi ^{e}$
, categorised into contaminated bubble, solid particle and clean bubble. Although the measured aspect ratios lie within a relatively narrow range (
$\chi ^e \lt 1.25$
), the predicted trends for
$\chi ^{cont}$
are consistent with the experimental data: the deformation of bubbles decreases (i.e. the aspect ratio approaches unity) as the contamination increases. However, in the nearly spherical regime (
$\chi ^e \lt 1.05$
), the model tends to underpredict the aspect ratio, suggesting a limitation in accuracy for weakly deformed bubbles. These results support the notion that surface contamination leads to reduced bubble deformation, while also indicating that further refinement of the model may be necessary in the low-deformation limit.
Comparison between the predicted aspect ratio
$\chi ^p$
and the experimentally measured aspect ratio
$\chi ^e$
. The diagonal line indicates perfect agreement. The predicted values are based on the correlation proposed by Kentheswaran et al. (Reference Kentheswaran, Antier, Dietrich and Lalanne2023), incorporating drag coefficient ratios and empirical aspect ratio models from Aoyama et al. (Reference Aoyama, Hayashi, Hosokawa and Tomiyama2016), Aoyama et al. (Reference Aoyama, Hayashi, Hosokawa and Tomiyama2018) and Chen et al. (Reference Chen, Hayashi, Hosokawa and Tomiyama2019).
As a side note, while Kentheswaran et al. (Reference Kentheswaran, Antier, Dietrich and Lalanne2023) provide a deformation-dependent drag coefficient ratio model using the Chen et al. (Reference Chen, Hayashi, Hosokawa and Tomiyama2019) model for clean bubble, we found that applying their model at an aspect ratio of unity yields a drag coefficient nearly 10
$\%$
different from the value given by Mei et al. (Reference Mei, Klausner and Lawrence1994). This discrepancy likely stems from the intermediate Reynolds number regime in which the model was developed. Consequently, while the model considered aspect ratio does not fully capture
$C_{D}^*\sim 0$
for nearly spherical bubbles (cap angles near
$0^\circ$
), it still offers a more realistic estimate than assuming a perfectly spherical shape. It should also be noted that using a drag coefficient model for solid spheres (Kentheswaran et al. Reference Kentheswaran, Antier, Dietrich and Lalanne2023) that explicitly accounts for aspect ratio tends to overestimate the drag by approximately 5
$\%$
when applied at an aspect ratio of 1, compared with the Cheng (Reference Cheng2009) model. When translated into
$C_{D}^*$
, this leads to a reduction of approximately 0.07, which should be taken into consideration when interpreting the results. In our study, since the observed aspect ratios remain below 1.25, we chose to base our drag analysis on the spherical approximation. This approach ensures consistency with drag correlations in the near-spherical regime, where deviations from unity aspect ratio are minor.
Appendix C. Effect of Triton X-100 concentration on surface tension with and without CNC
This section presents surface tension measurements of aqueous Triton X-100 solutions with and without CNC. Figure 13 shows the surface tension values as a function of Triton X-100 concentration. The measurements were conducted using the pendant drop method (Berry et al. Reference Berry, Neeson, Dagastine, Chan and Tabor2015), and the surface tension was recorded 5 seconds after droplet formation. At low concentrations, both w/ CNC and w/o CNC solutions exhibit surface tension values close to that of pure water (74 mN m−1). This suggests that, at the concentrations used in the main experiments, the contribution of CNC to surface tension is negligible. Furthermore, the surface tension curves of the two solutions nearly overlap over the entire concentration range tested, indicating that CNC does not significantly affect the adsorption behaviour of Triton X-100 under these conditions.
Effect of Triton X-100 concentration on surface tension with and without CNC.
However, at much higher CNC concentrations, CNC may strongly influence the surface tension and modify the Marangoni effect. Indeed, a slight reduction in bubble rise velocity, reminiscent of Marangoni-like effects, was observed. This phenomenon merits further investigation in future studies.
Open access
