2.2.1 Diffusive fluxes

Under the assumption of dilute species concentrations both in the liquid phase and on the interface, the diffusive fluxes are modelled via Fick’s law, i.e.

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}=-D\unicode[STIX]{x1D735}c\quad \text{in }\unicode[STIX]{x1D6FA}^{+}(t), & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}^{\unicode[STIX]{x1D6F4}}=-D^{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}c^{\unicode[STIX]{x1D6F4}}\quad \text{on }\unicode[STIX]{x1D6F4}(t). & \displaystyle\end{eqnarray}$$

Furthermore, homogeneous Neumann conditions (2.12) for the diffusive fluxes at the outer domain boundary are assumed, i.e.

(2.12) $$\begin{eqnarray}\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}=0\quad \text{on }\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}(t).\end{eqnarray}$$

2.2.2 Sorption process

To model the sorption process, two limiting situations can be considered: diffusion-controlled sorption (fast) and kinetically controlled sorption (slow); see Miller et al. (Reference Miller, Fainerman, Pradines, Kovalchuk, Kovalchuk, Aksenenko, Liggieri, Ravera, Loglio, Sharipova, Vysotsky, Vollhardt, Mucic, Wüstneck, Krägel, Javadi and Romsted2014). In the first case, the sorption process is much faster than the diffusive transport, while in the latter case the sorption process is slower than the diffusive transport, typically due to the presence of a kinetic barrier. Thus, the transfer rate $s^{\unicode[STIX]{x1D6F4}}$ will be determined in two different ways, while the transmission condition (2.9) always holds. In both cases, the effect of surfactant on the interfacial surface tension is described by the surface equation of state which in a general form reads

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D70E}_{0}=\unicode[STIX]{x1D6F1}(c^{\unicode[STIX]{x1D6F4}}).\end{eqnarray}$$

The function $\unicode[STIX]{x1D6F1}(c^{\unicode[STIX]{x1D6F4}})$ in (2.13) assumes a specific expression with respect to the sorption model employed; see Pesci et al. (Reference Pesci, Dieter-Kissling, Marschall, Bothe, Rahni, Karbaschi and Miller2015) for more details on the derivation of (2.13) and the full set of sorption models available in our library. For instance, in the Langmuir model, the surface tension equation of state reads

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}+RTc_{\infty }^{\unicode[STIX]{x1D6F4}}\text{ln}\left(1-\frac{c^{\unicode[STIX]{x1D6F4}}}{c_{\infty }^{\unicode[STIX]{x1D6F4}}}\right),\end{eqnarray}$$

where $R$ is the universal gas constant equal to $8.3144~\text{J}~(\text{mol}~\text{K})^{-1}$ , $T$ is the absolute system temperature in Kelvin and $c_{\infty }^{\unicode[STIX]{x1D6F4}}$ is the saturated surfactant concentration, i.e. the maximum number of adsorbed molecules per area. For the application case presented in § 4 it has been proven that a fast model is adequate to describe the sorption mechanism; see Aksenenko et al. (Reference Aksenenko, Makievski, Miller and Fainerman1998), Kovalchuk et al. (Reference Kovalchuk, Krägel, Makievski, Ravera, Liggieri, Loglio, Fainerman and Miller2004). For completeness we also report in § 2.2.4 the equations for the slow sorption mechanism as they are implemented in our solver.

2.2.3 Diffusion-controlled sorption

In the case of fast (as opposed to kinetically controlled transport) sorption, the ad- and desorption rates are locally in equilibrium, i.e.

(2.15) $$\begin{eqnarray}s^{ads}(c_{|\unicode[STIX]{x1D6F4}},c^{\unicode[STIX]{x1D6F4}})=s^{des}(c^{\unicode[STIX]{x1D6F4}}).\end{eqnarray}$$

This equality leads to an additional local relationship between $c^{\unicode[STIX]{x1D6F4}}$ and $c_{|\unicode[STIX]{x1D6F4}}$ , the so-called adsorption isotherm, which needs to be accounted for in the numerical solution. For instance, the Langmuir adsorption isotherm relates the surface and bulk surfactant concentrations by means of the Langmuir equilibrium constant $a$ , expressed in $\text{mol}~\text{m}^{-3}$ , and the saturated surface concentration:

(2.16) $$\begin{eqnarray}c^{\unicode[STIX]{x1D6F4}}=c_{\infty }^{\unicode[STIX]{x1D6F4}}\frac{c/a}{1+c/a}.\end{eqnarray}$$

2.2.4 Kinetically controlled sorption

In the case of kinetically controlled sorption, the source term at the interface is computed as

(2.17) $$\begin{eqnarray}s^{\unicode[STIX]{x1D6F4}}=s^{ads}(c_{|\unicode[STIX]{x1D6F4}},c^{\unicode[STIX]{x1D6F4}})-s^{des}(c^{\unicode[STIX]{x1D6F4}}),\end{eqnarray}$$

where $s^{ads}(c_{|\unicode[STIX]{x1D6F4}},c^{\unicode[STIX]{x1D6F4}})$ and $s^{des}(c^{\unicode[STIX]{x1D6F4}})$ describe the rate of ad- and desorption, respectively. Note that the rate of adsorption is a function of the bulk concentration near the interface and the concentration of the adsorbed species, while the desorption rate is usually assumed to be a function of the adsorbed species only. From (2.9) and the diffusive fluxes according to (2.10), a Neumann boundary condition for the bulk species equation is derived, namely

