1. Introduction
The phenomenon of gas being encapsulated in liquid drops is frequently observed when drops collide with solid surfaces (Pumphrey & Elmore Reference Pumphrey and Elmore1990; Rein Reference Rein1993; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012). It occurs in various scenarios, ranging from natural phenomena such as a raindrop falling to industrial applications such as the inkjet printing process (Singh et al. Reference Singh, Haverinen, Dhagat and Jabbour2010; Lohse Reference Lohse2022). After a drop impacts a solid surface, it often forms a stable configuration of sessile drops entrapping a small air bubble beneath their centre (see figure 1
$a$
) (Mehdi-Nejad, Mostaghimi & Chandra Reference Mehdi-Nejad, Mostaghimi and Chandra2003; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016). Other different configurations of stable compound sessile drops exist but are not commonly encountered in practice (Zhang et al. Reference Zhang, Gao, Li and Ding2021). Meanwhile, the impact of the drop generates oscillations of the bubble-containing (BC) sessile drop (Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005; Kern, Bostwick & Steen Reference Kern, Bostwick and Steen2021). The drop vibrations exhibit various interesting dynamic behaviours, including directional movement (Ding et al. Reference Ding, Zhu, Gao and Lu2018), Faraday-like waves (Vukasinovic, Smith & Glezer Reference Vukasinovic, Smith and Glezer2007) and liquid atomisation (James et al. Reference James, Vukasinovic, Smith and Glezer2003) of drops. However, previous studies of drop oscillations have not taken into account the presence of trapped bubbles at the bottom centre of sessile drops. This omission is a crucial gap in the existing research, as the presence of an inner bubble can affect the drop’s vibrational characteristics and induce a new class of oscillatory modes.
The oscillations of bubble-free sessile drops have been extensively studied within a theoretical framework of potential flow, assuming inviscid liquids and small amplitudes (linear) (Strani & Sabetta Reference Strani and Sabetta1984; Lyubimov, Lyubimova & Shklyaev Reference Lyubimov, Lyubimova and Shklyaev2004; Bostwick & Steen Reference Bostwick and Steen2009, Reference Bostwick and Steen2014; Bostwick & Steen Reference Bostwick and Steen2009; Ding & Bostwick Reference Ding and Bostwick2022b
; Sharma & Wilson Reference Sharma and Wilson2021; Zhang, Zhou & Ding Reference Zhang, Zhou and Ding2023; Zhang & Zhou Reference Zhang and Zhou2023). These works examined the effects of the contact angle
$\alpha$
, gravity and contact line (CL) mobility parameter
$\varLambda$
, on the frequency spectrum. A hemispherical drop with a free CL, referred to as the free semi-drop, often serves as a comparative prototype for stability analysis (Lyubimov, Lyubimova & Shklyaev Reference Lyubimov, Lyubimova and Shklyaev2006). This is because its spectrum, consisting of natural frequencies
$\lambda$
and mode shapes, inherits a subset of the spectrum of a spherical drop due to plane symmetry. This subset is described by the half–Rayleigh–Lamb (half–RL) spectrum (Lamb Reference Lamb1932; Lyubimov et al. Reference Lyubimov, Lyubimova and Shklyaev2006)
\begin{align} \lambda _{RL} ^2 = k\left ( {k - 1} \right )\left ( {k + 2} \right )\quad k,l \in \mathbb {N} ,\quad l\leq k, \quad k+l=\mathbb {N}_{{even}}, \end{align}
where
$k$
and
$l$
are the polar and azimuthal wavenumbers, respectively, and
$\mathbb {N}$
and
$\mathbb {N}_{{even}}$
denote natural and positive even numbers, respectively. Here, the frequency
$\lambda _{RL}$
is scaled by the capillary frequency
$\sqrt {\sigma /\varrho R^3}$
, where
$\sigma$
is the surface tension,
$\varrho$
the liquid density and
$R$
the drop radius. The most notable feature of the half–RL spectrum (1.1) is its spectral degeneracy, where modes with the same
$k$
share the same frequency. Another feature is that the Noether mode
$\{k=1,l=1\}$
of the free semi-drop has zero frequency due to its horizontal translational invariance. By pinning the CL in free semi-drops, Lyubimov et al. (Reference Lyubimov, Lyubimova and Shklyaev2004, Reference Lyubimov, Lyubimova and Shklyaev2006) first observed the breaking of spectral degeneracy. Bostwick & Steen (Reference Bostwick and Steen2014) later extended this to spherical-cap drops with varying contact angles
$\alpha$
using the Green’s function method, showing that varying
$\alpha$
from
$90^\circ$
also breaks spectral degeneracy. They also explored how the mobility parameter
$\varLambda$
, which relates contact angle deviation to CL velocity as
$\Delta \alpha = \varLambda U_{CL}$
, affects the spectrum. Recently, Zhang et al. (Reference Zhang, Zhou and Ding2023) used a boundary element method (BEM) model to study how gravity affects the spectrum of sessile drops with either a free or pinned CL, indicating the potential for exploring complex drop geometries (Viola & Verzicco Reference Viola and Verzicco2023). Using the same model, Zhang & Zhou (Reference Zhang and Zhou2023) examined the effects of reverse gravity in pendant drops. The above findings have been validated by numerical simulations (Basaran & DePaoli Reference Basaran and DePaoli1994; Sakakeeny & Ling Reference Sakakeeny and Ling2020, Reference Sakakeeny and Ling2021) and experiments (Noblin et al. Reference Noblin, Buguin and Brochard-Wyart2004, Reference Noblin, Buguin and Brochard-Wyart2005,Reference Noblin, Buguin and Brochard-Wyart2009; Chang et al. Reference Chang, Bostwick, Steen and Daniel2013; McCraney et al. Reference McCraney, Kern, Bostwick, Daniel and Steen2022; McCraney et al. Reference McCraney, Kern, Bostwick, Daniel and Steen2022).
Figure 1. (
$a$
) Diagram of a free BC semi-drop with inner
$R_2$
and outer
$R_1$
radii. (
$b$
–
$e$
) Cross-sections of minus and plus modes with
$\tau =0.75$
. (
$f\!$
,
$g$
) Frequency spectrum
$\lambda _\pm$
and amplitude ratios
$\xi _\pm$
for minus (blue solid) and plus (red dashed) modes, where curves with the same
$k_\pm$
and different
$l$
coincide into a single curve that exhibits spectral degeneracy. The blue dots in (
$f$
) correspond to cases of
$\tau = 0$
, given by the half–RL spectrum (1.1). (
$h$
–
$k$
) Scaling laws for the limiting cases of
$\tau \to 0$
and
$\tau \to 1$
, with green lines indicating asymptotic solutions (1.2) or (1.3). The values of
$k_\pm$
are 2, 3, 4, 5 and 6 along the arrows.
As previously mentioned, bubble entrapment at the centre of a drop’s base is widespread; however, its vibrational characteristics have not yet been thoroughly explored. Similar to free semi-drops (Lyubimov et al. Reference Lyubimov, Lyubimova and Shklyaev2006), one can infer that BC sessile drops with
$\alpha =90^\circ$
and free CLs (referred to as free BC semi-drops) inherit half (
$k+l=\mathbb {N}_{{even}}$
) of the spectrum of concentric spherical BC drops due to plane symmetry (see figure 1
$a$
). Spherical BC drops belong to a special class of compound drops with zero core density, which has been well studied in both concentric (Saffren et al. Reference Saffren, Elleman and Rhim1981; Bhattacharya Reference Bhattacharya2016; Shiryaev Reference Shiryaev2020) and eccentric (Lyubimov et al. Reference Lyubimov, Konovalov, Lyubimova and Egry2012; Sumanasekara & Bhattacharya Reference Sumanasekara and Bhattacharya2017) configurations. Compared with bubble-free drops, BC sessile drops exhibit new oscillatory modes dominated by the deformation of the inner bubbles (so-called bubble modes; see e.g. figure 1
$d,\!e$
). For free BC semi-drops, their modes, classified as bubble (new) and sloshing modes, inherit the following five properties from the modes of concentric spherical BC drops (Saffren et al. Reference Saffren, Elleman and Rhim1981; Bhattacharya Reference Bhattacharya2016) due to plane symmetry:
-
(i) The sloshing modes, also known as out-of-phase or minus modes, display vibrational directions that are opposite between the inner and outer surfaces, resulting in a phase difference of
$180^\circ$
(figure 1
$b$
,
$c$
). In these modes, the term ‘minus’ refers to the reduced coupling energy (i.e. frequency), analogous to coupled oscillators, which results in oscillators moving in opposite directions. This reduction in energy leads to the out-of-phase relationship. In contrast, the bubble modes, or plus modes, vibrate in phase between the two surfaces (figure 1
$d$
,
$e$
). -
(ii) For fixed
$k$
, plus modes always have higher frequencies than minus modes for the same size ratio
$\tau = R_2/R_1$
, where
$R_1$
and
$R_2$
are the radii of the drop and bubble, respectively. This higher frequency in plus modes is due to the greater coupling energy associated with in-phase motion. Additionally, in the limits
$\tau \to 0$
and
$\tau \to 1$
, there exist asymptotic solutions for the dimensionless frequencies (Bhattacharya Reference Bhattacharya2016)(1.2)Here, subscripts ‘
\begin{align} \lambda _-^a=\sqrt {k\left ( {k - 1} \right )\left ( {k + 2} \right )(1-\tau ^{(k+1)/2})} \quad \textrm {and} \quad \lambda _+^a=\sqrt {\frac {(k+2)(k^2-1)}{\tau ^3-\tau ^{(k+7)/2}}}. \end{align}
$-$
’ and ‘
$+$
’ denote minus and plus modes, respectively, and the characteristic frequency is
$\sqrt {\sigma /\varrho R_1^3}$
. As seen from (1.2), for
$k\gt 1$
modes, as
$\tau \to 0$
, frequencies
$\lambda _-$
of minus modes converge to the half-RL spectrum (1.1) (see the blue dots in figure 1
$f$
), while
$\lambda _+ \sim \tau ^{-3/2}\to \infty$
. In the other limit
$\tau \to 1$
,
$\lambda _-\sim (1-\tau )^{1/2}\to 0$
and
$\lambda _+ \sim (1-\tau )^{-1/2}\to \infty$
. Figure 1(
$h$
–
$j$
) illustrates the above scaling laws with green lines.
-
(iii) Compared with minus modes, plus modes always exhibit greater bubble deformation, characterised by the ratio
$\xi _+ \gt 1$
of inner to outer amplitudes (see figure 1
$g$
). Conversely, minus modes are characterised by
$\xi _- = 1/\xi _+ \lt 1$
. As
$\tau \to 0$
, asymptotic solutions
$\xi _\pm ^a$
have been derived (Bhattacharya Reference Bhattacharya2016)(1.3)This implies a scaling law
\begin{align} \xi _-^a=1/\xi _+^a=\frac {2k+1}{k+1}\tau ^{k+2}. \end{align}
$\xi _-=1/\xi _+\sim \tau ^{k+2}$
as
$\tau \to 0$
(figure 1
$k$
).
-
(iv) From an overview perspective, plus modes manifest U-shaped spectral lines, indicating a decrease followed by an increase in frequency as
$\tau$
increases, whereas the spectral lines of minus modes start flat and then decrease rapidly (figure 1
$f$
). -
(v) In addition to the aforementioned differences between plus and minus modes, both exhibit spectral degeneracy (Bostwick & Steen Reference Bostwick and Steen2014): modes with the same
$k_\pm$
have identical frequencies
$\lambda _\pm$
and amplitude ratios
$\xi _\pm$
(figure 1
$f\!$
,
$g$
).
In summary, free BC semi-drops inherit half of the spectrum of concentric spherical BC drops, offering a preliminary understanding of the vibrational characteristics of BC sessile drops. However, the model is limited to the case of contact angle
$\alpha =90^\circ$
and cannot consider the pinned CL. For bubble-free drops, Bostwick & Steen (Reference Bostwick and Steen2014) found that varying
$\alpha$
from
$90^\circ$
or CL pinning can lead to frequency shifts and the breaking of spectral degeneracy in free semi-drops. For BC sessile drops, varying
$\alpha$
or CL pinning introduces similar modifications to the spectrum of free BC semi-drops. Our work shows that these modifications include, not only frequency shifts and the breaking of spectral degeneracy, but also unexpected changes in the topological structure of spectral lines.
The compressibility of gas within internal bubbles is often neglected in BC droplets (e.g. Bhattacharya Reference Bhattacharya2016; Sumanasekara & Bhattacharya Reference Sumanasekara and Bhattacharya2017), as it introduces a distinct mode, known as the breathing or volume mode, which causes volume fluctuations or bubble collapse. In contrast, the previously discussed modes, referred to as shape modes, generally do not involve such volume changes. For spherical bubbles, shape and volume modes are governed by the RL (Rayleigh Reference Rayleigh1879; Lamb Reference Lamb1932) and Rayleigh–Plesset equations (Rayleigh Reference Rayleigh1917; Plesset Reference Plesset1949), respectively, and remain uncoupled in the linear regime, although nonlinear coupling may occur (see Feng & Leal (Reference Feng and Leal1997) for a review). In partially wetting bubbles (Ding & Bostwick Reference Ding and Bostwick2022a ), however, shape and volume modes become coupled at low static bubble pressures, resulting in an avoided crossing of their frequencies. At high bubble pressures, the first few shape modes remain unaffected, with the volume nearly conserved, implying that the bubble is incompressible. Our study focuses on the oscillations of BC sessile droplets at atmospheric pressure, where the volume mode occurs at a high frequency and has minimal impact on the first few shape modes, which are most easily excited. Nonetheless, we incorporate gas compressibility using the linearised ideal gas equation of state for physical fidelity.
In this paper, we study the vibrational characteristics of BC sessile drops with a free or pinned CL on a plane, assuming the liquid is inviscid and the oscillation amplitude is infinitesimal. Notably, our findings reveal the presence of avoided crossings in the spectral lines of minus and plus modes for BC sessile drops, extending beyond free semi-drops. The avoided crossing of two eigenvalues, also called a repulsion of eigenvalues, is related to symmetry breaking and has been observed in many fields, such as quantum systems (Faust et al. Reference Faust, Rieger, Seitner, Krenn, Kotthaus and Weig2012), oscillating bubbles (Ida Reference Ida2005) and Rayleigh–Bénard convection (Cliffe & Winters Reference Cliffe and Winters1986). Finally, we characterise these avoided crossings and discuss the potential occurrence of beat phenomena in coupled oscillators as a consequence (Thomson & Dahleh Reference Thomson and Dahleh1997), which is confirmed by our direct numerical simulations.
Figure 2. (
$a$
) Schematic diagram of a BC sessile drop on a plane in cylindrical coordinates
$(r,\varphi ,z)$
, where
$\varGamma _1$
and
$\varGamma _2$
represent the unperturbed surfaces of the drop and bubble, respectively, and
$\eta _{1,2}$
denote the corresponding perturbations. (
$b$
) Rendered mode shape with an azimuthal wavenumber
$l=4$
, showing
$n_+=1$
layer of perturbation on the bubble surface and
$n_-=7$
layers on the drop surface.
2. Mathematical formulation
Consider a BC sessile drop on a plane in the absence of gravity, where the unperturbed surfaces
$\varGamma _{1,2}$
of the drop and bubble form two encapsulated spherical caps (see figure 2
$a$
). We assume that the wettability of the solid surface is uniform, which ensures that the contact angle
$\alpha$
of the liquid remains constant, leading to the bubble’s contact angle being the complement of the drop’s (i.e.
$180^\circ - \alpha$
). Consequently, this configuration is governed by the outer CL radius
$R_1$
, inner CL radius
$R_2$
and contact angle
$\alpha$
. We define the characteristic length as
$L_*=R_1$
, which gives the characteristic pressure
$\sigma /L_*$
, time
$\sqrt {{\varrho L_*^3}/{\sigma }}$
, and potential
$\sqrt {{\sigma L_*}/{\varrho }}$
. Unless otherwise specified, all subsequent quantities are scaled by these characteristic values. Thus, the independent geometric parameters governing the drop configuration are the bubble size
$\tau =R_2/R_1$
and contact angle
$\alpha$
(figure 2
$a$
).
Under the assumptions of an ideal liquid (incompressible and inviscid) and by neglecting the densities of the ambient gas and the bubble, the linear oscillations of BC sessile drops can be analysed using an inviscid theoretical framework (Myshkis et al. Reference Myshkis, Babskii, Kopachevskii, Slobozhanin and Tyuptsov1987; Bostwick & Steen Reference Bostwick and Steen2014; Zhang et al. Reference Zhang, Zhou and Ding2023). This model is derived from potential flow theory and linearised free-surface governing equations using separation of variables (see, e.g. Zhang et al. Reference Zhang, Zhou and Ding2023). Unlike bubble-free drops, BC sessile drops require the consideration of gas compressibility (Ding & Bostwick Reference Ding and Bostwick2022a ), described by a polytropic process
\begin{align} \left (p_{g0}+p_g(t)\right ){\left (v_{g0}+v_g(t)\right )}^\kappa =p_{g0}{v_{g0}}^\kappa , \end{align}
where
$p_g$
is the pulsation part of the gas pressure in the bubble,
$v_g$
is the change in bubble volume,
$\kappa$
is the polytropic exponent and
$p_{g0}$
and
$v_{g0}$
denote equilibrium values of gas pressure and bubble volume, respectively. The polytropic process (2.1) can be incorporated into the free-surface governing equation for the inner surface
$\varGamma _{2}$
, where
${p}_g$
and
${v}_g$
are assumed to be infinitesimally small (see Appendix A for a derivation) (see also Ding & Bostwick Reference Ding and Bostwick2022a
).
It is worth noting that, although we account for the gas compressibility of the inner bubble, its high equilibrium pressure
$p_{g0}$
minimises this effect on the first few axisymmetric modes (Ding & Bostwick Reference Ding and Bostwick2022a
). This assumption is reasonable due to Earth’s high atmospheric pressure (
${\sim}2\times 10^3$
). For simplicity, we set an unrealistically large equilibrium pressure
$p_{g0}=10^{10}$
in our calculations. Under this condition, the differences between the results obtained at such high pressure and those at atmospheric pressure are negligible (see Appendix B), making the results applicable to practical situations. At this high pressure, the effect of gas compressibility on the frequency spectrum is negligible, allowing the inner bubble to be treated as incompressible. Note that non-axisymmetric modes are unaffected by gas compressibility because the bubble volume remains constant in those modes.
Aside from gas compressibility, deriving the governing equations for BC sessile drops is similar to those for bubble-free sessile drops (Lyubimov et al. Reference Lyubimov, Lyubimova and Shklyaev2006; Bostwick & Steen Reference Bostwick and Steen2014; Sharma & Wilson Reference Sharma and Wilson2021; Ding & Bostwick Reference Ding and Bostwick2022b
; Zhang et al. Reference Zhang, Zhou and Ding2023). Applying normal modes
$\psi (\boldsymbol {x}, t)=\phi (r, z) \mathrm {e}^{\textrm {i} \lambda t} \mathrm {e}^{\textrm {i} l \varphi }$
and
$\eta _{1,2}(s, \varphi , t)=y_{1,2}(s) \mathrm {e}^{\textrm {i} \lambda t} \mathrm {e}^{{\textrm {i}} l \varphi }$
in cylindrical coordinates
$(r,\varphi ,z)$
to the potential
$\psi$
(yielding velocity
$\boldsymbol {u} =\boldsymbol {\nabla } \psi$
) and perturbations
$\eta _{1,2}$
of the outer and inner surfaces (see figure 2
$a$
), we obtain a functional eigenvalue problem for oscillations of a BC sessile drop, expressed in axisymmetric coordinates
$(r,z)$
with a curvilinear coordinate
$s$
(see Appendix A for a detailed derivation)
\begin{align} \frac {1}{r}\frac {\partial }{{\partial r}}\left ( {r\frac {{\partial \phi }}{{\partial r}}} \right ) + \frac {{{\partial ^2}\phi }}{{\partial {z^2}}} - \frac {{{l^2}}}{{{r^2}}}\phi = 0\quad [D], \\[-24pt] \nonumber \end{align}
\begin{align} \frac {{\partial \phi }}{{\partial n}} = 0\quad [\partial {D^s}], \\[-2pt] \nonumber \end{align}
\begin{align} {\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )^{\prime \prime }} + \frac {{r^{\prime}}}{r}{\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )^\prime } + \left ( {2\sin ^2 \alpha - \frac {{{l^2}}}{{{r^2}}}} \right )\frac {{\partial \phi }}{{\partial n}} = - {\lambda ^2}\phi \quad [ \varGamma _{1}], \\[-24pt] \nonumber \end{align}
\begin{align} {\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )^{\prime \prime }} + \frac {{r^{\prime}}}{r}{\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )^\prime } + \left ( {\frac {2\sin ^2 \alpha }{\tau ^2}- \frac {{{l^2}}}{{{r^2}}}} \right )\frac {{\partial \phi }}{{\partial n}} -\frac {\kappa p_{g0}}{v_{g0}}\int _{\varGamma _2} {\frac {{\partial \phi }}{{\partial n}} {{\textrm {d}}}S }= - {\lambda ^2}\phi \quad [ \varGamma _{2}], \end{align}
\begin{align} \int _{\varGamma _1+\varGamma _2} {\frac {{\partial \phi }}{{\partial n}} {{\textrm {d}}}S } = 0\quad [\partial { D^l}], \end{align}
\begin{align} { {{{\pm \left ( {\frac {{\partial \phi }}{{\partial n}}} \right )}^\prime } + \chi _{1,2} \frac {{\partial \phi }}{{\partial n}} = 0} }\quad [\gamma _{1,2}]. \end{align}
Here,
$\phi$
and
$\partial \phi /\partial n$
are the potential and its normal derivative, respectively,
$\lambda$
is the scaled frequency and the prime denotes the derivative with respect to arc length
$s$
. The sign of the first term on the left-hand side of (2.2
$f$
) is positive for the outer CL
$\gamma _{1}$
and negative for the inner CL
$\gamma _{2}$
.
Equation (2.2
$a$
) is Laplace’s equation for potential flow in the bulk
$D$
, (2.2
$b$
) is the no-penetration condition on the solid surface, (2.2
$c$
) and (2.2
$d$
) are the free-surface governing equations for the drop surface
$\varGamma _1$
and bubble surface
$\varGamma _2$
, respectively, derived from the kinematic condition and normal pressure balance (Zhang et al. Reference Zhang, Zhou and Ding2023), (2.2
$e$
) is the volume conservation condition and (2.2
$f$
) describes the CL conditions.
In (2.2
$d$
), the term
$({\kappa p_{g0}}/{v_{g0}})\int _{\varGamma _2} {{{\partial \phi }}/{{\partial n}} {{\textrm {d}}}S }$
consists of a coefficient
${\kappa p_{g0}}/{v_{g0}}$
and an integral
$\int _{\varGamma _2} {{{\partial \phi }}/{{\partial n}} {{\textrm {d}}}S }$
, which represents the bubble’s volume change. This term characterises the effect of the gas compressibility. In the high-pressure limit
$p_{g0}\to \infty$
, the coefficient
${\kappa p_{g0}}/{v_{g0}}$
tends to infinity, acting similarly to a Lagrange multiplier that enforces the volume conservation condition
$\int _{\varGamma _2} {{\partial \phi }/{{\partial n}} {{\textrm {d}}}S }=0$
. Any deviation from this condition leads to a significantly large value for the associated term, causing the equation to approach infinity. Consequently, the system converges to a state where volume conservation is strictly satisfied. Thus, in the high-pressure limit, equation (2.2
$d$
) simplifies by removing the compressibility term entirely, as the volume conservation constraint is strictly enforced. The resulting equation can be viewed as the limiting case where bubble volume is perfectly conserved.
In (2.2
$f$
),
$\chi _{1,2}$
are boundary parameters dependent on the geometry at the CLs
$\gamma _{1,2}$
. For free CLs,
$\chi _1=-\cos \alpha$
and
$\chi _2=-{\cos \alpha }/{\tau }$
, while for pinned CLs,
$\chi _1$
and
$\chi _2$
are both set to
$+\infty$
(Zhang et al. Reference Zhang, Zhou and Ding2023). Substituting
$\chi _{1,2}$
into (2.2
$f$
) yields the free CL conditions
\begin{align} {{\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )}^\prime } -\cos \alpha \frac {{\partial \phi }}{{\partial n}} =0\quad [\gamma _{1}],\quad {{\left ( {\frac {{\partial \phi }}{{\partial n}}} \right )}^\prime } + \frac {\cos \alpha }{\tau } \frac {{\partial \phi }}{{\partial n}} = 0\quad [\gamma _{2}], \end{align}
and the pinned CL conditions
\begin{align} \frac {{\partial \phi }}{{\partial n}} = 0 \quad [\gamma _{1,2}]. \end{align}
Although classical spectral methods (Bhattacharya Reference Bhattacharya2016; Sumanasekara & Bhattacharya Reference Sumanasekara and Bhattacharya2017; Shiryaev Reference Shiryaev2020) have been successfully applied to spherical compound drops, extending these methods to solve problem (2.2) for BC sessile drops is challenging due to their complex geometry (see figure 2
$a$
). In this study, we use a BEM model, developed in our previous works (Zhang et al. Reference Zhang, Zhou and Ding2023; Zhang & Zhou Reference Zhang and Zhou2023), which handles arbitrary geometric boundaries, to numerically solve (2.2) for BC sessile drops. Minor extensions have been made to address the presence of two separate surfaces,
$\varGamma _1$
and
$\varGamma _2$
, and to incorporate the polytropic relation of gas compressibility in the free-surface governing equation (2.2
$d$
). The BEM model converts the functional eigenvalue problem (2.2) into a generalised matrix eigenvalue problem, with solutions providing the dimensionless frequencies
$\lambda$
and mode shapes
${\partial \phi }/{\partial n}$
. A brief procedure for solving (2.2) using the BEM model is outlined in Appendix C.
3. Mode identification and validation
Due to the presence of an inner bubble, BC sessile drops possess a new class of modes (i.e. plus modes) that are not present in bubble-free sessile drops (Bostwick & Steen Reference Bostwick and Steen2014). For bubble-free drops, it is known that modes are uniquely identified by the polar-azimuthal mode number pair
$\{k,l\}$
with
$k + l=\mathbb {N}_{{even}}$
(Bostwick & Steen Reference Bostwick and Steen2014). The mode numbers
$k$
and
$l$
determine the layer number of the mode shape (Chang et al. Reference Chang, Bostwick, Steen and Daniel2013)
\begin{align} n=(k-l)/2+1. \end{align}
For BC sessile drops, each mode has inner and outer surfaces with
$n_+$
and
$n_-$
layers, respectively (see figure 2
$a$
). Here, the subscripts ‘
$+$
’ and ‘
$-$
’ indicate inner and outer surfaces, respectively. By analysing the layer numbers
$n_\pm$
, we can calculate
$k_{\pm }$
using (3.1), which enables us to identify the modes of BC sessile drops with three mode numbers
$k_{+}$
,
$k_{-}$
and
$l$
. Note that plus and minus modes are dominated by deformations of the inner and outer surfaces, respectively (see property (iii) in § 1). This allows us to identify minus and plus modes using mode number pairs
$\{k_-,l\}_-$
and
$\{k_+,l\}_+$
, respectively. For example, the
$\{2,0\}_+$
mode indicates that the dominant deformation occurs on the inner surface with 2 layers (
$n_+=(k_+-l)/2+1=2$
).
Figure 3. Spectral lines of axisymmetric modes (
$l=0$
) for (
$a$
) free and (
$b$
) pinned BC semi-drops. In (
$a$
), solid and dashed lines represent minus and plus modes, respectively. In (
$b$
), spectral lines of plus and minus modes are discontinuous due to avoided crossings. Light blue and light pink arrows indicate spectral trends of minus and plus modes, respectively. Insets on the right show mode shapes at
$\tau =0.8$
, with blue for minus and red for plus modes.
For free BC semi-drops (i.e. with
$\alpha =90^\circ$
and free CLs), each mode has the same number of layers (
$n_- = n_+$
) on inner and outer surfaces, resulting in equal wavenumbers (
$k_- = k_+$
), as shown in figure 3(
$a$
). However, in pinned BC semi-drops (with pinned CLs), the layer numbers
$n_-$
and
$n_+$
usually differ (figure 6
$c,\!d$
). Despite this, the number pair
$\{k_\pm , l\}_\pm$
can still identify plus and minus modes, with spectral trends (light blue and light pink arrows in figure 3
$b$
) similar to spectral lines of minus and plus modes in free BC semi-drops (solid and dashed lines in figure 3
$a$
). In avoided crossings (figure 3
$b$
), the distinction between plus and minus modes becomes obscured, as will be discussed in § 4.1.2.
Figure 4. Comparison between the theoretical results (lines) for concentric spherical BC drops (Bhattacharya Reference Bhattacharya2016) and the BEM model results (symbols) for free BC semi-drops. The dimensionless frequencies
$\lambda _\pm$
in (
$a$
) and amplitude ratios
$\xi _\pm$
in (
$b$
) are normalised by their respective asymptotic solutions,
$\lambda ^a_\pm$
and
$\xi ^a_\pm$
.
As previously discussed, the frequency spectrum of free BC semi-drops aligns with that of concentric spherical BC drops for modes with
$k_\pm +l=\mathbb {N}_{{even}}$
due to plane symmetry (see figure 1
$a$
). Figure 4 shows excellent agreement between the theoretical results of Bhattacharya (Reference Bhattacharya2016) for concentric spherical BC drops and numerical results of our BEM model for free BC semi-drops in the
$\{2,0\}_\pm$
modes, thus validating our model. The observed spectral degeneracy for non-axisymmetric modes (see figure 1
$\!f$
) further confirms the accuracy of the BEM model, which will be discussed in the next section. In addition, the agreement of our model with the results obtained from numerical simulations of mode beating also supports our model (see § 4.3.2).
4. Results and discussion
In this section, we analyse the frequency spectra of BC sessile drops obtained by the BEM model. Specifically, we examine how CL pinning and varying
$\alpha$
affect the spectral lines of plus and minus modes. To illustrate this, we investigate the changes in the five properties of free BC semi-drops (see § 1) when varying
$\alpha$
from
$90^\circ$
or pinning the CL.
4.1. Bubble-containing semi-drops
We begin our discussion with free BC semi-drops because they inherit half the modes (
$k + l=\mathbb {N}_{{even}}$
) of concentric spherical BC drops (Bhattacharya Reference Bhattacharya2016), as discussed in § 1.
4.1.1. Breaking of spectral degeneracy
In free BC semi-drops, modes with the same
$k_\pm$
but different
$l$
have identical frequencies, a phenomenon known as spectral degeneracy (Bostwick & Steen Reference Bostwick and Steen2014). This degeneracy can be broken by varying
$\alpha$
from
$90^\circ$
, pinning the CL, or applying gravity for bubble-free sessile drops (Bostwick & Steen Reference Bostwick and Steen2014; Zhang et al. Reference Zhang, Zhou and Ding2023). Similarly, pinning the CL in free BC semi-drops also breaks spectral degeneracy.
Figure 5. Frequency spectrum for (
$a$
,
$b$
) free BC semi-drops and (
$c$
,
$d$
) pinned BC semi-drops: (
$a$
,
$c$
) even and (
$b$
,
$d$
) odd wavenumbers
$k_\pm$
and
$l$
. In (
$a,\!b$
), blue solid and red dashed lines indicate minus and plus modes, respectively. In (
$c,\!d$
), yellow circular areas highlight the avoided crossings of spectral lines.
Figure 5(
$a,\!b$
) shows the spectral lines of minus and plus modes in free BC semi-drops, with even
$l$
and odd
$l$
plotted separately. Spectral degeneracy manifests itself as: for a fixed
$k_-$
(
$k_+$
), there are
$\lfloor k_\pm /2\rfloor +1$
numbers of minus (plus) modes corresponding to the same spectral line, where
$\lfloor \rfloor$
denotes the floor function. For example, the spectral line with
$k_-=5$
(see figure 5
$b$
) corresponds to three minus modes with the same frequency:
$\{5,1\}_-$
,
$\{5,3\}_-$
and
$\{5,5\}_-$
. However, the
$\{1,1\}_\pm$
modes do not degenerate with other modes and always have zero frequency (figure 5
$b$
) due to horizontal translation invariance, similar to bubble-free drops (Zhang et al. Reference Zhang, Zhou and Ding2023).
In bubble-free drops (Bostwick & Steen Reference Bostwick and Steen2014), the breaking of spectral degeneracy by CL pinning causes spectral lines to split, with higher frequencies for higher
$l$
modes. We find that pinning the CL in free BC semi-drops also results in a similar pattern of spectral splitting. Figure 5(
$c,\!d$
) shows the frequency spectrum of pinned BC semi-drops. It is evident that CL pinning increases frequencies and causes spectral splitting, with higher
$l$
modes having higher frequencies. Additionally, near or avoided crossings of spectral lines for the same
$l$
are observed (highlighted in yellow in figure 5
$c,\!d$
), which will be discussed further in the next section.
4.1.2. Topological change of spectral lines: avoided crossing of frequencies
Figure 6. (
$a,\!b$
) Frequencies
$\lambda$
and amplitude ratios
$\xi$
for the minus mode
$\{6,0\}_-$
and plus mode
$\{2,0\}_+$
in pinned BC semi-drops with
$\tau \in [0.35,0.55]$
. The upper and lower spectral lines avoid crossing each other. The light blue and light red arrows indicate the spectral trends for the
$\{6,0\}_-$
and
$\{2,0\}_+$
modes, respectively. (
$c,\!d$
) Layer numbers
$n_{\pm }$
on the inner and outer surfaces for the (
$c$
) lower and (
$d$
) upper spectral lines. (
$e$
) Mode shapes at
$\tau =0.38$
, 0.45 and 0.52, where blue indicates the minus mode, red the plus mode and white the coupled modes in the avoided crossing.
The CL pinning not only increases frequencies but also changes the topology of spectral lines for free BC semi-drops, resulting in avoided crossings (see figure 5). Figure 6(
$a$
) shows that, before an avoided crossing, the lower (upper) spectral line corresponds to the minus (plus) mode, while after the crossing, these correspondences switch. This results in discontinuous spectral trends for minus and plus modes, such as the U-shaped trend for the
$\{6,0\}_-$
mode and the initially flat, then decreasing trend for the
$\{2,0\}_+$
mode in pinned BC semi-drops (see light red and light blue arrows in figures 3
$b$
and 6
$a$
).
From another perspective, a spectral line can be divided into regions with different growth rates with respect to
$\tau$
, separated by an avoided crossing. For instance, the lower spectral line in figure 6(
$a$
) shows an almost zero growth rate before the avoided crossing at
$\tau \approx 0.45$
, and then shifts to a region of rapid negative growth. A similar phenomenon is observed in coupled spherical caps (see figure 6 of Bostwick & Steen Reference Bostwick and Steen2013). Avoided crossings are common in coupled oscillator systems, such as the coupled spring oscillator (Novotny Reference Novotny2010). In our case, the inner and outer surfaces of BC sessile drops can be viewed as two sub-oscillators, thus forming a coupled oscillator system.
To understand the shifts in the plus and minus modes along a spectral line, we examine the mode shapes near an avoided crossing. Figure 6(
$b$
) shows the amplitude ratio
$\xi$
versus
$\tau$
for the two spectral lines in figure 6(
$a$
). As
$\tau$
increases, the amplitude ratio for the lower spectral line increases from the light blue region (
$\xi \lesssim 1/2$
) to the light red region (
$\xi \gtrsim 2$
), whereas the opposite is true for the upper spectral line. The middle region (
$1/2 \lesssim \xi \lesssim 2$
) corresponds to the avoided crossing around
$\tau = 0.45$
, where the deformations of inner and outer surfaces are comparable (see mode shapes with
$\tau =0.45$
in figure 6
$e$
).
Recall that, for minus modes,
$\xi \lt 1$
(outer surface amplitude greater), and for plus modes,
$\xi \gt 1$
(inner surface amplitude greater) (see property (iii) in § 1). We adjust this property slightly:
$\xi \gg 1$
indicates plus modes, and
$\xi \ll 1$
indicates minus modes. The mode number pair
$\{k_\pm , l\}_\pm$
is determined by the layer number
$n_\pm$
on the surface with greater deformation. Figure 6(
$c,\!d$
) shows the layer number
$n_\pm$
of modes along the upper and lower spectral lines. For free BC semi-drops,
$n_+=n_-$
always holds, consistent with concentric spherical BC drops (Bhattacharya Reference Bhattacharya2016). In contrast, for pinned BC semi-drops, the layer numbers on inner and outer surfaces sometimes differ, as shown in figure 6(
$c,\!d$
). Before the avoided crossing, modes on the lower spectral line exhibit greater deformation on the outer surface (
$\xi \ll 1$
) with
$n_-=4$
layers (i.e.
$k_-=6$
), and greater deformation on the inner surface (
$\xi \gg 1$
) with
$n_+=2$
layers (i.e.
$k_+=2$
) after the avoided crossing. Thus, it corresponds to the
$\{6,0\}_-$
mode before and the
$\{2,0\}_+$
mode after the avoided crossing. The upper spectral line shows the opposite behaviour. In the avoided crossing, both coupled modes have
$n_- = 4$
and
$n_+ = 2$
, corresponding to the spectral trends of
$\{6,0\}_-$
and
$\{2,0\}_+$
modes.
Most interestingly, the two modes in the avoided crossing cannot be classified as either plus or minus modes due to their comparable inner and outer surface deformations (
$\xi \approx 1$
). Additionally, their frequencies are very close but not identical, and they exhibit opposite directions of deformation on the inner and outer surfaces (figure 6
$e$
). In fact, these modes are coupled, similar to coupled oscillators (Novotny Reference Novotny2010), and are likely to exhibit a surprising beating phenomenon, discussed in detail in § 4.3.
In summary, we examine the topological changes between the spectra of free BC semi-drops and pinned BC semi-drops. We find that pinning the CL in free BC semi-drops causes the spectral lines to exhibit avoided crossings. Unlike in free BC semi-drops, the spectral lines of plus and minus modes for pinned BC semi-drops are no longer continuous but show spectral trends (see figure 3
$b$
). These trends still maintain the original spectral shapes: the spectral trends of plus modes are U-shaped, while those of minus modes are initially flat and then decrease rapidly. We also investigate mode transitions along spectral lines with an avoided crossing and find that plus and minus modes on continuous spectral lines transform into each other before and after the avoided crossing.
4.1.3. Velocity and pressure fields for modes near an avoided crossing
As discussed in § 4.1.2, the third zonal mode for pinned BC semi-drops can be classified as a minus mode
$\{6,0\}_-$
before the avoided crossing (around
$\tau =0.45$
), characterised by significant outer surface deformation with
$n_-=4$
layers. After the avoided crossing, it transitions into a plus mode
$\{2,0\}_+$
, where the inner surface deformation dominates, with
$n_+=2$
layers. In contrast, the fourth mode behaves in the opposite manner to the third mode. To further investigate the mode transitions between the minus and plus modes along the continuous spectral lines near the avoided crossing, we present the pressure and velocity fields for the modes near the avoided crossing (figure 7). Since the BEM model is established on the boundary, it provides only boundary values for the potential (pressure) and its normal derivative (velocity). Additional calculations are required to evaluate the pressure and velocity in the interior region (see step (vi) in Appendix C).
Figure 7. Velocity and potential (pressure) fields for zonal modes of pinned semi-drops with
$\tau =[0.38, 0.45, 0.52]$
corresponding to figure 6(
$a$
). The top row represents the fourth mode (higher frequencies), and the bottom row shows the third mode (lower frequencies). Plus modes
$\{2,0\}_+$
(upper-left and lower-right) and minus modes
$\{6,0\}_-$
(lower-left and upper-right) each have similar flow and pressure fields, while the middle modes at
$\tau =0.45$
show intermediate characteristics. In these contour plots, larger values are represented by purplish red, and smaller values by bright blue.
Figure 7 shows the velocity and pressure fields for typical modes near this avoided crossing (see figure 6
$a$
). For the third mode at
$\tau = 0.38$
(before the avoided crossing), the pressure field shows a steep gradient near the outer surface (denser contours in the bottom-left panel), resulting in greater velocity and deformation there. After the avoided crossing, at
$\tau =0.52$
, the pressure gradient shifts to the inner surface (bottom-right panel), causing a higher velocity near the inner surface. At the avoided crossing (
$\tau =0.45$
), the pressure gradient becomes more uniformly distributed, producing only slightly higher velocities at the inner surface. Similar pressure and velocity field characteristics appear for the fourth mode, but opposite to the third mode before and after the avoided crossing. As long as the modes are categorised as plus or minus, they exhibit similar flow patterns, even between the third and fourth modes. This suggests that the classification into plus and minus modes effectively captures the mode dynamics.
Figure 8 provides a more detailed view of the pressure and velocity transitions between plus and minus modes for
$\tau = [0.43, 0.44, 0.45, 0.46, 0.47]$
within the avoided crossing. For the third mode, the velocities near the vertices of the inner and outer surfaces are always opposite, while for the fourth mode, they are the same. As
$\tau$
increases, the maximum velocity (or deformation) in the third mode shifts from the outer to the inner surface, while in the fourth mode, it moves from the inner to the outer surface, illustrating the transition between minus and plus modes. These transitions highlight the complementary behaviour of the third and fourth modes, with opposing patterns before and after the avoided crossing. This clear distinction reinforces the effectiveness of the plus and minus classification in capturing the mode dynamics.
Figure 8. Velocity and pressure field plots, similar to figure 7, showing results for avoiding crossing at
$\tau = [0.43, 0.44, 0.45, 0.46, 0.47]$
. Panels (
$a$
) and (
$b$
) represent the third and fourth modes, respectively. As
$\tau$
increases, the maximum velocity of the third mode in (
$a$
) shifts from the outer surface to the inner surface, while the opposite occurs for the fourth mode in (
$b$
).
4.2. Frequency spectra of BC sessile drops on hydrophilic and hydrophobic surfaces
The previous section explores the frequency spectrum of BC semi-drops with
$\alpha =90^\circ$
. Here, we briefly examine the effects of varying
$\alpha$
on the frequency spectrum, reflecting realistic scenarios for aqueous drops on hydrophilic (
$\alpha \lt 90^\circ$
) or hydrophobic (
$\alpha \gt 90^\circ$
) surfaces. Note that, for
$\alpha \lt 90^\circ$
, the ratio
$\tau$
is limited to a maximum value of
$(1-\cos \alpha )/(1+\cos \alpha )$
, at which point the apexes of inner and outer surfaces contact each other. For
$\alpha \ge 90^\circ$
,
$\tau$
ranges between 0 and 1.
Figure 9 shows the
$l=0$
and
$l=1$
spectra of BC sessile drops for
$\alpha =30^\circ$
,
$60^\circ$
,
$120^\circ$
and
$150^\circ$
. Avoided crossings are present in all spectra, for both free and pinned CLs (highlighted circular areas in yellow). Thus, crossings of spectral line are specific to free BC semi-drops (see figure 3
$a$
), whereas avoided crossings are generally observed in all cases except for free BC semi-drops. In addition, varying
$\alpha$
affects the spectral trends of plus and minus modes. The term “spectral trend” is used instead of “spectral line” due to the presence of avoided crossings. Specifically, minus modes display a nearly horizontal spectral trend, while plus modes show a U-shaped trend with only the left half, without rising to infinity (figure 9). For lower-order minus modes without avoided crossings, their spectral lines are nearly horizontal. Furthermore, as
$\tau \to 0$
, minus modes converge to those of bubble-free sessile drops (Zhang et al. Reference Zhang, Zhou and Ding2023) (see blue dots at
$\tau =0$
in figure 9).
Figure 9. Frequency spectra of BC sessile drops with either a free or pinned CL for
$\alpha = 30^\circ$
,
$60^\circ$
,
$120^\circ$
and
$150^\circ$
and
$l = 0,1$
, where the horizontal axis is
$\tau$
and the vertical axis is the dimensionless frequency
$\lambda$
. In the
$l = 1$
spectra with free CLs (second column), red lines indicate zero-frequency rocking modes
$\{1,1\}_\pm$
. Light red arrows highlight the spectral trends of plus modes, and yellow circular areas mark avoided crossings. The blue dots at
$\tau = 0$
represent modes of bubble-free sessile drops (Zhang et al. Reference Zhang, Zhou and Ding2023).
In bubble-free sessile drops with free CLs, the Noether or rocking mode
$\{1, 1\}$
always has zero frequency due to the drops’ translational invariance, regardless of gravity (Zhang et al. Reference Zhang, Zhou and Ding2023; Zhang & Zhou Reference Zhang and Zhou2023). For drops with pinned CLs, however, this mode has a non-zero frequency. The rocking modes of BC sessile drops show similar behaviour: with free CLs, both plus and minus modes
$\{1, 1\}_\pm$
have zero frequency regardless of
$\alpha$
or
$\tau$
. With pinned CLs, the
$\{1,1\}_-$
mode has a small, nearly constant frequency and does not experience avoided crossings with plus modes, while the
$\{1,1\}_+$
mode undergoes avoided crossings with several minus modes and generally shows a spectral trend that decreases with increasing
$\tau$
.
Overall, varying
$\alpha$
from
$90^\circ$
causes spectra to show avoided crossings for both free and pinned CLs, similar to the effect of CL pinning in BC semi-drops. This suggests that avoided crossings are widespread in most cases, except for free BC semi-drops. Additionally, varying
$\alpha$
affects the shapes of spectral lines or trends for minus and plus modes. In pinned BC semi-drops, minus modes show flat trends that then decrease rapidly, while plus modes exhibit U-shaped trends. With varying
$\alpha$
, minus modes can display flat trends, and plus modes generally show a decreasing trend.
4.3. Beating of two coupled modes in an avoided crossing
4.3.1. Theoretical prediction of mode beating
Figure 10. Beating of two coupled modes in an avoided crossing: superposition of modes with (
$a$
) equal and (
$c$
) unequal (3 : 1) amplitudes on outer surfaces. In (
$a,\!c$
), red and blue solid lines represent the heights of apexes on outer and inner surfaces, respectively, with black dashed lines indicating envelope waves. Panels (
$b,\!d$
) show the superposition of mode shapes at maximum inner (
$t_1$
and
$t_3$
) and outer (
$t_2$
and
$t_4$
) amplitudes.
As described in § 4.1.2, two coupled modes in an avoided crossing have comparable amplitudes on inner and outer surfaces (
$\xi \approx 1$
) but in opposite directions (see figure 6
$e$
). These modes have natural frequencies that are very close but not identical, denoted as
$\lambda _1$
and
$\lambda _2$
(with
$\lambda _2 \gt \lambda _1$
). This small frequency difference allows both modes to be excited simultaneously. When both modes are excited with equal amplitude on the outer surfaces, mode beating occurs due to their superposition, as shown in figure 10(
$a$
).
Assuming the outer amplitudes of the coupled modes initially align, their outer amplitudes add together according to the law of superposition, while the inner amplitudes nearly cancel out due to their opposing directions. As the modes vibrate with a small frequency difference
$\delta \lambda =\lambda _2-\lambda _1$
, the phase difference between the outer surfaces evolves, causing them to become out of phase and cancel each other out, while the inner surfaces do the opposite (figure 10
$a$
). This cyclical process results in periodic variations in the modulated amplitude, a phenomenon known as mode beating, which is also seen in musical tones (Roberts Reference Roberts2016). This can be described as the superposition of the heights of the outer surface apexes of the coupled modes
\begin{align} {a_o} = a\cos ({\lambda _2}t) + a\cos ({\lambda _1}t) = 2a\underbrace {\cos \left ( {\frac {{{\lambda _2} + {\lambda _1}}}{2}t} \right )}_{{{carrier}}}\underbrace {\cos \left ( {\frac {{{\lambda _2} - {\lambda _1}}}{2}t} \right )}_{{{envelope}}}, \end{align}
where
$a_o$
is the displacement of the outer surface apex, and
$a$
is the outer amplitude of the modes. Equation (4.1) illustrates how two harmonic waves with equal amplitude
$a$
and slightly different frequencies combine into a carrier wave of frequency
$(\lambda _2+\lambda _1)/2$
, modulated by an envelope wave of frequency
$(\lambda _2-\lambda _1)/2$
.
To illustrate mode beating, we linearly superimpose modes
$\bar {\textrm {B}}$
and
$\textrm {B}$
in an avoided crossing (figure 6
$a$
) for a pinned BC semi-drop with
$\tau = 0.45$
. These modes have close frequencies
$\lambda _1 = 17.7359$
and
$\lambda _2 = 18.2605$
, and amplitude ratios
$\xi _1 = 1.9380$
and
$\xi _2 = 1.6977$
, where subscripts 1 and 2 correspond to modes
$\bar {\textrm {B}}$
and
$\textrm {B}$
, respectively. Figure 10(
$a$
) shows the time evolution of the linear superposition of modes
$\bar {\textrm {B}}$
and
$\textrm {B}$
with equal outer amplitudes of
$a=0.05$
. According to (4.1), the frequency of carrier waves (red and blue solid lines in figure 10
$a,\!c$
) is
$(\lambda _2+\lambda _1)/2=17.9982$
, while the frequency of envelope waves, or the beating frequency, is
$\lambda _b=(\lambda _2-\lambda _1)/2=0.2623$
. Thus, half of a beating period is
$\pi /\lambda _b=11.9771$
, corresponding to moments
$t_2$
and
$t_4$
in figure 10. Due to different amplitude ratios, the inner amplitudes do not fully cancel when the outer amplitudes are identical, as seen at moment
$t_2$
in figure 10(
$b$
). In practice, excited modes often have different amplitudes. Figure 10(
$c$
) shows the superposition of modes with unequal outer amplitudes, where the amplitudes are
$a = 0.005$
and
$a/3$
for modes
$\textrm {B}$
and
$\bar { {\textrm {B}}}$
, respectively. Despite this, the modulated amplitude still varies with the same beating frequency
$\lambda _b=0.2623$
. The mode shapes at moments
$t_3$
and
$t_4$
are shown in figure 10(
$d$
). For an intuitive understanding, the dynamic process of mode beating is available in supplementary movies 1 and 2. The next section conducts numerical simulations to confirm that mode beating can indeed occur.
4.3.2. Numerical simulation of mode beating
Although mode beating has been theoretically demonstrated, its practical occurrence remains uncertain. To confirm it, we perform direct numerical simulations using in-house developed codes to solve the Navier–Stokes equations for the natural oscillations of a pinned BC semi-drop with
$\tau = 0.45$
(figure 11
$a$
). The simulations employ a conservative level set method coupled with a continuous surface force model to accurately capture two-phase incompressible fluid interfaces (Ding, Spelt & Shu Reference Ding, Spelt and Shu2007; Chiu & Lin Reference Chiu and Lin2011). The CL is pinned at a given location by artificially setting the CL velocity to zero during calculations. The density and viscosity of the ambient fluid are set to 0.1 % of those of the drop to minimise the effects of the ambient fluid. We set the liquid viscosity
$\mu$
to achieve an Ohnesorge number
$Oh=\mu /\sqrt {\varrho \sigma L_*}=0.001$
.
Figure 11. Numerical simulation of natural oscillations of a pinned BC semi-drop with
$\tau =0.45$
. (
$a$
) Initial configuration at
$t = 0$
, with green dashed lines indicating equilibrium surfaces, and red and blue colours representing liquid and gas phases, respectively. (
$b,\!c$
) Time evolutions of outer (
$b$
) and inner (
$c$
) amplitudes. In (
$c$
), mode beating clearly appears with three half-periods
$T_1$
,
$T_2$
and
$T_3$
.
In numerical simulations, two coupled modes
$\textrm {B}$
and
$\bar {\textrm {B}}$
are excited simultaneously by setting one of their shapes as the initial configuration (figure 11
$a$
). This causes the pinned BC semi-drop to freely vibrate driven by restoring force due to surface tension. Time evolutions of outer
$a_o$
and inner
$a_i$
amplitudes are shown in figures 11(
$b$
) and 11(
$c$
), respectively. The inner surface exhibits clear mode beating with viscous decay, while the outer surface does not. Measured half-periods of mode beating are
$T_1=10.81$
,
$T_2=11.01$
and
$T_3=10.80$
, which are close to theoretical prediction of
$2\pi /(\lambda _2-\lambda _1)=11.9771$
. Minor discrepancies may be due to high sensitivity of the beating period to frequencies of the coupled modes. Nevertheless, periodic variations in modulated amplitude confirm the existence of mode beating (see figure 11
$c$
).
Although we use one coupled mode as the initial configuration for simulations, other modes, especially lower ones, are also excited. Table 1 lists the inviscid frequencies
$\lambda$
, amplitude ratios
$\xi$
and mode number pairs
$\{k_\pm , l\}_\pm$
for the first seven modes calculated by the BEM model. Note that the third and fourth modes, coupled in an avoided crossing, cannot be classified as plus or minus modes.
Table 1. Inviscid frequencies
$\lambda$
, amplitude ratios
$\xi$
, mode numbers
$\{k_\pm ,l\}_\pm$
and mode types for the first seven modes of a pinned semi-drop with
$\tau =0.45$
, listed by increasing frequency. The third and fourth modes in an avoided crossing cannot be classified as either plus or minus due to their comparable inner and outer amplitudes (
$\xi \approx 1$
).
To analyse excited modes in simulations, we use FFT (“fft” function in MATLAB) to convert amplitude evolutions of inner and outer surfaces in figure 11(
$b,\!c$
) to the frequency domain, as shown in figure 12. Results show that numerical frequencies agree well with theoretical predictions for minus modes (see red line protrusions in figure 12). However, frequencies of the third and fourth coupled modes and the seventh plus mode are slightly smaller than inviscid predictions, likely due to the high sensitivity of modes with larger inner amplitudes to viscosity, resulting in a greater viscous frequency shift for modes with
$\xi \gt 1$
. Additionally, the amplitude ratio
$\xi$
of the first mode,
$5.86571\times 10^{-4}/(9.0238\times 10^{-3})=0.0650$
, closely matches the theoretical value of 0.0639 (see table 1), further validating our model.
Finally, we explain why mode beating is not observed on the outer surface but is evident on the inner surface (figure 11
$b,\!c$
). This occurs because additional excited modes, especially minus modes, are primarily reflected on the outer surface (
$\xi \ll 1$
) and conceal the mode beating. The inner surface, which remains relatively undeformed, filters out these minus modes.
5. Summary and conclusions
We have investigated the natural oscillations of BC sessile drops with either a free or pinned CL on a plane, featuring a centre-bottom trapped bubble. Unlike bubble-free drops, BC sessile drops exhibit plus modes with greater inner surface deformation and minus modes with greater outer surface deformation. As
$\tau \to 0$
, minus modes converge to those of bubble-free drops (see figure 9). Due to plane symmetry, free BC semi-drops inherit half of the modes (
$k+l=\mathbb {N}_{{even}}$
) of concentric spherical BC drops (Bhattacharya Reference Bhattacharya2016), five properties of which are discussed in § 1. Plus and minus modes differ in the ratios
$\xi$
of inner to outer surface deformations and in their spectral line shapes (see figure 1).
To explore BC sessile drops beyond free BC semi-drops, we numerically solve the functional eigenvalue problem (2.2) using the BEM model, which provides frequencies
$\lambda$
and mode shapes. Modes are classified as minus for
$\xi \ll 1$
and plus for
$\xi \gg 1$
, with their mode number pairs
$\{k_\pm , l\}_\pm$
determined by the layer numbers
$n_\pm$
with larger deformations (i.e. (3.1)). The BEM model is validated against theoretical results of Bhattacharya (Reference Bhattacharya2016) (figure 4) and further confirmed by numerical simulations of mode beating (figure 12).
Pinning the CL in free BC semi-drops alters the topology of spectral lines. Specifically, crossings of spectral lines between minus and plus modes become avoided crossings (figure 3). Avoided crossings are common in coupled oscillator systems (Cliffe & Winters Reference Cliffe and Winters1986; Ida Reference Ida2005; Faust et al. Reference Faust, Rieger, Seitner, Krenn, Kotthaus and Weig2012). In BC sessile drops, the inner and outer surfaces act as a coupled oscillator, leading to avoided crossings. These avoided crossings break the spectral lines into spectral trends for plus and minus modes (figure 3
$b$
). Varying
$\alpha$
from
$90^\circ$
also results in avoided crossings for both free and pinned CLs. Thus, avoided crossings are typical in most cases, except for free BC semi-drops. Additionally, both plus and minus modes in free BC semi-drops exhibit spectral degeneracy, which can be broken by CL pinning or varying
$\alpha$
, similar to bubble-free drops (Bostwick & Steen Reference Bostwick and Steen2014).
In an avoided crossing, two coupled modes cannot be classified as plus or minus modes due to their comparable inner and outer surface deformations (
$\xi \approx 1$
). They share the same layer numbers
$n_\pm$
, corresponding to the spectral trends involved in the avoided crossing (figure 6). Due to their close frequencies
$\lambda _1$
and
$\lambda _2$
and opposite vibrational directions, mode beating occurs when both modes are excited simultaneously (figure 10), akin to the principle of binaural beats (Lane Reference Lane1925). Theoretical predictions for the half-period of mode beating,
$2\pi /(\lambda _2 - \lambda _1)$
, match the numerical simulations for a pinned BC semi-drop with
$\tau = 0.45$
, confirming the presence of mode beating (figure 11).
In practical applications, opaque drops often conceal bubbles, making their detection challenging. This study reveals that inner bubbles alter the vibrational properties of sessile drops, offering potential for predicting bubble size in opaque sessile drops. This is analogous to the work of Sumanasekara & Bhattacharya (Reference Sumanasekara and Bhattacharya2017), which examined how bubble position in spherical drops influences the frequency spectrum.
Supplementary movies.
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.31.
Funding.
We are grateful for the support of the National Natural Science Foundation of China (grant numbers 12102426, 12372263 and 12388101).
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Derivation of eigenproblem (2.2) with gas compressibility
Different from previous studies (e.g. Bostwick & Steen Reference Bostwick and Steen2014; Zhang et al. Reference Zhang, Zhou and Ding2023), this work integrates the bubble’s compressibility into the functional eigenvalue problem (2.2) by applying the polytropic relation (2.1) to the gas bubble (Ding & Bostwick Reference Ding and Bostwick2022a
). The derivation of eigenproblem (2.2) is based on three assumptions: (i) the liquid is inviscid, allowing the use of potential flow theory, (ii) the amplitude is infinitesimal, permitting linearisation of the governing equations, and (iii) the gas is ideal and obeys the adiabatic law (2.1) with a constant exponent
$\kappa$
. All subsequent derivations are presented in dimensionless form, scaled by the characteristic quantities introduced in § 2.
Disturbances induce a potential flow
$\boldsymbol {\nabla } \psi$
within the liquid domain
$D$
, which satisfies the Laplace equation as well as the no-penetration condition at the solid boundary
$\partial D^s$
and the kinematic condition at the liquid boundary
$\partial D^l$
, expressed as follows:
\begin{align} {\nabla ^2}\psi = 0\quad [D],\quad \frac {{\partial \psi }}{{\partial n}} = 0\quad [\partial D^s],\quad \frac {{\partial \psi }}{{\partial n}} = \frac {{\partial \eta }}{{\partial t}}\quad [\partial D^l]. \end{align}
Additionally, the pressure field
$p_l$
in the liquid is related to the potential through the linearised Bernoulli equation
\begin{align} \frac {{\partial \psi }}{{\partial t}} + {p_l} = 0. \end{align}
On the liquid surfaces
$\partial D^l=\varGamma _1\cup \varGamma _2$
, the perturbation
$\eta$
alters the curvature, generating a Laplace pressure
$P$
across the surfaces
\begin{align} P={\Delta _\varGamma }\eta + \left ( {{k_1}^2 + {k_1}^2} \right )\eta , \end{align}
where
$\Delta _\varGamma$
is the Laplace–Beltrami operator depending on the equilibrium surfaces
$\varGamma _{1,2}$
, and
$k_1$
and
$k_2$
are their principal curvatures. This Laplace pressure acts outward from the liquid region along the outer normal vector
$\boldsymbol {n}$
(see figure 2
$a$
).
In dynamic equilibrium, the liquid pressure at the inner and outer surfaces are balanced by the corresponding Laplace pressures on these surfaces. This can be expressed as
\begin{align} -p_l = P\quad [\varGamma _{1}],\quad p_g-p_l = P\quad [\varGamma _{2}], \end{align}
where
$p_g$
is the gas pressure resulting from perturbations, as defined in the adiabatic law (2.1). For convenience, the adiabatic law for the bubble is recalled here as
\begin{align} \left (p_{g0}+p_g(t)\right ){\left (v_{g0}+v_g(t)\right )}^\kappa =p_{g0}{v_{g0}}^\kappa . \end{align}
Since
${p}_g$
and
${v}_g$
are assumed to be infinitesimally small, linearising (A5) with normal modes
$p_g=\tilde {p}_g\mathrm {e}^{{\textrm {i}} \lambda t} \mathrm {e}^{{\textrm {i}} l \varphi }$
and
$v_g=\tilde {v}_g\mathrm {e}^{{\textrm {i}} \lambda t} \mathrm {e}^{{\textrm {i}} l \varphi }$
yields (Ding & Bostwick Reference Ding and Bostwick2022a
)
\begin{align} \tilde {p}_g =-\frac {\kappa {p_{g0}}}{v_{g0}}\tilde {v}_g, \end{align}
where the bubble volume change
$\tilde {v}_g$
is
\begin{align} {\tilde {v}_g} =-\frac {1}{{\textrm {i}}\lambda }\int _{{\varGamma _2}} {\frac {{\partial \phi }}{{\partial n}}{{\textrm {d}}}S} = \left \{ \begin{array}{l@{\quad}l} -\frac {1}{{\textrm {i}}\lambda }\int {\frac {{\partial \phi }}{{\partial n}}r{{\textrm {d}}}s}, & l=0, \\[3pt] 0, & l\geq 1. \end{array} \right . \end{align}
The minus sign in (A7) indicates that perturbation
$y_2$
(see unit normal
$\boldsymbol {n}$
in figure 2
$a$
) decreases the bubble volume. Note that in non-axisymmetric modes (
$l \geq 1$
), perturbations do not change the volume, so the gas compressibility has no effect on these modes.
Substituting (A3) into (A4) and comparing with (A2) yields
\begin{align} \frac {{\partial \psi }}{{\partial t}} = \ell \eta _1\quad [\varGamma _{1}],\quad p_g+\frac {{\partial \psi }}{{\partial t}} = \ell \eta _2\quad [\varGamma _{2}], \end{align}
where
$\ell$
is the operator defined as
$\ell \equiv \Delta _\varGamma +({{k_1}^2 + {k_2}^2})$
.
Applying normal modes
\begin{align} \psi (\boldsymbol {x}, t)=\phi (r, z) \mathrm {e}^{{\textrm {i}} \lambda t} \mathrm {e}^{{\textrm {i}} l \varphi },\quad \eta _{1,2}(s, \varphi , t)=y_{1,2}(s) \mathrm {e}^{{\textrm {i}} \lambda t} \mathrm {e}^{{\textrm {i}} l \varphi }, \end{align}
to (A1
$c$
) and (A8), and comparing the results, we obtain
\begin{align} \ell \frac {{\partial \phi }}{{\partial n}} = - {\lambda ^2}\phi \quad [\varGamma _{1}],\quad \ell \frac {{\partial \phi }}{{\partial n}}-{\textrm {i}}\lambda {\tilde {p}_g}{ =} - {\lambda ^2}\phi \quad [\varGamma _{2}]. \end{align}
In cylindrical coordinates
$(r,z)$
, the operator
$\ell \equiv \Delta _\varGamma +({{k_1}^2 + {k_2}^2})$
is expressed as
\begin{align} \ell (\cdot )=(\cdot )^{\prime\prime}+\frac {r^{\prime}}{r}(\cdot )^{\prime}+\left ({k_1}^2+{k_2}^2-\frac {l^2}{r^2}\right )(\cdot ), \end{align}
where the prime denotes the derivative with respect to arc length
$s$
, and
$k_1$
and
$k_2$
are the principal curvatures of surfaces
$\varGamma _{1,2}$
. These curvatures are given by
\begin{align} k_1=k_2=\sin \alpha \quad [\varGamma _{1}],\quad k_1=k_2=\sin \alpha /\tau \quad [\varGamma _{2}]. \end{align}
Finally, by applying the normal modes (A9
$a$
) to the Laplace equation (A1
$a$
) and the no-penetration condition (A1
$b$
), we derive (2.2
$a$
) and (2.2
$b$
). Substituting (A11) with (A12
$a$
) into (A10
$a$
), we obtain the governing equation (2.2
$c$
) for the outer surface. Similarly, substituting (A11) with (A12
$b$
) and (A6) with (A7) into (A10
$b$
) gives the governing equation (2.2
$d$
) for the inner surface. The volume conservation condition (2.2
$e$
) is directly derived, and the CL condition (2.2
$f$
) results from linearising Young’s equation (see Myshkis et al. Reference Myshkis, Babskii, Kopachevskii, Slobozhanin and Tyuptsov1987, p. 126).
Appendix B. Effect of bubble static pressure
$p_{g0}$
on axisymmetric modes
Although gas compressibility of the inner bubble is included in the governing equations (2.2), its high equilibrium pressure
$p_{g0}$
minimises its effect on the first few axisymmetric modes (Ding & Bostwick Reference Ding and Bostwick2022a
). Given Earth’s high atmospheric pressure, this assumption is reasonable. The bubble’s static pressure
$p_{g0}$
consists of the ambient pressure
$p_a$
and the static Laplace pressure across the two liquid surfaces.
Using water as an example, a droplet’s CL radius on Earth is typically a few millimetres (e.g.
$R_1=2\times 10^{-3}\, \rm m$
) with surface tension
$\sigma =7.28\times 10^{-2}\, \rm N\, m^-{^1}$
, giving a characteristic pressure of
$\sigma /R_1=36.4\, \rm N \, m^-{^1}$
. Under standard atmospheric conditions, the dimensionless background pressure is
$p_a\simeq 2.78\times 10^{3}$
. The static Laplace pressure, expressed as
$2\sin \alpha (1+1/\tau )$
, is negligible by comparison, typically ranging from
$ 10^0$
to
$ 10^2$
. Thus, the bubble’s static pressure is approximately the constant background pressure of
$2.78\times 10^{3}$
. For partially wetting bubbles, Ding & Bostwick (Reference Ding and Bostwick2022a
) showed that, when
$p_{g0}$
exceeds a critical threshold, the dependence of
$\lambda$
on
$p_{g0}$
reaches a high-frequency plateau region (see their figure 6). Beyond this threshold, further increases in
$p_{g0}$
have minimal effect on
$\lambda$
, indicating that gas compressibility no longer affects that mode.
Figure 13. Frequencies
$\textrm {Re}(\lambda )$
of axisymmetric modes against static bubble pressure
$p_{g0}$
for (
$a$
)
$\alpha =90^\circ$
,
$\tau =0.45$
, (
$b$
)
$\alpha =90^\circ$
,
$\tau =0.8$
, (
$c$
)
$\alpha =70^\circ$
,
$\tau =0.4$
and (
$d$
)
$\alpha =120^\circ$
,
$\tau =0.5$
, contrasting pinned (solid black) and free (dashed red) CLs. In (
$a,b$
), spectral lines for volume and shape modes cross with free CLs, while pinned CLs exhibit avoided crossings. Each spectral line shows two plateaus separated by a rising region at the avoided crossing, consistent with Ding & Bostwick (Reference Ding and Bostwick2022a
). Stars mark the critical pressure for volume mode instability. Crossings or avoided crossings occur between a volume mode (blue arrow) and shape modes (green arrow). In (
$a$
), the third and fourth coupled modes in table 1 appear at
$p_{g0} \gtrsim 20$
.
To further explore these pressure thresholds for BC sessile droplets, we compute curves of
$\lambda$
against
$p_{g0}$
for different configurations, contrasting free and pinned CLs, as shown in figure 13. In these calculations, the dimensionless bubble volume is
$v_{g0}=\pi (2-\cos \alpha )(1+\cos \alpha )^2{\tau }^3/(3\sin ^3\alpha )$
, and the polytropic exponent is
$\kappa =1.4$
for (2.2
$d$
).
At low pressure
$p_{g0}$
, the bubble inside the BC sessile droplet is easily compressed, resulting in a low-frequency volume mode (with only one layer of inner surface deformation). This mode is unstable with
$\lambda ^2\lt 0$
(thresholds marked by stars in figure 13), and
$\lambda ^2$
increases as
$p_{g0}$
rises (see the discontinuous spectral trend indicated by the light blue arrow in figure 13). The gas compressibility also affects the shape modes (which exhibit multiple layers of inner surface deformation), primarily at avoided crossings (see figure 13
$a$
). As shown in figure 13, the avoided crossings of the first few modes generally occur around
$p_{g0} = 10$
(the pressure threshold). Therefore, with Earth’s background pressure of
$2.78\times 10^{3}$
far exceeding this threshold, the first few axisymmetric modes are largely unaffected by gas compressibility, allowing the bubble to be considered incompressible at this point.
Appendix C. Solving (2.2) via a BEM model
Similar to the cases of ring-constrained drops (Bostwick & Steen Reference Bostwick and Steen2013; Wagoner et al. Reference Wagoner, Garg, Harris, Ramkrishna and Basaran2021) and compound drops (Bhattacharya Reference Bhattacharya2016; Sumanasekara & Bhattacharya Reference Sumanasekara and Bhattacharya2017; Shiryaev Reference Shiryaev2020), the eigenvalue problem (2.2) involves two separate free surfaces,
$\varGamma _1$
and
$\varGamma _2$
(see figure 2). However, extending classical spectral methods to BC sessile drops is challenging due to the lack of suitable basis functions.
In this work, we extend the BEM model developed by Zhang et al. Reference Zhang, Zhou and Ding(2023) to solve (2.2) for BC sessile drops. The BEM model can handle arbitrary shapes and has been effectively used for gravity-flattened (Zhang et al. Reference Zhang, Zhou and Ding2023) and elongated pendant drops (Zhang & Zhou Reference Zhang and Zhou2023). The key idea is to transform the Laplace equation (2.2
$a$
) in bulk
$D$
into a boundary integral equation on boundary
$\partial D$
, so as to construct a generalised matrix eigenvalue problem (Pozrikidis Reference Pozrikidis2002; Ebrahimian, Noorian & Haddadpour Reference Ebrahimian, Noorian and Haddadpour2013; Zhang et al. Reference Zhang, Zhou and Ding2023; Zhang et al. Reference Zhang, Zhou and Ding2023). We outline the essential steps for constructing the BEM model as follows:
Figure 14. Discretisation of the boundary
$\partial D=\varGamma _{1}+\varGamma _{2}+\partial D^s$
into uniform line elements, with collocation points positioned at the centre of each element. The element numbering
$i$
increases along the arrows, in the direction of increasing arc length
$s$
(figure 2
$a$
).
-
(i) We discretise the boundary
$\partial D$
into line elements (figure 14), and express the boundary integral equation corresponding to (2.2
$a$
) as a matrix equation (Zhang et al. Reference Zhang, Zhou and Ding2023)(C1)where
\begin{align} {\unicode{x1D646}}^{\,l} {\boldsymbol {\phi } ^ * } = {\unicode{x1D643}}^{\,l} \boldsymbol {\phi },\quad l=0,1,\ldots , \end{align}
$\boldsymbol {\phi }$
and
$\boldsymbol {\phi }^*$
are column vectors of the potential
$\phi$
and its normal derivative
${\partial \phi }/{\partial n}$
at the collocation points, respectively. In (C1), the influence matrices
$\unicode{x1D646}^{\,l}$
and
$\unicode{x1D643}^{\,l}$
, which depend only on the boundary shape
$\partial D$
, are determined by a standard BEM formulation (Pozrikidis Reference Pozrikidis2002; Zhang et al. Reference Zhang, Zhou and Ding2023) for each value of
$l$
.
-
(ii) By partitioning the boundary
$\partial D=\partial D^s+\partial D^l$
into the solid surface
$\partial D^s$
and the outer/inner surfaces
$\partial D^l=\varGamma _{1}+\varGamma _{2}$
, we can recast (C1) as(C2)where the subscripts ‘
\begin{align} \left [ \begin{array}{l@{\quad}l} \unicode{x1D646}_{1,1}^{\,l} & \unicode{x1D646}_{1,2}^{\,l}\\[3pt] \unicode{x1D646}_{2,1}^{\,l}&\unicode{x1D646}_{2,2}^{\,l} \end{array} \right ]\left \{ {\begin{array}{l} {\boldsymbol {\phi } _L^*}\\[3pt] {\boldsymbol {\phi } _S^*} \end{array}} \right \} = \left [ {\begin{array}{l@{\quad}l} \unicode{x1D643}_{1,1}^{\,l}&\unicode{x1D643}_{1,2}^{\,l}\\[3pt] \unicode{x1D643}_{2,1}^{\,l} & \unicode{x1D643}_{2,2}^{\,l} \end{array}} \right ]\left \{ {\begin{array}{*{20}{c}} {{\boldsymbol {\phi } _L}}\\[3pt] {{\boldsymbol {\phi } _S}} \end{array}} \right \}, \end{align}
$S$
’ and ‘
$L$
’ denote the quantities evaluated at the solid and liquid surfaces, respectively, and the numerals in the subscripts of the influence matrices
$\unicode{x1D646}$
and
$\unicode{x1D643}$
specify the positions of submatrices. Applying the no-penetration condition
$\boldsymbol {\phi }_S^*=\boldsymbol {0}$
to (C2) yields(C3)
\begin{align} \left [\unicode{x1D643}_{1,2}^{\,l}\left (\unicode{x1D643}_{2,2}^{\,l}\right )^{-1} \unicode{x1D646}_{2,1}^{\,l} - \unicode{x1D646}_{1,1}^{\,l}\right ] \boldsymbol {\phi }_L^*=\left [\unicode{x1D643}_{1,2}^{\,l}\left (\unicode{x1D643}_{2,2}^{\,l}\right )^{-1} \unicode{x1D643}_{2,1}^{\,l} - \unicode{x1D643}_{1,1}^{\,l}\right ] \boldsymbol {\phi }_L. \end{align}
-
(iii) In (C3), the quantities
$\boldsymbol {\phi }_L$
and
$\boldsymbol {\phi }_L^*$
on the liquid surface can be further distinguished as(C4)where subscripts ‘1’ and ‘2’ correspond to outer surface
\begin{align} \boldsymbol {\phi }_L=\left \{\begin{array}{l} \boldsymbol {\phi }_1 \\[3pt] \boldsymbol {\phi }_2 \end{array}\right \} \quad \text {and} \quad \boldsymbol {\phi }_L^*=\left \{\begin{array}{l} \boldsymbol {\phi }_1^* \\[3pt] \boldsymbol {\phi }_2^* \end{array}\right \}, \end{align}
$\varGamma _{1}$
and inner surface
$\varGamma _{2}$
, respectively. For
$\varGamma _{1}$
, the free-surface equation (2.2
$c$
) can be discretised using a central finite difference scheme and then incorporates the free/pinned CL condition (2.2
$f$
) through a ghost point method, yielding a system of linear equations (Zhang et al. Reference Zhang, Zhou and Ding2023)(C5)where
\begin{align} \widetilde {\unicode{x1D646}}^{\,l}_{1} \boldsymbol {\phi }_1^*=-\lambda ^2 \unicode{x1D644}_{1} \boldsymbol {\phi }_1, \end{align}
$\unicode{x1D644}_1$
denotes the identity matrix. The matrix
$\widetilde {\unicode{x1D646}}^{\,l}_{1}$
represents the discretisation of the operator
$\ell$
and is determined similarly to that for bubble-free sessile drops (Zhang et al. Reference Zhang, Zhou and Ding2023).
-
(iv) For
$\varGamma _{2}$
, the discretisation of the free-surface governing (2.2
$d$
) must account for the gas compressibility. For axisymmetric modes (
$l=0$
), the last term on the left-hand side of (2.2
$d$
) is(C6)which can be approximated as
\begin{align} \frac {\kappa p_{g0}}{v_{g0}}\int _{\varGamma _2} {\frac {{\partial \phi }}{{\partial n}} {{\textrm {d}}}S }= \frac {\kappa p_{g0}}{v_{g0}}\int {\frac {{\partial \phi }}{{\partial n}}r{{\textrm {d}}}s}, \end{align}
(C7)Here,
\begin{align} \frac {\kappa p_{g0}}{v_{g0}}\int {\frac {{\partial \phi }}{{\partial n}}r{{\textrm {d}}}s}\simeq \frac {2\pi \kappa p_{g0}}{v_{g0}} \boldsymbol {r}_B\boldsymbol {\phi }_2^*. \end{align}
$\boldsymbol {r}_B$
is a row vector consisting of the radii of collocation points on
$\varGamma _{2}$
. Then, the free-surface governing equation (2.2
$d$
) with (C7) can be discretised similarly to (C5) as(C8)with
\begin{align} \widetilde {\unicode{x1D646}}^{\,l}_{2} \boldsymbol {\phi }_2^* - \unicode{x1D646}_C \boldsymbol {\phi }_2^*=-\lambda ^2 \unicode{x1D644}_2 \boldsymbol {\phi }_2, \end{align}
(C9)For non-axisymmetric modes (
\begin{align} \unicode{x1D646}_C=\frac {2\pi \kappa p_{g0}}{v_{g0}}\left [ {\begin{array}{*{20}{c}} {{{\bf {r}}_B}}\\ \vdots \\ {{{\bf {r}}_B}} \end{array}} \right ]. \end{align}
$l\geq 1$
),
$\unicode{x1D646}_C=\boldsymbol {0}$
since gas compressibility has no effect.
-
(v) Combining the inner and outer surfaces, we assemble the matrix equations (C5) and (C8) into
(C10)with
\begin{align} \widetilde {\unicode{x1D646}}^{\,l} \boldsymbol {\phi }_L^*=-\lambda ^2 \unicode{x1D644} \boldsymbol {\phi }_L , \end{align}
(C11)Here,
\begin{align} \widetilde {\unicode{x1D646}}^{\,l}=\left [\begin{array}{c@{\quad}c} \widetilde {\unicode{x1D646}}^{\,l}_{1} & \boldsymbol {0} \\ \boldsymbol {0} & \widetilde {\unicode{x1D646}}^{\,l}_{2} - \unicode{x1D646}_C \end{array}\right ], \quad \unicode{x1D644} =\left [\begin{array}{c@{\quad}c} \unicode{x1D644}_1 & \boldsymbol {0} \\ \boldsymbol {0} & \unicode{x1D644}_2 \end{array}\right ]. \end{align}
$\unicode{x1D644}$
and
$\unicode{x1D644}_2$
denote identity matrices. By comparing equations (C3) and (C10), we obtain a matrix eigenvalue problem(C12)whose solutions give frequencies
\begin{align} \left [\unicode{x1D643}_{1,2}^{\,l}\left (\unicode{x1D643}_{2,2}^{\,l}\right )^{-1} \unicode{x1D643}_{2,1}^{\,l} - \unicode{x1D643}_{1,1}^{\,l}\right ]\tilde {\unicode{x1D646}}^{\,l} \boldsymbol {\phi }_L^*=\lambda ^2 \left [\unicode{x1D646}_{1,1}^{\,l} - \unicode{x1D643}_{1,2}^{\,l}\left (\unicode{x1D643}_{2,2}^{\,l}\right )^{-1} \unicode{x1D646}_{2,1}^{\,l}\right ] \boldsymbol {\phi }_L^*, \end{align}
$\lambda$
and mode shapes
$\boldsymbol {\phi }_L^*$
. So far, equation (C12) does not include the volume constraint (2.2
$e$
). Similar to gas compressibility, volume conservation is automatically maintained in non-axisymmetric modes, so (2.2
$e$
) does not need to be considered for
$l\geqslant 1$
. For
$l=0$
, the problem (C12) with volume constraint can be transformed into a modified unconstrained eigenvalue problem by projecting it into a constrained space (see Zhang et al. Reference Zhang, Zhou and Ding2023).
-
(vi) Once
$\boldsymbol {\phi }_L^*$
are determined on the liquid boundary
$\partial D^l$
, the potential
$\boldsymbol {\phi }_L$
on the same boundary can be easily computed using (C10). From this, we can then calculate the potential
$\boldsymbol {\phi }$
and its normal derivative
$\boldsymbol {\phi }^*$
across the entire boundary using (C2). To compute the potential in the interior, we use the integral representation of (2.2
$a$
) with the source point located in the interior (Pozrikidis Reference Pozrikidis2002)(C13)where
\begin{align} \phi ({\boldsymbol {x}_0}) = \int _{\partial D} {\left [{G^l}({\boldsymbol {x}},{\boldsymbol {x}_0})\frac {{\partial \phi }}{{\partial n}}({\boldsymbol {x}}) - \frac {{\partial {G^l}({\boldsymbol {x}},{\boldsymbol {x}_0})}}{{\partial n}}\phi (\boldsymbol {x})\right ]r{\textrm {d}}s(\boldsymbol {x})}, \end{align}
$\boldsymbol {x}_0$
is the interior source point,
$\boldsymbol {x}$
is a boundary point and
${G^l}({\boldsymbol {x}},{\boldsymbol {x}_0})$
is the Green’s function corresponding to (2.2
$a$
), depending on
$l$
(Zhang et al. Reference Zhang, Zhou and Ding2023). The potential
$\boldsymbol {\phi }$
at the interior point is then obtained by evaluating the integral (C13) using the boundary values
$\boldsymbol {\phi }$
and
$\boldsymbol {\phi }^*$
. To calculate the gradient of the potential, the gradient form of the integral representation is used (see Pozrikidis Reference Pozrikidis2002, pp. 21–22).
