1. Introduction
Bubbles oscillating in a fluid are used in numerous technological devices and processes, including microfluidic devices, chemical processes, emulsification, oil recovery, cleaning, drug delivery and flotation separation. In many of those processes, surfactants are introduced into the liquid. In addition to the technological applications, the analysis of the influence of surfactants on the frequency and damping of droplet and bubble oscillations provides tools for experimental measurements of rheological properties of an interface covered by a surfactant (Lu & Apfel Reference Lu and Apfel1990; Johnson & Stebe Reference Johnson and Stebe1994; Tian, Holt & Apfel Reference Tian, Holt and Apfel1995; Asaki & Marston Reference Asaki and Marston1997; Meier et al. Reference Meier, Greune, Mayboom and Hofmann2000; Freer, Wong & Radke Reference Freer, Wong and Radke2005; Ravera, Loglio & Kovalchuk Reference Ravera, Loglio and Kovalchuk2010; Lalanne, Masbernat & Risso Reference Lalanne, Masbernat and Risso2020).
The dynamics of the gas–liquid interface is strongly influenced by the presence of a surfactant (Levich Reference Levich1962). That influence is twofold. On one side, the surfactant locally diminishes the surface tension. Therefore, the non-uniformity of the surface concentration of the surfactant caused by an interfacial flow creates a tangential Marangoni stress that leads to an effective elasticity of the interface between the bubble and the surrounding liquid (the Gibbs elasticity). The mathematical description of the surfactant redistribution on a deformable surface taking into account advection, diffusion and the concentration change due to the surface deformation was formulated by Stone (Reference Stone1990) and Wong, Rumschitzki & Maldarelli (Reference Wong, Rumschitzki and Maldarelli1996). On another side, an interface containing a surfactant is subject to additional dissipation mechanisms known as the shear and dilatational interfacial viscosities (Scriven Reference Scriven1960).
While the radial bubble oscillations under the action of surfactants are relatively well explored (Glazman Reference Glazman1983; Karapantsios & Kostoglou Reference Karapaitsios and Kostoglou1999; Zhong & Ardekani Reference Zhong and Ardekani2022), the influence of surfactant on shape-changing oscillations has not yet been fully understood. In the absence of the surfactant, the shape-changing oscillations have been considered by Reid (Reference Reid1960) and Prosperetti (Reference Prosperetti1980) for both gas bubbles and liquid droplets in a similar way (because for that kind of oscillations the compressibility effects are of minor importance).
Estimates made for a millimetre air bubble in water show that the natural frequency of the quadrupole mode of capillary oscillations (Lamb Reference Lamb1932) of the bubble (approximately
$420\rm \:Hz$
) is an order of magnitude lower than the frequencies of monopole oscillations (Plesset Reference Plesset1949; Rayleigh Reference Rayleigh1879) of the bubble volume (approximately
$6.6\rm \:kHz$
). And with the growth of the bubble size, this difference increases. This should obviously increase the accuracy of the experiment aimed at observing precisely the oscillations of the shape of such a spherical interface. Moreover, the influence of surfactants further reduces, through the surface tension factor, the frequencies of capillary oscillations, having a relatively weak effect on volume oscillations.
The influence of the surfactant on the shape-changing oscillations of an interface between two fluids was the subject of the investigation by Miller & Scriven (Reference Miller and Scriven1968), Sparling & Sedlak (Reference Sparling and Sedlak1989), Lu & Apfel (Reference Lu and Apfel1991) and Tian et al. (Reference Tian, Holt and Apfel1995) (for a liquid–gas interface). A systematic investigation of small-amplitude shape oscillations of a liquid droplet in a gas in the presence of an insoluble surfactant on the interface was carried out by Lyubimov et al. (Reference Lyubimov, Konovalov, Lyubimova and Egry2011). The dispersion relation for the oscillation frequency that describes the transition from oscillatory to monotonic decay of disturbances caused by interfacial viscosities was derived.
The finite-amplitude oscillations of droplets and bubbles covered by a surfactant can be studied using numerical methods tracking the droplet and bubble shape change. Among them are the level-set method (Xu & Zhao Reference Xu and Zhao2003; Reusken & Zhang Reference Reusken and Zhang2013; de Langavant et al. Reference de Langavant, Guittet, Theillard, Temprano-Coleto and Gibou2017; Piedfert et al. Reference Piedfert, Lalanne, Masbernat and Risso2018), volume of fluid method (Drumright-Clarke & Renardy Reference Drumright-Clarke and Renardy2004; James & Lowengrub Reference James and Lowengrub2004) and the combination of the boundary element method with the finite volume method (Bazlekov, Anderson & Meijer Reference Bazlekov, Anderson and Meijer2004; Bazlekov, Anderson & Meijer Reference Bazlekov, Anderson and Meijer2006). A special Lagrangian–Eulerian method was suggested by Luo, Shang & Bai (Reference Luo, Shang and Bai2019).
The goal of the present paper is the analytical and numerical investigation of gas bubble oscillations within a viscous liquid in the presence of an insoluble surfactant. The finite viscosity and density of the gas in the bubble are taken into account. A novel numerical approach for the computation of finite-amplitude bubble oscillations is suggested. The physical formulation of the problem and the characteristic non-dimensional parameters are presented in § 2. In § 3, the mathematical formulation of the problem for small-amplitude oscillations is given. The results of the analysis of the dispersion relation that determines the frequency and damping of oscillations are discussed. Special attention is devoted to the case of small but finite viscosities of the liquid and the gas. The influence of the Gibbs elasticity, surface diffusion and interfacial viscosities is investigated in detail. In § 4, the results of numerical simulations of the full nonlinear system of equations are presented and compared with the predictions of the linear theory. Section 5 contains concluding remarks.
2. Problem formulation
We consider a gas bubble that is spherical at equilibrium with radius
$R_{0}$
, suspended in an unbounded liquid that remains motionless far from the bubble. It is assumed that the force of gravity is balanced by other forces, such as the force of acoustic radiation, therefore, the bubble does not rise. The influence of the bubble rising in a viscous liquid (Bozzano & Dente Reference Bozzano and Dente2009) on its natural frequencies and damping rates (Lalanne, Tanguy & Risso Reference Lalanne, Tanguy and Risso2013) is an important problem, but its consideration is beyond the scope of this paper. The action of gravity on the droplet shape is neglected, because the capillary constant is much larger than the radius of the bubble. The compressibility of the gas inside the bubble is neglected, which is valid when the condition
$\lambda _{2}=2\pi c_{2}/\omega \gg R_{0}$
is satisfied, where
$\lambda _{2}$
represents the acoustic wavelength in the gas calculated using the speed of sound
$c_{2}$
and the cyclic frequency of the system’s natural oscillations
$\omega$
.
The material parameters of the media include the liquid density
$\rho _{1}$
and its kinematic viscosity coefficient
$\nu _{1}$
, as well as the corresponding parameters for the gas
$\rho _{2}$
and
$\nu _{2}$
. Here and subsequently, the index
$j$
refers to either the external liquid (
$j=1$
) or the gas inside the bubble (
$j=2$
).
An insoluble surfactant is adsorbed at the gas–liquid interface, with the surface tension coefficient
$\gamma$
assumed to vary linearly with small deviations
$\varGamma '$
of the surfactant surface concentration from its equilibrium value
$\varGamma _{0}$
,
\begin{equation} \gamma =\gamma _{0}+\left(\frac{\mathrm{d}\gamma}{\mathrm{d}\Gamma}\right)_{\Gamma=\Gamma_{0}}{\varGamma '}, \end{equation}
where
$\gamma _{0}$
denotes the equilibrium value of the surface tension coefficient at
$\varGamma '=0$
(or
$\varGamma =\varGamma _{0}$
).
The motion of viscous incompressible media is described by the Navier–Stokes equations
\begin{equation} \frac {\partial \boldsymbol{u}_{\!j}}{\partial t}+\left (\boldsymbol{u}_{\!j}\boldsymbol{\cdot }\boldsymbol{\nabla }\right )\boldsymbol{u}_{\!j}=-\frac {\boldsymbol{\nabla }p_{\!j}}{\rho _{\!j}}+\nu _{\!j}\Delta \boldsymbol{u}_{\!j} \end{equation}
and the continuity equations
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\!j} = 0, \end{equation}
where
$t$
is the time, while
$\boldsymbol{u}_{\!j}$
and
$p_{\!j}$
are the velocity and pressure fields in the media, respectively.
The surfactant redistribution along the surface is described by the convective diffusion equation that accounts for the free surface deformation (Stone Reference Stone1990; Wong et al. Reference Wong, Rumschitzki and Maldarelli1996),
\begin{equation} \frac {\partial \varGamma }{\partial t}+\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\left (\varGamma \boldsymbol{u}_{s}\right )+\varGamma \boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{n}\left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\right )-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!s}\varGamma =D_{s}\boldsymbol{\nabla} _{\!s}^{2}\varGamma . \end{equation}
Here
$\varGamma$
is the surface surfactant concentration,
$D_{s}$
is the surface diffusion coefficient,
$\boldsymbol{\nabla} _{\!s}= (\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n} )\boldsymbol{\cdot }\boldsymbol{\nabla}$
is the surface gradient operator,
$\boldsymbol{u}_{s}= (\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n} )\boldsymbol{\cdot }\boldsymbol{u}$
, where
$\boldsymbol{u}$
is the interface velocity,
$\boldsymbol{I}$
is the identity matrix and
$\boldsymbol{n}$
is the external normal to the liquid boundary. The partial derivative of
$\varGamma$
with respect to time is calculated in the normal direction. The second term on the left-hand side describes surfactant convection along the surface, while the third and fourth terms account for changes in surface concentration due to local changes in boundary area.
In the depths of its volume, the liquid remains at rest and
\begin{equation} \boldsymbol{u}_{1}\to 0. \end{equation}
At the interface described by equation
$G ( \boldsymbol{r}, t ) = 0$
, where
$\boldsymbol{r}$
is the radius vector, the kinematic condition
\begin{equation} \frac {{\partial G}}{{\partial t}} + \boldsymbol{u}_{1} \boldsymbol{\cdot }\boldsymbol{\nabla }G = 0 \end{equation}
and continuity condition of velocity
\begin{equation} \boldsymbol{u}_{1}=\boldsymbol{u}_{2} \end{equation}
are set. The stress boundary conditions on the free surface were written as in Nadim (Reference Nadim1996),
\begin{equation} \boldsymbol{n}\boldsymbol{\cdot }\left (\boldsymbol{\varPi }_{\!1}-\boldsymbol{\varPi }_{\!2}\right )=-\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{\varPi }_{\!s}, \end{equation}
where
$\boldsymbol{\varPi }_{\!1}$
and
$\boldsymbol{\varPi }_{\!2}$
are the stress tensors in the liquid and gas phases, respectively, and
$\boldsymbol{\varPi }_{\!s}$
is the surface tensor describing the rheological properties of the interface according to the Boussinesq–Scriven model (Scriven Reference Scriven1960),
\begin{equation} \boldsymbol{\varPi }_{\!s}=\gamma \boldsymbol{I}_{\!s}+2\eta _{s}\boldsymbol{E}_{s}+\eta _{d}\boldsymbol{I}_{\!s}\left (\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{u}\right )\!, \end{equation}
where
$\eta _{s}$
and
$\eta _{d}$
are the shear and dilatational surface viscosities,
$\boldsymbol{I}_{\!s}= (\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n} )$
is the surface unit tensor and
$\boldsymbol{E}_{s}$
is the surface rate-of-strain tensor defined as
\begin{equation} \boldsymbol{E}_{s}=\frac {1}{2}\big \{ \boldsymbol{\nabla} _{\!s}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{I}_{\!s}+\boldsymbol{I}_{\!s}\boldsymbol{\cdot }\left (\boldsymbol{\nabla} _{\!s}\boldsymbol{u}\right )^{T}\big \} -\frac {1}{2}\boldsymbol{I}_{\!s}\left (\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{u}\right )\!. \end{equation}
To express the problem below in a dimensionless form, we choose the equilibrium bubble radius
$R_{0}$
as the unit of length, the capillary time
$\tau _{0}=\sqrt {(\rho _{1}+\rho _{2})R_{0}^{3}/\gamma _{0}}$
as the unit of time, the capillary pressure
$\gamma _{0}/R_{0}$
as the unit of pressure and the equilibrium surfactant concentration
$\varGamma _{0}$
as the unit of surface concentration.
The governing dimensionless parameters include the dimensionless densities of the media,
\begin{equation} \tilde {\rho }_{\!j}=\frac {\rho _{\!j}}{\rho _{1}+\rho _{2}}, \end{equation}
which satisfy the relation
$\tilde {\rho }_{1}+\tilde {\rho }_{2}=1$
. The dissipative parameters for the liquid and gas are defined as
\begin{equation} \delta _{1}=\frac {\nu _{1}\tau _{0}}{R_{0}^{2}},\quad \delta _{2}=\frac {\nu _{2}\tau _{0}}{R_{0}^{2}}. \end{equation}
Additionally, the parameters that characterize the shear and dilatational surface viscosities, the dimensionless Gibbs elasticity and the surfactant surface diffusion are given by
\begin{align}& h_{s}=\frac {\eta _{s}\tau _{0}}{\left (\rho _{1}+\rho _{2}\right )R_{0}^{3}},\quad h_{d}=\frac {\eta _{d}\tau _{0}}{\left (\rho _{1}+\rho _{2}\right )R_{0}^{3}}, \end{align}
\begin{align}&\quad k=\left (\frac {{\rm d}\gamma }{{\rm d}\varGamma }\right )_{\varGamma =\varGamma _{0}}\frac {\varGamma _{0}}{\gamma _{0}},\quad d_{s}=\frac {D_{s}\tau _{0}}{R_{0}^{2}}. \end{align}
The dissipative parameters
$\delta _{\!j}$
are defined as the inverse Reynolds numbers and specify the ratio of viscous forces to inertial forces. In the same way, the parameters of shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities give similar ratios, but for surface viscous stresses. As noted by Tian et al. (Reference Tian, Holt and Apfel1995), the dilatational surface viscosity is usually much larger than the shear surface viscosity (
$\eta _{d} \gg \eta _{s}$
) for a soluble surfactant (they can be compatible for insoluble surfactants). According to Karapantsios & Kostoglou (Reference Karapaitsios and Kostoglou1999) for a monolayer of stearic acid
$\eta _{d} = 1.5 \times 10^{ - 3}\,\rm N\: s\,m^{-1}$
. Then for the millimetre-sized air bubbles of interest to us in water,
$\delta _{1} \sim 0.005$
and
$h_{d} \sim 15$
. Note that the parameters of shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities increase with decreasing bubble radius
$R_{0}$
faster than the dissipative parameters
$\delta _{\!j}$
. Thus, for small gas bubbles, the role of surface viscosities increases. As mentioned in Lu & Apfel (Reference Lu and Apfel1990), the surface diffusion coefficient is
$D_{s} \sim 10^{ - 7}\,\rm m^{2}\,s^{-1}$
. Therefore, for millimetre bubbles
$d_{s} \sim 5 \times 10^{ - 4}$
.
It can be shown that, when the dependence of surface tension on surfactant concentration follows Henry’s linear law,
\begin{equation} \gamma =\gamma _{\textit{pure}}-R_{g}T\varGamma , \end{equation}
then
\begin{equation} k=-\textit{Ma}, \end{equation}
where the Marangoni number
\begin{equation} Ma=\frac {R_{g}T\varGamma _{0}}{\gamma _{0}} \end{equation}
is expressed in terms of the universal gas constant
$R_{g}$
, the absolute temperature
$T$
and the equilibrium surface tension coefficient
$\gamma _{0}$
, related to the surface tension coefficient in the absence of surfactant
$\gamma _{\textit{pure}}$
by
\begin{equation} \gamma _{0}=\gamma _{\textit{pure}}-R_{g}T\varGamma _{0}. \end{equation}
3. Small-amplitude oscillations of the system
3.1. Governing equations and boundary conditions
Let us first consider small-amplitude oscillations of the system. In this case, the flow induced by oscillations is described by linearized Navier–Stokes equations for viscous incompressible media,
\begin{equation} \frac {\partial \boldsymbol{u}_{\!j}}{\partial t}=-\frac {\boldsymbol{\nabla }p_{\!j}}{\tilde {\rho }_{\!j}}+\delta _{\!j}\Delta \boldsymbol{u}_{\!j}, \end{equation}
and the continuity equations (2.3).
The origin
$\mathrm{O}$
of the spherical coordinate system
$\{\mathrm{O};r,\vartheta ,\phi \}$
is placed at the centre of the equilibrium bubble (figure 1). It is known that, in the framework of the linearized problem (Miller & Scriven Reference Miller and Scriven1968), the natural frequencies and damping rates of the non-symmetric solution with a non-zero azimuthal number for the spherical harmonic coincide with those obtained for the axisymmetric case. Therefore, for the sake of simplicity, we further assume that the azimuthal component of the velocity is absent (
$u_{\phi } = 0$
) and the derivatives in the corresponding direction are equal to zero (
$\partial /\partial \phi \equiv 0$
).
Figure 1. Configuration and parameters of the problem.
Linearized boundary conditions are written at the unperturbed interface
$r=1$
. These are the kinematic boundary condition (3.2), the continuity conditions of the velocity field (3.3) and (3.4), and the balance conditions of the tangential (3.5) and normal (3.6) stresses, represented taking into account additional terms with the Marangoni force, as well as the shear and dilatational surface viscosities (Lu & Apfel Reference Lu and Apfel1991),
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \frac {\partial f}{\partial t}=u_{1r}, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad u_{1r}=u_{2r}, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad u_{1\vartheta }=u_{2\vartheta }, \end{align}
\begin{align}& \left [\rho \delta \left (\frac {1}{r}\frac {\partial u_{r}}{\partial \vartheta }+\frac {\partial u_{\vartheta }}{\partial r}-\frac {u_{\vartheta }}{r}\right )\right ]+\left (h_{s}+h_{d}\right )\frac {1}{r^{2}}\frac {\partial }{\partial \vartheta }\left (\frac {\partial u_{1\vartheta }}{\partial \vartheta }+u_{1\vartheta }\textrm {ctg}\vartheta \right ) \nonumber \\&\qquad\qquad\qquad\quad +2h_{s}\frac {u_{1\vartheta }}{r^{2}}+2h_{d}\frac {1}{r^{2}}\frac {\partial u_{1r}}{\partial \vartheta }+k\frac {1}{r}\frac {\partial {\varGamma '}}{\partial \vartheta }=0, \end{align}
\begin{align}&\qquad \left [p-2\rho \delta \frac {\partial u_{r}}{\partial r}\right ]+2h_{d}\frac {1}{r^{2}}\left (\frac {\partial u_{1\vartheta }}{\partial \vartheta }+u_{1\vartheta }\textrm {ctg}\vartheta +2u_{1r}\right )+k\frac {2}{r}{\varGamma '} \nonumber \\&\qquad\qquad\qquad\qquad\qquad\quad =2f+\frac {\partial ^{2}f}{\partial \vartheta ^{2}}+\textrm {ctg}\vartheta \frac {\partial f}{\partial \vartheta }. \end{align}
Here, square brackets denote the jump of the corresponding quantity when crossing from gas to liquid (a quantity related to the gas is subtracted from the quantity related to the liquid), while
$f$
and
$\varGamma '$
represent perturbations of the interface shape and the surfactant surface concentration, respectively, measured from the base state with the spherical bubble shape (
$r=1$
) and uniformly distributed surfactant. The transport of the latter is governed by the corresponding linearized equation
\begin{equation} \frac {\partial {\varGamma '}}{\partial t}+\frac {1}{r}\left (\frac {\partial u_{1\vartheta }}{\partial \vartheta }+u_{1\vartheta }\textrm {ctg}\vartheta \right )+2\frac {u_{1r}}{r}=d_{s}\frac {1}{r^{2}}\left (\frac {\partial ^{2}{\varGamma '}}{\partial \vartheta ^{2}}+\textrm {ctg}\vartheta \frac {\partial {\varGamma '}}{\partial \vartheta }\right )\!. \end{equation}
According to Chandrasekhar (Reference Chandrasekhar1959), the linear problem, (3.1) and (2.3), admits in a spherical coordinate system a solution of the toroidal type, axisymmetric in our case, corresponding to the deformation of the shape of the interface, and hence to its capillary oscillations. Under the condition of the solution being bounded at the centre of the bubble at
$r=0$
and decay at infinity at
$r=\infty$
, the velocity and pressure fields are found as
\begin{align}& u_{\!jr}=\frac {l(l+1)}{r^{2}}U_{\!j}(r)Y_{l}\left (\vartheta \right )\exp \left (-i\varOmega _{l}t\right )+\mathrm{c.c.}, \end{align}
\begin{align}&\,\,\,\, u_{\!j\vartheta }=\frac {1}{r}\frac {{\rm d}U_{\!j}}{{\rm d}r}\frac {\partial }{\partial \vartheta }Y_{l}\left (\vartheta \right )\exp \left (-i\varOmega _{l}t\right )+\mathrm{c.c.}, \end{align}
\begin{align}&\quad\,\, p_{\!j}=\varPi _{\!j}\left (r\right )Y_{l}\left (\vartheta \right )\exp \left (-i\varOmega _{l}t\right )+\mathrm{c.c.}, \end{align}
\begin{align}&\qquad\qquad\qquad \varPi _{1}=\frac {A_{1}}{r^{l+1}}, \end{align}
\begin{align}&\qquad\qquad\qquad \varPi _{2}=B_{1}r^{l}, \end{align}
\begin{align}&\quad\,\, U_{1}=\frac {i}{l}\frac {A_{1}}{\tilde {\rho }_{1}\varOmega _{l}}\frac {1}{r^{l}}+A_{2}r\frac {h_{l}^{\left (1\right )}\left (x_{1}r\right )}{h_{l}^{\left (1\right )}\left (x_{1}\right )}, \end{align}
\begin{align}&\,\,\,\, U_{2}=-\frac {i}{l+1}\frac {B_{1}}{\tilde {\rho }_{2}\varOmega _{l}}r^{l+1}+B_{2}r\frac {j_{l}\left (x_{2}r\right )}{j_{l}\left (x_{2}\right )}, \end{align}
where
$\varOmega _{l}$
is the dimensionless eigenfrequency for the harmonic with meridional number
$l$
, described by the spherical harmonic function
$Y_{l} (\vartheta )$
. Also,
\begin{equation} x_{1}=\sqrt {\frac {i\varOmega _{l}}{\delta _{1}}},\quad x_{2}=\sqrt {\frac {i\varOmega _{l}}{\delta _{2}}}, \end{equation}
and
$h_{l}^{(1)}$
,
$j_{l}$
are spherical Hankel and Bessel functions of the first kind, respectively.
It should be noted that, according to Prosperetti (Reference Prosperetti1980), the selection of bounded solutions (3.13) and (3.14) assumes the condition
$\mathrm{Re}\,\varOmega _{l} \geqslant 0$
for the real part of the eigenfrequency. Furthermore, by introducing a branch cut along the negative real axis in the complex plane to ensure an unambiguous definition of the square root, we guarantee that
$\mathrm{Im}\,x_{\!j} \geqslant 0$
. This approach maintains consistency in the mathematical treatment of the problem while ensuring physically meaningful solutions for the oscillation frequencies and damping characteristics of the system. The condition on the imaginary part of
$x_{\!j}$
is particularly important as it directly relates to the viscous dissipation mechanisms in both the liquid and gas phases, which fundamentally influence the bubble’s oscillation dynamics.
The deviation of the bubble surface from a spherical shape
$r=1$
and, accordingly, the disturbances of the surface concentration of the surfactant are represented as
\begin{align}& f=FY_{l}\left (\vartheta \right )\exp \left (-i\varOmega _{l}t\right )+\mathrm{c.c.}, \end{align}
\begin{align}& {\varGamma '}=\varLambda Y_{l}\left (\vartheta \right )\exp \left (-i\varOmega _{l}t\right )+\mathrm{c.c.} \end{align}
By substituting expressions (3.8)–(3.17) into relations (3.2)–(3.7) and performing appropriate transformations, we arrive at a system of linear algebraic equations for the amplitudes
$A_{1},A_{2},B_{1},B_{2}$
and
$F$
,
\begin{align}&\qquad\qquad\qquad\,\, \left (l+1\right )\frac {i}{\tilde {\rho }_{1}\varOmega _{l}}A_{1}+l\left (l+1\right )A_{2}+i\varOmega _{l}F=0, \end{align}
\begin{align}&\qquad\qquad\qquad\quad \frac {1}{l}\frac {i}{\tilde {\rho }_{1}\varOmega _{l}}A_{1}+A_{2}=-\frac {1}{l+1}\frac {i}{\tilde {\rho }_{2}\varOmega _{l}}B_{1}+B_{2}, \end{align}
\begin{align}&\qquad\qquad \frac {i}{\tilde {\rho }_{1}\varOmega _{l}}A_{1}-\left \{ Q\left (x_{1}\right )+2\right \} A_{2}=\frac {i}{\tilde {\rho }_{2}\varOmega _{l}}B_{1}-\left \{ P\left (x_{2}\right )+2\right \} B_{2}, \end{align}
\begin{align}& \left \{ M+2\left (l+1\right )H_{2}\right \} \frac {i}{\tilde {\rho }_{1}\varOmega _{l}}A_{1}-\left \{ i\tilde {\rho }_{1}\varOmega _{l}+MQ\left (x_{1}\right )-2\left (l-1\right )\left (l+2\right )H_{1}\right \} A_{2} \nonumber \\&\quad +\left \{ 2\left (l-1\right )\frac {i\delta _{2}}{\varOmega _{l}}B_{1}+i\tilde {\rho }_{2}\varOmega _{l}+2\tilde {\rho }_{2}\delta _{2}P\left (x_{2}\right )-2\left (l-1\right )\left (l+2\right )\tilde {\rho }_{2}\delta _{2}\right \} B_{2}=0, \end{align}
\begin{align}& \left \{ 1+2\left (l+1\right )\left (l+2\right )H_{2}\frac {i}{\tilde {\rho }_{1}\varOmega _{l}}\right \}\! A_{1}-2l\left (l+1\right )H_{2}Q\left (x_{1}\right )A_{2} -\left \{ 1+2l\left (l-1\right )\frac {i\delta _{2}}{\varOmega _{l}}\right \}\! B_{1}\nonumber \\&\quad +2l\left (l+1\right )\tilde {\rho }_{2}\delta _{2}P\left (x_{2}\right )B_{2}+\left (l-1\right )\left (l+2\right )F=0, \end{align}
where
\begin{align}&\quad M=l\left (l+1\right )\left (H_{2}-H_{1}\right )+2H_{1},\quad H_{1}=\tilde {\rho }_{1}\delta _{1}-h_{s},\nonumber \\&\qquad H_{2}=\tilde {\rho }_{1}\delta _{1}+h_{d}+\frac {k}{i\varOmega _{l}-l\left (l+1\right )d_{s}}, \end{align}
\begin{align}& Q\left (x_{1}\right )=l-1-\frac {x_{1}h_{l+1}^{\left (1\right )}\left (x_{1}\right )}{h_{l}^{\left (1\right )}\left (x_{1}\right )},\quad P\left (x_{2}\right )=l-1-\frac {x_{2}j_{l+1}\left (x_{2}\right )}{j_{l}\left (x_{2}\right )}. \end{align}
As in Prosperetti (Reference Prosperetti1980), the Hankel and Bessel function ratios in
$Q(x_{1})$
and
$P(x_{2})$
were computed using the following recurrence relations:
\begin{equation} Q\left (x_{1}\right )=l-1-\tilde {{H}}\left (l+1/2,x_{1}\right )\!,\quad P\left (x_{2}\right )=l-1-\frac {x_{2}^{2}}{\tilde {J}\left (l+3/2,x_{2}\right )}, \end{equation}
where
\begin{align} \tilde {{H}}\left (\nu ,x\right )&= \begin{cases} 1-ix, & \nu =1/2\\ 2\nu -\dfrac {x^{2}}{\tilde {{H}}\left (\nu -1,x\right )}, & \mathrm{otherwise}, \end{cases} \end{align}
\begin{align} \tilde {J}\left (\nu ,x\right )&= \begin{cases} x\cot x, & \nu =1/2\\ \dfrac {x^{2}}{2\left (\nu -1\right )-\tilde {J}\left (\nu -1,x\right )}, & \mathrm{otherwise}. \end{cases} \end{align}
The homogeneous system of linear algebraic equations (3.18)–(3.22) has non-trivial solutions when its determinant equals zero. This condition determines the complex eigenfrequencies of bubble oscillations
$\varOmega _{l}$
as functions of the problem parameters. If it is necessary to exclude the influence of a light internal medium, when
$\tilde {\rho }_{1}=1$
and
$\tilde {\rho }_{2}=0$
, then it is necessary to drop relations (3.19) and (3.20), and in the remaining equations (3.18), (3.21) and (3.22) set the amplitudes
$B_{1}$
and
$B_{2}$
equal to zero.
The eigenfrequencies were determined both analytically and numerically. In the latter case, the determinant of system (3.18)–(3.22) or the reduced system was found using Gaussian elimination with pivoting. The requirement of vanishing real and imaginary parts of the complex determinant leads to a system of two dispersion relations. Its solutions for the real part (oscillation frequency) and imaginary part (damping rate) of the eigenfrequency were determined numerically using the Newton–Raphson method. According to the remarks (Prosperetti Reference Prosperetti1980) on the spherical Hankel function, the presence of an external medium makes the number of possible solutions of the system of transcendental algebraic equations finite (this is infinite for the configuration of a drop in a vacuum when
$\tilde{\rho}_{2}=1$
and
$\tilde{\rho}_{1}=0$
). Among numerical solutions, those showing the least damping were selected. The obtained results are presented in the following section.
3.2. Results
3.2.1. Estimates for the case of small bulk viscosity
It is evident that small viscosities of the liquid and gas should not significantly affect the capillary oscillations of the system. Therefore, the main part of the eigenfrequency equals the Rayleigh frequency (Lamb Reference Lamb1932),
\begin{equation} \varOmega _{l0}=\sqrt {\frac {l(l^{2}-1)(l+2)}{l\tilde {\rho }_{1}+(l+1)\tilde {\rho }_{2}}}. \end{equation}
However, the correction
$\Delta \varOmega _{l}=\varOmega _{l}-\varOmega _{l0}$
may become noticeable at finite values of the surface viscosity parameters
$h_{s}$
and
$h_{d}$
, as well as in several other situations discussed below.
In the case when
$\delta _{\!j}\to 0$
and
$\varOmega _{l}\approx \varOmega _{l0}$
, the following asymptotic formulae can be used:
\begin{align} Q\left (x_{1}\right )&\approx ix_{1}+\ldots , \end{align}
\begin{align} P\left (x_{2}\right )&\approx -ix_{2}+\ldots . \end{align}
Analysis of the system (3.18)–(3.22) revealed several limiting cases for which approximate estimates of the correction
$\Delta \varOmega _{l}$
can be obtained:
-
(i) first, with the dominance of Gibbs elasticity or dilatational surface viscosity, if
$h_{s}\ll 1$
,
$h_{d}\sim 1$
or
$|k| / (\varOmega _{l0})_{\tilde {\rho }_{2}=0}\sim 1$
and
$d_{s}\lesssim 1$
, then(3.31)
\begin{equation} \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}\approx -\frac {1+i}{2\sqrt {2}}\frac {\left (l+2\right )^{2}}{l}\sqrt {\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}\sqrt {\delta _{1}}-\frac {2\left (l-1\right )\left (l+2\right )}{l}ih_{s}; \end{equation}
-
(ii) second, if shear surface viscosity dominates and
$h_{s}\sim 1$
,
$h_{d}\ll 1$
,
$|k| / (\varOmega _{l0})_{\tilde {\rho }_{2}=0}\ll 1$
and
$d_{s}\ll 1$
, then(3.32)
\begin{equation} \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}\approx -\frac {1+i}{2\sqrt {2}}l\sqrt {\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}\sqrt {\delta _{1}}-\frac {2\left (l+1\right )}{\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}k-2\left (l+1\right )ih_{d}; \end{equation}
-
(iii) third, if the complicating factors are small as
$h_{s},h_{d},|k| / (\varOmega _{l0})_{\tilde {\rho }_{2}=0},d_{s}\ll \sqrt {\delta _{1}}$
, then(3.33)
\begin{align} \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}&\approx -\left (2l+1\right )\left (l+2\right )i\delta _{1}-\frac {1}{2}\frac {l+2}{l-1}\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0} k \nonumber \\&\quad -\frac {1}{2}l\left (l-1\right )\left (l+2\right )ih_{s}-\frac {1}{2}\left (l+1\right )\left (l+2\right )^{2}ih_{d}. \end{align}
In all three expressions presented above, the dynamic effect of gas in the bubble was neglected (
$\tilde {\rho }_{2}=0$
). Here, the real part of the correction
$\Delta \varOmega _{l}$
indicates the shift in the eigenfrequency of oscillations, while the imaginary part relates to the damping rate taken with the opposite sign (see (3.8)–(3.10), (3.16) and (3.17)).
The expression for the Rayleigh frequency for the light bubble configuration is used above,
\begin{equation} \left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}=\sqrt {(l^{2}-1)(l+2)}, \end{equation}
with which the correction
$ (\Delta \varOmega _{l} )_{\tilde {\rho }_{2}=0}$
should be summed up in the three approximations (3.31)–(3.33).
If the gas density in the bubble is considered non-zero, the contributions (3.31)–(3.33) should be corrected as follows:
-
(i) for the first case,
(3.35)
\begin{equation} \Delta \varOmega _{l}\approx \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}-\frac {1+i}{2\sqrt {2}}\frac {\left (l-1\right )^{2}}{l}\sqrt {\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}\:\tilde {\rho }_{2}\sqrt {\delta _{2}}; \end{equation}
-
(ii) for the second case,
(3.36)
\begin{equation} \Delta \varOmega _{l}\approx \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}-\frac {1+i}{\sqrt {2}}\frac {\left (l+1\right )^{2}}{l^{3/2}}\frac {1}{\sqrt {\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}}\:\tilde {\rho }_{2}\sqrt {\delta _{2}}; \end{equation}
-
(iii) for the third case,
(3.37)
\begin{equation} \Delta \varOmega _{l}\approx \left (\Delta \varOmega _{l}\right )_{\tilde {\rho }_{2}=0}-\frac {1+i}{2\sqrt {2}}\frac {\left (2l+1\right )^{2}}{l}\sqrt {\left (\varOmega _{l0}\right )_{\tilde {\rho }_{2}=0}}\:\tilde {\rho }_{2}\sqrt {\delta _{2}}. \end{equation}
And they are already summed up with the full Rayleigh frequency
$\varOmega _{l0}$
(see (3.28)). From expressions (3.35) and (3.36), it follows that the influence of gas is small due to the parameter
$\tilde {\rho }_{2}\approx 0$
compared with the liquid contribution from (3.31) and (3.32) and can be neglected. In expressions (3.33) and (3.37), the linear contribution of liquid viscosity must be compared with the root contribution of gas viscosity.
It is known that when deriving approximation relations for the correction to the natural frequency and the damping rate for the case of an interface between two media, it is important to know the order of magnitude of the ratios of their densities and viscosities. This is emphasized, for example, in Lu & Apfel (Reference Lu and Apfel1991), where the corresponding asymptotic formulae are obtained, but only for comparable dynamic properties of the media. In the present work, the internal medium was initially considered a gas with negligible density and viscosity. This allows us to exclude it from consideration at the stage of deriving the dispersion relation, to which the technique of expansions in series in a small parameter was then applied. Note (Miller & Scriven Reference Miller and Scriven1968) that the presence of surfactants promotes a root damping law in viscosity, which in its absence is typical for the interface between two liquids.
The asymptotic formulae in the present paper, (3.31) and (3.32), were found to be consistent with the asymptotic expressions available from Lu & Apfel (Reference Lu and Apfel1991), in which the density and viscosity of the internal medium formally tend to zero. Additionally, the contribution of the shear surface viscosity was included in (3.31), and the contribution of the dilatational surface viscosity was included in (3.32).
The case of (3.33) in the present work should be mentioned separately. Here, the parameters of the Gibbs elasticity, as well as the shear and dilatational surface viscosities, were considered to be simultaneously small, which ensures a linear contribution to the correction to the natural frequency in the viscosity of the liquid. At the same time, for the case of two media with comparable dynamic properties considered in Lu & Apfel (Reference Lu and Apfel1991), the contributions of their viscosities were root-valued even with such a weak influence of the surfactant. In terms of the contributions of the surface viscosities, (3.33) corresponds to the expression from Lu & Apfel (Reference Lu and Apfel1991).
If it is necessary to take into account the viscous contribution of the light gas in the bubble, it is necessary to return to the general dispersion relation for the two media. Such contributions (see (3.35)–(3.37)) turn out to be proportional to the square root of the gas viscosity and to the relative gas density. Qualitatively, they coincide with the asymptotic expressions available from Lu & Apfel (Reference Lu and Apfel1991). Complete quantitative correspondence is possible for (3.35) and (3.37), but not for (3.36). A comparison of the above approximation relations and numerical calculation data is shown in figure 2(a,b).
Figure 2. Dependences of the frequency shift (a) and damping rate (b) on the dissipative parameters of the media, when
$\delta _{1}=\delta _{2}$
. The relative density of the light phase is
$\tilde {\rho }_{2}=0.01$
, the meridional number is
$l=2$
and the Gibbs elasticity parameter is
$k=0$
. The surface viscosity parameters are
$h_{s}=1$
and
$h_{d}=0$
. The surface diffusion parameter is
$d_{s}=0$
. The solid line presents the numerical calculation. The dashed line shows the summary result of the asymptotic formulae (3.32) and (3.36). The long-dash line presents the corresponding approximation relation in Lu & Apfel (Reference Lu and Apfel1991).
As seen from (3.31)–(3.33), the transition to a square root dependence on liquid viscosity for the damping law, caused by surfactants, is typically accompanied by a decrease in the eigenfrequency. This trend is counteracted by the contribution proportional to the dimensionless Gibbs elasticity
$k\lt 0$
, which, conversely, can lead to an increase in frequency. The small shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities associated with surfactants increase the damping rate but either does not affect or weakly affects the frequency.
If we assume that the surface viscosity parameters
$h_{s}$
and
$h_{d}$
and the dimensionless Gibbs elasticity
$k$
remain finite, while the surfactant surface diffusion parameter
$d_{s}=0$
, then setting the determinant of system (3.18)–(3.22) to zero in the leading order with respect to
$\delta _{\!j}$
yields a simplified dispersion relation,
\begin{equation} \left (a\varOmega _{l}+1\right )\big(\varOmega _{l}^{2}-\varOmega _{l0}^2\big) +b\varOmega _{l}^{2}+c\varOmega _{l}=0. \end{equation}
Here
\begin{align}& a=\frac {l\left (l+1\right )\left (h_{s}+h_{d}\right )-2h_{s}}{l\left (l+1\right )k}i,\quad b=-\frac {4\left (l-1\right )\left (l+2\right )}{l\tilde {\rho }_{1}+\left (l+1\right )\tilde {\rho }_{2}}\frac {h_{s}h_{d}}{k}, \end{align}
\begin{align}&\qquad\qquad\qquad\quad c=\frac {4\left (l-1\right )\left (l+2\right )h_{s}}{l\tilde {\rho }_{1}+\left (l+1\right )\tilde {\rho }_{2}}i. \end{align}
A particular case of (3.38) for the configuration of a liquid droplet suspended in a vacuum, where
$\tilde {\rho }_{2}=1$
and
$\tilde {\rho }_{1}=0$
, can be found in Lyubimov et al. (Reference Lyubimov, Konovalov, Lyubimova and Egry2011).
Figure 3. The region of non-oscillatory solutions of (3.38) for the light bubble.
Figure 4. Dependences of the oscillation frequency (a) and damping rate (b) on the shear surface viscosity parameter. The dissipative parameter is
$\delta _{1}=0.001$
, the dilatation surface viscosity parameter is
$h_{d}=2$
, and the surface diffusion parameter is
$d_{s}=0$
. The dashed and long-dash lines show the results of the asymptotic formula (3.38). Solid lines `1’ and `2’ present two numerical solutions.
Figure 5. Dependences of the oscillation frequency (a) and damping rate (b) on the shear surface viscosity parameter. The dissipative parameter is
$\delta _{1}=0.001$
, the dilatation surface viscosity parameter is
$h_{d}=0.1$
and the surface diffusion parameter is
$d_{s}=0$
. The dashed and long-dash lines show the results of the asymptotic formula (3.38). Solid lines `1’ and `2’ present two numerical solutions.
Figure 3 shows a shaded region in the surface viscosity parameter space (
$h_{s}-h_{d}$
) where all solutions of (3.38) are non-oscillatory for the case of light bubble, when
$\tilde {\rho }_{1}=1$
and
$\tilde {\rho }_{2}=0$
. Here, the meridional number is taken as
$l=2$
, and the dimensionless Gibbs elasticity equals
$k=-1$
. The numerical results presented in figures 4(a,b) and 5(a,b) demonstrate that non-oscillatory damping is unattainable. For a fixed meridional number, no more than two solutions of the dispersion relation can be distinguished, which are approximately described by (3.38). The solution presented by lines `1’ in the figures can be called capillary, and the solution shown by lines `2’ is associated with the Gibbs elasticity and is absent for the parameter
$k=0$
. Such solutions are characterized by a non-zero real part of the eigenfrequency, which always tends to zero with increasing surface dissipation for the `elastic’ solution 2 and can remain finite for the capillary solution 1 (see figure 5
a). Also, the solutions can alternate in the magnitude of their damping decrement. So that the solution related with the Gibbs elasticity may become the least rapidly decaying.
As a result of the approximation, the ratio of the spherical Hankel functions of the first kind disappears from the simplified dispersion relation (3.38). This ratio is complex for the real argument. The latter would correspond to the non-oscillatory regime. According to Prosperetti (Reference Prosperetti1980), the presence of this ratio excludes the possibility that non-oscillatory solutions would be exact solutions of the general dispersion relation for a configuration with an external media, where the general expression for the velocity includes the indicated Hankel function.
3.2.2. Numerical data for cases of small and finite viscosities
Let us analyse the influence of Gibbs elasticity on the shift of the eigenfrequency (figure 6
a) and the damping rate (figure 6
b). The quadrupole mode with meridional number
$l=2$
was considered. A sufficiently small dissipative parameter for the liquid
$\delta _{1}=0.001$
was taken. For the gas in the bubble,
$\tilde {\rho }_{2}=0$
was assumed. Surface viscosity and diffusion were excluded from this consideration (
$d_{s},h_{s},h_{d}=0$
). Calculations show that the capillary mode remains the least rapidly damping. The dashed lines in the figures show estimate (3.31) at finite values of the Gibbs elasticity parameter
$k$
and small liquid viscosity
$\delta _{1}\to 0$
. The long-dash lines show a similar estimate (3.33) for small
$k$
. It can be stated that the asymptotic formulae obtained above have a sufficiently wide range of applicability.
Figure 6. Dependences of the frequency shift (a) and damping rate (b) of the capillary mode on the Gibbs elasticity parameter. The dissipative parameter is
$\delta _{1}=0.001$
.
The maximum damping rate, as well as the real part of the eigenfrequency, is observed at sufficiently small absolute values of parameter
$k$
(figure 6
b). Then, as
$|k|$
increases, the surfactant effect weakens and then reaches saturation. However, at finite
$|k|$
, the damping rate remains greater, and the real part of the eigenfrequency smaller, than in the case without surfactant.
Figure 7(a,b) show the dependences of the eigenfrequency and damping rate on the Gibbs elasticity parameter
$k$
for different values of the surface diffusion parameter
$d_{s}$
. It can be seen that as the latter increases, when the surface concentration of surfactant levels out faster along the interface, the influence of Gibbs elasticity on the characteristics of free oscillations of the system decreases. It can be expected that the same effect would be produced by redistribution of surfactants in the volume of liquid.
Figure 7. Dependences of the oscillation frequency (a) and damping rate (b) of the capillary mode on the Gibbs elasticity parameter. The surface diffusion parameter takes values
$d_{s}=0$
(curve 1),
$d_{s}=1$
(2),
$d_{s}=10$
(3) and
$d_{s}=100$
(4).
Figure 8. Dependences of the oscillation frequency (a) and damping rate (b) on the Gibbs elasticity parameter. The dissipative parameter is
$\delta _{1}=1/300$
. Solid lines `1’ and `2’ represent two different modes.
As shown in figure 8(a,b), constructed at
$d_{s}=0$
, the situation becomes more complicated at a slightly larger value of the dissipative parameter
$\delta _{1}=1/300$
. In this case, a competition arises between two solutions of the dispersion relation, obtained at a fixed meridional number
$l=2$
and showing alternately the least damping, which are presented in the figures by solid lines. The capillary mode shown by lines `1’ in the figures becomes more rapidly damping, as the value
$|k|$
increases, than the mode related to Gibbs elasticity, shown by lines `2’. The latter solution can be creeping, that is, it can demonstrate a small eigenfrequency and weak damping. Figure 8(a) shows the existence of a finite gap between the eigenfrequencies of the two solutions, even when their damping rates coincide (figure 8
b).
Let us analyse the case of a finite value of the dissipative parameter
$\delta _{1}=1$
(figure 9
a,b). In this case, the mode related to Gibbs elasticity, shown by lines `2’ in the figures, always has the least damping compared with the capillary mode shown by lines `1’. Also, as can be seen from figure 9(a), the real part of the eigenfrequency for the two modes behaves differently as
$|k|$
increases.
Figure 9. Dependences of the oscillation frequency (a) and damping rate (b) on the Gibbs elasticity parameter. The dissipative parameter is
$\delta _{1}=1$
. Solid lines `1’ and `2’ represent two different modes.
Figure 10. Dependences of the oscillation frequency (a) and damping rate (b) on the shear surface viscosity parameter
$h_{s}$
for the `elastic’ mode. The Gibbs elasticity parameter is
$k=-1$
, and the dissipative parameter is
$\delta _{1}=1/300$
. The dilatational surface viscosity parameter takes values
$h_{d}=0$
(curve 1),
$h_{d}=0.1$
(2) and
$h_{d}=1$
(3).
Figure 11. Dependences of the oscillation frequency (a) and damping rate (b) on the dilatational surface viscosity parameter
$h_{d}$
for the `elastic’ mode. The Gibbs elasticity parameter is
$k=-1$
, and the dissipative parameter is
$\delta _{1}=1/300$
. The shear surface viscosity parameter takes values
$h_{s}=0$
(curve 1),
$h_{s}=0.1$
(2) and
$h_{s}=1$
(3).
Figures 10(a,b) and 11(a,b) show the dependences of the eigenfrequency and damping rate (for the least rapidly damping mode that interests us) on the parameters of shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities. Here, the Gibbs elasticity parameter is taken as
$k=-1$
, the dissipative parameter equals to
$\delta _{1}=1/300$
, and the surface diffusion factor is
$d_{s}=0$
. With this choice of parameters, figures 12(a,b) and 13(a,b) demonstrate that the mode related to the Gibbs elasticity is the least rapidly damping. Regarding the damping rate (figure 10
b and 11
b), the same type of dependence on surface viscosities is found as on parameter
$k$
, where reaching a maximum is replaced by weakening and saturation of the effect. The real part of the eigenfrequency (figures 10
a and 11
a) can either decrease or increase with increasing surface viscosities.
Figure 12. Dependences of the oscillation frequency (a) and damping rate (b) on the shear surface viscosity parameter
$h_{s}$
. The Gibbs elasticity parameter is
$k=-1$
, the dissipative parameter is
$\delta _{1}=1/300$
and the dilatation surface viscosity parameter is
$h_{d}=0$
. Solid lines `1’ and `2’ represent two different modes.
Figure 13. Dependences of the oscillation frequency (a) and damping rate (b) on the dilatation surface viscosity parameter
$h_{d}$
. The Gibbs elasticity parameter is
$k=-1$
, the dissipative parameter is
$\delta _{1}=1/300$
and the shear surface viscosity parameter is
$h_{s}=0$
. Solid lines `1’ and `2’ represent two different modes.
Figure 14. Dependences of the oscillation frequency (a) and damping rate (b) on the shear surface viscosity parameter
$h_{s}$
for the least rapidly damping mode. The Gibbs elasticity parameter is
$k=-0.05$
, and the dissipative parameter is
$\delta _{1}=0.001$
. The dilatation surface viscosity parameter takes values
$h_{d}=0$
(curve 1),
$h_{d}=0.1$
(2) and
$h_{d}=1$
(3).
Figure 15. Dependences of the oscillation frequency (a) and damping rate (b) on the dilatational surface viscosity parameter
$h_{d}$
for the least rapidly damping mode. The Gibbs elasticity parameter is
$k=-0.05$
, and the dissipative parameter is
$\delta _{1}=0.001$
. The shear surface viscosity parameter takes values
$h_{s}=0$
(curve 1),
$h_{s}=0.1$
(2) and
$h_{s}=1$
(3).
Figures 14(a,b) and 15(a,b) were constructed for the dissipative parameter
$\delta _{1}=0.001$
, the surface diffusion factor
$d_{s}=0$
and the Gibbs elasticity parameter
$k=-0.05$
. The latter approximately corresponds to the maximum effect of Gibbs elasticity according to figure 6(a,b). Qualitatively, the same results are demonstrated as in the case of the previous drawings.
It can be said that there are two regimes. In the first regime, the damping rate, as well as the real part of correction to the eigenfrequency, change quite rapidly with an increase in some parameter associated with the surfactant, be it the absolute value of the Gibbs elasticity (figure 6 a,b), or the shear surface viscosity (figure 10 a,b) or the dilatational surface viscosity (figure 11 a,b). The reason is that the Marangoni force, as well as the surface viscous stresses, affects the pressure field, which is included in the condition of the balance of normal stresses, which leads to the change of the eigenfrequency of the system. With a further increase in the parameters, the saturation regime appears, the transition to which occurs, as a rule, through passing the maximum. Here, the velocity field on the interface changes in such a way as to reduce the Marangoni force or the surface viscous stresses. Thus, the shear surface viscosity factor is weakened if the tangential component of the velocity on the interface tends to zero. However, this gives the interface the properties of a membrane, which leads to an increase in vorticity near the boundary, which turns out to be significantly greater than in the case of a free surface without any surfactant. Such a solution through the term containing viscous stresses in the condition of the balance of normal stresses also leads to a change of the eigenfrequency.
4. Direct numerical simulation based on full nonlinear equations
4.1. Mathematical model
This section is devoted to direct numerical simulation of natural oscillations of a gas bubble in a viscous liquid with insoluble surfactant at the gas–liquid interface. The numerical simulation was performed for a domain representing a circular cylinder with the bubble located along its axis. The distance from the sidewalls of the cylinder to the bubble was chosen sufficiently large to exclude their influence on the flow near the bubble (approximately
$15$
radii). A computational experiment showed that the value of the velocity components at a distance of 12 radii does not exceed
$10^{-6}$
. Calculations were carried out within an axisymmetric approach, in which case the computational domain represents a rectangle (figure 16). The numerical solution of the problem under consideration in a full three-dimensional formulation, allowing for asymmetric oscillation modes, is a rather labour-intensive task. The main difficulty is associated with the development of a computational algorithm with explicit identification of the shape of the interphase boundary (tracking method), which will allow for dynamic boundary conditions on it. The `volume of fluid’ and `level set’ algorithms (capture method), which are widely used in computational fluid dynamics, are not capable of solving problems with the boundary conditions used in this work.
Figure 16. The computational domain.
The motion of the liquid phase was described by the Navier–Stokes (2.2) and continuity (2.3) equations, which are projected onto the axes of the cylindrical coordinate system
$\{r, z\}$
(figure 16).
The gas in the bubble was considered ideal and described by the Mendeleev–Clapeyron equation of state
\begin{equation} pV=R_{g}T, \end{equation}
where
$R_{g}\:(\rm J\:mol^{-1}\,K^{-1})$
is the universal gas constant,
$T\:(K)$
is the temperature and
$V\:(m^{3})$
is the volume of the gas bubble. The viscosity of the gas in the bubble is neglected. The boundary conditions at the interface include (2.6)–(2.10), written in projections on the axes of the Cartesian coordinate system
$\left \{n, \tau \right \}$
(figure 16), normally associated with the free surface. As a result, the conditions of balance of normal and tangential stresses on the free surface in dimensionless form are written as
\begin{align}& p_{2}-p_{1}+2\delta _{1}\frac {\partial u_{n}}{\partial n} =-\overline {\gamma }\left (K_{1}+K_{2}\right )-\left \{K_{1}\left (h_{s}+h_{d}\right )-K_{2}\left (h_{\tau }-h_{d}\right )\right \}\left (\frac {\partial u_{\tau }}{\partial \tau }+u_{n}K_{1}\right ) \nonumber \\&\qquad\qquad\qquad\qquad\,\, +\left \{K_{1}\left (h_{s}-h_{d}\right )-K_{2}\left (h_{s}+h_{d}\right )\right \}\frac {u_{r}}{r}, \end{align}
\begin{align}& \delta _{1}\left (\frac {\partial u_{\tau }}{\partial n}+\frac {\partial u_{n}}{\partial \tau }-u_{\tau }K_{1}\right ) =\frac {\partial }{\partial \tau }\left \{\overline {\gamma }+\left (h_{s}+h_{d}\right )\left (\frac {\partial u_{\tau }}{\partial \tau }+u_{n}K_{1}\right )-\left (h_{s}-h_{d}\right )\frac {u_{r}}{r}\right \} \nonumber \\&\qquad\qquad\qquad\qquad\qquad\quad\,\, +2h_{s}\left (\frac {\partial u_{\tau }}{\partial \tau }+u_{n}K_{1}-\frac {u_{r}}{r}\right )\sqrt {\frac {1}{r^{2}}-K_{2}^{2}}, \end{align}
where
$K_{1}$
and
$K_{2}$
are the negative principal curvatures of the concave surface, and
$\overline {\gamma } =\gamma /\gamma _{0}$
is a dimensionless coefficient of surface tension. The transport of surfactants along the boundary is described by (2.4). The tension coefficient is calculated in accordance with the Henry isotherm, which assumes a linear dependence on the concentration (2.15). We do not use the more general Langmuir law, since we believe that even in the numerical solution, the deviations of the surface tension coefficient from its equilibrium value are small and obey a linear dependence on the surface concentration of the surfactant.
On the walls of the cylinder, no-slip conditions are specified, and on the axis of symmetry, symmetry conditions are specified. At the initial moment of time, the bubble has the shape of a flattened ellipsoid and is in a liquid at rest. The dimensionless semiaxes of the ellipsoid are equal
$(1-\varDelta )$
and
$(1-\varDelta )^{-0.5}$
, where
$\varDelta$
is a geometric parameter. The surfactant is distributed uniformly along the boundary with a concentration of
$\varGamma _{0}$
.
4.2. Numerical method
The numerical solution of the problem was conducted using an original Lagrangian–Eulerian computational technique. The solution domain was covered with a fixed staggered grid with non-uniform spacing. The grid spacing was refined in the vicinity of the free surface and gradually increased with distance from it. The phase interface was represented as an ordered set of marker particles that served as nodes of a moving computational grid.
The equations of motion (2.2) and continuity (2.3) were discretized using the finite volume method, with the SIMPLE (semi-implicit method for pressure-linked equations) algorithm (Patankar Reference Patankar1980) employed for pressure–velocity coupling. Three types of control volumes were identified in the domain: those not containing the liquid phase where computations were not performed; complete control volumes where standard calculations were carried out; and partially filled control volumes through which the free surface passed. Calculations in the latter required the construction of irregular control volumes and the use of values of the computed variables from the interface.
The computation of characteristics on the free surface was performed in two stages. At the first stage, a local Cartesian coordinate system associated with the boundary was introduced at each point of the surface. The equations for stresses (4.2) and (4.3) and the continuity (2.3) were projected onto its axes, and their finite-difference analogues were obtained using a finite-difference approach. The discretized continuity equation and the balance of tangential stresses were used to determine the normal and tangential components of velocity, while the equation for the balance of normal stresses was used to determine the pressure. At the second stage, the convective diffusion (2.4) was discretized, and the surface concentration was calculated. This implementation of direct numerical simulations is valid only for zero viscosity of the dispersed phase, but the conditions for normal and shear stresses at the interphase boundary are explicitly satisfied, taking into account dependence of gas pressure on bubble volume.
The movement of marker particles was carried out in accordance with the finite-difference analogues of the kinematic boundary conditions. The free surface is approximated using a cubic spline, the coefficients of which are calculated according to the method described in Zav’ialov et al. (Reference Zav’ialov, Kvasov and Miroshnichenko1980). In this case, the polar coordinate system is used. When calculating the curvature of the surface, the first and second derivatives are determined using a spline. During the deformation of the boundary over time, areas with local clustering or rarefaction of markers could form on it, which adversely affected the stability of the numerical calculation. In this regard, a procedure for redistributing markers was periodically performed using the mentioned spline with an interval of 0.25 dimensionless time units.
Figure 17. Temporal dependences of the velocity of the south pole (a), the total amount of surfactant on the surface (b) and the distribution of surface concentration (c) at t = 0.5. Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
,
$h_{s}=0.5$
and
$h_{d}=0.5$
.
$\varDelta =0.005$
(1, 2, 3, 4 are
$h=1/20$
,
$1/40$
,
$1/80$
,
$1/160)$
.
Figure 18. Temporal dependences of the velocity of the south pole (a), the total amount of surfactant on the surface (b) and the distribution of surface concentration (c) at t = 0.5. Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
,
$h_{s}=0$
, and
$h_{d}=0$
.
$\varDelta =0.005$
(1, 2, 3, 4 are
$h=1/20$
,
$1/40$
,
$1/80$
,
$1/160)$
.
Tests were performed to verify the approximation convergence on a sequence of grids. The dependences of the velocity of the south pole on time, calculated on different grids, for the case with and without surface viscosity are shown in figures 17(a) and 18(a), respectively. Table 1 shows the values of time and velocity at the first two extremum points, demonstrating the convergence of the calculation algorithm. Here,
$h$
is the minimum grid spacing in the vicinity of the bubble. It can be seen that the proposed computational technique demonstrates an order of accuracy not lower than the first. All further calculations were performed on a grid with spacing
$h=1/80$
.
Table 1. Values of time and velocity at the first two extremum points for different grid spacings.
Figures 17(b) and 18(b) show a graph of the temporal behaviour of the total amount of surfactant, i.e. the integral of the surfactant concentration
$\varGamma _1$
over the bubble surface. This demonstrates the fulfilment of the conservation law. The characteristic
$\varGamma _{1}$
is calculated as a surface integral of the concentration
$\varGamma$
over the interface, computed using the rectangle method. The step-like changes in the characteristics on coarse grids occur at the step of performing the procedure for redistributing the markers. Figures 17(c) and 18(c) illustrate the distribution of the surface concentration of surfactant along the free surface. The coordinate
$\varphi$
is the angle by which the polar axis must be rotated anticlockwise in order to reach the current point of the free boundary. The polar axis originates from a point located midway between the north and south poles and is directed vertically downward.
4.3. Calculation results
At the initial moment of time, the free surface of the bubble begins to deform under the action of surface tension forces and the Marangoni effect. Subsequently, the motion takes on a decaying oscillatory character, illustrated in figure 19, which shows graphs of the change in axial velocity and surface concentration at the north pole over time. It can be seen that with increasing values of surface viscosities, both the damping rate and frequency increase. Thus, the first maximum of speed is achieved at time 0.625, 0.535, 0.47 for the cases in figures 19 (a), 19(b) and 19(c), respectively. During periods when the velocity
$V_{n}$
has positive values, the bubble deforms from a flattened state (A) to one elongated along the axis of symmetry (D). The change in concentration for the selected similarity numbers is mainly associated with local changes in surface area and convective transport.
Figure 19. Dependence of axial velocity and concentration
$\varGamma$
of the north pole on time. Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
and
$\varDelta =0.005$
. Here (a)
$h_{s}=h_{d}=0$
, (b)
$0.2$
and (c) are
$0.5$
.
When deforming from a flattened state (A) to a state close to spherical (B), the area of the selected region near the north pole increases, and the concentration consequently decreases. During deformation from the spherical state (B) to an ellipsoid elongated along the axis (C), the area continues to increase, while the concentration changes its trend to increasing due to enhanced convective transport. On the negative sections of
$V_{n}$
(C–D–E), the picture is exactly the opposite. At the same time, a slight phase difference between velocity and concentration is observed. Non-zero Marangoni stresses make an additional contribution to the convective transport of surfactants along the surface, since they are taken into account through the boundary conditions for the jump in tangential stresses at the interface. The resulting phase difference is a consequence of the mutual influence of convective transport and local changes in surface area. In the absence of surface viscosity, the velocity lags behind the concentration, while with increasing viscosity coefficients
$h_{s}$
and
$h_{d}$
the situation changes, and the velocity begins to lead the concentration.
Figure 20 shows the distributions of calculated characteristics on the free surface for different values of surface viscosity at times corresponding to spherical surface shapes (point B, figure 19
a). It can be seen that with increasing values of
$h_{s}$
and
$h_{d}$
, the absolute values of the tangential velocity component decrease, which leads to a weakening of the convective transport mechanism of surfactant along the surface. Thus, at zero surface viscosity, the concentration changes due to area change and convective transport practically balance each other at point B. With increasing surface viscosity, the term associated with local area changes varies insignificantly, while convective terms decrease, leading to a phase lag. The velocity distributions at the time when the surface has a maximally elongated shape along the axis of symmetry (point C, figure 19
a) are close to zero.
Figure 20. Distributions of tangential velocity (a), normal velocity (b) and concentration (c) along the free boundary. Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
and
$\varDelta =0.005$
(solid, dot–dash, dashed lines are for
$h_{s}=h_{d}=0$
,
$0.2$
,
$0.5)$
.
Figure 21. Distribution of radial (a), axial (b) velocity components, pressure (c) and streamlines (d) at
$t=0.6$
.
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
,
$h_{s}=h_{d}=0$
and
$\varDelta =0.005$
.
The flow kinematics near the bubble at the moment when the surface shape is close to spherical (point B, figure 19) is demonstrated in figure 21. The flow is characterized by the formation of two symmetric circulation zones. Calculations showed that accounting for surface viscosity in the considered range does not qualitatively affect the flow in the bulk phase.
Figure 22 demonstrates the effect of the initial deviation of the shape from a sphere on the damping dynamics. It can be seen that changing the parameter
$\varDelta$
in the range from
$0.005$
to
$0.1$
does not significantly affect the values of the damping rate and oscillation period. The difference in the integral values of the surfactant at the boundary is explained by the difference in the areas of the initial surface shape. As a result, the radii of the equilibrium spherical shapes of the bubble will differ depending on
$\varDelta$
, since they will have different concentrations of the surfactant and, as a consequence, different surface tension coefficients.
Figure 22. Dependences of velocity (a), axial coordinate (b) of the north pole and amount of surfactant (c) on time for different initial surface shapes. Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
and
$h_{s}=h_{d}=0.2$
(1, 2, 3, 4, 5 are for
$\varDelta =0.005$
,
$0.01$
,
$0.02$
,
$0.05$
,
$0.1$
).
Figure 23. Dependence of natural frequency (a) and damping rate (b) on the Gibbs elasticity parameter in the absence of surface viscosity (points are for numerical results, solid and dashed lines are for linear theory). Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
and
$\varDelta =0.005$
.
Figure 24. Dependence of natural frequency (a) and damping rate (b) on the shear surface viscosity parameter
$h_{s}$
(points are for numerical results, lines are for linear theory). Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
,
$h_{d}=0$
and
$\varDelta =0.005$
.
Figure 25. Dependence of natural frequency (a) and damping rate (b) on the dilatational surface viscosity parameter
$h_{d}$
(points are for numerical results, lines are for linear theory). Here
$\delta _{1}=1/300$
,
$d_{s}=1.176\times 10^{-6}$
,
$k=-1$
and
$\varDelta =0.005$
(solid, dashed lines are for
$h_{s}=0$
,
$1)$
.
Figure 23(a,b) present a comparison of numerical modelling results with results obtained within the framework of linear theory for the case of absence of surface viscosity in terms of natural frequency and damping rate characteristics. The solid and dashed lines represent two solutions of the linear theory, corresponding alternately to the least rapidly damping mode regime. Figures 24(a,b) and 25(a,b) demonstrate a similar comparison taking into account non-zero values of shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities. It is evident that the numerical solution tends to the least rapidly damped mode of natural oscillations, which is the mode related to the Gibbs elasticity (compare the parameters used here and data of figures 12
b and 13
b). The damping rate was determined from the extremum points on the curve of the dependence of the north pole coordinate on time. For cases with large damping rate values (greater than
$1.5$
), the oscillations damped quickly enough, making it difficult to extract a sufficient number of points for approximation. The discrepancies between the theoretical and numerical results become significant when taking into account
$h_{s}$
. The observed discrepancies with increasing surface viscosity
$h_{s}$
are associated with the presence of a nonlinear term in the equation of balance of tangential stresses (the last term on the right-hand side of (4.3)) in the formulation for the numerical problem. For the case
$h_{s} = 0$
, the results agree qualitatively and quantitatively.
5. Conclusions
Taking into account the combined influence of bulk viscosities of the liquid and gas, shear and dilatational surface viscosities, Marangoni force and surfactant surface diffusion, we have investigated small capillary oscillations of a gas bubble in a viscous liquid with insoluble surfactant at the gas–liquid interface. The key dimensionless parameters of the problem are the Gibbs elasticity,
$k$
, the coefficients of surface viscosities,
$h_{s}$
and
$h_{d}$
, and the surface diffusion,
$d$
.
For small oscillation amplitudes, the problem was reduced to solving linearized Navier–Stokes equations and continuity equations in spherical coordinates. The boundary conditions on the bubble surface included the balance of normal and tangential stresses, the kinematic condition for interface deformation and the surfactant transport equation. Solutions of toroidal type were sought in the form of spherical harmonics.
The dispersion relation for the complex oscillation frequency was obtained, taking into account the influence of surface viscosities and Gibbs elasticity. Analytical estimates were obtained for the frequency shift and damping rate in limiting cases of small viscosities, where surface effects dominate. It was shown that the gas contribution to dissipation is negligible due to its much lower density compared with the liquid.
Gibbs elasticity significantly affects the dynamics of the gas bubble. As the absolute value of the elasticity parameter
$k$
increases, the oscillation frequency decreases while the damping rate grows. Maximum dissipation is observed at moderate absolute values of
$k$
, after which the effect is weakening and saturates. Surface diffusion weakens the influence of Gibbs elasticity by levelling the surfactant concentration and reducing surface tension gradients.
The shear
$h_{s}$
and dilatational
$h_{d}$
surface viscosities increase the damping rate but may either increase or decrease the oscillation frequency.
The analysis showed that non-oscillatory damping (zero real part of frequency) is unattainable – even in the limit of large viscosities, a small non-zero frequency persists.
For cases with noticeable viscous dissipation in the liquid, the coexistence of two modes was discovered: the capillary mode (with high frequency and damping) and a creeping ‘elastic’ mode (with low frequency and weak damping). The latter is associated with Gibbs elasticity and may become dominant as the absolute value of parameter
$k$
increases due to its weak damping.
Using an original Lagrangian–Eulerian method, direct numerical simulation of natural oscillations of a gas bubble in liquid with insoluble surfactant at the interface was performed, based on full nonlinear equations. The obtained dependences of damping rate and frequency on Gibbs elasticity and dimensionless parameters of shear and dilatational viscosities agree well with the results obtained within the linear theory. The distributions of kinematic characteristics and surfactant surface concentration along the free boundary and in its vicinity were determined.
The results of this study provide a comprehensive understanding of the complex interplay between surface rheology, surfactant transport and bulk viscous effects in bubble dynamics, while establishing the limits of applicability of linear theory through detailed numerical validation. The developed computational approach enables accurate modelling of finite-amplitude oscillations with arbitrary surface rheological properties.
A potential application of the obtained results can be the determination of rheological characteristics of surfactant-covered surfaces, which are very important for planning experimental studies on the behaviour of liquids containing bubbles, by measurement of the frequency and the damping rate of shape-changing bubble oscillations. Also, the present paper can be considered as a starting point for solving the broader problem of the quantitative description of technological processes involving such liquids, among them the flotation separation, which is the main method of mineral processing. For that goal, the analysis has to be extended by the investigation of the influence of bubble rising, which creates a non-uniform distribution of the surfactant on the bubble surface. It is expected that bubble oscillations can have a significant effect on the interaction of solid particles with gas bubbles in that process.
Acknowledgements
The work was supported by the Russian Science Foundation (grant no. 24-11-00269).
Declaration of interests
The authors report no conflict of interest.