(2.18) $$\begin{eqnarray}(\unicode[STIX]{x1D735}c)_{|\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=-s^{\unicode[STIX]{x1D6F4}}/D.\end{eqnarray}$$

3.2.1 Equation discretization

A finite volume method is applied to discretize the species transport equation in the liquid phase. In this case, the transported quantity is the surfactant molar concentration $c$ . The transport equation in integral form can be derived from (2.6). Applying Fick’s law (2.10) to describe the diffusive fluxes it reads

(3.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int _{V(t)}c\,\text{d}V+\int _{\unicode[STIX]{x2202}V(t)}(c(\boldsymbol{v}-\boldsymbol{w})-D\unicode[STIX]{x1D735}c)\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S=0.\end{eqnarray}$$

The fully discretized transport equation for the control volume $V_{P}$ then reads

(3.2) $$\begin{eqnarray}\frac{3c_{P}^{n}V_{P}^{n}-4c_{P}^{o}V_{P}^{o}+c_{P}^{oo}V_{P}^{oo}}{\unicode[STIX]{x0394}t}+\mathop{\sum }_{f}\unicode[STIX]{x1D719}_{f}c_{f}^{n}=\mathop{\sum }_{f}D_{f}(\unicode[STIX]{x1D735}c)_{f}^{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{f},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{f}=\unicode[STIX]{x1D64E}_{f}\boldsymbol{\cdot }(\boldsymbol{u}-\boldsymbol{w})_{f}$ is the face flux. We denote the discrete velocity as $\boldsymbol{u}$ to distinguish between the discrete and the continuous quantities. The superscripts $n,o$ and $oo$ represent values evaluated at the new time instance $t^{n}$ and the two previous time instance $t^{o}=t^{n}-\unicode[STIX]{x0394}t$ and $t^{oo}=t^{o}-\unicode[STIX]{x0394}t$ . The discretized concentration field is defined in the cell centres $P$ as $c_{P}$ . Then, as required by the discretization of the diffusive and convective terms, the quantities $(\unicode[STIX]{x1D735}c)_{f}$ and $c_{f}$ have to be approximated at the face centres.

The surfactant transport on the interface can also be derived from (2.6). Applying the finite area method, the local discretized form of the equation is obtained for the control surface $S_{P}$ ,

(3.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{3(c_{P}^{S})^{n}(S_{P})^{n}-4(c_{P}^{S})^{o}(S_{P})^{o}+(c_{P}^{S})^{oo}(S_{P})^{oo}}{\unicode[STIX]{x0394}t}+\mathop{\sum }_{e}(\unicode[STIX]{x1D719}_{e}^{S}c_{e}^{S,n})\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{e}D_{e}^{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D735}_{S}c^{S})_{e}^{n}\boldsymbol{\cdot }(\boldsymbol{m}_{e}L_{e})+s_{P}^{S}S_{P}\end{eqnarray}$$

with the relative edge flux $\unicode[STIX]{x1D719}_{e}^{S}=(\boldsymbol{m}_{e}L_{e})\boldsymbol{\cdot }(\boldsymbol{u}-\boldsymbol{w})_{\Vert }$ and $\unicode[STIX]{x1D735}_{S}$ representing the discrete counterpart of the surface gradient operator $\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}$ . The quantity $c^{S}$ denotes the discretized counterpart of the continuous quantity $c^{\unicode[STIX]{x1D6F4}}$ in the face centre $c_{P}^{S}$ or interpolated on the edge centre $c_{e}^{S}$ . The terms $s_{P}^{S}$ is the discrete source term which can be split in explicit $s_{P,exp}^{S}$ and implicit $s_{P,imp}^{S}$ parts, respectively. In case of fast sorption processes the source term appears only in an explicit form, thus the splitting is not necessary.

The diffusion terms (bulk and surface transport) can be decomposed into orthogonal and non-orthogonal contributions, treating the first one implicitly and the second one explicitly; see Tuković & Jasak (Reference Tuković and Jasak2008).

Equations in the bulk and on the interface are solved with a preconditioned bi-conjugate gradient (PBiCG) linear solver, with a diagonal incomplete-LU preconditioner (DILU) and tolerance $1\times 10^{-12}$ .

3.2.2 Sorption process

The coupling between bulk and interfacial surfactant, transport is achieved applying a Dirichlet (fast sorption) or a Neumann (slow sorption) boundary condition to the diffusive term in (3.2) and the respective constitutive equation for the source term in (3.3) derived from the sorption model. In our code, a sorption model library is available (Pesci et al. Reference Pesci, Dieter-Kissling, Marschall, Bothe, Rahni, Karbaschi and Miller2015) where multiple sets of models, both fast and slow sorption models, are implemented. Depending on the chosen model, the solver will automatically use the respective boundary conditions and source term. As in this work we use a fast sorption model, slow sorption is not treated in this section, but its numerical treatment can be found in Pesci et al. (Reference Pesci, Dieter-Kissling, Marschall, Bothe, Rahni, Karbaschi and Miller2015, Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017).

For diffusion-controlled (fast) sorption processes, the source term for the surface concentration equation is computed from the transmission condition (2.9) as

(3.4) $$\begin{eqnarray}s^{\unicode[STIX]{x1D6F4}}=\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=-D(\unicode[STIX]{x1D735}c\boldsymbol{\cdot }\boldsymbol{n})_{\unicode[STIX]{x1D6F4}}=:s_{fast}^{\unicode[STIX]{x1D6F4}}.\end{eqnarray}$$

Then the discretized surface transport (3.3) is solved to obtain the new surface concentration field of the surfactant species. Since the adsorption isotherm $c^{\unicode[STIX]{x1D6F4}}=f(c_{|\unicode[STIX]{x1D6F4}})$ is known, e.g. equation (2.16), the value of

(3.5) $$\begin{eqnarray}c_{|\unicode[STIX]{x1D6F4}}=f^{-1}(c^{\unicode[STIX]{x1D6F4}})\quad \text{at }\unicode[STIX]{x1D6F4}\end{eqnarray}$$

is taken as a Dirichlet boundary condition for the discretized surfactant bulk (3.2). After solving the interfacial and bulk surfactant transport equations, the surface tension $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}(c^{\unicode[STIX]{x1D6F4}})$ is updated according to the chosen sorption model.

3.2.3 Implicit SGS model

Consider the species, in this case surfactant, transport problem in the liquid phase. The species transport along the bubble interface is mainly governed by two transport processes, namely advection in the streamwise direction and diffusion in the interface normal direction. To characterize the flow the following dimensionless numbers are used. The Reynolds number defined as $Re=uL/\unicode[STIX]{x1D708}$ , the Péclet number as $Pe=Re\cdot Sc=uL/D$ and the Schmidt number as $Sc=\unicode[STIX]{x1D708}/D$ , where $u$ is the velocity, $L$ the characteristic length, $\unicode[STIX]{x1D708}$ the kinematic viscosity and $D$ the molecular diffusivity. For the bubbles under investigation (Reynolds number $Re\approx 10^{2}$ and Péclet number $Pe\approx 10^{7}$ ), the species transport is dominated by advection, leading to a very thin concentration boundary layer around the bubble ( $\unicode[STIX]{x1D6FF}_{s}\approx 10^{-6}~\text{m}$ ). Thus, a fully resolved three-dimensional (3-D) DNS for the species transport is not feasible due to the high computational costs. In previous studies, for instance in Cuenot, Magnaudet & Spennato (Reference Cuenot, Magnaudet and Spennato1997), this issue was faced using a very fine grid on the axisymmetric case with bubbles at steady state, i.e. with a non-deformable interface. Moreover, the hydrodynamics is solved only in the liquid phase. This approach is not suitable for the study of the initial transient of the bubble rise and the effect of surfactant on it. An effective solution to the thin species boundary layer problem is the use of a subgrid-scale (SGS) model, a by now standard approach in mass transfer problems (Bothe & Fleckenstein Reference Bothe and Fleckenstein2013) to approximate the surfactant boundary layer in the vicinity of the bubble. The main idea behind the SGS model is to employ an appropriate model function to compute the numerical (SGS) fluxes on all cell faces of an interface cell. These SGS fluxes are used to correct the numerical fluxes to accurately predict the species transport close to the interface, even if the concentration boundary layer is fully embedded in a single cell layer. Our approach is based on the latest development of the SGS model presented in Weiner & Bothe (Reference Weiner and Bothe2017), although here the transport equation is coupled to the sorption process at the interface and solved implicitly to improve the numerical stability and to allow for larger time steps. In Weiner & Bothe (Reference Weiner and Bothe2017) the transport equations are solved explicitly with a direct modification of diffusive fluxes and concentration values at the required faces. Since our solution is implicit, i.e. the fluxes contain the unknown variable $(c_{f})^{n},(\unicode[STIX]{x1D735}c)_{f}^{n}$ , we modify the diffusion coefficient and the advective term as described in the § B.1. It has been shown by Weiner & Bothe (Reference Weiner and Bothe2017) that the SGS model can reduce the mesh resolution requirements near the interface by a factor of ten or more.

Applying the SGS model to the bulk surfactant transport results in the following discretized transport equation (from (3.2)) solved with locally modified diffusion coefficients and advection flux field:

(3.6) $$\begin{eqnarray}\frac{3c_{P}^{n}V_{P}^{n}-4c_{P}^{o}V_{P}^{o}+c_{P}^{oo}V_{P}^{oo}}{\unicode[STIX]{x0394}t}+\mathop{\sum }_{f}\unicode[STIX]{x1D719}_{f}^{SGS}c_{f}=\mathop{\sum }_{f}D_{f}^{SGS}\unicode[STIX]{x1D64E}_{f}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}c)_{f}.\end{eqnarray}$$

