1 Introduction
Surface active agents, so-called surfactants, are present in most multiphase contactors, either as contaminants or added on purpose to change the way how phases interact. In froth flotation, for example, a so-called frother is used to separate hydrophobic from hydrophilic particles. The frother is surface active and renders the particles in question hydrophobic. The particles can then attach to air bubbles, which rise to the surface of the floatation cell and form a froth that can be removed. The efficiency of flotation cells is determined by the probability of bubble–particle collisions, and therefore by the interaction of gas, liquid, particles and frother. The example of froth flotation demonstrates how complex a system involving surfactants can be. But also systems as simple as a single air bubble rising in tap water may be determined by the presence of surfactants. Experiments have shown that bubbles rising in purified water can reach terminal velocities that are two times higher than in tap water; see Clift, Grace & Weber (Reference Clift, Grace and Weber1978, p. 172, figure 7.3). This demands that the used substance system must be well determined in order to obtain reliable and reproducible results. The most challenging but also most astonishing property of a surfactant is that even traces of it, which modify cohesion forces on a molecular level, can cause a tremendous change in the macroscopically observed, sometimes metre-sized, flow patterns.
Levich’s Physicochemical Hydrodynamics (Levich Reference Levich1962) is one of the first textbooks containing a theoretical treatment of surface forces resulting from an inhomogeneous distribution of a surface active substance on the interface of a rising bubble, and it also describes in much greater detail, for the interested reader, some of the basic concepts outlined hereafter. Bubbles rising in a pure liquid are characterized by a mobile interface, meaning that the fluid elements forming the gas–liquid interface are movable and can be exchanged or displaced. Therefore, the velocity gradients present in the liquid around a rising bubble are smaller than those around a solid body, and less energy is dissipated in the liquid. Consequently, under the same driving force, bubbles rise faster than solid particles. If impurities are present in the surrounding liquid, however, the observed rise velocity varies somewhere between that of particles with a fully mobile and fully immobile or rigid interface. This observation gave rise to the idea of a partially immobilized interface, which is useful to derive simplified models to account for the influence of surfactants, but which can be misleading sometimes. It is important to clarify that the inhomogeneous surfactant distribution causes additional surface specific forces which in turn change the flow pattern around a rising bubble. The surfactant itself cannot render a fluid particle (partially) rigid.
In this work a substance is called surface active if its molecules, present in the liquid bulk phase, accumulate at the gas–liquid interface and lower the surface tension. The process of accumulation is characterized by two steps (see Chang & Franses (Reference Chang and Franses1995, § 4 and the reference therein)): (i) the exchange of molecules between a surface and a subsurface layer, which is only a few molecule diameters in width, and (ii) the transfer of molecules from the bulk liquid into the subsurface layer. The first step is called adsorption and the latter (bulk) mass transfer. In this work we consider only cases of diffusion-controlled adsorption, meaning that the diffusive transport of surfactant molecules from the bulk into the subsurface layer is much slower than their adsorption such that the surfactant concentrations in the surface and subsurface layer are always locally in equilibrium. Because the interface of a rising bubble is mobile and constantly entrained by the surrounding bulk liquid, the adsorbed surfactant is transported to the rear of the bubble, where it accumulates. As a consequence, there is a region in the rear part with high surfactant concentration and lowered surface tension, while the upper part stays almost uncontaminated and the surface tension is unchanged. In the transition zone between contaminated and uncontaminated interface segments, strong gradients of surfactant concentration and surface tension result. These surface tension gradients lead to additional, so-called Marangoni forces, acting from points of low towards points of high surface tension. These tangential interface forces have to be balanced by shear forces in the liquid phase. The arising viscous forces act against the Marangoni forces from the top to the bottom and, hence, add to the overall drag force.
The described mechanisms and experimental observations led Davis & Acrivos (Reference Davis and Acrivos1966) to propose a mathematical model which incorporates the idea of a ‘stagnant cap’. The interface is divided at a certain polar angle into two rotationally symmetric segments, one fully covered with surfactant and one completely clean. The contaminated cap is stagnant, meaning that the velocity at the interface is zero in a reference frame moving with the bubble centre, and the shear stress at the cap is equal to the surface tension gradient. The clean bubble front instead is characterized by zero shear stress. The dividing angle is often referred to as stagnant cap angle. Such a clear separating circle is a strong idealization, assuming that the transition zone from fully contaminated to uncontaminated surface is small compared to the bubble size. A variety of theoretical and numerical studies based on the stagnant cap concept have appeared in the last decades, e.g. He, Maldarelli & Dagan (Reference He, Maldarelli and Dagan1991), Fdhila & Duineveld (Reference Fdhila and Duineveld1996), Liao & McLaughlin (Reference Liao and McLaughlin2000), Zhang & Finch (Reference Zhang and Finch2001), Dukhin et al. (Reference Dukhin, Kovalchuk, Gochev, Lotfi, Krzan, Malysa and Miller2015, Reference Dukhin, Lotfi, Kovalchuck, Bastani and Miller2016). One drawback of stagnant cap based models is that dynamic effects cannot be easily included, especially when the assumption of rotational symmetry is violated, as occurs in most applications. In fact, experiments show that the bubble motion is highly transient, especially after the bubble release. Sam, Gomez & Finch (Reference Sam, Gomez and Finch1996) describe the typical transient rise of single bubbles under the influence of different surface active agents (frothers) as a three stage process that was then observed several times in experiments. After release, the bubble accelerates until a maximum terminal velocity is reached; in the second stage, the rise velocity starts to reduce until, given sufficient time, a plateau is reached. The constant plateau velocity defines the third stage. Interestingly, the first and second stages depend on the liquid bulk concentration of the surfactant, but the plateau velocity in the third stage seems to be fully determined by the surfactant type alone. Furthermore, the authors observed in their experiments that all investigated bubbles (bubble diameter
$d_{b}<3~\text{mm}$
), after an initial deformation to an ellipsoidal shape, were almost spherical at the top of the column. Also, an influence of the frother concentration on the bubble path was reported: for bubbles showing path instability, the oscillation frequency decreased from the bottom to the top of the column with increasing frother concentration. Even in the case of large bubbles, the path at the column top was rectilinear. Since the work of Mougin & Magnaudet (Reference Mougin and Magnaudet2002), it is known that helical and zig-zagging trajectories of bubbles in the spherical and ellipsoidal regimes are associated with pairs of rotating or symmetric vortices in the bubble wake. Sometimes during the initial acceleration, a transition from zig-zag to helical paths can be observed. The reverse transition, from helical to zig-zag, was only reported recently by Tagawa, Takagi & Matsumoto (Reference Tagawa, Takagi and Matsumoto2014) for contaminated systems. The authors infer that a similar transition between different wake structures may happen. A strong surfactant influence on wake structure, path and shape were also visualized and comprehensively studied by Huang & Saito (Reference Huang and Saito2017a
,Reference Huang and Saito
b
). The possible impact of Marangoni forces on lift and drag was deduced from the bubble motion. All previously mentioned experimental results contribute to a partial understanding of, and description for, processes occurring on the reactor scale, for instance why the gas hold-up in flotation cells increases from the bottom to the top. However, to fully understand the transient behaviour of contaminated systems, complementary local field information of surfactant concentration, velocity and pressure at the interface and in the liquid bulk is necessary, which is currently only accessible via direct numerical simulations (DNS). Early numerical studies assuming rotational symmetry, Fdhila & Duineveld (Reference Fdhila and Duineveld1996), Liao & McLaughlin (Reference Liao and McLaughlin2000), were only able to find a qualitative agreement with the previously described experimental observations, presumably because of too many limiting assumptions in the mathematical model. However, more sophisticated, fully three-dimensional DNS solving the coupled problems of two-phase hydrodynamics, and surfactant transport in the bulk and on the interface, Tasoglu, Demirci & Muradoglu (Reference Tasoglu, Demirci and Muradoglu2008), Tuković & Jasak (Reference Tuković and Jasak2008), could also only partially reproduce and explain the typical three stage process. As we will show in the following sections this is mainly due to the studied parameter range. The authors study Péclet numbers (
$Pe$
) below
$10^{3}$
(calculated with the kinematic viscosity of the bulk liquid and the molecular diffusivity of the dissolved surfactant in the bulk). For typical systems instead,
$Pe$
ranges from
$10^{4}$
to
$10^{7}$
. The Péclet number is a measure of the ratio of convective to diffusive transport of a dilute species. High values of
$Pe$
are associated with thin boundary layers forming along the bubble surface, which determine the surfactant transfer, and hence, the ad- and desorption. The boundary layer width is approximately three to four orders of magnitude smaller than the bubble size (Weiner & Bothe Reference Weiner and Bothe2017) which is why it is extremely demanding to resolve them in a DNS.
In this work we use an arbitrary Lagrangian–Eulerian (ALE) interface-tracking approach (Muzaferija & Perić Reference Muzaferija and Perić1997; Jasak & Tuković Reference Jasak and Tuković2006; Tuković & Jasak Reference Tuković and Jasak2008, Reference Tuković and Jasak2012) combined with a recently introduced subgrid-scale model methodology (Weiner & Bothe Reference Weiner and Bothe2017) for the surfactant transfer, which allows us to study realistic systems and to find a good agreement with experimental results. The results for a single rising bubble influenced by different amounts of soluble surfactant are discussed. We present local and global quantities which explain how the surfactant distribution in the bulk and on the interface is related to the macroscopically observed bubble motion, and examine thoroughly different contributions to the overall drag and lift forces. It is the authors’ intention to provide detailed information which could lead to better scale-reduced models accounting for the influence of contamination in bubbly flows.
2 Mathematical model
The mathematical model for two-phase flows employs a sharp interface representation, meaning that the interface is represented as a surface of zero thickness with unknown time-dependent shape and location. Consider a fluid domain
$\unicode[STIX]{x1D6FA}$
containing two immiscible fluids, separated by a deformable interface. The interface,
$\unicode[STIX]{x1D6F4}(t)$
, separates the domain into two sub-domains,
$\unicode[STIX]{x1D6FA}^{+}(t)$
and
$\unicode[STIX]{x1D6FA}^{-}(t)$
, corresponding to the two bulk phases. The presence of surfactant in the denser phase and on the interface is taken into account. Under the hypothesis of incompressible Newtonian fluids, isothermal conditions and absence of phase change and chemical reactions, the governing equations are based on the conservation of mass, momentum and surfactant molar mass. For the latter, the additional assumption of negligible inertia of the adsorbed surfactant on the interface is fundamental.
2.1 Hydrodynamics
The velocity and the pressure field are obtained from the standard two-phase Navier–Stokes equations for incompressible Newtonian fluids. In local formulation, the continuity equation and the momentum balance in the bulk phases
$\unicode[STIX]{x1D6FA}^{\pm }(t)$
read
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}\boldsymbol{v})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{v}\otimes \boldsymbol{v})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}^{visc}+\unicode[STIX]{x1D70C}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{v}$
is the barycentric velocity,
$p$
the pressure,
$\unicode[STIX]{x1D70C}$
the density,
$\unicode[STIX]{x1D64E}^{visc}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{v}+(\unicode[STIX]{x1D735}\boldsymbol{v})^{\mathsf{T}})$
the viscous stress tensor with
$\unicode[STIX]{x1D707}$
being the dynamic viscosity and
$\boldsymbol{g}$
the acceleration due to gravity. The two bulk phases, separated by the moving interface
$\unicode[STIX]{x1D6F4}(t)$
, are coupled via transmission (or jump) conditions at the interface:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x27E6}\boldsymbol{v}\unicode[STIX]{x27E7}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=\boldsymbol{v}^{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x27E6}p\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}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{v}^{\unicode[STIX]{x1D6F4}}$
is the interface velocity with
$\boldsymbol{v}^{\unicode[STIX]{x1D6F4}}=\boldsymbol{v}_{|\unicode[STIX]{x1D6F4}}$
(the notation
$\cdot _{|\unicode[STIX]{x1D6F4}}$
denotes the trace of a quantity defined in
$\unicode[STIX]{x1D6FA}^{\pm }$
on the interface), and
$\unicode[STIX]{x1D705}$
the surface curvature defined as
$\unicode[STIX]{x1D705}=-\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$
, with
$\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }$
representing the surface divergence. Note that the surface gradient of a quantity
$\unicode[STIX]{x1D719}(\boldsymbol{x})$
is defined as:
$\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D719}(\boldsymbol{x})=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x})-\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}(\boldsymbol{x})(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}(\boldsymbol{x}))$
at
$\boldsymbol{x}\in \unicode[STIX]{x1D6F4}$
, where
$\unicode[STIX]{x1D719}$
is extended to a neighbourhood of
$\unicode[STIX]{x1D6F4}$
as a differentiable function. Then, the surface divergence of a vector
$\boldsymbol{f}$
is defined as
$(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{f})(\boldsymbol{x})=\text{tr}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{f})(\boldsymbol{x})$
. The symbol
$\unicode[STIX]{x1D70E}$
denotes the surface tension coefficient. In contaminated systems, the surface tension coefficient depends on the local concentration of surfactant on the interface
$\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}(c^{\unicode[STIX]{x1D6F4}})$
. The notation
$\unicode[STIX]{x27E6}\cdot \unicode[STIX]{x27E7}$
stands for the jump of a physical quantity, e.g.
$\unicode[STIX]{x1D719}$
, across the interface. The jump of
$\unicode[STIX]{x1D719}$
is defined as
$\unicode[STIX]{x27E6}\unicode[STIX]{x1D719}\unicode[STIX]{x27E7}(t,\boldsymbol{x})=\lim _{h\rightarrow 0+}(\unicode[STIX]{x1D719}(t,\boldsymbol{x}+h\boldsymbol{n}_{\unicode[STIX]{x1D6F4}})-\unicode[STIX]{x1D719}(t,\boldsymbol{x}-h\boldsymbol{n}_{\unicode[STIX]{x1D6F4}})),$
$\boldsymbol{x}\in \unicode[STIX]{x1D6F4}(t)$
. The system of equations governing the hydrodynamic problem is completed by appropriate initial and boundary conditions.
Figure 1. Domain representation for two-phase flows system.
2.2 Surfactant transport
The core part of the mathematical model consists of the surfactant transport equations in the liquid phase and on the interface for moving domains. Let
$V(t)$
be a control volume moving with velocity
$\boldsymbol{w}$
inside the fluid domain
$\unicode[STIX]{x1D6FA}$
. The boundary of the control volume is denoted by
$\unicode[STIX]{x2202}V(t)$
, with
$\boldsymbol{n}$
being the outer unit normal to
$V(t)$
. The intersection between the interface and the control volume
$\unicode[STIX]{x1D6F4}(t)\cap V(t)$
is denoted as
$S(t)$
, with the boundary curve
$\unicode[STIX]{x2202}S(t)$
and the outer unit normal
$\boldsymbol{m}\bot \boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$
to
$\unicode[STIX]{x2202}S(t)$
; see figure 1. The integral balance of surfactant molar mass for a moving control volume
$V(t)$
in absence of chemical reactions or any other source term (Bothe, Prüss & Simonett Reference Bothe, Prüss, Simonett, Escher and Chipot2005) reads
$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left[\int _{V(t)}c\,\text{d}V+\int _{S(t)}c^{\unicode[STIX]{x1D6F4}}\,\text{d}S\right] & = & \displaystyle -\int _{\unicode[STIX]{x2202}V(t)}c(\boldsymbol{v}-\boldsymbol{w})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S-\int _{\unicode[STIX]{x2202}V(t)}\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S+\nonumber\\ \displaystyle & & \displaystyle -\,\int _{\unicode[STIX]{x2202}S(t)}c^{\unicode[STIX]{x1D6F4}}(\boldsymbol{v}^{\unicode[STIX]{x1D6F4}}-\boldsymbol{w})\boldsymbol{\cdot }\boldsymbol{m}\,\text{d}l-\int _{\unicode[STIX]{x2202}S(t)}\boldsymbol{j}^{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{m}\,\text{d}l,\end{eqnarray}$$
where
$c$
is the molar concentration of surfactant in the bulk (
$\text{mol}~\text{m}^{-3}$
),
$c^{\unicode[STIX]{x1D6F4}}$
is the surface molar concentration of surfactant on the interface (
$\text{mol}~\text{m}^{-2}$
) and
$\boldsymbol{j}$
and
$\boldsymbol{j}^{\unicode[STIX]{x1D6F4}}$
are the diffusive fluxes in the bulk phase and on the interface, respectively. In local formulation the equations for surfactant transport in the bulk phase and on the interface read
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}c+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(c\boldsymbol{v}+\boldsymbol{j})=0\quad \text{in }\unicode[STIX]{x1D6FA}\setminus \unicode[STIX]{x1D6F4}(t), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}^{\unicode[STIX]{x1D6F4}}c^{\unicode[STIX]{x1D6F4}}+\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }(c^{\unicode[STIX]{x1D6F4}}\boldsymbol{v}^{\unicode[STIX]{x1D6F4}}+\boldsymbol{j}^{\unicode[STIX]{x1D6F4}})=s^{\unicode[STIX]{x1D6F4}}\quad \text{on }\unicode[STIX]{x1D6F4}(t), & \displaystyle\end{eqnarray}$$
where the sorption term
$s^{\unicode[STIX]{x1D6F4}}$
satisfies
$$\begin{eqnarray}s^{\unicode[STIX]{x1D6F4}}+\unicode[STIX]{x27E6}\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x27E7}=0\quad \text{on }\unicode[STIX]{x1D6F4}(t),\end{eqnarray}$$
and (2.8) is a dynamic boundary condition for (2.7). The transport equations in the bulk and on the interface are completed by appropriate initial and boundary conditions.
The system of equations (2.7)–(2.9) is not closed, i.e. additional relations are needed to determine the diffusive fluxes and the source terms as functions of the primitive variables. The derivation of the transport equations can be found, for instance, in Stone (Reference Stone1990), Bothe et al. (Reference Bothe, Prüss, Simonett, Escher and Chipot2005), Alke & Bothe (Reference Alke and Bothe2009).
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.
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}=-D\unicode[STIX]{x1D735}c\quad \text{in }\unicode[STIX]{x1D6FA}^{+}(t), & \displaystyle\end{eqnarray}$$
$$\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.
$$\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
$$\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
$$\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.
$$\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:
$$\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
$$\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
$$\begin{eqnarray}(\unicode[STIX]{x1D735}c)_{|\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=-s^{\unicode[STIX]{x1D6F4}}/D.\end{eqnarray}$$
3 Numerical model
The solution procedure is based on the arbitrary Lagrangian–Eulerian (ALE) interface-tracking method, originally presented by Hirt, Amsden & Cook (Reference Hirt, Amsden and Cook1974), later further developed by Muzaferija & Perić (Reference Muzaferija and Perić1997) and extended by Tuković & Jasak (Reference Tuković and Jasak2012). Collocated finite volume/finite area methods are applied to solve the transport equations on unstructured meshes of general topology with moving mesh support. The interface is represented by a computational surface mesh (boundary mesh) advected in a semi-Lagrangian manner under the enforcement of jump conditions at the interface, whereas the volume mesh is updated through automatic mesh motion with Laplacian smoothing in order to preserve a high mesh quality. The interface divides the computational domain into two disconnected sub-domains. The coupling between the two is enforced by the boundary conditions for pressure and velocity at
$\unicode[STIX]{x1D6F4}(t)$
derived from the jump conditions (2.3) to (2.5). The governing equations are discretized in time using a second-order backward scheme known also as Gear’s method; see Ferziger & Perić (Reference Ferziger and Perić1996). The two fluid domains
$\unicode[STIX]{x1D6FA}^{\pm }(t)$
are discretized by a finite number of convex polyhedral control volumes
$V_{P}$
. The centroid of the control volume is denoted by
$P$
and the one of the neighbouring cell by
$N$
. The cell faces
$f$
are of polygonal shape with area
$S_{f}$
and area normal vector
$\unicode[STIX]{x1D64E}_{f}$
. In analogy to the volume discretization, the interface
$\unicode[STIX]{x1D6F4}(t)$
is subdivided into polygonal control areas (the computational surface mesh can be seen as the boundary of the volume mesh, that is the faces approximating the interface belong to the boundary cells of the volume mesh). The centre of a control area is again denoted by
$P$
and the neighbouring one by
$N$
. The two control areas are separated by the edge
$e$
, characterized by the edge vector
$\boldsymbol{e}$
, length
$L_{e}$
and bi-normal
$\boldsymbol{m}_{e}$
(perpendicular to both
$\boldsymbol{e}$
and the edge normal vector
$\boldsymbol{n}_{e}=(\boldsymbol{n}_{1}+\boldsymbol{n}_{2})/2$
).
3.1 Hydrodynamics and mesh motion
The pressure–velocity coupling is solved by applying the iterative pressure implicit with splitting of operators (PISO) algorithm (Issa Reference Issa1986). A modified version of the Rhie–Chow interpolation suggested in Tuković & Jasak (Reference Tuković and Jasak2012) is employed to prevent a decoupling of pressure and velocity. A detailed description of the flow field solution and the mesh motion can be found in Tuković & Jasak (Reference Tuković and Jasak2012), Pesci et al. (Reference Pesci, Dieter-Kissling, Marschall, Bothe, Rahni, Karbaschi and Miller2015).
3.2 Surfactant transport
In our system only one surfactant species is considered, while in Dieter-Kissling, Marschall & Bothe (Reference Dieter-Kissling, Marschall and Bothe2015a ,Reference Dieter-Kissling, Marschall and Bothe b ) the methodology and the results for multicomponent surfactant systems in free-surface flows were presented. For cases where a fully resolved DNS for the species transport is not feasible due to high computational costs and numerical stability issues, a subgrid-scale model is employed.
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
$$\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
$$\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}$
,
$$\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
$$\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
$$\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:
$$\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.
$$\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}$$
3.3 Validation
The validation of the pure hydrodynamics has been conducted comparing with the experimental data by Duineveld (Reference Duineveld1995) for single bubbles rising in pure water and can be found in Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017). There, rise velocity and aspect ratio for bubbles with radii ranging between 0.75 and 1.0 mm were considered and found in very good agreement with the experimental data; see Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017, p. 418, figure 15.13). The validation for the sorption source term for fast and slow sorption processes can also be found in Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017). The validation of the implicit SGS model can be found in appendix C.
4 Results and discussion
A single air bubble rising in aqueous solution contaminated by surfactant is considered. For this prototypical problem, a direct comparison with experimental results is possible. The experimental data and a short description of the corresponding set-up can be found in Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017). More details on the experimental set-up are presented in Krzan & Malysa (Reference Krzan and Malysa2002), Krzan, Zawala & Malysa (Reference Krzan, Zawala and Malysa2007). Briefly, a digital camera was used to record the bubble motion at various distances from the orifice. Four to eight images of the bubble were obtained for each camera position illuminating the region of interest with a strobe frequency from 100 to 200 Hz. The higher frequency was used for the initial acceleration stage. From the distances between the subsequent positions of the bubble and knowing the strobe frequency, the local bubble velocity is computed. The measurement at each camera position was performed at least three times and mean local velocity values were calculated.
Table 1. Fluid properties.
Table 2. Surfactant (
$\text{C}_{12}\text{DMPO}$
) properties, fast Langmuir adsorption model parameters.
4.1 Simulations set-up
The material properties used in the simulations are reported in tables 1 and 2. The bubble diameter is
$d_{B}=1.45~\text{mm}$
. The initial shape of the bubble is a sphere positioned in the centre of a spherical domain with radius 20 times the bubble radius. The computational domain is divided into two sub-domains, one representing the gas phase and the other one representing the liquid phase. The meshes used for the simulations consist of polyhedral cells in the gas phase and prismatic cells with polyhedral base in the liquid phase, as can be seen in figure 2. The calculation is performed in a moving reference frame (MRF) that follows the bubble centre during its rise, while the interface is deformable. The presence of a non-inertial reference frame located in the centre of the bubble is taken into account via a correction in the momentum equation (
$\unicode[STIX]{x1D70C}\boldsymbol{a}_{MRF}$
added to the momentum equation) and the velocity boundary condition at the outer domain boundary,
$\boldsymbol{v}_{out}=-\boldsymbol{v}_{MRF}$
(the boundary condition inletOutlet available in OpenFOAM is used. The inlet velocity is set to
$-\boldsymbol{v}_{MRF}$
, at the outlet a zeroGradient condition is set). A (constant) time step
$\unicode[STIX]{x0394}t\approx O(10^{-6})~\text{s}$
is chosen to fulfil the criterion for the interface numerical stability (Tuković & Jasak Reference Tuković and Jasak2012), i.e.
$\unicode[STIX]{x0394}t<\sqrt{\unicode[STIX]{x1D70C}^{+}\text{min}(l_{PN})^{3}/(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E})}$
with
$\text{min}(l_{PN})$
being the minimum distance between two face centres on the interface. The surfactant used in the experiments is the non-ionic dodecyl-dimethyl-phosphine-oxide (
$\text{C}_{12}\text{DMPO}$
). Its sorption process is modelled via the fast Langmuir sorption model. For the simulations the bubble shape is initialized as a sphere with zero initial velocity. To model the surfactant transport in the bulk phase in the vicinity of the interface, the SGS model described in § 3.2.3 and appendix B is used. From the available experimental data we consider the clean case and three different initial surfactant concentrations as a reference. The surface diffusivity
$D^{\unicode[STIX]{x1D6F4}}$
is only an estimate since it is not possible to accurately measure it. Nevertheless, a parameter study with
$D^{\unicode[STIX]{x1D6F4}}$
varying in the range of
$[10^{-6}\cdots 10^{-9}]~\text{m}^{2}~\text{s}^{-1}$
confirmed that its variation has only a minor effect on the sorption dynamics and rise velocity because the transport is advection dominated.
Figure 2. Three-dimensional computational domain for a rising bubble. Inner, outer and surface (dark grey on the right) meshes.
The selected experimental results from Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017) are given in figure 3 and they will be the base for our discussion of the simulation results. According to Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017), the average accuracy of the experimental data is
$\pm 5\,\%$
. Three different initial concentrations in the liquid phase (the surfactant concentration in the gas phase is set to zero) are considered, a relatively small one,
$c_{0,1}=2\times 10^{-3}~\text{mol}~\text{m}^{-3}$
, an intermediate one,
$c_{0,2}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$
, and a relatively high one,
$c_{0,3}=5\times 10^{-2}~\text{mol}~\text{m}^{-3}$
. To these three initial surfactant concentrations correspond Marangoni numbers
$Ma$
of 34, 49 and 70, which express the ratio between surface tension and viscous forces.
$Ma$
is computed as
$$\begin{eqnarray}Ma=\frac{RTc_{\infty }^{\unicode[STIX]{x1D6F4}}}{\unicode[STIX]{x1D707}^{+}U_{max}},\end{eqnarray}$$
where
$U_{max}$
is the peak rise velocity reached by the bubble. Moreover, the respective surface equilibrium concentrations computed from the Langmuir isotherm (2.16) are
$c_{eq,1}^{\unicode[STIX]{x1D6F4}}=1.2175\times 10^{-6}~\text{mol}~\text{m}^{-2}$
,
$c_{eq,2}^{\unicode[STIX]{x1D6F4}}=2.596\times 10^{-6}~\text{mol}~\text{m}^{-2}$
and
$c_{eq,3}^{\unicode[STIX]{x1D6F4}}=3.801\times 10^{-6}~\text{mol}~\text{m}^{-2}$
. In figure 3, the well-known velocity profile of rising bubbles under the effects of surfactants can be observed. The bubble rising in clean water (crosses), thus with a fully mobile surface, after an initial acceleration reaches a constant velocity that is the terminal velocity. The same can be observed for bubbles rising in highly contaminated solutions (filled circles). After an initial acceleration, the bubble velocity reaches a constant value, although it is much lower than the velocity for a mobile surface. At intermediate concentrations (empty circles, triangles) there is still an acceleration phase, but after reaching the peak velocity the bubble decelerates. The bubbles keep decelerating until they reach a quasi-steady terminal velocity which is similar to the case with very high contamination.
Figure 3. Experimental bubble centre rise velocities in rise direction
$y$
. Data from Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017).
In applications involving bubbly flows, it is fundamental to correctly reproduce the initial transient stage of the bubble rise, because it determines the position of the bubble and perhaps also how it will interact with other bubbles. Thus in §§ 4.2–4.4 the attention is focused on correctly reproducing the transient velocity profiles.
4.2 Discussion on under-resolved species boundary layers
The simulation results for the clean case have already been compared to the experimental ones in Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017) showing a very good agreement. These results are reproduced in § 4.5.1 with additional information about the bubble path.
Figure 4. Comparison between simulations without (black lines) and with SGS model (grey line);
$c_{0}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$
; simulated time
$t=1~\text{s}$
.
The surfactant transport problem is a typical case with highly nonlinear concentration profiles at the interface in a very thin boundary layer. Thus a standard linear interpolation from the cell centres to the face centres leads to over- or underestimated diffusive and convective fluxes normal to the interface, resulting in an unphysically thick boundary layer. Only thanks to the application of the SGS model described in § 3.2.3 it becomes possible to study cases with real diffusion coefficients for the surfactant in the liquid phase. The usage of physical diffusivities is imperative to get the correct transient velocity since it not only affects the surfactant bulk transport but also the sorption mechanism itself, as described in § 3 and in Pesci et al. (Reference Pesci, Marschall, Ulaganathan, Kairaliyeva, Miller, Bothe, Bothe and Reusken2017). A comparison between the standard interpolation and the flux correction by the SGS model is given in figure 4. The results there refer to the intermediate surfactant bulk concentration
$c_{0,2}=8\times 10^{-3}~\text{mol}~\text{m}^{-3}$
. The first set of simulations is run without SGS modelling to test the sensitivity to different diffusivities with a fixed mesh resolution (first cell thickness
$l\approx 16~\unicode[STIX]{x03BC}\text{m}$
). For a realistic diffusivity, the rise velocity is overpredicted; see figure 4. On the other hand, increased diffusion coefficients result in thicker species boundary layers that can be resolved by this mesh, but at the same time they speed up the adsorption process and, consequently, the rise velocity approaches the steady-state value too quickly.
Figure 4 depicts also the velocity profile obtained with the SGS approach and the physical diffusivity. The initial transient velocity is reproduced much better, but the velocity peak is still overestimated. This difference can be explained considering the bubble formation and detachment time in the experiments. As it is known from experimental works, e.g. Krzan et al. (Reference Krzan, Zawala and Malysa2007), Małysa et al. (Reference Małysa, Zawala, Krzan, Krasowska, Miller and Liggieri2011), Ulaganathan (Reference Ulaganathan2016), the initial transient velocity depends strongly on the time of bubble formation and release. During the bubble formation process, the newly generated bubble surface is exposed to the contaminated solution. Thus, when the bubble detaches from the capillary, its interface holds already a certain amount of surfactant. This relatively small (not above
$10\,\%$
of
$c_{eq}^{\unicode[STIX]{x1D6F4}}$
) initial surface contamination influences the peak rise velocity. From the experiments, the adsorption time for detaching bubble is known to be approximately 1.6 s, hence, during this time there would be a diffusion of surfactant towards the growing bubble surface. The surface coverage at release is a function of time and bulk surfactant concentration, and it can be estimated as
$$\begin{eqnarray}c_{0}^{\unicode[STIX]{x1D6F4}}(t)=\frac{1}{3}\left(2c_{0}\sqrt{\frac{3Dt}{7\unicode[STIX]{x03C0}}}\right),\end{eqnarray}$$
a formula taken from Dukhin, Kretzschmar & Miller (Reference Dukhin, Kretzschmar and Miller1995, pp. 118–119). A summary of the estimated surface coverages at detachment is reported in table 3. Within our simulation set-up, different detachment times can be investigated varying the initial surfactant surface concentration. Before presenting these results, a mesh sensitivity study of the full problem with SGS modelling is necessary. Note that for the simulations corresponding to figures 4 and 5 the initial surface concentration was set to zero,
$c^{\unicode[STIX]{x1D6F4}}(t=0)=0~\text{mol}~\text{m}^{-2}$
.
Figure 5. Rise velocity for three initial surfactant bulk concentrations. Results for two mesh resolutions (continuous lines – fine mesh; dotted lines – coarse mesh); simulated time up to
$t=0.6~\text{s}$
.
Table 3. Initial surface coverage estimates at release time
$t_{rel}=1.6~\text{s}$
with
$D=5\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$
.
4.3 Mesh sensitivity study
To study the dependence of the numerical results with respect to the mesh resolution, simulations with different initial bulk concentrations and zero initial surface coverage are performed on two different meshes, a fine one (
${\approx}$
320 000 cells) with a first layer thickness of
$l\approx 8~\unicode[STIX]{x03BC}\text{m}$
and 3700 faces on the interface, and a coarser one (
${\approx}$
160 000 cells) with a first layer thickness of
$l\approx 16~\unicode[STIX]{x03BC}\text{m}$
and 2400 faces on the interface. As can be noticed from figure 5 the biggest difference between fine and coarse mesh is encountered in the decelerating phase for the smallest initial bulk concentration. In fact, for higher
$c_{0}$
, the bubble rises slower, thus the Reynolds number is smaller and consequently the hydrodynamic boundary layer thicker. A thicker hydrodynamic boundary layer is then well resolved by a coarser mesh, too. Even though there is a small difference between the coarse and the fine mesh results, for the simulations that are reported below we decided to use the coarser mesh because of the required computational time. Only for the least contaminated bubble, the bubble path is reported both for the coarse and the fine meshes, since the helical path was more pronounced in the latter case. The fine and coarse cases ran in parallel (MPI) on three and five cores, respectively, with the interface (liquid side) and its counterpart (gas side) on the same processor. The computations took between thirty and sixty days to reach 1 s of simulated physical time.
4.4 Bubble shape and path under the influence of surfactant
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 Effects of surfactant on the onset of path instability
The simulation results which agree best with the experimental ones from figure 6 (grey curves) are selected for the rest of the discussion and reported in figure 11. All the simulations run until reaching
$t=1~\text{s}$
of physical time. In figure 11, also the estimated velocity from the correlation for fully contaminated systems proposed by Tomiyama et al. (Reference Tomiyama, Kataoka, Zun and Sakaguchi1998, equation (33)) is plotted. All the simulation results, including the least contaminated case, are in very good agreement with this estimated velocity at quasi-steady state. A further indicator of agreement between experimental and numerical results at quasi-steady state may be the comparison of the standard deviation of the rise velocities for
$0.1~\text{m}<y<0.16~\text{m}$
. In fact, the numerical results show pronounced oscillations that are not clearly visible from the experiments. The values for the standard deviation reported in table 4 show a similar trend, i.e. oscillations decrease with increasing bulk concentration. The magnitude is in agreement between simulations and experiments, too. For completeness, also the frequency of the horizontal velocity and the vortex shedding from the rear part of a rigid sphere are computed as reported in Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014) from Tsuge & Hibino (Reference Tsuge and Hibino1971) and Kim & Pearlstein (Reference Kim and Pearlstein1990), respectively, and compared to the simulation results. The frequency of the bubble horizontal velocity is computed according to Tsuge & Hibino (Reference Tsuge and Hibino1971) as
$f=(u_{b}/d_{e})0.1c_{D}^{0.734}$
, where
$u_{b}$
is the averaged quasi-steady velocity and
$d_{e}$
is the bubble equivalent diameter. The frequency of the vortex shedding from the rear part of a rigid sphere is
$f_{v}=\unicode[STIX]{x1D714}\unicode[STIX]{x1D708}_{l}Re/\unicode[STIX]{x03C0}d_{e}^{2}$
(from Kim & Pearlstein Reference Kim and Pearlstein1990), where
$\unicode[STIX]{x1D714}$
is taken equal to 0.30 as in Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014), Kim & Pearlstein (Reference Kim and Pearlstein1990) and
$Re$
is computed based on
$u_{b}$
. The drag coefficient
$c_{D}$
is computed equating the drag to the buoyancy force as in Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014, equation (2.7)), thus
$c_{D}=4d_{e}g/3u_{b}^{2}$
, where
$g$
is the gravitational acceleration. From table 5 it can be seen that the oscillation frequencies of the vertical velocity are approximately twice the horizontal ones, as expected from de Vries, Biesheuvel & van Wijngaarden (Reference de Vries, Biesheuvel and van Wijngaarden2002), Mougin & Magnaudet (Reference Mougin and Magnaudet2006), Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014). Moreover, the least and intermediate contaminated cases show a good agreement between the numerical and literature results. The most contaminated case is not relevant for this comparison since the velocity oscillations are not as regular as the other two cases.
Figure 11. Bubble rise velocity, influence of the initial bulk concentration with pre-contaminated surface.
Table 4. Standard deviation from the mean velocity value at quasi-steady state,
$0.1~\text{m}<y<0.16~\text{m}$
.
Table 5. Oscillation frequencies of the velocity components compared to the frequencies
$f$
and
$f_{v}$
reported in Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014) from Tsuge & Hibino (Reference Tsuge and Hibino1971) and Kim & Pearlstein (Reference Kim and Pearlstein1990), respectively.
Figure 12. Bubble rise velocity, influence of the initial bulk concentration with
$u_{x^{\prime }}=\sqrt{u_{x}^{2}+u_{z}^{2}}$
.
Figure 13. Study on the effects of the initial bulk concentration on the bubble path.
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:
$$\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,
$$\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
$$\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
$$\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}$$
$$\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.
$$\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
$$\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:
$$\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.
4.6 Local velocity and surface fields under the influence of surfactant
Figure 25 shows the velocity field in the liquid phase close to the bubble, while on the bubble surface the local Marangoni force vectors are depicted for the three initial concentrations at different time instances. At
$t=0.072~\text{s}$
the bubble rising in the most contaminated solution has already reached a surfactant distribution characteristic of the steady state. In the lower hemisphere, where the surfactant concentration is the highest and uniformly distributed, the Marangoni force is almost zero, while the surface coverage is not yet the equilibrium one. In fact, the surfactant species is still adsorbed, see figure 15. For the other two initial bulk concentrations, a longer initial transient stage is visible. The surface coverage is much smaller at the beginning of the rise, while much higher and more confined Marangoni stresses are visible. For the cases on the left and in the middle of figure 25, it is clearly visible that the line where the flow detaches corresponds to the region where the Marangoni forces are the highest. As the bubbles are rising, more and more surfactant is adsorbed and the region where the Marangoni stresses are present moves towards the upper hemisphere. The bubble in the middle, at
$t=0.9~\text{s}$
has reached a similar state as the most contaminated bubble in terms of Marangoni stresses and terminal velocity, even though the surface coverage is approximately 60 % less; see figure 15. It is reasonable to predict that the least contaminated bubble, if simulated for a longer time, would reach a similar state as the other two bubbles, but with an even lower surface coverage.
Figure 25. Velocity vectors (bulk) and Marangoni forces (interface) at different time instances for
$Ma=34,49,70$
.
Figure 26. From left to right, surfactant distribution on
$\unicode[STIX]{x1D6F4}$
, interface velocity field and sorption source term at different time instances for the intermediate bulk concentration
$Ma=49$
.
To have a better understanding of the variation of the Marangoni forces and their local distribution, one can consider the adsorption, advection and diffusion processes on the interface; see figures 25 and 26. Three different stages during the bubble rise can be identified.
(i) After being released, the bubble undergoes a strong acceleration due to the buoyancy force. The surface coverage is low and uniform and, therefore, the interface is fully mobile. A thin concentration boundary layer forms at the interface and the adsorption rates are the highest. The first stage may be very short, depending on the initial surface and bulk concentrations.
(ii) Due to the high mobility of the interface, the surfactant is quickly advected to the rear part of the bubble. As a consequence, the surface coverage becomes less uniform and surface tension gradients that are strong enough to locally reduce the tangential interface velocity in the rear part arise. The flow detaches, and vortices are shed. The interface below the detachment ring is almost stagnant, and the adsorption rates are small because the concentration difference with respect to the bulk decreases and no new surfactant is transported there by convection. The front of the bubble is still mobile and the adsorbed surfactant is quickly transported towards the cap. As a consequence, the transition from a very small to a very high contamination happens in a small belt above the ‘stagnant cap’ zone. Here the highest surface tension gradient and hence Marangoni forces are observed.
(iii) The transition from the second to the third stage happens on a larger time scale than between the first two stages. The convective surfactant transport in the bubble front slowly decreases. This happens, on the one hand, because the bubble decelerates (for small Marangoni numbers), and on the other hand due to the decreasing overall mobility of the interface. The narrow transition zone with high concentration gradients widens and the surfactant distribution in the front becomes approximately linear. Consequently, the resulting Marangoni forces have a smaller magnitude but act almost uniformly on the entire upper hemisphere. The integral tangential viscous force due to the Marangoni stresses is, therefore, higher than in stage two.
To see a further transition to a fourth stage, a much longer physical time would have to be simulated since also the adsorption steadily decreases. Such an investigation shall be part of future studies.
5 Conclusion and outlook
The focus of the current work is on the dynamics of single bubbles rising in a contaminated solution with surfactant. Within this study, it has been possible to investigate realistic length and time scales thanks to a subgrid-scale model, and the available experimental results for the rising bubble case could be reproduced well. The necessity of a subgrid-scale model has been proven via specific test cases involving thin species boundary layers. Note that the same methodology that allowed us to simulate realistic surfactant systems can be applied to mass transfer problems to eventually study the effect of surfactant on mass transfer.
We firstly investigated the influence of the initial surface coverage on the rise velocity. In fact, in the experiments, there is a certain detachment time including the bubble formation until the release. In this time adsorption mechanisms are already occurring, such that the bubble is pre-contaminated at release. The results show that the initial transient stage is very sensitive to the initial surface concentration. With a parameter study varying the initial surface contamination, we could find the initial surface coverage corresponding to the experiments, a value that was not known a priori. For very high bulk concentrations, we demonstrated that a lower initial surface contamination than the one suggested by the theory (equation (4.2)) was already sufficient to obtain the correct bubble transient velocity. This information is fundamental in view of application cases because from the initial stage depends, for instance, the position of the bubble in a channel or column.
The focus then moved on to study the influence of the initial bulk concentration on the rise velocity and bubble dynamics. From the simulation results, global and local quantities can be evaluated. The bubble path depends both on the initial surface and bulk contaminations. For the least contaminated case, a transition from helical to a zig-zag path is observed, as in the experimental work by Tagawa et al. (Reference Tagawa, Takagi and Matsumoto2014). It has also been found that the quasi-steady-state velocity can be reached without an equilibrium of ad- and desorption. Moreover, the transfer of surfactant in the sub-layer in a steady-state regime for the bubble rise velocity is close to the mass transfer at a solid particle. The local vorticity fields have been used to characterize the flow type in the vicinity of the bubble to understand the formation of vortices in the bubble wake.
The forces acting on the bubble surface have been studied considering their contribution to lift and drag forces. The dynamic pressure force, being the major contributor to the lift force, is responsible for the deviation from a rectilinear path. The steady-state terminal velocity is a consequence of the overall drag force. In fact, 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 due to the reduced mobility of the interface. These two counter-acting effects lead to an approximately constant overall drag force. In other Reynolds regimes, for example for very small bubbles such as those considered in Takemura (Reference Takemura2005) that rise along a straight path even if contaminated, these mechanisms could perhaps be different.
From the local distribution of the Marangoni forces, it has been shown that the detachment of the flow from the bubble surface occurs where the Marangoni stresses are the highest. The quasi-steady-state situation corresponds to a more uniform distribution of the Marangoni forces on the upper hemisphere of the bubble surface. These findings are relevant for deriving simplified models such as an improved stagnant cap model. In fact, one should refer to the quasi-steady state, not in terms of ‘fully contaminated’ surface, but regarding a certain Marangoni stress distribution. The latter depends on the surfactant distribution on the interface and, above a certain threshold, not on the amount of surfactant on
$\unicode[STIX]{x1D6F4}$
. This implies that at steady state the surface concentration is not necessarily equal to the equilibrium concentration.
Considering the local adsorption, advection and diffusion processes at the interface, three different stages during the bubble rise have been identified. A first stage where the adsorption rates are the highest, a second stage where the transport at the front of the bubble is advection dominated while in the rear part it is diffusion dominated, and a third stage with a uniform distribution of the Marangoni stresses in the upper hemisphere of the bubble. A further transition to a fourth stage is foreseeable, but a much longer physical time would have to be simulated since also the adsorption steadily decreases. Such an investigation shall be part of future studies.
Acknowledgements
We kindly acknowledge the financial support by the German Research Foundation (DFG) within the Priority Program SPP1740 ‘Reactive Bubbly Flows’, Project BO1879/13-2, and the Collaborative Research Center 1194 ‘Interaction between Transport and Wetting Processes’, Project B02. We also thank Dr. habil. Reinhard Miller and his group at the Max Plank Institute of Colloids and Interfaces for providing the experimental data and the anonymous reviewers for their insightful comments. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt.
Figure 27. Overview of the algorithm to solve the full problem: hydrodynamics with mesh motion, surfactant transport and sorption.
Appendix A. Algorithm overview
In figure 27, a schematic overview of the numerical solution procedure is depicted.
Appendix B. Implicit SGS model for advection-dominated problems
B.1 Implicit SGS model description
The SGS model for advection-dominated transport is based on a simplified 2-D problem formulation of the species advection–diffusion equation (2.7). Consider the species transport in the vicinity of a bubble surface. Close to the interface
$\unicode[STIX]{x1D6F4}$
, a situation as sketched in figure 28 is encountered.
Figure 28. Simplified 2-D model for species transport close to the bubble surface, figure based on Weiner & Bothe (Reference Weiner and Bothe2017).
For high Péclet numbers, constant species concentration in the gas phase (the diffusivity in the gas phase is much higher than that of the liquid phase) and a fully developed and quasi-stationary boundary layer, equation (2.7) can be reduced to
$$\begin{eqnarray}v\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}=D\frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x^{2}}\quad \text{for }x\geqslant 0\text{ and }y\geqslant 0\end{eqnarray}$$
with the boundary conditions
$$\begin{eqnarray}c(x,y=0)=c_{\infty },\quad c(x\rightarrow \infty ,y>0)=c_{\infty },\quad c(x=0,y>0)=c_{|\unicode[STIX]{x1D6F4}}.\end{eqnarray}$$
This problem has an analytical solution, describing the species distribution normal to the interface for a given boundary layer thickness
$\unicode[STIX]{x1D6FF}(y)$
,
$$\begin{eqnarray}c(x,y)=c_{|\unicode[STIX]{x1D6F4}}+(c_{\infty }-c_{|\unicode[STIX]{x1D6F4}})\;\text{erf}\left(\frac{x}{\unicode[STIX]{x1D6FF}(y)}\right)\end{eqnarray}$$
with
$\unicode[STIX]{x1D6FF}(y)=\sqrt{4Dy/v}$
. The physical profile derived from the local substitute problem is adopted to compute the fluxes over the faces in the interface cells. The free model parameter
$\unicode[STIX]{x1D6FF}$
is computed iteratively to be consistent with the cell centred concentration value. The computation of the SGS model parameter is reported in § B.2.
Consider now the discretized species (surfactant) transport equation in the liquid phase (3.2), and reported here in a condensed form,
$$\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}F_{f}^{A}=\mathop{\sum }_{f}F_{f}^{D},\end{eqnarray}$$
where
$F_{f}^{A}=\unicode[STIX]{x1D719}_{f}c_{f}^{n}$
and
$F_{f}^{D}=D_{f}(\unicode[STIX]{x1D735}c)_{f}^{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64E}_{f}$
are the advective and diffusive species fluxes, respectively. Recall from § 2 that this equation is completed by the initial condition
$$\begin{eqnarray}c(t=0,\boldsymbol{x})=c_{0},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA}^{+}(t=0)\end{eqnarray}$$
and the Dirichlet boundary condition imposed at the bubble surface
$\unicode[STIX]{x1D6F4}(t)$
in case of fast sorption (as outlined in § 2.2.2 and (3.5)), i.e.
$$\begin{eqnarray}c(t,\boldsymbol{x})=f^{-1}(c^{\unicode[STIX]{x1D6F4}}(t)),\quad \boldsymbol{x}\in \unicode[STIX]{x1D6F4}(t).\end{eqnarray}$$
When applying the SGS model, the goal is to correctly represent the species distribution around the interface, even if the concentration boundary layer is completely contained in the first cell layer (i.e. when the DNS cannot resolve the boundary layer). To achieve this, a correction of the diffusive and advective species fluxes is introduced on the first cell faces normal to
$\unicode[STIX]{x1D6F4}$
to counteract the otherwise overestimated numerical fluxes.
Figure 29. Two-dimensional sketch for the SGS model with enlarged view of the region near the interface.
$\tilde{\unicode[STIX]{x1D6FA}}^{\pm }(t),\tilde{\unicode[STIX]{x1D6F4}}(t)$
are the discretized counterparts of
$\unicode[STIX]{x1D6FA}^{\pm }(t),\unicode[STIX]{x1D6F4}(t)$
.
B.1.1 Diffusion
The diffusive species fluxes
$F_{f}^{D}$
at the faces belonging to
$\unicode[STIX]{x1D6F4}$
,
$f_{i}^{\unicode[STIX]{x1D6F4}}$
and at the first cell faces opposite to
$\unicode[STIX]{x1D6F4}$
,
$f_{i}^{\unicode[STIX]{x1D6F4},o}$
, are considered; see figure 29 for the notation. We compute the desired numerical diffusive fluxes at the relevant faces
$f_{i}^{\ast }=f_{i}^{\unicode[STIX]{x1D6F4}},f_{i}^{\unicode[STIX]{x1D6F4},o}$
as
$$\begin{eqnarray}F_{f_{i}^{\ast }}^{D,num}=-D_{f_{i}^{\ast }}S_{f_{i}^{\ast }}(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{num},\end{eqnarray}$$
where
$D_{f_{i}^{\ast }}$
is a corrected diffusion coefficient to counteract the numerical effects of the under-resolved species boundary layer. To derive an expression for
$D_{f_{i}^{\ast }}$
we use the diffusive fluxes coming from the SGS modelling
$$\begin{eqnarray}F_{f_{i}^{\ast }}^{D,SGS}=-DS_{f_{i}^{\ast }}(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{SGS},\end{eqnarray}$$
where
$D$
is the molecular diffusivity and
$(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{SGS}$
is provided by the SGS model; see § B.2 for the analytical expression. Our goal is to compute
$D_{f_{i}^{\ast }}$
such that the numerical diffusive fluxes, coming from the standard discretization, equal the SGS fluxes,
$$\begin{eqnarray}F_{f_{i}^{\ast }}^{D,num}\overset{!}{=}F_{f_{i}^{\ast }}^{D,SGS}.\end{eqnarray}$$
Thus, we impose
$$\begin{eqnarray}D_{f_{i}^{\ast }}(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{num}\overset{!}{=}D(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{SGS},\end{eqnarray}$$
to get an expression for the modified diffusion coefficients to be substituted in the discretized transport equation,
$$\begin{eqnarray}D_{f_{i}^{\ast }}=D\frac{(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{SGS}}{(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{num}}.\end{eqnarray}$$
To simplify the notation, below we will address
$D_{f_{i}^{\ast }}$
as
$D^{SGS}$
, where
$D^{SGS}$
contains the modified local values from the SGS model in the required faces. For the other faces the standard molecular diffusivity is kept. In case the estimated boundary layer thickness is more than 1000 times larger than the first cell width, the SGS correction is not applied to avoid non-physical diffusive fluxes; see § B.2.
B.1.2 Advection
The SGS correction of the advective species fluxes
$F_{f}^{A}$
is necessary only at the first cell faces opposite to
$\unicode[STIX]{x1D6F4}$
,
$f_{i}^{\unicode[STIX]{x1D6F4},o}$
because the velocity normal to the interface in a moving reference frame is zero. Our aim would be to correct directly the concentrations with the prescribed value from the SGS model
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}$
. However, this cannot be done within an implicit framework, thus we correct the advective fluxes to match the prescribed SGS concentration. The numerical fluxes are computed as
$$\begin{eqnarray}F_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{A,num}=c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}},\end{eqnarray}$$
where
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}$
is the concentration value interpolated to the face centre and
$\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}$
is a modified advective flux.
The species fluxes computed with the SGS face value are
$$\begin{eqnarray}F_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{A,SGS}=c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}\end{eqnarray}$$
where
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}$
is provided by the SGS model. Enforcing the SGS fluxes to be equal to the numerical ones
$$\begin{eqnarray}F_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{A,num}\overset{!}{=}F_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{A,SGS},\end{eqnarray}$$
we get the equality
$$\begin{eqnarray}c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}\overset{!}{=}c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}\end{eqnarray}$$
from which we compute the corrected advective fluxes
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}=\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}\frac{c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}}{c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}}.\end{eqnarray}$$
Also for the advective term, to simplify the notation, we will address
$\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}$
as
$\unicode[STIX]{x1D719}^{SGS}$
, where
$\unicode[STIX]{x1D719}^{SGS}$
contains the modified local values from the SGS model in the required faces. For the other faces the original numerical fluxes are kept. Note that
$\unicode[STIX]{x1D719}^{SGS}=\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}~c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}/c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}$
. Thus, if
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}$
and
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}$
are interpolated with the same scheme, the modification of the advective term at the interested faces translates into enforcing the
$c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}$
; in fact
$\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}(c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}/c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num})c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}=\unicode[STIX]{x1D719}_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{SGS}$
. This also assures that our method remains conservative.
The advection correction via the SGS model is applied only if the concentration profile in the first three cell layers close to the interface is monotonic, see § B.3 for more details on exception handling. This condition is fundamental to avoid non-physical (unbounded) concentrations; see Versteeg & Malalasekera (Reference Versteeg and Malalasekera1995--2007).
B.2 Algorithm for the SGS model parameter calculation
In this section, the main steps to compute the SGS model parameter
$\unicode[STIX]{x1D6FF}$
are explained. We adopt an iterative approach, as described in Weiner & Bothe (Reference Weiner and Bothe2017), to find the model parameter
$\unicode[STIX]{x1D6FF}$
that fulfils
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D702}}_{C}=\frac{\bar{c}-c_{|\unicode[STIX]{x1D6F4}}}{c_{\infty }-c_{|\unicode[STIX]{x1D6F4}}}\overset{!}{=}\frac{1}{V}\int _{V}\unicode[STIX]{x1D702}(x/\unicode[STIX]{x1D6FF})\,\text{d}V=\unicode[STIX]{x1D702}_{SGS},\end{eqnarray}$$
where
$\bar{\unicode[STIX]{x1D702}}_{C}$
is the volume averaged cell-centred value coming from the finite volume discretization, which has to be equal to the volume average computed with the SGS model. Above,
$\unicode[STIX]{x1D702}$
is given as
$$\begin{eqnarray}\unicode[STIX]{x1D702}(x,y)=\frac{c(x,y)-c_{|\unicode[STIX]{x1D6F4}}}{c_{\infty }-c_{|\unicode[STIX]{x1D6F4}}}=\text{erf}(x/\unicode[STIX]{x1D6FF}(y))\end{eqnarray}$$
according to (B 3). The quantity
$\bar{c}$
is the average concentration in an interface cell (
$c_{i}$
in figure 29),
$c_{|\unicode[STIX]{x1D6F4}}$
is the bulk concentration at the interface (
$c_{f_{i}^{\unicode[STIX]{x1D6F4}}}$
in figure 29). The iterative solution based on (B 17) requires the evaluation of the volume integral. Here only the main steps from Weiner & Bothe (Reference Weiner and Bothe2017) are reported. The iterative algorithm is based on the work of Ahn & Shashkov (Reference Ahn, Shashkov, Brewer and Marcum2008) and uses a combined Newton–Bisection method to search for
$\unicode[STIX]{x1D6FF}$
which converges very quickly, usually after three iterations. The maximum number of iterations is set to 10. As initial guess for
$\unicode[STIX]{x1D6FF}_{0}$
the first two terms of a series expansion for the inverse error function are taken, that is
$\unicode[STIX]{x1D6FF}_{0}=(l/2)/(0.5\unicode[STIX]{x03C0}(\unicode[STIX]{x1D702}_{c}+\unicode[STIX]{x03C0}/12\unicode[STIX]{x1D702}_{c}^{3}))$
, with
$l$
being the first cell thickness. Bounding values for
$\unicode[STIX]{x1D6FF}$
are taken equal to
$\unicode[STIX]{x1D6FF}_{min}=1\times 10^{-15}$
and
$\unicode[STIX]{x1D6FF}_{max}=10\unicode[STIX]{x1D6FF}_{0}$
. The convergence tolerance is set to
$\text{tol}=1\times 10^{-9}$
.
In each time step, there is an initialization step for the required parameters. The result of the iterative procedure will be a vector containing all the
$\unicode[STIX]{x1D6FF}$
values (for all the interface cells). The algorithm is displayed as pseudo-code in Algorithm 1. Note that the formula to compute the residual has been corrected with respect to Weiner & Bothe (Reference Weiner and Bothe2017).
Exception handling. Before the iterative procedure is started a check that the values of
$\unicode[STIX]{x1D702}_{c}$
are between 0 and 1 is done. If the maximum number of iterations is reached without a converged value for
$\unicode[STIX]{x1D6FF}$
or if the computed
$\unicode[STIX]{x1D6FF}$
is larger than the first cell thickness by a factor of 1000, then
$\unicode[STIX]{x1D6FF}$
is set to
$-1$
and the SGS correction will not be applied at the corresponding face.
B.3 Correction of diffusive and advective fluxes within the SGS modelling
After the iterative computation of the model parameter
$\unicode[STIX]{x1D6FF}$
, the SGS correction is applied to the diffusive and advective fluxes as explained in § 3.2.3. The various steps for the flux correction are reported in Algorithm 2.
Exception handling. The diffusive and advective fluxes are corrected only if the iterative procedure to compute
$\unicode[STIX]{x1D6FF}$
converged, that is
$\unicode[STIX]{x1D6FF}>0$
. Furthermore, a check that the gradient and the concentration close to the interface are non-zero is included,
$|(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{num}|>10^{-15}$
and
$|c_{f_{i}^{\unicode[STIX]{x1D6F4},o}}^{num}|>10^{-15}$
. If these checks fail, the standard discretization is used.
An additional exception handling is implemented specifically for the correction of the diffusive fluxes at the second layer of faces
$f_{i}^{\unicode[STIX]{x1D6F4},o}$
. The SGS correction is applied only if the ratio between the SGS gradient and the numerical one is smaller than unity,
$|(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{SGS}/(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\ast }}^{num}|<1$
. If the correction factor is larger than one, the SGS model application is not necessary and the diffusivity will not be corrected at the respective face.
The last exception regards the correction of the advective fluxes. The SGS model correction is applied only if the concentration profile within the first three cells close to the interface is monotonic. If we number the cell centres from the interface outwards as
$c_{1},c_{2},c_{3}$
, then the SGS correction is applied only if
$(c_{1}-c_{2})(c_{2}-c_{3})>0$
.
Appendix C. Validation of the SGS model for species transfer
To validate the solution of the species transfer problem with SGS modelling, four test cases with increasing complexity are presented. The local Sherwood number
$Sh_{loc}$
is used for comparison with the reference solution.
Figure 30. SGS 2-D model problem set-up.
C.1 Two-dimensional model problem
This test case refers directly to the simplified problem formulation on which the SGS model is based. The implementation of the SGS model has been validated against the analytical solution taken from Weiner & Bothe (Reference Weiner and Bothe2017) and reported in § B.1. The problem set-up under investigation is sketched in figure 30. All the simplifying assumptions of the model problem are fulfilled if the computational domain size is large enough. The distance between the interface and the boundaries in the
$x$
-direction is approximately 50 times the maximum species boundary layer thickness, to ensure that the presence of a finite domain is negligible. The presence of the gas phase is modelled via the boundary condition for the species concentration at
$\unicode[STIX]{x1D6F4}$
. The boundary and initial conditions can be found in figure 30. Four different mesh resolutions are considered from 5 to
$40~\unicode[STIX]{x03BC}\text{m}$
. As we are interested in advection-dominated problems, a high Péclet number of
$Pe=10^{5}$
is chosen. The local Sherwood number is computed as
$$\begin{eqnarray}Sh_{loc}(y_{i})=(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\unicode[STIX]{x1D6F4}}}\frac{L_{y}}{c_{i|\unicode[STIX]{x1D6F4}}-c_{\infty }}\end{eqnarray}$$
with the normal derivative at the interface
$(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\unicode[STIX]{x1D6F4}}}$
, the concentration in the boundary cell centre
$c_{i|\unicode[STIX]{x1D6F4}}$
and the species concentration far away from the interface
$c_{\infty }$
. Without the SGS model the gradient is computed as
$(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\unicode[STIX]{x1D6F4}}}=(c_{i|\unicode[STIX]{x1D6F4}}-c_{f_{i}^{\unicode[STIX]{x1D6F4}}})/d_{i}$
, where
$d_{i}$
is the distance between the boundary face centre and the boundary cell centre, and
$c_{f_{i}^{\unicode[STIX]{x1D6F4}}}$
is the concentration at the interface face; otherwise
$(\unicode[STIX]{x2202}_{n}c)_{f_{i}^{\unicode[STIX]{x1D6F4}}}^{SGS}$
is used.
Figure 31 depicts the comparison between the analytical solution and the numerical results obtained with and without the SGS model. When the problem is solved with linear interpolation, the relatively coarse meshes are not able to predict the solution precisely. The finest mesh (
$5~\unicode[STIX]{x03BC}\text{m}$
) provides a good approximation of the local Sherwood number except for the region close to the inlet. All the cases where the SGS model is applied are in very good agreement with the reference solution. The enlarged view in figure 31 shows also mesh convergence for the SGS model results.
Figure 31. Local Sherwood number for the 1-D model problem.
C.2 Spherical bubbles at small Reynolds number
A spherical bubble at small Reynolds number is considered. For this case a semi-analytical solution of the species transport equation is possible. The velocity field is based on the solution of Satapathy & Smith (Reference Satapathy and Smith1960) (spherical particle of radius
$r_{b}$
rising in a larger sphere
$R$
). On top of this velocity field, the species transport equation can be solved numerically using a very high grid resolution (cell thickness
$l\approx 0.06~\unicode[STIX]{x03BC}\text{m}$
close to the interface). Four different molecular diffusivities are considered corresponding to Schmidt numbers of
$Sc=10^{4},10^{5},10^{6},10^{7}$
, where
$Sc$
is the ratio between viscous and molecular diffusion
$\unicode[STIX]{x1D708}/D$
. The bubble radius is
$r_{b}=1~\text{mm}$
and the Reynolds number is set to
$Re=0.56$
. The local Sherwood number
$Sh_{loc}(\unicode[STIX]{x1D703}_{i})$
is computed as in (C 1), where
$\unicode[STIX]{x1D703}_{i}$
is the polar angle, i.e. the angle following a streamline on the bubble surface from the top
$(\unicode[STIX]{x1D703}=0)$
to the bottom
$(\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0})$
. The bubble equivalent diameter
$d_{eq}$
is taken as reference length.
C.2.1 Axisymmetric species transfer with given velocity field
The species transport is solved on top of the velocity field provided by the solution of Satapathy and Smith for the different Schmidt numbers. The results obtained with the SGS model are compared to the mesh independent direct numerical solution. The set-up for these simulations is depicted in figure 32. The fluid properties for the liquid side (identified with a
$^{+}$
) can be found in table 6. Four different mesh resolutions are considered with a cell thickness
$l$
close to the interface ranging from 5 to
$40~\unicode[STIX]{x03BC}\text{m}$
. The four different diffusion coefficients are
$10^{-8},10^{-9},10^{-10}$
and
$10^{-11}~\text{m}^{2}~\text{s}^{-1}$
. The species concentration at the interface
$\unicode[STIX]{x1D6F4}$
is set to
$c_{|\unicode[STIX]{x1D6F4}}=1~\text{mol}~\text{m}^{-3}$
, while the initial bulk concentration in
$\unicode[STIX]{x1D6FA}$
is set to
$c_{0}=c_{\infty }=0$
.
Figure 32. Domain used to solve the species transport with the given analytical velocity field.
Figure 33. Local Sherwood number for the species transfer problem with given Satapathy–Smith velocity profile.
In figure 33 an overview of the results obtained by applying the SGS model compared to the reference solutions is reported. Figure 33 shows a very good agreement between the numerical results using the SGS modelling and the respective references for each tested Schmidt number. This test case shows that the two coarsest meshes (
$l=40,20~\unicode[STIX]{x03BC}\text{m}$
) are not fully capable of properly resolving the species transport for the highest Schmidt number, under-predicting the Sherwood number in the upper part of the bubble. Such behaviour has to be considered in the application case set-up with surfactant transport and sorption, mainly in the choice of the mesh resolution.
Figure 34. Local Sherwood number for the species transfer problem with given Satapathy–Smith velocity profile. Black symbols: with SGS modelling, grey symbols: linear interpolation.
Table 6. Fluid properties for the Satapathy–Smith case.
For completeness, in figure 34 the comparison between the cases with and without SGS modelling is reported. The results obtained applying the SGS model are coloured in black, while the ones obtained with a linear interpolation method are grey. Already for
$Sc=10^{5}$
the standard discretization is inadequate to correctly describe the species transfer close to the interface for the given mesh resolutions. This comparison confirms again that with the SGS model one can save several mesh refinement levels.
C.2.2 Species transfer with computed velocity field
The species transport problem from a rising bubble is considered. The full 3-D problem, hydrodynamics and species transfer, is solved within the interface-tracking framework, see the algorithm in figure 27. The case set-up follows the one described in § 4.1. The interface consists of polyhedral faces with an edge length of approximately
$50~\unicode[STIX]{x03BC}\text{m}$
and first cell layer thicknesses of
$l=12~\unicode[STIX]{x03BC}\text{m}$
and
$l=25~\unicode[STIX]{x03BC}\text{m}$
. The initial shape of the bubble is a sphere of radius
$r_{b}=1~\text{mm}$
. The bubble is positioned in the centre of a spherical domain of radius
$10r_{b}$
. The fact that the interface is deformable is not relevant to the Satapathy–Smith case, because due to the choice of the fluid properties, the bubble does not deform significantly.
The initial and boundary conditions for the transferred species are the same as for the semi-analytical solution. The fluids properties are given in table 6. For this test case, the smallest and the highest Schmidt numbers are considered, i.e.
$Sc=10^{4},10^{7}$
. As a reference, the semi-analytical solution presented in the former paragraph is used. The calculated velocity profile in the interface-tracking framework slightly differs from the Satapathy–Smith solution (less than 1.2 %, see Weber (Reference Weber2016, § 4.1.2)), because the latter is based on a Stokes flow. This small difference can have some impact on the concentration profile close to the interface. In figure 35 the results in terms of Sherwood number for the 3-D case are reported. As can be seen from the two graphs there is a good agreement between the reference solution and the numerical one employing the SGS model. As anticipated, the reference solution is computed based on the Satapathy–Smith velocity profile, thus, since we are dealing with highly nonlinear functions (species concentration close to
$\unicode[STIX]{x1D6F4}$
), small deviations in the velocity field could be enough to produce the observed discrepancies in the results. In figure 35 also the results without the SGS model are plotted. For small Schmidt numbers, figure 35(a), the standard discretization provides results in good agreement with the reference solution, while the Sherwood numbers resulting from the SGS modelling show a sensitivity to the mesh resolution. On the other hand, for high Schmidt numbers and the given mesh resolution, figure 35(b), the standard discretization provides underestimated Sherwood numbers, while the ones obtained with the SGS model are in good agreement with the reference.
Figure 35. Local Sherwood numbers for the species transfer problem with Satapathy–Smith set-up. Black symbols: with SGS modelling, grey symbols: linear interpolation.
C.3 Two-dimensional deformable bubbles at higher Reynolds number
As a final step to validate the SGS model, simulations of a 2-D bubble rising in contaminated water are performed with and without SGS modelling for the surfactant transport. This setting aims to demonstrate that the SGS model predicts the surfactant transfer well under dynamic conditions, e.g. when the bubble deforms, accelerates or decelerates or when the flow detaches and vortices form. As can be seen from figure 37, the results where the SGS modelling was employed are matching the mesh independent results obtained with standard interpolation.
For these tests, the intermediate initial concentration is used, i.e.
$c_{0}=0.008~\text{mol}~\text{m}^{-3}$
, with
$c_{0}^{\unicode[STIX]{x1D6F4}}=0~\text{mol}~\text{m}^{-2}$
. Different bulk diffusivities,
$D=5\times 10^{-7},5\times 10^{-8},5\times 10^{-9},5\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$
, and mesh resolutions, first cell thicknesses
$l_{i}=16,12,8,3,1.7,1.2~\unicode[STIX]{x03BC}\text{m}$
, are considered. The changes between
$l_{i}=1.7~\unicode[STIX]{x03BC}\text{m}$
and
$l_{i}=1.2~\unicode[STIX]{x03BC}\text{m}$
in rise velocity and surfactant transport are always less than 1.15 %. Therefore we consider the results on the finest mesh employing standard discretization as mesh independent and use them as reference solution (lines in the plots). The results for mesh resolutions with first cell thicknesses equal to 16, 8, 3 and
$1.2~\unicode[STIX]{x03BC}\text{m}$
are selected for the plots below. The results for higher surfactant bulk diffusivities are not depicted in figure 37 because they look qualitatively similar. In fact, if the mesh resolution is sufficient, then all the results lay on the reference curve.
Figure 37(a) shows that the rise velocities obtained applying the SGS modelling are all in agreement with the reference. On the other hand, the results obtained with standard interpolation follow a very different trend. Only the
$3~\unicode[STIX]{x03BC}\text{m}$
mesh gets close to the reference. Not only the rise velocities are in good agreement with the reference, but also the total amount of surfactant on the interface, as shown in figure 37(b) for different diffusion coefficients. As can be seen from the graph, the results obtained with the SGS modelling are all laying on the reference curves, while for
$D=5\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$
, only the
$3~\unicode[STIX]{x03BC}\text{m}$
case with standard interpolation tends to the correct result. So far we considered global quantities for comparison. Further confirmation that the SGS modelling is performing well and corresponding to the standard interpolation results is given by the local Sherwood numbers for the different diffusivities at
$t=0.2~\text{s}$
, see figure 37(c). Here it can be seen that all the cases where the SGS model has been used deliver a very good approximation of the local Sherwood number. Moreover, the shape of the local Sherwood number profile reflects the flow field around the bubble, see figure 36.
Figure 36. Flow field around and inside the rising bubble. The bubble surface is coloured by the surfactant concentration;
$t=0.2~\text{s}$
.
Figure 37. Simulation results for a 2-D bubble rising in contaminated water; comparison between cases with and without SGS modelling. Black symbols: with SGS modelling, grey symbols: linear interpolation.
Appendix D. Computation of the forces acting on the interface
The starting point for the derivation of the expression to compute the different forces acting on the interface is the interfacial jump condition (2.5), namely
$$\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}$$
Equation (D 1) can be decomposed into normal and tangential components to the interface, in order to better understand the meaning of each force acting on the bubble. With some notions of tensor calculus and knowing the definition of the surface operators, the projection of the jump condition in the direction normal
$(\boldsymbol{n}_{\unicode[STIX]{x1D6F4}})$
to the interface reads
$$\begin{eqnarray}\unicode[STIX]{x27E6}p_{tot}\unicode[STIX]{x27E7}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}+2\unicode[STIX]{x27E6}\unicode[STIX]{x1D707}\unicode[STIX]{x27E7}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}},\end{eqnarray}$$
while in the direction tangential
$(\boldsymbol{t}_{\unicode[STIX]{x1D6F4}})$
to
$\unicode[STIX]{x1D6F4}$
we obtain
$$\begin{eqnarray}-\unicode[STIX]{x27E6}\unicode[STIX]{x1D707}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}\unicode[STIX]{x27E7}-\unicode[STIX]{x27E6}\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{v})\boldsymbol{\cdot }\boldsymbol{n}\unicode[STIX]{x27E7}-\unicode[STIX]{x27E6}\unicode[STIX]{x1D707}\unicode[STIX]{x27E7}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})=\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D70E}.\end{eqnarray}$$
As before, if we indicate with
$A$
the liquid side and with
$B$
the gas side, we can specify all the terms in the jump brackets as follows,
$$\begin{eqnarray}p_{tot,B}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}-p_{tot,A}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}+2\unicode[STIX]{x1D707}_{B}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}-2\unicode[STIX]{x1D707}_{A}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}},\end{eqnarray}$$
for the normal direction, and
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle -\unicode[STIX]{x1D707}_{B}[(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v})_{B}-(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{v})_{B}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}+\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})]\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\unicode[STIX]{x1D707}_{A}[(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v})_{A}+(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{v})_{A}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}+\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})]=\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D70E},\end{eqnarray}$$
for the tangential direction. Note that the interface normal is
$\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$
with
$\boldsymbol{n}_{A}=\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$
, while
$\boldsymbol{n}_{B}=-\boldsymbol{n}_{\unicode[STIX]{x1D6F4}}$
. Each term in (D 4) and (D 5), when multiplied by the face area will give a force contribution.
(I) Marangoni force
(i) area specific force at face
$i\in \unicode[STIX]{x1D6F4}$
(D 6)$$\begin{eqnarray}\boldsymbol{f}_{i}^{ma}=\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\unicode[STIX]{x1D70E}_{i}\end{eqnarray}$$
(ii) resultant force on
$\unicode[STIX]{x1D6F4}$
(D 7)where$$\begin{eqnarray}\boldsymbol{F}^{ma}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{i}^{ma}S_{f_{i}},\end{eqnarray}$$
$N_{f}$
is the number of faces on the interface and
$S_{f_{i}}$
the face area.
(II) Capillary pressure force
(i) area specific force at face
$i\in \unicode[STIX]{x1D6F4}$
(D 8)$$\begin{eqnarray}\boldsymbol{f}_{i}^{ca}=\unicode[STIX]{x1D70E}_{i}k_{i}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}\end{eqnarray}$$
(ii) resultant force on
$\unicode[STIX]{x1D6F4}$
(D 9)$$\begin{eqnarray}\boldsymbol{F}^{ca}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{i}^{ca}S_{f_{i}}.\end{eqnarray}$$
(III) Total pressure force jump
(i) area specific force at face
$i\in \unicode[STIX]{x1D6F4}$
(D 10)$$\begin{eqnarray}\boldsymbol{f}_{i}^{p_{tot}}=(p_{tot,B_{i}}-p_{tot,A_{i}})\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}\end{eqnarray}$$
(ii) resultant force on
$\unicode[STIX]{x1D6F4}$
(D 11)$$\begin{eqnarray}\boldsymbol{F}^{p_{tot}}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{i}^{p_{tot}}S_{f_{i}}.\end{eqnarray}$$
(IV) Dynamic pressure force jump
(i) area specific force at face
$i\in \unicode[STIX]{x1D6F4}$
(D 12)where the dynamic pressure is computed as$$\begin{eqnarray}\boldsymbol{f}_{i}^{p_{dyn}}=(p_{dyn,B_{i}}-p_{dyn,A_{i}})\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}},\end{eqnarray}$$
$p_{dyn}=p_{tot}-p_{hydro}$
, with the hydrostatic pressure
$p_{hydro}=\unicode[STIX]{x1D70C}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x}_{f_{i}}$
(ii) resultant force on
$\unicode[STIX]{x1D6F4}$
(D 13)$$\begin{eqnarray}\boldsymbol{F}^{p_{dyn}}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{i}^{p_{dyn}}S_{f_{i}}.\end{eqnarray}$$
(V) Normal viscous force
(i) area specific forces at face
$i\in \unicode[STIX]{x1D6F4}$
(D 14)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\bot ,B_{i}}^{visc}=2\unicode[STIX]{x1D707}_{B}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})_{i}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}} & \displaystyle\end{eqnarray}$$
(D 15)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\bot ,A_{i}}^{visc}=-2\unicode[STIX]{x1D707}_{A}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})_{i}\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}} & \displaystyle\end{eqnarray}$$
(D 16)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\bot ,i}^{visc}=\boldsymbol{f}_{\bot ,B_{i}}^{visc}+\boldsymbol{f}_{\bot ,A_{i}}^{visc} & \displaystyle\end{eqnarray}$$
(ii) resultant forces on
$\unicode[STIX]{x1D6F4}$
(D 17)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\bot ,B}^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\bot ,B_{i}}^{visc}S_{f_{i}} & \displaystyle\end{eqnarray}$$
(D 18)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\bot ,A}^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\bot ,A_{i}}^{visc}S_{f_{i}} & \displaystyle\end{eqnarray}$$
(D 19)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\bot }^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\bot ,i}^{visc}S_{f_{i}}. & \displaystyle\end{eqnarray}$$
(VI) Tangential viscous force
(i) area specific forces at face
$i\in \unicode[STIX]{x1D6F4}$
(D 20)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\Vert ,B_{i}}^{visc}=\unicode[STIX]{x1D707}_{B}[-(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v})_{B_{i}}+(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{v})_{B_{i}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}-\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})_{i}] & \displaystyle\end{eqnarray}$$
(D 21)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\Vert ,A_{i}}^{visc}=\unicode[STIX]{x1D707}_{A}[(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v})_{A_{i}}+(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{v})_{A_{i}}\boldsymbol{\cdot }\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}+\boldsymbol{n}_{\unicode[STIX]{x1D6F4}_{i}}(\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D6F4}}\boldsymbol{\cdot }\boldsymbol{v})_{i}] & \displaystyle\end{eqnarray}$$
(D 22)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{\Vert ,i}^{visc}=\boldsymbol{f}_{\Vert ,B_{i}}^{visc}+\boldsymbol{f}_{\Vert ,A_{i}}^{visc} & \displaystyle\end{eqnarray}$$
(ii) resultant forces on
$\unicode[STIX]{x1D6F4}$
(D 23)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\Vert ,B}^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\Vert ,B_{i}}^{visc}S_{f_{i}} & \displaystyle\end{eqnarray}$$
(D 24)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\Vert ,A}^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\Vert ,A_{i}}^{visc}S_{f_{i}} & \displaystyle\end{eqnarray}$$
(D 25)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{\Vert }^{visc}=\mathop{\sum }_{i}^{N_{f}}\boldsymbol{f}_{\Vert ,i}^{visc}S_{f_{i}}. & \displaystyle\end{eqnarray}$$
If we write the jump condition in terms of global forces then we obtain the following expression
$$\begin{eqnarray}\boldsymbol{F}^{p_{tot}}+\boldsymbol{F}^{visc}=\boldsymbol{F}^{ca}+\boldsymbol{F}^{ma}\end{eqnarray}$$
that can serve as a check of the fulfilment of the jump condition at the interface at the end of the simulation.
