1. Introduction
Bubble and droplet coalescence is relevant in a wide range of industrial processes including injection printing (Lohse Reference Lohse2022), chemical reactors (Schlüter et al. Reference Schlüter, Herres-Pawlis, Nieken, Tuttlies and Bothe2021) and wastewater treatment (Terasaka et al. Reference Terasaka, Hirabayashi, Nishino, Fujioka and Kobayashi2011). When coalescence takes place on a solid substrate, the release of excess potential energy due to different total surface areas of the original bubbles or droplets and the merged one can cause the latter to jump off the substrate. The phenomenon was first documented for droplets on super-hydrophobic surfaces (Boreyko & Chen Reference Boreyko and Chen2009) and more recently for bubbles attached to solid surfaces (Soto et al. Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018). Evidently, the phenomenon is highly relevant in applications where bubble or droplet removal from a solid surface is essential for the continuation of a process, e.g. in drop-wise condensation (Miljkovic et al. Reference Miljkovic, Enright, Nam, Lopez, Dou, Sack and Wang2013), self-cleaning surfaces (Wisdom et al. Reference Wisdom, Watson, Qu, Liu, Watson and Chen2013), boiling (Zhou et al. Reference Zhou, Qi, Zhang, Wei, Yang and Cao2020) and water splitting (Zhang et al. Reference Zhang, Wang, Feng, Zhen, Liu, Cao, Zhao and Guo2024).
Recent interest in bubble removal from solid surfaces has been fuelled by its application in gas-evolving electrochemical processes, specifically green hydrogen production (Kadyk et al. Reference Kadyk, Bruce and Eikerling2016; Zhao et al. Reference Zhao, Ren and Luo2019; Bashkatov et al. Reference Bashkatov, Hossain, Mutschke, Yang, Rox, Weidinger and Eckert2022; Sepahi et al. Reference Sepahi, Verzicco, Lohse and Krug2024), where gas bubbles attached to an electrode can impede the reaction sites and deteriorate the energy efficiency of the electrolysis process (Vogt & Balzer Reference Vogt and Balzer2005; Angulo et al. Reference Angulo, van der Linde, Gardeniers, Modestino and Fernández Rivas2020). Indeed, in recent years, considerable effort has been devoted to the passive (e.g. Brussieux et al. Reference Brussieux, Viers, Roustan and Rakib2011; Lake et al. Reference Lake, Soto and Varanasi2022) or active (e.g. Koza et al. Reference Koza, Mühlenhoff, Żabiński, Nikrityuk, Eckert, Uhlemann, Gebert, Weier, Schultz and Odenbach2011; Bolaños-Jiménez et al. Reference Bolaños-Jiménez, Rossi, Rivas, Kähler and Marin2017) methods of promoting bubble detachment; this calls for a deep understanding of the underlying physics.
Bubble detachment off solid surfaces has been the subject of several studies in the past, in which the effects of Marangoni forces (Lubetkin Reference Lubetkin2002; Bashkatov et al. Reference Bashkatov, Babich, Hossain, Yang, Mutschke and Eckert2023; Park et al. Reference Park, Liu, Demirkır, van der Heijden, Lohse, Krug and Koper2023) and electric or magnetic fields (Fernandez et al. Reference Fernandez, Diao, Dunne and Coey2010; Baczyzmalski et al. Reference Baczyzmalski, Karnbach, Mutschke, Yang, Eckert, Uhlemann and Cierpka2017; Hossain et al. Reference Hossain, Bashkatov, Yang, Mutschke and Eckert2022) on the size and release frequency of the bubbles from an electrode have been examined, among others. Certain studies focused more directly on the effect of bubble coalescence on detachment. For instance, Fernandez et al. (Reference Fernandez, Maurer, Martine, Coey and Mobius2014) reported that the coalescence efficiency of micro-bubbles formed on a micro-electrode can significantly alter the bubble detachment patterns. Bashkatov et al. (Reference Bashkatov, Hossain, Mutschke, Yang, Rox, Weidinger and Eckert2022) showed that coalescence-induced jump off can cause bubble removal from a solid surface even in zero-gravity condition. Later on, Bashkatov et al. (Reference Bashkatov, Park, Demirkır, Wood, Koper, Lohse and Krug2024) demonstrated that bubble coalescence can lead to substantially earlier bubble departure and higher current densities as compared with a purely buoyancy-driven bubble removal. Interestingly, these authors also reported possible re-attachment of the already detached bubbles due to coalescence with those still attached to the surface.
When two coalescing bubbles merge, the released surface potential energy is mainly consumed by the sum of viscous dissipation and the ’adhesion energy’, which is the term used to describe the energy required to overcome the work of surface tension force at the three-phase contact line. The latter contribution is obviously dependent on the contact angle and it trends zero as the contact angle approaches 180
$^\circ$
(perfect wetting condition). According to this energy balance, Lv et al. (Reference Lv, Peñas, Le The, Eijkel, van den Berg, Zhang and Lohse2021) used scale analysis to derive a criterion for occurrence of coalescence-driven jump off as a function of the radii of the coalescing bubbles and physical fluid properties. While Lv et al. (Reference Lv, Peñas, Le The, Eijkel, van den Berg, Zhang and Lohse2021) considered a case of negligible adhesion energy, their work was complemented by Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025) who added this effect using additional experimental and numerical data. As a result, the criterion proposed by Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025) introduces two asymptotic detachment behaviours; at contact angles close to 180
$^\circ$
, the viscous dissipation is dominant and its relative magnitude in terms of the Ohnesorge number (
$\textit{Oh}=({\mu _l}/({(\rho _l R\sigma )^{1/2}})$
) determines occurrence of jump-off,
$\mu _l$
,
$\rho _l$
,
$\sigma$
and
$R$
being the liquid viscosity, liquid density, surface tension and bubble radius, respectively. The other asymptotic behaviour, which is not the subject of the present work, occurs at finite contact angles where adhesion energy is the main factor determining the jump-off. Central to the analysis presented by Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025) is a proposed scaling of viscous dissipation and the length scale
$\lambda$
over which it takes place.
The mechanism behind coalescence-induced jump off and its scaling laws have been studied in details mostly for droplets. Boreyko & Chen (Reference Boreyko and Chen2009), who first observed the phenomenon, also provided an estimate of the jump velocity using a scale analysis based on the assumption that the entire released potential energy is converted to the translational kinetic energy of the merged droplet. Wang et al. (Reference Wang, Yang and Zhao2011) followed up by further elaborating the energy conservation analysis to develop an analytical equation for the jump velocity. Liu et al. (Reference Liu, Sun, Sun and Ai2012) derived a dynamic equation for the shape of droplets after coalescence. These works set the stage for further experimental investigations, e.g. by Liang & Keblinski (Reference Liang and Keblinski2015), Wang et al. (Reference Wang, Liang, Jiang, Zheng, Lan and Ma2018), Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018) and Zhang et al. (Reference Zhang, Wang, Feng, Zhen, Liu, Cao, Zhao and Guo2024), of the subject matter. While these studies provided significant insight into the coalescence dynamics, detailed understanding of the flow field and energy budget requires access to direct numerical simulations (DNS) as well. Liu et al. (Reference Liu, Ghigliotti, Feng and Chen2014) conducted one such simulation campaign for droplets of same size at a contact angle of 0
$^\circ$
(which is analogous to 180
$^\circ$
in case of bubbles). These authors studied a wide range of
$ \textit{Oh}$
which enables observing the transition from the state at which the viscous and inertial forces are dominant (at high
$ \textit{Oh}$
) to the one at which capillary and inertial forces are dominant the (at low
$ \textit{Oh}$
). Their results showed, among others, that contrary to the assumption made by earlier researchers, only a small fraction of the released surface energy is converted to the translational kinetic energy of the merged droplet and the rest is converted to oscillatory kinetic energy before being dissipated through viscous forces. Other numerical studies have focused on questions such as the symmetry-breaking effects of the substrate (Chen & Lian Reference Chen and Lian2018) or the effect of the contact angle (Nam et al. Reference Nam, Kim and Shin2013; Cheng et al. Reference Cheng, Xu and Sui2016), contact line dynamics (Chen et al. Reference Chen, Lu and Tryggvason2019) and droplet size ratio (Wang & Ming Reference Wang and Ming2019). A recent review of the subject can be found in Liu et al. (Reference Liu, Zhao, Zheng, Cheng, Zhang and Tee2021).
Compared with droplets, only a few dedicated numerical studies have been reported for the coalescence-driven jump of bubbles. One should note that droplet and bubble coalescence differ in at least one fundamental aspect; while viscous dissipation is mostly contained within a droplet, it takes place outside a bubble, which rules out the radius as a convenient length scale for dissipation in the latter case. Iwata et al. (Reference Iwata, Zhang, Lu, Gong, Du and Wang2022) investigated coalescence of hydrogen bubbles in water focusing on developing a regime map for bubble jump off based on the equilibrium contact angle. They reported a critical contact angle of approximately 150
$^\circ$
(30
$^\circ$
in the original reference due to a different definition) under which the coalescence cannot lead to the jump off for any bubble size. Zhao et al. (Reference Zhao, Hu, Cheng, Huang and Gong2022) studied a similar problem but for a single coalescing bubble size of
$R=230$
μm and included the contact angle hysteresis. These numerical works focused on the conditions under which jump off is determined by the adhesion energy and, as a result, the contact angle. As described above, such behaviour prevails at lower
$ \textit{Oh}$
, where inertio-capillary forces are dominant. In the present study, however, our interest is towards jump off when it is determined by viscous dissipation. To isolate this effect, we assume perfect wetting, which fully eliminates the adhesion energy effect. Subsequently, we study the coalescence at a wide range of
$ \textit{Oh}$
spanning the viscous-dissipation-controlled jump-off threshold. Unlike the previous works, in which developing jump-off maps have been the main goal, we strive to study the flow phenomena including energetic aspects of bubble coalescence as well as the scaling laws as viscous effects start to dominate the problem. In doing so we study the interface evolution, velocity and energy content of the system across nearly two orders of magnitude of
$\textit{Oh}$
, i.e.
$0.02\leqslant \textit{Oh}\leqslant 0.86$
at different diameter ratios.
The paper is organised in the following way. In § 2 we describe the numerical methodology and the exact definition of the problem. Section 3 contains the results and the discussion, where the evolution of the interface, the kinematics and the energetics of coalescence are discussed in three consecutive sub-sections. Finally, in § 4 we summarise the main outcomes of this work.
2. Methodology
2.1. Numerical solution
We consider coalescence of two bubbles of radii
$R_1$
and
$R_2$
deposited on a flat substrate with an initial contact angle of 180
$^\circ$
(as measured in the gas phase). The ‘coalescence’ is initiated by allowing an overlap of 0.01
$R_1$
at time
$t=0\,\text{s}$
. As shown in Appendix B the results are not affected by this choice of the initial overlap. We employ the Basilisk (basilisk.fr) flow solver to solve the incompressible Navier Stokes equations with surface tension
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}
\begin{equation} \rho \left (\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\boldsymbol{u} \right )= -\boldsymbol{\nabla }\!p+\boldsymbol{\nabla }\boldsymbol{\cdot }\big[\mu \big(\boldsymbol{\nabla }\!\boldsymbol{u}+\boldsymbol{\nabla }\!\boldsymbol{u}^T\big )\big ] +\sigma \kappa \boldsymbol{n}\delta _s, \end{equation}
where
$\rho$
,
$\boldsymbol{u}$
,
$p$
,
$\mu$
,
$\sigma$
,
$\kappa$
,
$\boldsymbol{n}$
and
$\delta _s$
denote density, velocity vector, pressure, viscosity, surface tension, interface curvature, interface normal vector and Dirac delta function at the interface. In our notation
$(x,y,z)$
denote the coordinates along the coalescence direction, normal to the substrate and the third direction normal to both, respectively, and
$(u,v,w)$
are the corresponding velocity components. The interface location is determined based on the volume fraction field
$\chi$
governed by a pure advection transport equation
\begin{equation} \frac {\partial \chi }{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\chi =0. \end{equation}
The entire domain is treated as a single fluid in which the physical properties in each control volume are an average of the liquid values (
$\rho _l,\mu _l$
) and gas values (
$\rho _g,\mu _g$
) weighted by
$\chi$
, whose value is 0 and 1 in the control volumes occupied by the gas and liquid, respectively. Basilisk employs a numerical scheme similar to Gerris (Popinet Reference Popinet2003, Reference Popinet2009); the Navier Stokes equations are integrated using the projection method (Chorin Reference Chorin1969) with second-order spatial discretisation specifically using the Bell–Collela–Glatz unsplit upwind scheme for the advection term in the momentum equation. A geometric volume of fluid (VoF) method based on piecewise linear interface calculation (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Tryggvason et al. Reference Tryggvason, Scardovelli and Zaleski2011) is used for the interface reconstruction and flux calculation of
$\chi$
. The interface curvature and interface normal vector are calculated using the well-balanced height-function method described in Popinet (Reference Popinet2018).
The problem is solved in a cubic domain, as shown in figure 1(b), with symmetry boundary conditions at all boundaries except for the substrate where a no-slip boundary condition is applied for velocity and fully wet boundary conditions for the volume fraction.
Numerically, the fully wet boundary condition implies that the ‘ghost’ cells are always fully filled with liquid; therefore, the cells adjacent to the substrate may be considered as the interfacial cells, without formation of a three-phase contact line.
A multi-level Cartesian grid with adaptive refinement based on wavelet decomposition of
$\boldsymbol{u}$
and
$\chi$
(Van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018) is used to locally refine the grid. The minimum mesh size
$\varDelta_{\textit{min}}$
and grid refinement thresholds are chosen based on a comprehensive mesh convergence study described in Appendix A. Accordingly, the minimum mesh size used in the simulations is between approximately
$\varDelta_{\textit{min}} = 0.003R_1$
for the cases with the largest
$ \textit{Oh}$
and
$\varDelta_{\textit{min}}= 0.023R_1$
for the cases with the smallest
$ \textit{Oh}$
. The grid refinement threshold (Van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018) is set at
$10^{-3}$
for
$\chi$
and
$10^{-3}$
or
$10^{-4}$
for velocity depending on
$ \textit{Oh}$
. The grid details are summarised in table 1.
Physical and numerical configuration of the different simulations carried out. Here,
$\epsilon _v$
and
$\epsilon _c$
denote the grid refinement thresholds.
Table 1. Long description
A table with nine rows and seven columns. The columns are labeled Colour, Code, R1/R2, Oh, εv, εc, and Δmin/R1. The table presents different configurations for simulations with varying values of R1/R2 and Oh. Each row corresponds to a specific configuration with associated values for εv, εc, and Δmin/R1. The configurations are grouped by different values of R1/R2 (1.0, 1.3, 1.5) and Oh (0.019, 0.027, 0.038, 0.045, 0.085, 0.271, 0.858). The table includes specific values for each configuration, such as 1.0_019 with Oh 0.019, εv 0.0058, εc 10^-3, and Δmin/R1 0.0058.
Illustration of initialisation of coalescing bubbles, with zoom on the area where the bubbles overlap (a) and computational domain (b). Note that the origin of coordinate is on the substrate and vertically aligned with the centre of volume of the bubbles at the moment of coalescence.
Two main simplifying assumptions are adopted; gravity is neglected since the
$ \textit{Oh}$
range of interest leads to bubble diameters much smaller than the capillary length for most real fluid pairs, and the flat substrate is assumed to be perfectly hydrophilic, meaning that no pinning forces occur between the substrate and the bubbles.
2.2. Problem overview
Based on a dimensional analysis, it can be shown that the problem is governed by four independent non-dimensional parameters: the Ohnesorge number
$\textit{Oh}=({\mu _l}/({(\rho _l \sigma R_m)^{1/2}})$
, radius ratio
$R_1/R_2$
, density ratio
$\rho _l / \rho _g$
and viscosity ratio
$\mu _l/\mu _g$
, where the subscripts
${\scriptscriptstyle \Box }_l$
and
${\scriptscriptstyle \Box }_g$
represent the properties of liquid and gas, respectively,
${\scriptscriptstyle \Box }_1$
and
${\scriptscriptstyle \Box }_2$
signify the two bubbles and
$R_m=(R_1^3+R_2^3)^{1/3}$
is the radius of the merged bubble at equilibrium. The density ratio and viscosity ratios are kept constant at
$\rho _l / \rho _g=740$
and
$\mu _l/\mu _g=56$
, which roughly correspond to the properties of oxygen bubbles in 30 % hydrogen peroxide
$\text{H}_2\text{O}_2$
inspired by the experiments of Lv et al. (Reference Lv, Peñas, Le The, Eijkel, van den Berg, Zhang and Lohse2021). Our study covers a range of Ohnesorge numbers between
$0.019$
and
$0.858$
at three different radius ratios of
$1$
,
$1.3$
and
$1.5$
. Different configurations are listed in table 1. Note that the variation of
$ \textit{Oh}$
in the present simulations is achieved by varying the radius at fixed values of properties (
$\rho _l=1110$
kg m
$^{-3}$
,
$\mu _l=0.00113$
kg m
$^{-1}$
s
$^{-1}$
,
$\sigma =0.071$
kg s
$^{-2}$
). With these property values, the smallest
$ \textit{Oh}$
value in the present simulations corresponds to a radius of 44 μm.
Two different sets of scales are examined in this study. The time and velocity scales arising from density and surface tension are referred to as the ’inertio-capillary’ scales
\begin{equation} \begin{aligned} t_{ic} = \sqrt {\frac {\rho _lR_m^3}{\sigma }}, \quad & u_{ic} = \sqrt {\frac {\sigma }{\rho _l R_m}}. \end{aligned} \end{equation}
The ones arising from viscosity and surface tension
\begin{equation} \begin{aligned} t_{vc} =\frac {R_m\mu _l}{\sigma }, \quad u_{vc} = {\frac {\sigma } {\mu _l} }, \end{aligned} \end{equation}
are referred to as the ’visco-capillary’ scales. Note that the ratios of the above time and velocity scales are
$\textit{Oh}^{-1}$
and
$ \textit{Oh}$
, respectively. In the following, we frequently report the dimensionless values of time and velocity in terms of
$t/t_{ic}$
and
$\boldsymbol {u}/u_{ic}$
.
Bubble interface shapes on the
$z=0$
plane for different values of
$ \textit{Oh}$
at
$R_1/R_2 = 1$
. Colour code is identical to table 1.
3. Results and discussion
3.1. Evolution of the interface
Figure 2 shows the evolution of bubble interface at
$R_1/R_2=1$
and different non-dimensional times after the coalescence for the entire studied range of
$ \textit{Oh}$
. It is observed that, apart from the two largest values of
$ \textit{Oh}$
(
$ \textit{Oh} \geqslant 0.271$
), the interface shape is not significantly affected by
$ \textit{Oh}$
as long as the time is normalised by
$t_{ic}$
. The near independence from
$ \textit{Oh}$
, and hence viscosity, is expected at small
$ \textit{Oh}$
. In this range of
$ \textit{Oh}$
, the merged bubble goes through a sequence of phases after coalescence. Initially, a ’neck’ forms at the location of coalescence, which later evolves into capillary waves propagating towards the opposite apex of the two original bubbles. When the waves on the lower side of the bubbles reach the points of contact with the substrate (which occurs between the snapshots at
$t/t_{ic}=0.28$
and
$0.42$
) this leads to a detachment from the substrate since there is no pinning force in the current simulations. We refer to this incident as a ‘symmetric detachment’, which should be distinguished from the ‘jump-off’ as it is merely a result of the near-symmetric propagation of the capillary waves. At
$t/t_{ic}\approx 0.57$
the waves converge at the opposite apex of the bubbles, leading to a ‘lemon-like’ shape. After this point, the bottom of the merged bubble starts pushing against the substrate. As a result, an increasing asymmetry is observed around the
$xz$
-plane, and eventually, the merged bubble can jump off the substrate. A notable observation in figure 2 is that, while in most of the snapshots only the largest two
$ \textit{Oh}$
do not conform to others, this is not the case from the moment of lemon-like shape to approximately
$t/t_{ic}=1$
, indicating that the viscous effects possibly extend to a wider range of
$ \textit{Oh}$
in this period of bubble evolution.
At
$ \textit{Oh} \geqslant 0.271$
, where viscous effects are clearly present, the shape of interface and the sequence of events are different. Here, the wave front seems to be smoothened by viscous effects over time. As a result no clear lemon-like shape can be identified for these cases. Moreover, the evolution time is larger when normalised with
$t_{ic}$
, hinting at transition to an alternative visco-capillary time scale. We also observe no symmetric detachment and jump-off event in this range due to strong damping of the capillary waves. The case with the intermediate
$ \textit{Oh}$
of 0.085 shows some similarity with these two cases in that it does not show a clear lemon-like shape and jump-off.
Figure 3 shows the evolution of bubble interfaces at
$R_1/R_2=1.5$
. Similar to at
$R_1/R_2=1$
, for the two cases at the largest
$ \textit{Oh}$
, the bubble shape as a function of
$t/t_{ic}$
does not conform with the others. The difference is, however, more gradual with an increase of
$ \textit{Oh}$
. Moreover, due to the initial asymmetry around both the
$yz$
- and
$xz$
-planes, the capillary waves reach the apex at different times and a clear lemon-like shape is not identified. Furthermore, in all cases, except the one at
$ \textit{Oh} =0.858$
, the smaller original bubble is permanently detached during the simulation time and the merged bubble has a single area of contact with the substrate until it jumps – if it does. The detachment of the smaller original bubbles always occurs at
$t/t_{ic}\approx 0.3$
.
In order to elucidate the role of the substrate in the bubble jump-off, we conduct additional simulations for bubble coalescence without a substrate (referred to as free coalescence hereafter) for the diameter ratio of 1. Figure 4 depicts the bubble shape evolution overlaid by the wall-normal velocity field comparing a coalescence with and without a substrate. Here, the cases at the largest and the smallest studied values of
$ \textit{Oh}$
and
$R_1/R_2=1$
are shown at three moments corresponding roughly to the symmetric detachment, lemon-like shape and jump-off moments for the case with
$ \textit{Oh} =0.019$
. At
$ \textit{Oh} =0.019$
(the case in which the bubble jump off is observed) not only the bubble shape but also the velocity field remains relatively symmetric vertically (i.e. around a
$y$
-normal plane) as long as the bottom of the merged bubble has not touched the substrate.
Bubble interface shapes on the
$z=0$
plane for different values of
$ \textit{Oh}$
at
$R_1/R_2 = 1.5$
. Colour code is identical to table 1.
At
$t/t_{ic}=0.57$
, when the substrate-free bubble has crossed the imaginary substrate (that is the plane of
$y=0$
), it can be clearly observed that the vertical symmetry starts to be disturbed in the presence of substrate. The shape and the velocity field are, however, still relatively similar for simulations with and without a substrate, at least for the upper half of the bubble. It is at the moment of jump off (
$t/t_{ic}\approx 1.27$
) that the breach of vertical symmetry and similarity is most pronounced. Here, the substrate clearly pushes the bulk of the merged bubble upwards leading to a total loss of symmetry, which is essential for the jump-off. The vertical velocity at the top of the bubbles with and without substrate have different signs at this moment because the substrate-free bubble is undergoing a retraction after reaching its maximum vertical extent, leading to a negative
$y$
-velocity at the top. As shown in the second row of figure 4 similar patterns are observed at
$ \textit{Oh} =0.858$
when it comes to the symmetry-breaking effect of the substrate despite the obvious differences in the evolution time scales and bubble shapes. Moreover, at
$ \textit{Oh} =0.858$
no jump off is observed as will be discussed later.
Bubble interface at the
$z=0$
plane overlaid by the normalised wall-normal velocity at different times at
$ \textit{Oh}=0.019$
(top row) and
$ \textit{Oh}=0.858$
(bottom row) and
$R_1/R_2 = 1$
. In each plot, the left half shows coalescence without a substrate and the right half coalescence on a substrate.
In figure 5, the bubble volume crossing the
$y=0$
plane after a ’free coalescence’ (referred to as the displaced volume hereafter) is plotted against the normalised time at the entire range of
$ \textit{Oh}$
. It can readily be seen that the bubble does not cross the imaginary substrate up to a certain time (around
$t/t_{ic}=0.4$
for most cases), which is consistent with the observation that the vertical symmetry is nearly maintained in the early phases after coalescence on a substrate. In figure 5, the moment of jump off for the corresponding case with substrate is marked with a circle symbol for the values of
$ \textit{Oh}$
that the jump off is observed. Moreover, the jump off of a bubble on a substrate occurs shortly after the displaced volume is at its peak, meaning that the substrate-free bubble has started shrinking vertically. It is also observed in figure 5 that, as the value of
$ \textit{Oh}$
increases, the maximum displaced volume decreases due to the damping effect of viscosity on the bubble oscillation. The displaced volume keeps oscillating around its equilibrium value until asymptotically reaching it. At the largest value of
$ \textit{Oh}$
, however, the curve shows no maximum, meaning that the oscillation is overdamped. In this case, the maximum displaced volume is its final asymptotic value.
Normalised volume of the merged bubble in a substrate-free coalescence crossing an imaginary plane at
$y=0$
at different values of
$ \textit{Oh}$
. Colour code is identical to table 1. The horizontal dotted line indicates the normalised volume of the spherical cap below the imaginary substrate at the final equilibrium shape, a sphere of radius
$R_m$
.
In order to study the oscillation patterns of the bubbles on the substrate, the cases at the largest and the smallest values of
$ \textit{Oh}$
have been run for an extended time up to
$t/t_{ic}=10$
, and the variations of some bubble dimensions are shown in figure 6. We plot the maximum extension of the bubble in
$x$
-direction (
$d_x$
) as well as the neck diameters, defined as the maximum extension of the bubble in the
$y$
- (
$d_y$
) and
$z$
-directions (
$d_z$
) on the plane of
$x=0$
. The same dimensions are plotted using dashed lines for the corresponding free coalescence. At
$ \textit{Oh}=0.019$
, it is observed that the moment of jump off coincides with the maximum extension of the bubble in the
$y$
- and
$z$
-directions and slightly after its minimum extension in the
$x$
-direction (i.e. a nearly oblate ellipsoidal shape). As mentioned above this is shortly after the corresponding ’free bubble’ has reached its maximum vertical extension, or its most oblate shape. After the jump-off, the bubble freely oscillates with an amplitude decaying due to the viscous effects. Notably, the oscillation pattern of the free and confined bubbles are very similar in the
$z$
-direction and relatively similar in the
$x$
-direction with a phase difference originating from the time when the confined bubble is still in contact with the wall. The difference is more pronounced in
$y$
-direction, in which the effect of wall confinement is directly exerted. Here, the amplitude of oscillation is significantly suppressed in the presence of the substrate. At
$ \textit{Oh}=0.858$
, both free and confined bubbles show overdamped behaviours. However, the presence of the substrate and the resulting additional viscous dissipation seems to cause a stronger damping effect on the motion of the bubble.
Evolution of maximum extension of the merged bubble in the
$x$
-direction (a) and the neck diameters in the
$y$
-direction (b) and
$z$
-direction (c) against normalised time at
$ \textit{Oh}=0.019$
(red) and
$ \textit{Oh}=0.858$
(blue) and
$R_1/R_2 = 1$
. Solid and dashed lines indicate coalescence on a substrate and coalescence without a substrate, respectively. For the case with smaller
$ \textit{Oh}$
, the moments of symmetric detachment, lemon-like shape and jump off are indicated by square, triangle and sphere symbols, respectively.
Normalised duration of different events after coalescence up to the moment of jump off (a), and normalised time period of bubble oscillation for the first period after jump off (b) against
$ \textit{Oh}$
. The horizontal dashed lines indicate
$\tau _{n=2}/2$
and
$\tau _{n=2}$
in the left and right panels, respectively.
We further analyse the motion of bubble interface in figure 7(a) where the normalised time of different events, namely symmetric detachment, lemon-like shape and jump-off, are plotted. Note that the first two events are only defined for the diameter ratio of unity, and that not all the events occur at larger Ohnesorge numbers. As expected, the normalised times are either independent or weak function of
$ \textit{Oh}$
; the mild increase with
$ \textit{Oh}$
is clearly the effect of viscous damping. Particularly, the time of jump off is constant up to
$t/t_{ic}\approx 0.04$
–0.05). We added the time of reaching the lemon-like shape from the experiments of Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018) at a smaller
$ \textit{Oh}$
than all those studied here, which conforms well with the trend of the present data. As discussed by those authors, the evolution of the bubble at their studied value of
$ \textit{Oh}$
perfectly matches an inviscid solution. Finally, in figure 7(b) we pay attention to the oscillation period of the bubble after jump-off,
$T_{\textit{osc.}}$
Precisely, the period between the first two minima of
$d_z$
(as defined under figure 6) is plotted in this figure. Similar to the time of jump-off, the oscillation period scales well with
$t_{ic}$
and increases with an increase in the value of
$R_1/R_2$
.
In figure 7 we add horizontal dashed lines corresponding to the oscillation period of an ellipsoidal bubble in an inviscid liquid to serve as a reference. The ellipsoidal mode is the second oscillation mode whose period for a radius
$R_m$
can be calculated to be
$\tau _{n=2}=1.81t_{ic}$
(Rayleigh Reference Rayleigh1879; Lamb Reference Lamb1924). Essentially, here, we compare the current problem (figure 8
a–c) with a prolate ellipsoidal bubble starting at the interface and separating from the substrate when it reaches an oblate shape (figures 8
d and 8
e). For
$R_1/R_2=1$
our results predict a jump-off time of between
$1.32t_{ic}$
and
$1.37t_{ic}$
, which is meaningfully longer than
$\tau _{n=2}/2=0.9t_{ic}$
, as indicated in figure 7(a). After the bubble has jumped off the substrate, however, there is a close similarity between the period of its oscillation and that of an ellipsoidal bubble, as shown in figure 7(b). The minor difference is attributed to the viscous effects, which become more pronounced at larger values of
$ \textit{Oh}$
. This similarity is an indication that the higher oscillation modes are already mostly disappeared at the moment of jump-off. This point will be further elaborated in § 3.4.
Illustration of initial shape, lemon-like shape and jump-off shape (from
$a$
to
$c$
). Illustration of perfect ellipsoidal bubble oscillation (
$d,e$
). For the latter the given time period is based on the inviscid problem.
3.2. Kinematics of coalescence
Figure 9 shows the evolution of the vertical velocity of the bubble centre of mass for all studied cases. At
$R_1/R_2 = 1$
, the velocity increases sharply between the time of symmetrical detachment and lemon-like shape, at which the acceleration visibly drops. The maximum vertical velocity occurs at a point after the lemon-like shape but before jump-off. As expected, after the jump-off, when the bubble enters the free oscillation regime, the velocity drops as a result of the drag force exerted by the surrounding liquid. It can be observed that the velocity of all bubbles at
$ \textit{Oh}\leqslant 0.045$
collapse up to
$t\approx 2t_{ic}$
if both time and velocity are normalised in inertio-capillary units. This is also the case at
$R_1/R_2=1.3$
and 1.5, but with some scatter around the velocity peaks. At the larger diameter ratios, however, a somewhat different peaking behaviour is observed. Here, a pronounced first local maximum occurs once the deformation reaches the apex (the analogues of lemon-like shape) of the smaller bubble. The moment of jump off is also observed to retard further behind the maximum vertical velocity of the bubble.
Normalised vertical velocity of the merged bubble centre of mass against normalised time for all simulated cases. Colour code is identical to table 1. The moment of jump off is indicated by circle symbols for the cases where jump off occurs.
Figure 9 also reveals that, starting at a certain
$ \textit{Oh}$
value, the curves start to deviate from a common behaviour. To further elaborate this observation, in figure 10, the maximum centre of mass velocity normalised by
$u_{ic}$
is plotted against
$ \textit{Oh}$
for all studied cases. Two additional simulations have been conducted at
$ \textit{Oh}=0.495$
and 1.213 for
$R_1/R_2=1$
to complete the picture. The result in the figure show two distinct asymptotic behaviours; at
$ \textit{Oh}$
below roughly 0.05,
$V^{\textit{CM}}_{\textit{max}}/u_{ic}$
is independent of
$ \textit{Oh}$
. At
$ \textit{Oh}$
above roughly 0.5, it shows an
$\textit{Oh}^{-1}$
-dependence, which can be translated to
$V^{\textit{CM}}_{\textit{max}}/u_{vc}$
, being independent of
$ \textit{Oh}$
. Both observations are physically justified; at small
$ \textit{Oh}$
, the jump velocity scales with
$u_{ic}$
as
$ \textit{Oh}$
varies, indicating no viscous influence. In a dimensional sense, it translates to a
$V^{\textit{CM}}_{\textit{max}}\propto R_m^{-1/2}$
scaling when the radius is the only varying quantity. As
$ \textit{Oh}$
approaches unity, viscous forces dominate the inertial ones so the jump velocity scales purely with
$u_{vc}$
, which in a dimensional sense yields
$V^{\textit{CM}}_{\textit{max}}$
, being independent of
$R_m$
.
Normalised maximum vertical velocity of bubble centre of mass against
$ \textit{Oh}$
for all studied cases.
Finally, it is useful to inspect the kinematics of coalescence beyond the bubble centre of mass velocity. To this end, the contours of velocity magnitude overlaid by velocity vectors are shown in figure 11 for the largest and smallest values of
$ \textit{Oh}$
. We show only the velocity in the liquid phase as it dominates the viscous dissipation, which is the main topic of discussion in the rest of the paper. At
$ \textit{Oh}=0.019$
, the location of highest velocity clearly follows the capillary wave front that propagates the deformation. At the moment of lemon-like shape, the velocity magnitude shows a strong maximum at the apex, where capillary waves from different sides converge. As the moment of jump off is approached, there is a strong downward velocity near the bottom of the bubble, which ensures the continuity as the bulk of bubble is moving upwards at this point. Remarkably, the maximum velocity around the moment of lemon-like shape is significantly larger than the rest of the bubble life time (note the different colour scales). It is also notable that, at least up to the moment of lemon-like shape, the largest liquid velocities in the system are normal to the bubble interface following the movement of the interface itself. This is somewhat different from the notion of an interface-parallel shear flow. At
$ \textit{Oh}=0.858$
, where the viscous effects dominate over inertia, a decrease in the magnitude of velocity and a steady increase in the thickness of the layer affected by the movement of bubble interface are observed. We also observe that the maximum velocity is either at the top of the merged bubble (initial stage) or at the apex.
Contours of absolute velocity in the liquid phase overlaid by velocity vectors at the planes of
$z=0$
(top row) and
$y=R_1$
(bottom row). The left and right halves of each plot correspond to
$ \textit{Oh}=0.858$
and
$ \textit{Oh}=0.019$
, respectively. Since the velocity magnitude is given by the colour scales, the size of vectors are adjusted separately in each plot for the best clarity.
Figure 11. Long description
The image contains six contour plots with velocity vectors overlaid, arranged in two rows and three columns. Each plot represents the velocity magnitude in the liquid phase at different times, with the left and right halves of each plot corresponding to different Ohnesorge numbers (Oh). The top row shows the velocity contours and vectors in the y-x plane, while the bottom row shows them in the z-x plane. The color scale indicates the velocity magnitude, with separate adjustments made to the size of the vectors in each plot for clarity. The axes are labeled with normalized units, and the time progression is indicated at the top of each column.
Variation of potential energy, kinetic energy and accumulated viscous dissipation at
$ \textit{Oh}=0.019$
(red) and
$ \textit{Oh}=0.858$
(blue) and
$R_1/R_2 = 1$
against normalised time. The sum of the three contributions is also plotted for both cases as well as a horizontal line at unity.
3.3. Energetics of coalescence
Figure 12 shows the energy budget of the system for the largest and the smallest values of
$ \textit{Oh}$
normalised by total excess surface potential energy
$E_{\textit{tot}}=4\pi \sigma [(R_1^2+R_2^2)- R_m^2]$
. Since the system is initially at rest and the gravitational potential energy is neglected, energy conservation requires that
$E_{\textit{tot}}$
is equal to the sum of three components: the surface potential energy difference with the equilibrium value (simply referred to as potential energy hereafter)
\begin{equation} E_P= \sigma (A_s-4\pi R_m^2), \end{equation}
where
$A_s$
is the total interface surface area, the kinetic energy
\begin{equation} E_K= \int _\varOmega \frac {\rho }{2}\|\boldsymbol {u} \|^2\texttt {d}\mathcal{V} , \end{equation}
where
$\int _\varOmega \square \texttt {d}\mathcal{V}$
indicates volume integration over the entire domain, and the total accumulated viscous dissipation
\begin{equation} E_\mu =\int _0^t\dot {E}_\mu \texttt {d}t^\prime =\int _0^t\int _\varOmega \phi \texttt {d}\mathcal{V}\texttt {d}t^\prime, \end{equation}
where
$\dot {E_\mu }$
is the instantaneous dissipation rate and
$\phi = \mu {(\partial u_i}/{\partial x_{\!j}}) (({\partial u_i}/{\partial x_{\!j}}) + ({\partial u_{\!j}}/{\partial x_i} ))$
(
$i,j$
are summation indices in Cartesian coordinates) is the local volumetric dissipation rate. Sum of the three budget terms is also plotted in figure 12. We observe a less than 2 % drop in the computed sum, which represents the numerical error and a minor energy flux at the domain boundaries. The moments of symmetric detachment, lemon-like shape and jump off for the smaller Ohnesorge number are specified by the rectangle, triangle and diamond symbols in the figure.
As observed in figure 12, the potential and kinetic energies of the system vary in time while their sum is reduced due to viscous dissipation. As expected, the two cases shown in the figure display very different behaviours when it comes to the evolution of potential and kinetic energies. At
$ \textit{Oh}=0.019$
, the exchange between kinetic and potential energies leads to alternating maxima and minima, whose values decay in time. The first local maximum of kinetic energy occurs shortly after the moment of lemon-like shape. After this point, the kinetic energy decreases monotonically until the jump-off, but nearly the entirety of this loss is consumed by viscous dissipation since in this interval the surface potential energy barely changes. Only shortly before the jump off a peak in
$E_K$
emerges. The moment of jump off itself corresponds to the first local minimum of
$E_K$
. Following the jump-off, the merged bubble keeps oscillating while the motion is being asymptotically damped by viscous dissipation, as also indicated in the previous subsections. The case with
$ \textit{Oh}=0.858$
clearly shows the characteristics of an over-damped oscillatory system, with the potential energy monotonically decreasing in time. In this specific case, nearly the entire released potential energy is consumed by dissipation with very small rise in kinetic energy.
Figure 13 depicts evolutions of the three energy components for all cases run in this study. Generally, both kinetic and potential energies show similar trends down to
$ \textit{Oh}=0.045$
. At
$ \textit{Oh}=0.085$
, a more pronounced deviation starts to emerge while for the cases at
$ \textit{Oh}\geqslant 0.271$
the behaviour is visibly different. It is furthermore observed that, as the value of
$R_1/R_2$
increases, a second local peak in
$E_k$
appears at larger values of
$ \textit{Oh}$
. The double peak behaviour is linked to the fact that capillary waves take different amounts of time to reach the apex of the two original bubbles leading to a more gradual conversion of potential energy leading to two weaker peaks instead of one strong one. Finally, it is observed that, as
$R_1/R_2$
grows, the moment of jump off does not coincide with the first minimum (maximum) of kinetic (potential) energy, as it does at
$R_1/R_2=1$
, and happens with an increasing time delay.
(a–c) Variation of kinetic energy (top row), (d–f) potential energy (middle row) and (g–i) accumulated dissipated energy (bottom row) at
$R_1/R_2=1$
(left column),
$R_1/R_2=1.3$
(middle column) and
$R_1/R_2=1.5$
(right column) against normalised time. Colour code is identical to table 1.
Variation of normalised dissipation rate at
$R_1/R_2=1$
(a),
$R_1/R_2=1.3$
(b) and
$R_1/R_2=1.5$
(c) against normalised time. Results in the range
$0.019\leqslant Oh \leqslant 0.085$
are shown. Triangle and circle symbols indicate the moments of lemon-like shape and jump-off. Colour code is identical to table 1.
Figure 14. Long description
Three line graphs depict the variation of normalized dissipation rate against normalized time. Panel A: The line graph shows the normalized dissipation rate on the vertical axis and normalized time on the horizontal axis. The data is represented by multiple lines in different colors, with triangle and circle symbols indicating specific moments. Panel B: Similar to Panel A, this line graph also shows the normalized dissipation rate against normalized time with different colored lines and symbols. Panel C: This line graph follows the same pattern as Panels A and B, displaying the normalized dissipation rate against normalized time with various colored lines and symbols.
Notably, accumulated viscous dissipation curves displayed in the last row of figure 13 show a similar growth rate in all cases up to, and including,
$ \textit{Oh} =0.271$
. We quantify this by calculating the duration in which half of the total energy is dissipated. The duration varies between
$0.56t_{ic}$
and
$0.66t_{ic}$
in the range
$0.019\leqslant$
$ \textit{Oh} \leqslant 0.271$
for
$R_1/R_2=1$
(and varies less for the other diameter ratios). In this sense,
$t_{ic}$
can be considered as the time scale at which the released surface energy is dissipated for the studied range. At this point, one can use a scaling argument to estimate the length scale of viscous dissipation
$\tilde {\lambda }$
– that is the length at which the variation of velocity in the vicinity of the bubble takes place. Taking
$u_{ic}$
as the velocity scale follows
$E_\mu \sim (\mu _l{u_{ic}^2}/{\tilde {\lambda }^2} ) (R_m^2 \tilde {\lambda } )t_{ic}$
, where the first and second brackets on the right-hand side represent the scales of
$\phi$
and the volume at which it takes place. That the total dissipation should also scale with the excess surface energy (
$\sim \sigma R_m^2$
) requires
\begin{equation} \tilde {\lambda }\sim R_m\textit{Oh} \end{equation}
for the scales to remain consistent. Note that
$\tilde {\lambda }$
in (3.4) is equivalent of the diffusion length scale based on the inertio-capillary velocity (
$\nu _l/u_{ic}=R_mOh$
). It must be born in mind that
$\tilde {\lambda }$
is simply the global length scale of dissipation if
$t_{ic}$
and
$u_{ic}$
are the relevant time and velocity scales and it is not necessarily locally relevant. The tilde
$\tilde {\Box }$
sign in our notation is included to emphasise this nature.
With
$\tilde {\lambda }$
determined from (3.4) the scale of
$\dot {E}_\mu$
will be
$\mu _l{(( R_mu_{ic}^2})/({{\textit{Oh}}}))=(({\mu _l \sigma })/{\rho _l}) ({1}/{\textit{Oh}})$
. Notably, the viscosity is cancelled out at this scale due to the division by
$ \textit{Oh}$
, indicating that the dissipation rate in the system is independent of the value of viscosity, which is consistent with the notion of an inertio-capillary regime. We plot the dissipation rate normalised by
$(({\mu _l \sigma })/{\rho _l})({1}/{\textit{Oh}})$
in figure 14 for the values of
$ \textit{Oh}$
up to 0.085. It is observed that, up to the moment of jump-off, the normalised dissipation rate is of the order of unity, and moreover, a satisfactory collapse of the curves is achieved for almost the entire range. Only the largest
$ \textit{Oh}$
(
${=}0.085$
) shows some deviation.
It is also observed in figure 14 that the normalised dissipation rate varies (and generally decreases) in time, which is related to the fact that the instantaneous velocity of the interface and the thickness of the boundary layer around it both vary in time (note that, as discussed before,
$\tilde {\lambda }$
is the ‘global’ length scale of dissipation given a time scale
$t_{ic}$
and not a local boundary layer thickness). Figure 14 suggests that the dissipation rate is at its highest immediately after the coalescence and at the moment of lemon-like shape, where the boundary layer thickness and the interface velocity are at the smallest and largest, respectively. We observe that the collapse of different curves is particularly good between the maximum corresponding to the lemon-like shape and the minimum shortly before the jump-off, the latter being clearly linked to the minimum of
$E_K$
observed in figure 13.
While figure 14 shows the variation in time of viscous dissipation integrated in space, it is insightful to investigate how the dissipation is spatially distributed. To this end, figure 15 depicts normalised volumetric dissipation rate for three values of
$ \textit{Oh}$
corresponding to the smallest and the largest collapsing curves in figure 14 (
$ \textit{Oh}=0.019$
and
$0.045$
) as well as the one at which the deviation is more pronounced (
$ \textit{Oh}=0.085$
). We show four moments in time: before the lemon-like shape, approximately at the lemon-like shape, after the lemon-like shape and the approximate moment of jump-off. Because the outlined scaling analysis indicates that the volumetric dissipation rate scales with
$\mu _l( {u_{ic}^2}/{\tilde {\lambda }^2})\sim \mu _l({u_{ic}^2}/({(R_m\textit{Oh})^2}))$
, we plot
$\phi$
normalised by
$\mu _l({u_{ic}^2}/({(R_m\textit{Oh})^2}))$
in this figure. The maxima of the colour scales in this figure are roughly adjusted based on the dissipation rate at different moments.
Normalised volumetric dissipation rate
$\phi ^*=\phi (( { (R_m\textit{Oh})^2})/({\mu _l u_{ic}^2}))$
at
$ \textit{Oh} =0.085$
(left),
$ \textit{Oh} =0.045$
(middle) and
$ \textit{Oh} =0.019$
(right). Panels (a–d) correspond to different times. Note the logarithmic colour scale.
Figure 15 demonstrates that viscous dissipation is typically very concentrated at certain regions adjacent to the interface. Prior to the moment of lemon-like shape (
$t/t_{ic}=0.28$
), dissipation is maximised around the capillary wave front, particularly at the bottom of the bubble where the interface meets the wall. At this stage, both the magnitude and thickness of the high-dissipation region increase with viscosity (or Ohnesorge number) as also the normalised
$\dot {E}_\mu$
curves in figure 14 suggest. As the lemon-like shape appears (
$t/t_{ic}\approx 0.57)$
, the dissipation is highly concentrated around the apex. It was previously shown that the velocity magnitude in the liquid phase reaches its highest value at this point, which is consistent with the dissipation pattern. When comparing different values of
$ \textit{Oh}$
at this moment, the curvature of interface decreases at the apex, resulting in less intensive and less concentrated dissipation. Note that the peak value of normalised
$\dot {E}_\mu$
in figure 14 decreases with an increasing
$ \textit{Oh}$
, which is in line with this pattern. The last rows in figure 15 shows the results after the moment of lemon-like shape, where dissipation is mainly limited to the lower half of the merged bubble and wall shear seems to be the dominant effect. The local maxima of dissipation are, however, smaller than the earlier stages of bubble life time where the fast-moving capillary wave front create concentrated regions of high viscous dissipation.
Before closing the presentation of energetic aspects, it is worthwhile examining the energy of the system at jump off in more detail. To this end, we define the jump energy
$E_{\!j}$
as the sum of kinetic and potential energies (
$E_K+E_P$
) at the moment of jump-off. In other words, the jump energy is the part of the total released surface energy that remains in the system to be dissipated after the bubble jumps. Obviously,
$E_{\!j}$
must have a finite value for the jump off to happen. It must be born in mind that, if the bubble is attached to the substrate with a contact angle of smaller than 180
$^\circ$
(which is the case in reality), part of
$E_{\!j}$
must be spent to overcome adhesion, therefore, a higher
$E_{\!j}$
means that a higher adhesion energy barrier can be overcome. Figure 16 depicts the computed values of
$E_{\!j}/E_{\textit{tot}}$
in a
$R_1/R_2$
–
$ \textit{Oh}$
map for all simulated cases. The cases where jump off is not observed are indicated by cross symbols in this map. A clear trend is revealed in this figure; the jump energy increases as both
$ \textit{Oh}$
and
$R_1/R_2$
decrease. Considering that the jump energy is linked to a higher likelihood of jump off for ’pinned’ bubbles, this trend matches the experimental observation of Lv et al. (Reference Lv, Peñas, Le The, Eijkel, van den Berg, Zhang and Lohse2021), which is itself aligned with a scale analysis conducted by the same authors invoking the balance of surface potential energy and viscous dissipation.
Remaining energy of the system (
$E_{\!j}=E_K+E_P$
) at the moment of jump off for different simulated cases. Where jump off is not observed, a cross symbol is used. Red lines of different shades show iso-lines of
$E_{\!j}/E_{\textit{tot}}=0.17$
, 0.13, 0.09 and 0.05 (from left to right) interpolated based on the existing data. The jump-off criterion suggested by Demirkırr et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025) is shown by a dashed line in the map.
Figure 16 shows that, above a certain value of
$ \textit{Oh}$
, roughly corresponding to the onset of the visco-capillary scaling regime, bubbles cease to jump off after coalescence in the studied range of
$R_1/R_2$
. In the absence of gravity, this translates to all the released energy being asymptomatically dissipated towards the equilibrium spherical shape without the centre of mass of the merged bubble passing
$y=R_m$
. In the present study, we do not continue the simulations towards the equilibrium as the large additional computational cost would not be justified. In figure 16, we also depict the jump-off criterion developed by Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025). It can be seen that the slope of this line in the map, which is essentially a result of the scale analysis conducted by these authors, is between vertical and that of the present
$E_{\!j}$
iso-lines. Also the line lies to the left (smaller
$ \textit{Oh}$
values) of where our simulations predict the jump-off threshold. We remind the reader that the above criterion is based on a scale analysis and the line plotted here corresponds to the predicted asymptotic behaviour at the limit of small adhesion energy. Considering that this criterion is tuned with data at finite contact angles (i.e. with pinning), a minor mismatch with the current DNS results is not unexpected. It is in general expected that a finite contact angle would globally displace the jump-off threshold towards lower
$ \textit{Oh}$
values, therefore matching with the present observation. Furthermore, Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025) report the necessity of increasing the effective Ohnesorge number for asymmetric coalescence events with a corrective term function of
$R_1/R_2$
, recognising the tendency of the jump-off threshold to ’bend’ towards lower
$ \textit{Oh}$
numbers at higher values of the radius ratio, in agreement with our simulations. Interestingly, Lv et al. (Reference Lv, Peñas, Le The, Eijkel, van den Berg, Zhang and Lohse2021) attribute this bending effect especially to the presence of a finite contact angle, suggesting that the pinning force is more relevant at higher
$R_1/R_2$
.
3.4. Discussion
Before concluding, it is worthwhile discussing a few key insights that link several of the results presented in the previous three subsections. A motivation of the present work has been to quantify the range of
$ \textit{Oh}$
at which a transition from the inertio-capillary to visco-capillary force dominance occurs. The most clear evidence of such transition can be sought in figure 10 where, at adequately large (small) values of
$ \textit{Oh}$
, a nearly perfect scaling of bubble velocity with the visco-capillary (inertio-capillary) velocity is observed. The region at which neither perfect scaling is observed (hereafter the transition region) is at
$0.05\leqslant$
$ \textit{Oh}$
$\leqslant 0.5$
while the most visible transition takes place between the values of roughly 0.1 and 0.3. The transition interval is only slightly affected by the radius ratio in the studied range. The start of deviation from inertio-capillary scaling at
$ \textit{Oh}\approx 0.05$
is in agreement with other observations. In figures 2 and 3, the normalised bubble shapes collapse nicely when depicted as a function of
$t/t_{ic}$
up to
$ \textit{Oh}=0.045$
and the first case that shows meaningful deviation is the one at
$ \textit{Oh}=0.085$
. Moreover, as shown in figure 14, when plotting the normalised viscous dissipation rate against
$t/t_{ic}$
the case at
$ \textit{Oh}=0.085$
is the first one showing a significant deviation. To complete the picture, there seems to be a clear link between this regime transition and jump off. In this regard, our results indicate that, as
$ \textit{Oh}$
values grow at a fixed radius ratio, the bubble stops to jump off around the onset of transition away from the inertio-capillary regime (see figure 16).
At this point, it can also be insightful to draw a parallel with the results on coalescence of droplets on superhydrophobic surface reported by Liu et al. (Reference Liu, Ghigliotti, Feng and Chen2014). These authors reported the inertio-capillary regime to hold until approximately
$ \textit{Oh}\approx 0.1$
. While these authors did not explicitly examine the visco-capillary scaling, they indicated full viscous behaviour to be reached at approximately
$ \textit{Oh}=0.3$
. These values are similar to the core transition region we find for bubble coalescence underlying a close analogy. The fact that viscous dissipation for bubbles and droplets occurs mainly on the opposite sides of the interface does not seem to be a major determining factor in this discussion, which can be justified by the fact that at
$ \textit{Oh}\approx 0.1$
the viscous length scale
$\tilde \lambda$
is approximately
$0.1R_m$
, meaning that the viscous region inside a droplet is not significantly modified by the confined geometry of the droplet (for a extended discussion on dissipation length scales inside droplets readers are referred to the recent work of Sanjay et al. (Reference Sanjay, Zhang, Lv and Lohse2025)).
Our results (e.g. figure 10) also indicate a rather ‘gradual’ transition to the visco-capillary regime. This is arguably due to the non-uniform variation of viscous dissipation rate in both time and space causing the dissipation to become significant at certain regions even at relatively small
$ \textit{Oh}$
values while in other regions this happens at larger
$ \textit{Oh}$
values. In particular, our results highlight the moment the merged bubble shape becomes lemon-like as a result of the capillary waves meeting at its opposite apex (while a perfect lemon-like shape is only seen at
$R_1/R_2$
= 1, a similar dynamics prevails at other ratios when the waves reach at the apex of the larger bubble). The results in figure 14 indicate that the total dissipation rate peaks at the moment of lemon-like shape, featuring it as a major dissipative event during the bubble life time. Spatially, the dissipation rate is very local at this moment and mainly concentrated around the apex of the bubble (see figure 15, where large interface velocities resulting from the convergence of the capillary waves lead to a large dissipation rate at the two apexes of the bubble). It is also notable that the normalised bubble shapes deviate from the inertio-capillary behaviour at smaller values of
$ \textit{Oh}$
immediately after the moment of lemon-like shape and near the apex, as shown in figures 2 and 3, which indicates stronger viscous effect in this specific time and region.
In relation to the high local viscous dissipation at the moment of lemon-like shape, one naturally expects such an event to intensify the decay of faster capillary waves propagating on the bubble interface. Such an expectation indeed matches the data in figure 7, according to which the time of reaching the lemon-like shape is much shorter than the time between that moment and that of jump off. For
$R_1/R_2=1$
, these time intervals are
$0.53t_{ic}$
and
$0.78t_{ic}$
, respectively. As a reference, Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018) reported a time interval of
$0.47t_{ic}$
for reaching the lemon-like shape (
$0.67t_{ic}$
in the original work as those authors defined
$t_{ic}$
based on
$R_1$
), which they calculated to be equivalent to the angular speed of the 20th inviscid bubble oscillation mode (the angular propagation speed of the nth mode is determined from
$\omega ^2_{n} = {(n-1)(n+1)(n+2)}/{t^2_{ic}}$
following Rayleigh (Reference Rayleigh1879) and Lamb (Reference Lamb1924)). Using the same approach, the present time of lemon-like shape corresponds to an effective propagation speed of the 15th inviscid mode. Meanwhile, the interval of
$0.78t_{ic}$
between the moments of lemon-like shape and jump off is very close to half that of an ellipsoidal oscillation (see figures 8
d and 8
e), which corresponds to the second inviscid mode. While the above analysis based on the effective propagation mode is of a speculative nature without a full modal analysis, which is outside of the present scope, the clear slow down of the propagation speed and a possible elimination of the higher modes after the moment of lemon-like shape is at least qualitatively in line with high viscous dissipation at this moment.
4. Conclusions
DNS are performed using the VoF implementation in the Basilisk flow solver (basilisk.fr) for two bubbles coalescing on a flat surface with a contact angle of 180
$^\circ$
. The problem is investigated at three different diameter ratios between 1 and 1.5 and the Ohnesorge number is varied between 0.019 and 0.858. Results show the jump velocity of the merged bubble to scale perfectly with the inertio-capillary velocity up to
$ \textit{Oh}\approx 0.05$
and with the visco-capillary velocity above
$ \textit{Oh}\approx 0.5$
. This ‘transition’ is most pronounced at roughly 0.1
$\lt$
$ \textit{Oh}$
$\lt$
0.3. Our results for the time evolution of the normalised bubble shape and dissipation rate also confirm the beginning of deviation from inertio-capillary scaling at around
$ \textit{Oh}$
= 0.05. Notably, the time evolution of viscous dissipation rates collapses below this value of
$ \textit{Oh}$
if the length scale of viscous dissipation
$\tilde \lambda =R_m$
$ \textit{Oh}$
is used for the normalisation. We also report a threshold of
$ \textit{Oh}=0.05$
–0.1, above which no jump off is observed, which matches well with the onset of the viscous effects.
The maximum vertical velocity of the bubble is found to be between
$0.19u_{ic}$
(for
$R_1/R_2=1$
) and
$0.25u_{ic}$
(for
$R_1/R_2=1.5$
) for the studied cases falling within the ‘inertio-capillary regime’. This is similar to the values of
$\approx 0.2u_{ic}$
for coalescence of droplets on super-hydrophobic surfaces reported by Liu et al. (Reference Liu, Ghigliotti, Feng and Chen2014). The observed start of deviation from the inertio-capillary regime in the present work is also similar to the one reported by those authors. This similarity can be attributed to the fact that, at the range of
$ \textit{Oh}$
corresponding to the inertio-capillary regime, the viscous dissipation length scale
$\tilde \lambda$
is at least an order of magnitude smaller than the bubble/droplet radius, which renders the effect of the confined droplet geometry less relevant.
At a certain time after the two bubbles merge, the capillary waves originating from the point of coalescence meet at the opposite apex of the merged bubble, resulting in a lemon-like shape. We note the significance of this event as the rate of viscous dissipation peaks at this moment. One should note that the spatial distribution of viscous dissipation rate is also highly non-uniform. In particular, the region of high dissipation is observed adjacent to the capillary wave front, and particularly around the apex of the bubble at the moment of lemon-like shape, where the largest local dissipation rate through the bubble life time is measured. After this moment, the dissipation rate generally drops, and the region of high dissipation moves near the substrate as a result of the high shear rate at the wall.
Furthermore, we study the evolution of surface potential energy, kinetic energy and accumulated viscous dissipation in the system. It is observed that, globally, the dissipation takes place at very similar rates up to
$ \textit{Oh}\approx 0.271$
with a time constant comparable to
$t_ {ic}$
, the inertio-capillary time scale based on
$R_m$
. Quantifying the non-dissipated energy of the system at the moment of jump off, a clear trend is observed: less energy is available to the system right before the jump off as either the diameter ratio or the Ohnesorge number increases. In the studied range of
$R_1/R_2$
, we do not observe a jump off similar to the cases at larger
$ \textit{Oh}$
starting at a threshold in the
$ \textit{Oh}$
range of
$0.045$
–
$0.085$
. This is somewhat smaller than the range of 0.03–0.04 that the criterion proposed by Demirkır et al. (Reference Demirkır, Yang, Bashkatov, Sanjay, Lohse and Krug2025), based on a scale analysis, predicts. The difference is attributed to the fundamental difference in the methodologies.
Overall, the present results shed light on certain aspects of the bubble coalescence problem on a flat substrate, at the limit of small adhesion and spanning the inertio-capillary and visco-capillary scaling regimes. One must, however, note that many of the simplifications do not hold in real-world applications; the substrate is not necessarily flat, the coalescing bubbles are not necessarily isolated in an otherwise quiescent liquid and the contact line dynamics also affects the process. Adding such effects as the Marangoni flows, electrostatic forces and, at larger bubble sizes, the gravity, the present work can only be regarded as a departure point for more elaborate numerical studies with systematically increasing complexity.
Acknowledgements
The authors thank Professor D. Lohse for the insightful discussions on the scaling and energy analysis. The simulations in this work were conducted on the resources provided by Danish e-Infrastructure Cooperation (project DeiC-AU-N2-2025120).
Funding
This work was supported by Innovation Fund Denmark through the Grand Solution research project LC-H2 (2077-00021B).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Grid study
The study of the mesh size
$\varDelta_{\textit{min}}/R_1$
has been carried out at varying
$R_1/R_2$
and
$ \textit{Oh}$
. In figure 17 the velocity of the centre of mass of the bubble is shown for different mesh sizes and at two values of
$ \textit{Oh}$
. The collapse of the lines proves that mesh independence is reached for
$ \textit{Oh}= 0.012$
at
$\varDelta_{\textit{min}}/R_1 = 0.0058$
, and for
$ \textit{Oh}= 0.079$
at
$\varDelta_{\textit{min}}/R_1 = 0.0234$
.
Furthermore, we report that quantities such as bubble shape and the work of viscous forces do not require the same mesh refinement as quantities such as velocity and the kinetic and potential energies. For example, in figure 18, where the evolution of the bubble shape is reported, we show that the shape outlines collapse at
$\varDelta_{\textit{min}}/R_1 = 0.0117$
for
$ \textit{Oh}= 0.012$
and at
$\varDelta_{\textit{min}}/R_1 = 0.0469$
for
$ \textit{Oh}= 0.079$
. As a result, coarser meshes can be used for the study of quantities merely related to the bubble interface geometry, translating to longer simulated physical time for the same amount of computational cost.
Evolution of the velocity of the centre of mass of the coalescing bubbles at varying mesh sizes
$\varDelta_{\textit{min}}/R_1$
for
$ \textit{Oh}= 0.012$
(a) and
$ \textit{Oh}= 0.079$
(b).
Bubble interface outline at varying
$\varDelta_{\textit{min}}/R_1$
reported at different
$t/t_{ic}$
for
$ \textit{Oh}= 0.012$
(a) and
$ \textit{Oh}= 0.079$
(b).
In figure 19 we show the different meshes used to produce the results presented in figures 17 and 18.
Overall, different minimum mesh size
$\varDelta_{\textit{min}}/R_1$
(ranging from
$0.0234$
to
$0.0058$
) can be employed to reach mesh independence at different values of
$ \textit{Oh}$
. Importantly, in all cases, for the final mesh we selected, energy is conserved within less than
$2\,\%$
error. Notably, minimum mesh size increases loosely with
$\textit{Oh}^2$
which translates to a constant mesh size if normalised with the intrinsic length scale resulting from viscosity, surface tension and density, i.e.
$\mu _l^2/(\rho _l\sigma )$
.
Mesh refinement contours at varying
$\varDelta_{\textit{min}}/R_1$
reported at
$t/t_{ic}= 0.57$
for
$ \textit{Oh}= 0.012$
(a) and
$ \textit{Oh} = 0.079$
(b).
Evolution of the velocity of the centre of mass of the coalescing bubbles at varying tolerances
$\epsilon _u$
for
$ \textit{Oh} = 0.012$
(a) and
$ \textit{Oh} = 0.079$
(b).
Regarding the variations of
$R_1/R_2$
, the main impact on the accuracy of the simulation is that
$\varDelta_{\textit{min}}/R_1$
will no longer be equal to
$\varDelta_{\textit{min}}/R_2$
, as
$R_2$
decreases. For this reason, after a series of preliminary tests, we have decided to limit
$R_1/R_2$
to
$1.5$
, to offer a meaningful range of values without significantly reducing
$\varDelta_{\textit{min}}/R_2$
.
To ensure that the results are not only independent of the minimum mesh size but the mesh adaptation pattern, in figure 20 we show the overlapping of the vertical velocity of the centre of mass of the bubbles at varying
$\epsilon _u$
, which is the most important parameter governing the mesh adaptation with respect to capturing the velocity gradients and therefore the viscous dissipation. Figure 21 shows the actual meshes at changing
$\epsilon _u$
. We notice that, in the case of symmetric bubbles, the highly refined region extends for a shorter distance (typically less than
$1R_1$
) until jump off, while in the case of non-symmetric bubbles the smaller bubble leaves a relatively large velocity wake.
Mesh size contours at
$\epsilon _v= 10^{-1}$
(left),
$\epsilon _v= 10^{-2}$
and
$\epsilon _v= 10^{-3}$
(c), reported for
$ \textit{Oh} = 0.012$
(a,b,c) and
$ \textit{Oh} = 0.079$
(d,e, f).
Appendix B. Influence of initial bubble overlap
As explained in the body of the paper, the two bubbles are initialised with an initial overlap to initiate coalescence in the simulation. The effect of the initial shape on the development of coalescence is studied in this section for bubbles of equal radius as a function of mesh size
$\varDelta_{\textit{min}}/R_1$
and within a range of initial overlapping
$r_0$
as a percentage of
$R_1$
.
In figures 22 and 23 we show the shape outline and the vertical velocity of the centre of mass of coalescing bubbles at
$ \textit{Oh} =0.0079$
(corresponding to the reference experimental work by Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018)) at three different values of
$r_0/R_1$
, namely 0.5, 1 and 1.2. A clear collapse of the results indicates independency of the results from the overlap size in the studied range.
Figure 24 shows the impact of
$r_0$
on the neck radius during the early stages of coalescence compared with the experimental data from Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018). For the very early times immediately after the coalescence (
$t/t_{ic} \lt 0.1$
), minor discrepancies are observed. Nonetheless, for
$t/t_{ic} \gt 0.1$
all the reported simulations i.e.
$r_0 = 0.5\,\%R_1, 1\,\%R_1$
and
$1.2\,\%R_1$
, are in good agreement with experimental data, showing no significant influence by the way coalescence is initiated in the simulations.
As a conclusion, we chose to utilise an initial overlap of
$r_0 = 1\,\%R_1$
for all the simulations carried on in our study, as it showed faster run time than the thinner ones.
Bubble interface at varying initial overlap
$r_0$
reported for different
$t/t_{ic}$
. Mesh size is
$\varDelta_{\textit{min}}/R_1 = 0.0117$
and
$ \textit{Oh} =0.0079$
.
Evolution of the velocity of the centre of mass of the coalescing bubbles at varying initial overlap
$r_0$
. Mesh size is
$\varDelta_{\textit{min}}/R_1 = 0.0117$
and
$ \textit{Oh} =0.0079$
.
Evolution of the normalised neck radius
$r_y/R$
over time for a range of initial overlap
$r_0$
values. The red circles represent experimental data from Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018) for bubbles of equal size having
$R_1 = R_2 = 250$
μm translating to
$ \textit{Oh} =0.0079$
. Mesh size is
$\varDelta_{\textit{min}}/R_1 = 0.0117$
.
Coalescence of carbon-dioxide bubbles of
$300$
μm radius immersed in water. The time difference between subsequent experimental snapshots from (a) to (b) is
$75$
μs while it is
$150$
μs for snapshots (b) and (c) and
$225$
μs between snapshots (c) and (d) and it is adjusted by a factor 1.2 in the simulations. Blue shades show the present numerical framework and black shades the experimental results from Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018).
Evolution of the normalised neck
$r_y/R_1$
for the bubble in figure 25 over time. Mesh size
$\varDelta_{\textit{min}}/R_1 = 0.0117$
.
Appendix C. Comparison with experiment
There have been no experimental studies for bubble coalescence on perfectly wet substrates due to the obvious physical constraints. However, experimental results from Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018) obtained at very large contact angles (measured in the gas phase) provide a possibility for direct comparison with experimental data. In figure 25 we report the evolution of the shape over time in a frame by frame comparison with experimental results from Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018). We used the same physical properties. As can be seen in the figure, our framework correctly captures the shape of the bubble at different sequences. The only noticeable difference is that at
$t = 225$
μs the simulated coalescing bubbles detach from the substrate before the bottom of the merged bubble touches the substrate at approximately
$t = 675$
μs (for a deep discussion on symmetrical detachment we refer the reader to figure 4 and the related considerations). In the experiment such a symmetric detachment is not observed as the experiment ensures that the bubbles are pinned at super-hydrophobic holes on the substrate. Apart from this deviation, the bubble evolution is captured well.
In figure 26, we report the evolution of the normalised neck
$r_y/R_1$
between the bubbles for our simulation together with the experimental data by Soto et al. (Reference Soto, Maddalena, Fraters, van der Meer and Lohse2018). The neck between the bubbles is captured accurately when compared with the experiment. One should note that the experiments cited above have been run in real-world conditions, while in the numerical simulation the gravitational effect is neglected. This is, however, not expected to cause a significant effect since the size of the bubbles in the experiments is around an order of magnitude smaller than the corresponding capillary length
$l_c= \sqrt {({\sigma }/({\Delta \rho g}))} \approx 2.7\,\text{mm}$
. There is moreover the effect of the three-phase contact line, which, as explained above, seems to be minor for the current experiments.
ϵv
ϵc
z=0
Oh
R1/R2=1
z=0
Oh
R1/R2=1.5
z=0
Oh=0.019
Oh=0.858
R1/R2=1
y=0
Oh
Rm
x
y
z
Oh=0.019
Oh=0.858
R1/R2=1
Oh
Oh
τn=2/2
τn=2
a
c
d,e
Oh
z=0
y=R1
Oh=0.858
Oh=0.019
Oh=0.019
Oh=0.858
R1/R2=1
R1/R2=1
R1/R2=1.3
R1/R2=1.5
R1/R2=1
R1/R2=1.3
R1/R2=1.5
0.019⩽Oh⩽0.085
ϕ∗=ϕ(((RmOh)2)/(μluic2))
Oh=0.085
Oh=0.045
Oh=0.019
Ej=EK+EP
Ej/Etot=0.17
Δmin/R1
Oh=0.012
Oh=0.079
Δmin/R1
t/tic
Oh=0.012
Oh=0.079
Δmin/R1
t/tic=0.57
Oh=0.012
Oh=0.079
ϵu
Oh=0.012
Oh=0.079
ϵv=10−1
ϵv=10−2
ϵv=10−3
Oh=0.012
Oh=0.079
r0
t/tic
Δmin/R1=0.0117
Oh=0.0079
r0
Δmin/R1=0.0117
Oh=0.0079
ry/R
r0
R1=R2=250
Oh=0.0079
Δmin/R1=0.0117
300
75
150
225
ry/R1
Δmin/R1=0.0117
Open access