The derivation of $\unicode[STIX]{x1D719}^{SGS}$ and $D^{SGS}$ is reported in § B.1.

3.2.4 SGS model and fast sorption

The inverse expression of the adsorption isotherm (3.5) serves as a Dirichlet boundary condition for the bulk transport. The bulk transport is coupled to the surface balance via the source term (3.4). Also, for computing the source term, we apply the locally corrected SGS diffusion coefficients, i.e.

(3.7) $$\begin{eqnarray}s_{f_{i}^{\unicode[STIX]{x1D6F4}}}^{\unicode[STIX]{x1D6F4}}=-D_{f_{i}^{\unicode[STIX]{x1D6F4}}}^{SGS}(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\unicode[STIX]{x1D6F4}}}^{num}.\end{eqnarray}$$

4.4.1 Initial surface coverage

We vary the detachment time via pre-contaminating the bubble surface, while the initial shape deformation at detachment is neglected. Since (4.2) provides only an estimate of the initial surface coverage at release, we found it appropriate to conduct a parameter study varying $c_{0}^{\unicode[STIX]{x1D6F4}}$ for the different bulk concentrations to obtain a more precise value of the initial surface contamination.

Figure 6.

Study on the effects of the initial surface coverage on the rise velocity; simulated time $t=1~\text{s}$ .

Figure 6(a) shows that for a small initial bulk surfactant concentration the surface coverage at detachment must have been almost zero (estimated value ${\approx}1\,\%c_{eq,1}^{\unicode[STIX]{x1D6F4}}$ ), since the simulation results for $c_{0}^{\unicode[STIX]{x1D6F4}}=0~\text{mol}~\text{m}^{-2}$ are the closest to the experimental ones. After reaching the peak velocity, the bubble starts to decelerate until the rise velocity oscillates around its steady-state value. The most noticeable difference between the experimental and the numerical results for the case in figure 6(a) is that in the simulation the bubble decelerates sooner than in the experiments. This discrepancy can result from small perturbations occurring at different times for simulations and experiments. In fact, the case studied is strongly sensitive to the onset of path instability. Perturbations triggering path instabilities are caused by different mechanisms in experiments and simulations. In experiments, perturbations could derive for instance from initial shape deformations. In numerical simulations, such perturbations can be numerical errors which are highly dependent on the mesh topology. Moreover, the discrepancy between experiments and simulations is only more pronounced for the least contaminated case which is also the case with the highest oscillations in the experimental data; see figure 13(a) and table 5. For intermediate and high initial surfactant bulk concentrations, the presence of initial surface contamination is evident; see figure 6(b,c). The higher the initial bulk concentration, the more contaminated the bubble surface at release and the lower the velocity peaks. Figure 6(b) shows that the best agreement between numerical and experimental results is obtained with an initial contamination of approximately $2\,\%c_{eq,2}^{\unicode[STIX]{x1D6F4}}$ which is in agreement with the estimated value in table 3. For the highest initial bulk concentration, see figure 6(c), a very good agreement with the experimental results is found already for $c_{0}^{\unicode[STIX]{x1D6F4}}\approx 5\,\%c_{eq,3}^{\unicode[STIX]{x1D6F4}}$ , which is a smaller value than the predicted one by (4.2). In fact, with a further increase of the initial surface contamination above the $5\,\%c_{eq,3}^{\unicode[STIX]{x1D6F4}}$ , the rise velocity profile almost does not change any more.

It is also interesting to note from figure 6 that after the initial transition period, all the bubble rise velocity values present small amplitude oscillations around a similar mean velocity value.

4.4.2 Effects of the initial surface coverage on bubble shape and path

In this section, the effects of the detachment time, or better the initial surface coverage for the simulations, are investigated in terms of bubble shapes and paths. Consider the intermediate initial bulk concentration, $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ , that is the case shown in figure 6(b). The velocity profiles for the different initial surface coverages are plotted again in figure 7 but over time. In figure 7, five time instances are marked where the bubble shape and the surface coverage are then compared and studied in figure 8. In figure 8, from the bottom to the top, the five bubbles are shown in their rise at the selected time instances (every column shows one of the bubbles rising), while from left to right the initial surface concentration increases (see the surface coverage at $t=0~\text{s}$ ). The bubble surfaces are coloured by the local surfactant surface concentration. From figures 7 and 8 it is clearly visible that increasing $c_{0}^{\unicode[STIX]{x1D6F4}}$ results in a less deformed interface and a slower bubble. In fact, for $c_{0}^{\unicode[STIX]{x1D6F4}}=0\,\%c_{eq}^{\unicode[STIX]{x1D6F4}},2.15\,\%c_{eq}^{\unicode[STIX]{x1D6F4}}$ and $3.6\,\%c_{eq}^{\unicode[STIX]{x1D6F4}}$ , respectively, the bubble surface is still deforming and reaches its maximum aspect ratio ( $AR=1.27,1.1,1.06$ , respectively) with the peak velocity. During the deceleration phase the bubbles go back to a more spherical shape; see $t=0.066~\text{s}$ . For the two cases on the right of figure 8 with the highest initial surface coverage, the amount of surfactant on the interface is high enough to result in an almost not deformed interface ( $AR=1.04$ ). These bubbles accelerate until reaching the quasi-steady-state velocity and their shape remains spherical.

If we consider the latest time ( $t=0.4~\text{s}$ ) in figure 8, the bubbles have a similar velocity although, surprisingly, they do not have the same surface coverage. Moreover, with different $c_{0}^{\unicode[STIX]{x1D6F4}}$ (and/or different initial bulk concentrations $c_{0}$ , see figure 5) we obtain similar terminal velocities, but with a different final surface coverage that is not yet the equilibrium value, $c_{eq}^{\unicode[STIX]{x1D6F4}}$ , and not even close to it. To confirm this, we show in figure 9 the total amount of surfactant on the interface with respect to time. Here it can be seen that even at $t=1~\text{s}$ the total amount of surfactant on the interface is less than 30 % of the equilibrium value. For the smallest initial surface concentration, the total amount of surfactant on the interface grows more rapidly than in the other cases. This behaviour can be explained by the fact that Péclet and Reynolds numbers are higher for smaller $c_{0}^{\unicode[STIX]{x1D6F4}}$ . Also, the concentration difference between bulk and interface is larger (for a given bulk concentration and varying the initial surface concentration). This results in stronger advective transport, thus thinner concentration boundary layers. Instead, from $t\approx 0.6~\text{s}$ , when the bubbles have approximately the same terminal velocity, the total amount of surfactant on the interface grows similarly for each bubble.

In figure 10 the respective bubble paths are depicted. From the top view in figure 10(a) it can be observed that all the bubbles follow a zig-zag path, but the onset of path instability occurs later for less contaminated surfaces, as shown by the path front view, figure 10(b).

Figure 7.

Influence of the detachment time for the initial bulk concentration, $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ .

Figure 8.

Influence of the detachment time on the bubble shape and local surface coverage for the initial bulk concentration $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ .

Figure 9.

Influence of the detachment time on the total amount of surfactant on the interface divided by the respective equilibrium value for the initial bulk concentration $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ .

Figure 10.

Effects of the initial surface coverage on the bubble path for the initial bulk concentration $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ .

4.5.1 Bulk concentration

The velocity components along the rise direction $y$ and in the $x$ – $z$ plane for the three bubbles under investigation are reported in figure 12. The respective bubble paths are given in figure 13 (top view $x$ – $z$ in 13 a and lateral view $x^{\prime }$ – $y$ in 13 b, where $x^{\prime }=\sqrt{x^{2}+z^{2}}$ ).

Even though the bubbles reach a similar terminal velocity, their lateral velocity components and paths show a significant difference. The bubble rising in the least contaminated aqueous solution ( $c_{0}^{1}$ , $Ma=34$ ) follows first a helical path until it starts to oscillate around its terminal velocity ( $t\approx 0.35~\text{s}$ ) and then turns into a zig-zag path. The amplitude of this zig-zag path is around one bubble diameter. While the shift from zig-zag to a helical path was already observed for clean bubbles (Cano-Lozano et al. Reference Cano-Lozano, Martínez-Bazán, Magnaudet and Tchoufag2016) the transition from helical to zig-zag trajectory occurs only in presence of the surfactant and was first reported by Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014). Our simulation results can serve as a further confirmation of this phenomenon. For the intermediate surfactant bulk concentration ( $c_{0}^{2}$ , $c_{0}^{\unicode[STIX]{x1D6F4}}=2.15\,\%c_{eq}^{\unicode[STIX]{x1D6F4}}$ , $Ma=49$ ), see also figure 10, after the initial transient stage when the bubble accelerates and then decelerates towards its quasi-steady state, the bubble follows a zig-zag path (starting from $t\approx 0.11~\text{s}$ ) with an amplitude around 0.7 bubble diameters.

The bubble rising in the most contaminated solution ( $c_{0}^{3}$ , $c_{0}^{\unicode[STIX]{x1D6F4}}=5\,\%c_{eq}^{\unicode[STIX]{x1D6F4}}$ , $Ma=70$ ), after the initial acceleration, at $t\approx 0.22~\text{s}$ starts to follow a zig-zag path, but with a pronounced drift towards one side. Lateral migration is a known effect both from experimental and numerical works (de Vries et al. Reference de Vries, Biesheuvel and van Wijngaarden2002; Albert et al. Reference Albert, Kromer, Robertson and Bothe2015) for bubbles close to the path instability regime. Our own studies have confirmed this trend, too. For small bubbles rising in clean systems, the lateral drift is almost zero, while for larger bubbles (but not yet path unstable) a significant migration can be observed. The lateral migration can be observed also for the bubble under investigation ( $d=1.45~\text{mm}$ ) rising in clean water, as reported in figure 14. In fact, looking at the lateral components of the rise velocity (figure 14 a) it can be noticed that they are non-zero. This causes the drift visualized in the top view of the bubble path, see figure 14(b). We assume that in our set-up the instabilities are triggered by the unstructured nature of the computational mesh.

Figure 14.

Bubble rising in clean water, evidence of the later drift.

Figure 15.

Temporal evolution of the total amount of surfactant on the interface divided by the respective equilibrium values for the three selected initial surface and bulk concentrations.

Figure 16.

Global Sherwood number referred to surfactant transfer. The surface area variation is less than 3 %.

The temporal evolution of the total amount of surfactant on the interface is depicted in figure 15. It is remarkable that for all studied cases the surface coverage is much smaller than the respective equilibrium concentration. Nonetheless, the quasi-steady state terminal velocity is reached. This finding is relevant because it shows that the steady-state velocity can be reached without an equilibrium between ad- and desorption and without the bubble being ‘fully contaminated’. This situation will also have a large impact on the mass transfer processes in contaminated systems. From the slopes of the depicted curves in figure 15, it becomes visible that the bubble rising in the most contaminated liquid is adsorbing the surfactant much quicker than in the other two cases. The bubbles rising in low and medium contaminated liquid follow a similar trend, even though the bulk concentration in the latter case is approximately four times higher. There are mainly three effects causing this behaviour: (i) With increasing surfactant bulk concentration the initial concentration difference between interface and bulk increases, and therefore also the driving force for mass transfer is higher. (ii) The first effect is mitigated because at the same time the bubble accumulates surfactant quicker. (iii) Since the surfactant distribution on the interface is coupled with the bubble hydrodynamics via the Marangoni forces, the shape of the surfactant boundary layer changes. In general, an increasing amount of surfactant will slow down the bubble, and therefore decrease the advective transport which in turn decreases the driving force for mass transfer. The last effect may be expressed as the dimensionless surfactant gradient at the sub-layer in terms of the global Sherwood number, figure 16. In the initial state, when the bubble is formed in the experiment, or at the very beginning of the numerical simulation, the bubble is stagnant, and a surfactant boundary layer forms very quickly at the liquid–gas interface, driven by pure diffusion. This process is not depicted in figure 16, since the concentration difference is the highest and the boundary layer formation happens on a time scale much smaller than the course of the bubble rise from the initial release up to the quasi-steady state, i.e. $O(1)~\text{s}$ . When the bubble starts to rise it accelerates and the initial boundary layer becomes thinner due to the strong convective transport. The cleaner the system, the higher the maximum rise velocity, and hence the more pronounced this effect will be. After the initial increase in the acceleration phase, the Sherwood number decreases rapidly as the bubble decelerates. When the bubble velocity reaches a quasi-steady state, the Sherwood number for the cases with low and medium contamination keeps decreasing, but at a much slower rate. This is because the Marangoni forces are constantly increasing with increasing surface contamination. The Marangoni forces, in turn, influence the shape of the hydrodynamic boundary layer, and therefore also of the surfactant boundary layer. For the most contaminated bubble, a further increase of surfactant on the interface does not lead to an increase of the Marangoni forces. A more detailed view of all forces acting on the bubble will be given in the §§ 4.5.3 and 4.6. In figure 16 the correlations for mass transfer problems based on the boundary layer theory from Lochiel & Calderbank (Reference Lochiel and Calderbank1964) are plotted, too. Two limiting situations are considered, that is a fully mobile interface (Lochiel & Calderbank Reference Lochiel and Calderbank1964, equation (58)) and a solid particle (Lochiel & Calderbank Reference Lochiel and Calderbank1964, equation (86)). It is very interesting to notice that the global Sherwood number computed for the adsorbed surfactant tends to a value very close to the predicted one for solid particles. Moreover, for the least contaminated case, the global Sherwood number at the beginning of the rise is comparable with the one of a clean bubble. Note that there are only two reference lines given, based on the Reynolds number of the clean case, $Re=544$ , and the average Reynolds number for the contaminated cases, $Re=235$ . We did not want to place much emphasis on the mass transfer similarity to solid particles since the physical effects leading to a comparable quantitative outcome in both cases are actually very different.

So far, we have described what we could observe from the simulation results in terms of rise velocity, surface coverage and path. Nevertheless, to really disclose the bubble dynamics, a study of the local flow field in the proximity of the interface and the forces acting on the bubble surface, in particular, the local and global Marangoni forces generated by a non-uniform surface tension and their interplay with deformable interfaces, viscous and pressure forces is performed below.

Figure 17.

Vorticity contour plot ( $\unicode[STIX]{x1D714}_{y}=\pm 20~~\text{s}^{-1}$ ) at different time instances, $c_{0}=2\times 10^{-3}~\text{mol}~\text{m}^{-3}$ , $c_{0}^{\unicode[STIX]{x1D6F4}}=0$ .

4.5.2 Vorticity

The flow type around the bubble may be characterized by the vorticity $(\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\otimes \boldsymbol{u})$ contour plots in the rise direction reported here at various time instances for the three different initial surfactant bulk concentrations; see figures 17–19. Common to all the cases is the strong vorticity production already very close to the interface due to the presence of Marangoni forces. This behaviour related to the surfactant presence is not encountered for path unstable bubbles rising in clean water; see for instance the vorticity distribution in Mougin & Magnaudet (Reference Mougin and Magnaudet2006, figures 8 and 9). Moreover, at the end of each period, that is when the bubble completes a full turn (from $t_{1}$ to $t_{5}$ in figure 17 for example), the streamwise vorticity does not vanish. In the least contaminated case, the bubble follows first a helical and then a zig-zag path. This behaviour is confirmed by the vorticity contour plots in figure 17. The figures from 17(c) to 17(l) refer to time instances when the bubble path is helical, while the figures from 17(m) to 17(q) refer to the zig-zag trajectory. As already observed by other authors, Ellingsen & Risso (Reference Ellingsen and Risso2001), Mougin & Magnaudet (Reference Mougin and Magnaudet2006), Cano-Lozano et al. (Reference Cano-Lozano, Martínez-Bazán, Magnaudet and Tchoufag2016), along with the helical trajectory, the vortical structure is formed by two counter-rotating vortices of opposite sign that produce a bubble inclination in both $x$ and $z$ directions. The two vorticity regions are wrapping around each other without any symmetry plane. On the other hand, when the bubble exhibits a zig-zag trajectory, the inclination changes only in one direction. In this case, the wake structure consists of two counter-rotating vortices with a symmetry plane. Common to both trajectories, at each cycle (from one velocity peak to another which corresponds from one side to the other of the path in the $x^{\prime }$ – $y$ view) the two vortices interchange their signs. Due to the high mobility of the interface in the initial stage, the bubble reaches a high terminal velocity and deforms. After the onset of the path instability, the trajectory is helical. With increasing surface contamination, a symmetry between the wake vortices is established and the trajectory changes from helical to zig-zag. Interestingly this happens when the rise velocity is already very close to its quasi-steady value. We, therefore, conclude that not only the pure deceleration but also the indirect influence of the Marangoni forces on the flow pattern around the bubble cause the observed transition. A similar zig-zag trajectory can be observed for the bubble in figure 18. Also in this case two counter-rotating vortices with a symmetry plane are present. A different behaviour is observed for the most contaminated case; see figure 19. The bubble follows a zig-zag trajectory, but the motion is accompanied by a lateral migration. The vortical structure is composed by two counter-rotating vortices with a symmetry plane, but the duration of each half-cycle is not constant any more, as it was for the cases in figures 17 and 18, due to the drift. Considering figure 19 from $t_{1}$ to $t_{3}$ , the vorticity production is much higher than from $t_{4}$ to $t_{7}$ . This means that a bigger portion of fluid around the interface is influenced by bubble motion. Instead, at the sample times $t_{5}$ and $t_{6}$ the vorticity production is much less, thus the fluid around the bubble will be less perturbed and the drift towards the left side lasts longer. At $t_{7}$ the same conditions as in $t_{1}$ are restored. It seems to be a superimposition of clean case migration and contaminated case oscillation. A possible explanation will be given in § 4.5.3.

Figure 18.

Vorticity contour plot ( $\unicode[STIX]{x1D714}_{y}=\pm 40~~\text{s}^{-1}$ ) at different time instances, $c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$ , $c_{0}^{\unicode[STIX]{x1D6F4}}=2\,\%c_{eq,2}^{\unicode[STIX]{x1D6F4}}$ .

Figure 19.

Vorticity contour plot ( $\unicode[STIX]{x1D714}_{y}=\pm 10~~\text{s}^{-1}$ ) at different time instances, $c_{0}=5\times 10^{-2}~\text{mol}~\text{m}^{-3}$ , $c_{0}^{\unicode[STIX]{x1D6F4}}=5\,\%c_{eq,3}^{\unicode[STIX]{x1D6F4}}$ .

4.5.3 Forces acting on the interface

Several experimental works have derived correlations for the global lift and drag coefficients of single rising bubbles, e.g. Tomiyama et al. (Reference Tomiyama, Kataoka, Zun and Sakaguchi1998). In our work, we focus on the local forces acting on the interface and how they influence the integral lift and drag forces. The interfacial jump condition (2.5) is considered in order to evaluate the forces acting on the interface:

(4.3) $$\begin{eqnarray}\unicode[STIX]{x27E6}p_{tot}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D64E}^{visc}\unicode[STIX]{x27E7}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}+\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D70E},\end{eqnarray}$$

where $p_{tot}$ is the total pressure, the sum of dynamic and hydrostatic contributions. Within the algorithm, equations (2.1) to (2.5) are solved for the modified pressure, or dynamic pressure $p^{dyn}$ as we will refer to it, that is the total pressure minus the hydrostatic contribution,

(4.4) $$\begin{eqnarray}p^{dyn}=p^{tot}-p^{hydro}\end{eqnarray}$$

with $p^{hydro}:=\unicode[STIX]{x1D70C}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x}$ . This means that in (2.2) the gravity term disappears and the transmission condition (2.5) has to be adapted according to the relation (4.4), too. For clarity, we recall that $\boldsymbol{f}^{ma}=\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D70E}$ is the area specific Marangoni force, while $\boldsymbol{f}^{ca}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$ is the area specific capillary pressure force. Equation (4.3) at each interface element reads

(4.5) $$\begin{eqnarray}\boldsymbol{f}_{B}^{p_{tot}}-\boldsymbol{f}_{A}^{p_{tot}}-\boldsymbol{f}_{B}^{visc}+\boldsymbol{f}_{A}^{visc}=\boldsymbol{f}^{ca}+\boldsymbol{f}^{ma},\end{eqnarray}$$

where $\boldsymbol{f}^{\ast }$ are the area specific forces $\boldsymbol{f}^{\ast }=\boldsymbol{f}^{\ast }(\boldsymbol{x}^{\unicode[STIX]{x1D6F4}},t)$ (the superscript ‘*’ stands for ‘ $p_{tot}$ ’, ‘ $visc$ ’, ‘ $ca$ ’ or ‘ $ma$ ’), $A$ represents the liquid phase and $B$ the gas phase. The symbols $\boldsymbol{f}^{p^{tot}}$ and $\boldsymbol{f}^{visc}$ indicate the total pressure and viscous forces, respectively. Comparing the magnitude of the forces between the sides $A$ and $B$ it can be noticed that $\boldsymbol{f}_{B}^{\ast }$ is always at least one order of magnitude smaller than the respective force from the $A$ side, thus in the following analysis it will be neglected.

The local force balance at the interface (4.5) is projected in the normal and tangential directions to the interface. For the liquid side ( $A$ , dropped from here onwards) the two balances read

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle -\boldsymbol{f}^{p_{tot}}+\boldsymbol{f}_{\bot }^{visc}=\boldsymbol{f}^{ca}\quad \text{normal to }\unicode[STIX]{x1D6F4}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\Vert }^{visc}=\boldsymbol{f}^{ma}\quad \text{tangential to }\unicode[STIX]{x1D6F4}. & \displaystyle\end{eqnarray}$$

The total pressure force can be further decomposed into the hydrostatic and the dynamic contributions, i.e.

(4.8) $$\begin{eqnarray}\boldsymbol{f}^{p_{tot}}=\boldsymbol{f}^{p_{hydro}}+\boldsymbol{f}^{p_{dyn}}.\end{eqnarray}$$

Integrating the area specific forces $\boldsymbol{f}^{\ast }(\boldsymbol{x}^{\unicode[STIX]{x1D6F4}},t)$ over the interface, we get the resultant force $\boldsymbol{F}^{\ast }(t)$ on $\unicode[STIX]{x1D6F4}$ as

(4.9) $$\begin{eqnarray}\boldsymbol{F}^{\ast }(t)=\int _{\unicode[STIX]{x1D6F4}}\boldsymbol{f}^{\ast }(\boldsymbol{x}^{\unicode[STIX]{x1D6F4}},t)\,\text{d}A.\end{eqnarray}$$

Thus, the following forces are acting on the bubble surface: the hydrostatic pressure force $\boldsymbol{F}^{p_{hydro}}$ , the dynamic pressure force $\boldsymbol{F}^{p_{dyn}}$ , normal and tangential viscous forces $\boldsymbol{F}_{\bot }^{visc}$ , $\boldsymbol{F}_{\Vert }^{visc}$ , the Marangoni force $\boldsymbol{F}^{ma}$ and the capillary pressure force $\boldsymbol{F}^{ma}$ . The hydrostatic pressure force is approximately constant over time, so we do not analyse it. As can be observed from (4.6) and (4.7) the tangential viscous force is balanced by the Marangoni force. Thus we can just consider one of them, say $\boldsymbol{F}_{\Vert }^{visc}$ . For the same reason, we drop the capillary pressure force as it is equal in magnitude to the sum of total pressure force and normal viscous force. We are left with three integral forces, $\boldsymbol{F}_{\Vert }^{visc}$ , $\boldsymbol{F}^{p_{dyn}}$ and $\boldsymbol{F}_{\bot }^{visc}$ , that are decisive for understanding the bubble dynamics. Each force may be written as the sum of contributions parallel and perpendicular to the bubble velocity vector. The parallel component we refer to as drag and the remaining component as the lift force:

(4.10) $$\begin{eqnarray}\boldsymbol{F}^{\ast }(t)=\boldsymbol{F}_{Lift}^{\ast }+\boldsymbol{F}_{Drag}^{\ast },\end{eqnarray}$$

Figure 20.

Schematic representation of the lift and drag directions.

Figure 21.

Integral lift force contributions, influence of the initial bulk concentration.

Figure 22.

Integral drag force contributions, influence of the initial bulk concentration.

Figure 23.

Lift (grey) and drag (black) due to tangential viscous forces along the path. Note that the lift force is depicted ten times larger than the drag force.

as depicted in figure 20. The drag force governs the bubble acceleration/deceleration and the lift force the bubble’s change in direction. Figures 21 and 22 show the contributions from the three integral forces mentioned above to lift and drag. The different line types correspond to the various initial bulk concentrations. In order to have a common reference, the magnitude of the forces has been made non-dimensional with respect to the buoyancy force. As can be noticed from figure 21, the major contribution to the lift force is from the dynamic pressure force (up to 50 % of the buoyancy force). The tangential viscous force contribution to the lift does not exceed 5 %, while the normal viscous force contribution is below 1 %. Considering the lift contribution of the dynamic pressure and the bubbles’ paths in figure 13, one can see that a wider trajectory corresponds to a higher lift force (in terms of helical or zig-zag radius); the lower the Marangoni number, the higher the dynamic pressure force and the wider the path. We can see that the lateral motion is mainly driven by the dynamic pressure force. Whether or not the Marangoni forces/tangential viscous forces decrease the lateral motion directly will be clarified in § 4.6. From the plot of the force magnitude, we cannot draw any conclusion on the direction of the bubble motion. For instance, it is not possible to deduce from this plot when the least contaminated bubble ( $Ma=34$ ) changes its trajectory from helical to zig-zag. These aspects will be investigated later in this section; see figures 23 and 24. Consider now the force contributions to the drag force, see figure 22. As for the lift, the main contribution comes from the dynamic pressure force, although for the drag, tangential and normal viscous forces cannot be neglected. In the first graph in figure 22, the contribution of the tangential viscous force to the drag is reported. Increasing Marangoni numbers, i.e. higher surfactant concentrations, lead to a higher drag contribution of $\boldsymbol{F}_{\Vert }^{visc}$ . When the bubble reaches the quasi-steady state, after approximately 0.4 s, the tangential viscous force (as the Marangoni force) is still slowly increasing. We believe that this is due to the fact that the equilibrium value of the interfacial concentration has not yet been reached, and thus surfactant is still accumulating on the interface, changing its properties and consequently the Marangoni force. On the other hand, it can be seen from figure 22 that the drag contribution of the normal viscous force decreases with time. At the beginning of the bubble rise, there is a stronger change of the velocity normal to the interface, resulting in higher viscous stresses. In fact, the drag due to viscous forces is the highest for the lowest Marangoni number. For increasing Marangoni numbers, this contribution becomes more and more negligible; see for instance the line corresponding to $Ma=70$ . To conclude the analysis on the drag force, consider the dynamic pressure contribution to it in figure 22(c). During the initial part of the acceleration phase at the beginning of the rise, the dynamic pressure force contributions reach values comparable to the gravitational force, being the highest for the least contaminated bubble, that is the one with the highest rise velocity. After this initial phase, the contribution of the dynamic pressure force to the drag drops and oscillates at approximately 60 % of the buoyancy force. As pointed out in the previous section, all studied surfactant bulk concentrations lead to a similar quasi-steady terminal velocity, even though ad- and desorption are not in equilibrium and the total surface coverage varies significantly. The steady-state terminal velocity is a consequence of the overall drag force. For higher surfactant bulk concentrations, the viscous drag force increases due to higher surface tension gradients. At the same time, the dynamic pressure force decreases as a result of the decreasing mobility of the interface. These two counteracting effects lead to an approximately constant overall drag force. In figures 23 and 24 the integral force contributions to the lift and drag from the tangential viscous force and the dynamic pressure force are depicted as vectors along the bubble path. In the two figures, the coordinate $x^{\prime }$ correspond to the direction along which each bubble is translating in a horizontal plane. From these plots, one can clearly deduce how the forces are changing the bubble trajectory. The main contribution to the lift comes from the dynamic pressure; see figure 21. Thus, the deviation from a rectilinear path is mainly caused by the dynamic pressure force and not directly by the tangential viscous force (in response to the Marangoni force). Yet, with increasing contamination, the lateral motion of the bubble decreases, and this effect may be caused by a non-axisymmetric (with respect to the rise velocity vector) distribution of the surfactant on the interface. As can be seen in figure 23, the Marangoni effect is actually adding to the lift. However, the reduction of the dynamic pressure is much stronger, and consequently, the overall lift is reduced. Regarding the drag component, the dynamic pressure force is still the dominating contribution, but the tangential viscous force contributes in comparable amounts to the drag.

Figure 24.

Lift (grey) and drag (black) due to dynamic pressure forces along the path. Note that the lift force is depicted two times larger than the drag force.

Even though the dynamic pressure force is the dominating component, locally the flow field is governed by the Marangoni stresses. A study of the local fields is performed in § 4.6.