1. Introduction
Ostwald ripening is an important multiphase process ubiquitously occurring in nature, daily life and industrial applications, such as CO
$_{2}$
storage (Blunt Reference Blunt2022), the emulsion system (Doan-Nguyen et al. Reference Doan‐Nguyen, Jiang, Koynov, Landfester and Crespy2021) and bubble evolution (Chen et al. Reference Chen, Chen, Xu, Hu, Zhan, Li, Zhuang and Li2022a
). The main mechanism of Ostwald ripening lies in the dissolution of smaller particles or droplets or bubbles into the solution, and their migration towards larger ones, ultimately leading to thermodynamic equilibrium of the system with minimum free energy. Ostwald ripening is typically classified into two categories, diffusion-controlled and reaction-controlled. The former is primarily limited by slow diffusion, while the latter is determined by the low interface reaction rate exhibiting a growth law proportional to time raised to the power
$1/2$
(Watanabe, Inaoka & Ito Reference Watanabe, Inaoka and Ito2016). Zhu, Zhang & Lohse (Reference Zhu, Zhang and Lohse2018) studied multi-bubble ripening in the diffusion regime by considering the attractive or repulsive effects of adjacent bubbles, similar to the evaporation dynamics of adjacent droplets (Hu et al. Reference Hu, Wang, Man and Doi2017). Watanabe et al. (Reference Watanabe, Suzuki, Inaoka and Ito2014, Reference Watanabe, Inaoka and Ito2016) elucidated the influence of temperature on ripening mechanisms, demonstrating a transition from reaction-controlled to diffusion-controlled ripening, with reaction-controlled ripening dominating in bubble systems at low temperatures. In moderately wet foams, Lambert et al. (Reference Lambert, Cantat, Delannay, Mokso, Cloetens, Glazier and Graner2007) revealed through experiments that the bubble ripening process deviates from a purely diffusive mechanism, exhibiting partial reaction-controlled characteristics, underscoring the critical role of reaction kinetics under wet conditions. Through theoretical and numerical analyses, Tjhung, Nardini & Cates (Reference Tjhung, Nardini and Cates2018) revealed how reaction-controlled inverse Ostwald ripening in active matter systems drives microphase separation of bubbles and clusters, exhibiting unique dynamics of self-organised criticality and non-equilibrium phase separation. These works show that although diffusion-controlled Ostwald ripening is widely observed, interface-controlled ripening is also a physically relevant mechanism occurring when interfacial kinetics becomes rate-limiting. Understanding this less-explored reaction/interface-controlled regime for attached bubbles, which has a basis in reality under certain conditions, is the primary focus of our theoretical and numerical investigation in the current work.
From a practical standpoint, the effective prevention or control of ripening continues to pose a major challenge across both scientific research and industrial practice. Recently, various methods have been proposed to control Ostwald ripening. For example, Guo et al. (Reference Guo, Zhang, Wang, Zhang, Xu and Sun2024) showed that by modifying the composition or adjusting the interfacial properties of the two phases, ripening in emulsions can be effectively suppressed. In the work of Doan-Nguyen et al. (Reference Doan‐Nguyen, Jiang, Koynov, Landfester and Crespy2021), silica shells were generated via sol–gel reactions to restrict the oil cores formed by Ostwald ripening, enabling the fabrication of ultra-small nanocapsules with diameters as small as 6 nm. Additionally, Wang et al. (Reference Wang, Zhou, Wang, Hu, Li, Lin, Gao, Zhang and Wang2021b ) demonstrated that the adsorption of sufficient surfactants onto nanobubbles could extend the lifespan of nanobubble clusters. These studies indicate that by tuning the properties of the continuous phase, the dispersed phase, and the interfaces, Ostwald ripening can be effectively mitigated. Although these proposed methods may work well for ripening in an unconfined space, they may become less effective for systems with heterogeneous surfaces or complex structures, such as bubble ripening on rough or chemically treated substrates or in porous media (Xu, Bonnecaze & Balhoff Reference Xu, Bonnecaze and Balhoff2017; Wang et al. Reference Wang, Mehmani and Xu2021a ; Mehmani & Xu Reference Mehmani and Xu2022).
Owing to its inherent complexity, understanding the mechanisms of ripening on heterogeneous substrates and subsequently controlling this process has attracted significant attention from researchers. For instance, Zhu et al. (Reference Zhu, Zhang and Lohse2018) explored the diffusive dynamics of multiple surface nanobubbles, employing numerical simulations to elucidate their stability and ripening behaviour under the influence of contact line pinning and gas oversaturation. It demonstrated that pinned nanobubbles resist Ostwald ripening, while unpinned ones exhibit competitive growth, with the ripening process decelerating as bubble distances increase over time. Zhang et al. (Reference Zhang, Chen, Zhang and Zhang2020) confirmed through lattice Boltzmann modelling (LBM) that heterogeneous substrates or solid nanoparticles can act as scaffolds to suppress Ostwald ripening of droplets, with the scaffolds serving to pin the droplets. The above studies illustrate that by regulating the properties of a solid substrate, such as surface roughness, chemical homogeneity and adsorption energy, the migration pathways and dynamic behaviours of particles or droplets can be effectively altered, thereby suppressing or reversing the ripening process. However, up to now, an in-depth theoretical or modelling analysis of how the surface heterogeneity induced contact line pinning (under the constant contact radius, i.e. CCR condition) affects the ripening dynamics remains limited.
In this study, inspired by the work of Zhu et al. (Reference Zhu, Zhang and Lohse2018), we conduct a combined theoretical and numerical investigation to further investigate the ripening process of bubbles under the reaction-controlled regime, when the contact lines are pinned (under CCR) due to a heterogeneous substrate with high contact angle hysteresis. Our work not only contrasts with but also complements the diffusion-controlled mechanism studied by Zhu et al. (Reference Zhu, Zhang and Lohse2018). Our investigation is organised as follows. In § 2, we first theoretically derive the governing equation of the ripening process, then analyse the diverse ripening processes not limited to classical Ostwald ripening. Interestingly, egalitarianism of bubble curvature radii is found for various ripening processes, based on the free energy analysis. In § 3, we numerically investigated the diverse ripening processes and their final states, by applying Shan–Chen pseudo-potential two-phase LBM (Shan Reference Shan2008; Chen et al. Reference Chen, Kang, Mu, He and Tao2014; Qin et al. Reference Qin, Del Carro, Mazloomi Moqaddam, Kang, Brunschwiler, Derome and Carmeliet2019), due to its ease in automatically capturing the evolution of a liquid–gas interface. The simulation results are highly consistent with the theoretical analysis, providing valuable insights into the study of ripening on the solid substrate with CCR. Section 4 concludes the current work.
Dual bubbles ripening on the heterogeneous substrate with CCR under six different conditions.
2. Theoretical derivations and analyses
2.1. Derivation of the bubble growth rate with CCR
2.1.1. Two-dimensional case
In this subsubsection, we theoretically derive the governing equations of dual bubbles ripening on the solid surface with CCR but varied contact angle. The theoretical derivations and numerical modelling are basically conducted in two dimensions, since this can reflect the core physical mechanisms under certain geometrical constraints. For instance, Huang et al. (Reference Huang, Su, Yang, Li, Chen, Li, Zhou, Li and Song2017) proposed a universal patterning method by controlling the evolution of two-dimensional (2-D) liquid foams, where inverse Ostwald ripening, regulated by interfacial Laplace pressure, was found to dominate bubble evolution under geometric constraints of micropillars, aligning with the interface-limited characteristics of reaction-controlled ripening. Knüsel et al. (Reference Knüsel, Riedinger, Rossinelli, Ott, Mule and Norris2020) investigated 2-D Ostwald ripening in quasi-2-D CdS (cadmium sulphide, a key semiconductor material) nanoplatelets, observing thickness transitions and material transfer. Under the 2-D condition, the bubbles are placed on the solid substrate with different contact radii
$r_{B1},r_{B2}$
and contact angles
$\theta _{1},\theta _{2}$
, respectively, as shown in figure 1. The curvature radius
$R_B$
of the bubble can be expressed by the contact radius
$r_B$
and contact angle
$\theta$
as
\begin{align} R_B=\frac {r_B}{\sin \theta },\quad \theta \lt \frac {\pi }{2} \ \text{or} \ \theta \gt \frac {\pi }{2}. \end{align}
The volume and interfacial area of the 2-D bubbles can be given as
\begin{align} V_{\textit{cap}}=\left (\theta -\sin \theta \cos \theta \right )R_B^2,\quad S_{\textit{cap}}=2 \theta R_B. \end{align}
In reaction-controlled ripening (which refers to the phase-transition speed being much slower than that of mass transferring) in two dimensions, the kinetic equation can be written as
\begin{align} \frac {{{\textrm {d}}{V_{{{cap }}}}}}{{{\textrm {d}}t}} = {S_{{{cap }}}} \times k \big({C_P} - {C_{{R_B}}}\big), \end{align}
where the left-hand side represents the growth or shrinkage of the bubbles, and the right-hand side indicates the inlet or outlet of the solute (dissolved gas). Here,
$C_P$
is the actual solute concentration in the solution close to the bubble interface, and the concentration difference
$C_P-C_{R_B}$
drives the ripening kinetics, while
$C_{R_B}$
is the solute concentration on the surface of the bubble. The parameter
$k$
denotes the rate of dissolution or phase transition in the free space. Equation (2.3) illustrates that the ripening speed is determined by the concentration difference
$C_P-C_{R_B}$
. In this work, we focus on the reaction-controlled ripening as mentioned in the Introduction, i.e. the diffusion process is much faster than the interfacial reaction process (Chen et al. Reference Chen, Feng, Hu, Du and Wang2022b
). In this way, the actual concentration
$C_P$
is simplified as the average concentration
$C_m$
of the solute, i.e.
$C_P=C_m$
. To obtain the solute concentration on the surface of the bubble
$C_{R_B}$
, the linearised Gibbs–Thompson equation for the solute system is utilised, i.e. (Lifshitz & Slyozov Reference Lifshitz and Slyozov1961; Voorhees Reference Voorhees1985)
\begin{align} C_{R_B}=\left .C(r)\right |_{r=R_B}=C_{\infty }+\frac {\lambda }{R_B}, \quad \lambda =\frac {\sigma V_m C_{\infty }}{R_g T}, \end{align}
where
$C(r)$
represents the solute concentration in the solution at a distance
$r$
from the centre of the bubble,
$C_{\infty }$
refers to the solute concentration in the solution when the interface of gas and solution is planar,
$\lambda$
is the concentration variable associated with the critical radius, and
$\sigma$
,
$V_m$
,
$R_g$
and
$T$
represent surface tension, molar volume of the gas phase, gas constant and thermodynamic temperature, respectively. Similarly, the average concentration
$C_m$
can be expressed as
${C_m} = {C_\infty } + ({\lambda }/{{{R_C}}})$
, where
$R_c$
is the critical radius for the ripening process. By plugging (2.4) and
$C_P={C_m} = {C_\infty } + ({\lambda }/{{{R_C}}})$
into (2.3), we obtain
\begin{align} \frac {{{\textrm {d}}{V_{{{cap }}}}}}{{{\textrm {d}}t}} = {S_{{{cap }}}} \times k\lambda \left ( {\frac {1}{{{R_c}}} - \frac {1}{{{R_B}}}} \right )\!. \end{align}
By plugging (2.2) into (2.5), the left-hand side of the equation is derived as
\begin{align} \frac {{{\rm d}{V_{{{cap }}}}}}{{{\rm d}t}} = 2{R_B}(\theta - \sin \theta \cos \theta )\frac {{{\rm d}{R_B}}}{{{\rm d}t}} + 2R_B^2\sin ^2\theta\, \frac {{{\rm d}\theta }}{{{\rm d}t}}. \end{align}
By substituting (2.6) into (2.5) and dividing by
$2R_B$
, we can obtain
\begin{align} (\theta -\sin \theta \cos \theta ) \dot {R_B}+\sin ^2 \theta\, R_B \dot {\theta }=\theta \times k \lambda \left (\frac {1}{R_c}-\frac {1}{R_B}\right )\!. \end{align}
For the bubbles ripening with CCR, by utilising (2.1), we have
\begin{align} \dot {R_B}=\frac {{\rm d}}{{\rm d} t}\left (\frac {r_B}{\sin \theta }\right )= - \frac {{{R_B}}}{{\tan \theta }}\dot \theta . \end{align}
By plugging (2.8) into (2.7), and with some algebra, we can obtain the kinetic equation of dual bubbles ripening on the heterogeneous substrate under CCR as
\begin{align} \dot {R_B}=\frac {\theta }{\theta -\sin \theta \cos \theta -\sin ^2 \theta \tan \theta } \times k \lambda \left (\frac {1}{R_C}-\frac {1}{R_B}\right )=k_R \lambda \left (\frac {1}{R_c}-\frac {1}{R_B}\right )\!. \end{align}
In (2.9), we define the bubble growth rate
$k_R$
for bubbles ripening under CCR, which is a function of parameter
$k$
and contact angle
$\theta$
, i.e.
\begin{align} \quad k_{R,\text{2-D}}=\frac {\theta }{\theta -\sin \theta \cos \theta -\sin ^2 \theta \tan \theta }\times k. \end{align}
As introduced above, the parameter
$k$
is the fixed growth rate of bubble ripening in the free space, without the existence of a solid substrate. The obtaining of
$k\lambda$
in free space is given in Appendix A. In other words, the bubble growth rate with CCR (
$k_R$
) is determined by the instantaneous contact angle
$\theta$
and the growth rate
$k$
in free space. Since the Ostwald ripening in free space is already well studied in the literature, in this work we focus on the difference of the ripening on the heterogeneous substrate with CCR by investigating the normalised bubble growth rate
$k_R/k$
. As shown in figure 2,
$k_R/k$
increases monotonically with the contact angle
$\theta$
. We also observe an important ‘transition contact angle’, i.e.
${\theta _t} = 90^\circ$
, when the bubble growth rate
$k_R/k$
changes its sign from negative to positive (or the opposite). This transition is very interesting, since it can lead to instantaneous bubble growth to shrinkage (or the opposite) during the ripening process, making it deviate from the classical Ostwald ripening. For comparison, the bubble growth rate
$k_{\theta }$
for dual bubbles ripening on homogeneous flat surfaces with constant contact angle (CCA) is analysed in Appendix B. We can see that the transition only occurs under the CCR condition.
2.1.2. Three-dimensional case
To compare the difference of bubble growth rate in two and three dimensions, we further conduct the theoretical derivation for reaction-controlled ripening of dual bubbles under the CCR mode in three dimensions. Similarly to the 2-D case, the three-dimensional (3-D) spherical-cap-shaped bubbles are placed on the solid substrate with different contact radii
$r_{B1}, r_{B2}$
and contact angles
$\theta _1, \theta _2$
. The curvature radius
$R_B$
of the bubble can be expressed as
$R_B = {r_B}/{\sin \theta }$
, with
$\theta \lt {\pi }/{2}$
or
$\theta \gt {\pi }/{2}$
. The volume and gas–liquid interfacial area of a spherical cap in 3-D are given by
\begin{align} V_{\textit{cap}} = \frac {\pi R_B^3}{3} \big(2 - 3 \cos \theta + \cos ^3 \theta \big), \quad S_{\textit{cap}} = 2 \pi R_B^2 (1 - \cos \theta ). \end{align}
Following the procedures of (2.3)–(2.9), we obtain the effective growth rate
$k_R$
under CCR in the 3-D case as
\begin{align} k_{R,\text{3-D}} = \frac {2 \cos \theta }{ \cos \theta - 1} \times k. \end{align}
Similar to (2.9), we can obtain the kinetic equation of dual bubbles ripening under CCR in 3-D as
\begin{align} \dot {R_B}=\frac {2 \cos \theta }{ \cos \theta - 1} \times k \lambda \left (\frac {1}{R_C}-\frac {1}{R_B}\right )=k_{R,\text{3-D}} \lambda \left (\frac {1}{R_c}-\frac {1}{R_B}\right )\!. \end{align}
For comparison, the final expressions for the effective growth rate
$k_R$
in 2-D and 3-D are written as
\begin{align} k_{R,\text{3-D}} = \frac {2 \cos \theta }{ \cos \theta - 1} \times k, \quad k_{R,\text{2-D}}=\frac {\theta }{\theta -\sin \theta \cos \theta -\sin ^2 \theta \tan \theta }\times k . \end{align}
Equation (2.14) shows that in the 3-D case, the normalised growth rate
$k_R/k$
also varies with the contact angle
$\theta$
. Similar to the 2-D case, a transition point exists at
$\theta _t = 90^\circ$
, where the bubble growth rate changes sign, i.e. from shrinkage to growth or vice versa. This behaviour leads to non-classical ripening characteristics, distinct from Ostwald ripening. In figure 2, the variation trends of
$k_R/k$
with the contact angle
$\theta$
in 2-D and 3-D cases are compared. The trends of the two curves are almost identical, despite small differences in values. This indicates that it is reasonable to preliminarily study the ripening of double bubbles under the CCR mode in the 2-D case. In the next subsection, the critical radius and the final evolution equations during dual bubble ripening under CCR are derived for both the 2-D and 3-D cases.
The comparison of normalised growth rate
$k_R/k$
against contact angle
$\theta$
for bubbles ripening with CCR between 2-D and 3-D cases.
2.2. Derivation of the evolution equation for bubble curvature radii under CCR mode
2.2.1. Two-dimensional case
Since we consider the dual bubbles ripening in this work, we first derive the critical radius
$R_c$
as well as the explicit bubble radius evolution equations of a dual-bubble system. The continuity equation for the bubble distribution function
$\mathcal{F} (R_B, t )$
is written as
\begin{align} \frac {\partial \mathcal{F}}{\partial t}+\frac {\partial }{\partial R_B}\left (\mathcal{F} \frac {\mathrm{d} R_B}{\mathrm{\,d} t}\right )=0. \end{align}
To determine the critical radius
$ R_c$
, we use the volume conservation principle. Since the bubble volumes are constrained as spherical caps, and mass exchange occurs between bubbles, the total gas volume remains constant:
\begin{align} \frac {{\rm d}}{{\rm d}t} \int V_{\textit{cap}}(R_B)\, F(R_B, t) \, {\rm d}R_B = 0, \end{align}
where
$ F(R_B, t)$
is the bubble size distribution function. The volume conservation condition for the double-bubble system is
\begin{align} \frac {{\rm d}V_{\textit{cap}1}}{{\rm d}t} + \frac {{\rm d}V_{\textit{cap}2}}{{\rm d}t} = 0. \end{align}
Substituting (2.2), we have
\begin{align} 2 R_{1}\, f(\theta _1) \times\dot {R}_{1} + 2 R_{2}\, f(\theta _2) \times\dot {R}_{2} = 0, \end{align}
where
$ R_1$
and
$ R_2$
are the curvature radii of the two bubbles. To simplify, we define the geometric factors
$ f(\theta _1)$
and
$ f(\theta _2)$
:
\begin{align} f(\theta _1) = \theta _1 - \sin \theta _1 \cos \theta _1 - \sin ^2 \theta _1 \tan \theta _1, \\[-28pt] \nonumber \end{align}
\begin{align} f(\theta _2) = \theta _2 - \sin \theta _2 \cos \theta _2 - \sin ^2 \theta _2 \tan \theta _2. \\[6pt] \nonumber \end{align}
Substituting (2.9) into (2.18), we obtain
\begin{align} R_{1} \theta _1 \left ( \frac {1}{R_c} - \frac {1}{R_{1}} \right ) + R_{2} \theta _2 \left ( \frac {1}{R_c} - \frac {1}{R_{2}} \right ) = 0. \end{align}
By solving the above equation, we obtain the final expression for the critical radius as
\begin{align} R_c = \frac {R_1 \theta _1 + R_2 \theta _2}{\theta _1 + \theta _2}, \end{align}
where
$R_1$
and
$R_2$
are the curvature radii of the two bubbles. By combining (2.9) and (2.22), we can obtain the kinematic equations for the 2-D dual-bubble ripening process as
\begin{align} \dot {R_1}=k_{R,\text{2-D}}\lambda \left (\frac {1}{R_c}-\frac {1}{R_1}\right )\!,\quad \dot {R_2}=k_{R,\text{2-D}}\lambda \left (\frac {1}{R_c}-\frac {1}{R_2}\right )\!,\quad R_c=\frac {R_1 \theta _1 + R_2 \theta _2}{\theta _1 + \theta _2}. \end{align}
2.2.2. Three-dimensional case
Similarly, we derive the critical radius of double bubbles ripening under the CCR mode in the 3-D case. The analysis is still based on the volume conservation of double bubbles. Substituting (2.11) into (2.17), we have
\begin{align} \pi R_{1}^2\, F(\theta _1) \times\dot {R}_{1} + \pi R_{2}^2\, F(\theta _2) \times\dot {R}_{2} = 0, \end{align}
where
$ R_1$
and
$ R_2$
are the curvature radii of the two bubbles. To simplify, we define the geometric factor
$ F(\theta )$
:
\begin{align} F(\theta _1) = \big(2 - 3 \cos \theta _1 + \cos ^3 \theta _1 \big) - \frac {1}{3} \big(3 \sin \theta _1 \tan \theta _1 - 3 \cos \theta _1 \sin ^2 \theta _1 \big), \\[-28pt] \nonumber \end{align}
\begin{align} F(\theta _2) = \big(2 - 3 \cos \theta _2 + \cos ^3 \theta _2 \big) - \frac {1}{3} \big(3 \sin \theta _2 \tan \theta _2 - 3 \cos \theta _2 \sin ^2 \theta _2 \big). \\[6pt] \nonumber \end{align}
Substituting (2.13) into (2.24), we obtain
\begin{align} R_{1}^2 (1 - \cos \theta _1) \left ( \frac {1}{R_c} - \frac {1}{R_{1}} \right ) + R_{2}^2 (1 - \cos \theta _2) \left ( \frac {1}{R_c} - \frac {1}{R_{2}} \right ) = 0. \end{align}
By solving this equation, we finally obtain the expression for the critical radius under 3-D conditions as
\begin{align} R_c = \frac {R_1^2 (1 - \cos \theta _1) + R_2^2 (1 - \cos \theta _2)}{R_1 (1 - \cos \theta _1) + R_2 (1 - \cos \theta _2)}. \end{align}
The kinematic equations for the 3-D dual-bubble ripening process are written as
\begin{align} \dot {R_1} & =k_{R,\text{3-D}}\lambda \left (\frac {1}{R_c}-\frac {1}{R_1}\right )\!,\quad \dot {R_2}=k_{R,\text{3-D}}\lambda \left (\frac {1}{R_c}-\frac {1}{R_2}\right )\!, \nonumber \\ R_c & = \frac {R_1^2 (1 - \cos \theta _1) + R_2^2 (1 - \cos \theta _2)}{R_1 (1 - \cos \theta _1) + R_2 (1 - \cos \theta _2)}. \end{align}
In both 2-D and 3-D cases,
$R_c$
satisfies
$\min (R_1, R_2 ) \lt R_c \lt \max (R_1, R_2 )$
. Furthermore, when the radius ratio satisfies
$ R_1/R_2 \in [1, 2]$
, for all the contact angles
$ \theta _1, \theta _2 \in [0^\circ , 180^\circ ]$
, the ratio between the critical radii in the 2-D and 3-D cases lies within the interval
$[0.8749, 1.0863]$
. This indicates that the critical radii in the 2-D and 3-D cases are relatively close within a certain range. Figure 3 presents the ratio of 2-D to 3-D critical radii
$R_c$
for different contact angles of a two-bubble system with
$R_1/R_2 = 2$
, highlighting that the dimensional effect is very limited on the critical radii.
Ratio of the critical radii
$R_c$
in two and three dimensions for different contact angles of two bubbles when
$R_1/R_2 = 2$
.
2.3. Ripening process under diverse bubble conditions
In § 2.2, both (2.23) and (2.29) illustrate that the bubble growth or shrinkage (
$\dot {R_B}\gt 0$
or
$\dot {R_B}\lt 0$
) is the mutual effect of the growth rate (
$k_R$
) and the relative bubble radius compared to the critical radius (
${1}/{R_c}-{1}/{R_B}$
). As introduced in § 2.1, the growth rate
$k_R$
is related to the bubble contact angle
$\theta$
. To thoroughly investigate the dual-bubble ripening on the heterogeneous substrate with CCR, we summarise all the possible conditions of the bubbles with different contact angles and contact radii, as illustrated in figure 1. To better describe and analyse the evolution process of the bubble curvature radius
$R_B$
, we define the following three terminologies to represent the corresponding ripening processes.
-
(i) Ostwald ripening (O-R): the
$R_B$
of the bubble with larger
$R_B$
increases, while the
$R_B$
of the bubble with smaller
$R_B$
reduces. -
(ii) Reversed ripening (R-R): the
$R_B$
of the bubble with larger
$R_B$
reduces, while the
$R_B$
of the bubble with smaller curvature increases. -
(iii) Consistent ripening (C-R): the
$R_B$
values of the two bubbles increase or reduce at the same time.
We note that here we focus on the instantaneous ripening status under the six conditions initially considered in figure 1. However, the ripening status may transit from one to another when the contact angle of the bubble crosses
$90^\circ$
or the relative difference of bubble radii reciprocals
${1}/{R_c}-{1}/{R_B}$
changes the sign. With the three defined processes, the ripening processes of dual bubbles under the above-categorised six conditions are analysed as follows.
(a,b) For the conditions in figure 1(a,b), the R-R regime occurs:
\begin{align} \begin{array}{l} \kern100pt {\theta _1} \lt {\theta _2} \lt \dfrac {\pi }{2}\quad \Rightarrow\quad {k_{R1}} \lt {k_{R2}} \lt 0,\\[3pt](\textit{a})\quad{R_1} \lt {R_c} \lt {R_2}\quad \Rightarrow\quad \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \gt 0,\ \dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \lt 0,\\[12pt](\textit{b})\quad {R_1} \gt {R_c} \gt {R_2}\quad \Rightarrow\quad \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \lt 0,\ \dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \gt 0. \end{array} \end{align}
(c,d) For the conditions in figure 1(c,d), the O-R regime occurs:
\begin{align} \begin{array}{l} \kern100pt {\theta _2} \gt {\theta _1} \gt \dfrac {\pi }{2}\quad \Rightarrow\quad {k_{R2}} \gt {k_{R1}} \gt 0,\\[3pt](\textit{c})\quad {R_1} \lt {R_c} \lt {R_2}\quad \Rightarrow\quad \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \lt 0,\ \dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \gt 0,\\[12pt](\textit{d})\quad{R_1} \gt {R_c} \gt {R_2}\quad \Rightarrow\quad \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \gt 0,\ \dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \lt 0. \end{array} \end{align}
(e, f ) For the conditions in figure 1(e,f), the C-R regime occurs:
\begin{align} \begin{array}{l} \kern100pt {\theta _1} \gt \dfrac {\pi }{2} \gt {\theta _2}\quad \Rightarrow\quad {k_{R1}} \gt 0 \gt {k_{R2}} ,\\[3pt](\textit{e})\quad {R_1} \lt {R_c} \lt {R_2} \Rightarrow \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \lt 0,\dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \lt 0,\\[12pt](\textit{f})\quad{R_1} \gt {R_c} \gt {R_2} \Rightarrow \dfrac {{{\rm d}{R_1}}}{{{\rm d}t}} = {k_{R1}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_1}}}} \right ) \gt 0,\dfrac {{{\rm d}{R_2}}}{{{\rm d}t}} = {k_{R2}}\left ( {\dfrac {1}{{{R_c}}} - \dfrac {1}{{{R_2}}}} \right ) \gt 0. \end{array} \end{align}
The theoretical analysis under the 3-D condition yields results consistent with those in the 2-D case. From the above theoretical analysis, we can see that the bubbles ripening on heterogeneous substrates with CCR may experience three different processes of Ostwald ripening, reversed ripening or consistent ripening, depending on the different bubble conditions. Namely, the classical theory of Ostwald ripening in the free surface is not sufficient to categorise the complex processes of bubbles ripening with CCR. The surprising observations lead to a more challenging problem, i.e. how could we unify the three different ripening processes considering all the six different bubble conditions?
2.4. The egalitarianism of bubble curvature radii
To unify the various ripening processes under diverse conditions, we propose two analysis methods, from the aspects of driving chemical potential and minimisation of free energy. From the first aspect, we understand that ripening is driven by the chemical potential difference between the two bubbles. Namely, the gas in the bubble with high chemical potential dissolves into the solution and is then transported to the bubble with lower chemical potential. The chemical potential can be represented as the gas pressure
$p_g=p_s+p_L$
in the bubble, where
$p_s$
and
$p_L=\sigma /R_B$
are the pressure in the solution and the Laplace pressure due to bubble curvature, respectively. For the dual-bubble system, the solution pressure
$p_s$
is uniform and the surface tension
$\sigma$
is fixed. Therefore, the pressure difference of the bubbles can only come from the different curvature radii
$R_B$
. In the current study of dual bubbles ripening on the heterogeneous substrate with CCR, the contact angles vary with time, leading to simultaneous variations of
$R_B$
. It is straightforward that the final equilibrium state can only be reached at two conditions: (i) when
$R_B$
values of the dual bubbles become equal and the driving chemical potential vanishes, i.e. the egalitarianism of bubble curvature radii is obtained finally; (ii) when one of the bubbles grows to the maximum at the expense of the extinction of the other bubble. Here, we explain why the second condition would not occur. The dimensionless interfacial area (
$S_{cap}/r_B$
) against volume (
$V_{cap}/r_B^2$
) is shown in figure 4. This curve indicates how the interfacial free energy varies with the bubble volume under CCR but varied contact angle. For a dual-bubble system, if one bubble (A) vanishes, then it approximates to the coordinate of (
$V_{cap}/r_B^2=0, S_{cap}/r_B=2$
) with slope zero, indicating that the change of interfacial free energy is zero. However, the volume of the other bubble (B) has to increase since the total volume is conserved, making the free energy increase at a certain rate (positive slope). In this way, the total free energy of the dual-bubble system increases, which violates the free energy minimisation law. Therefore, it is impossible that one bubble finally vanishes; instead, the
$R_B$
values of the two bubbles become equal.
The relation between dimensionless interfacial area and volume of a bubble with CCR but varied contact angle.
From the second aspect, we analyse the free energy evolution of the dual-bubble system. Under the condition of heterogeneous substrate with CCR, the bubble states are characterised by their contact angles
$\theta _1$
,
$\theta _2$
and contact radii
$r_1$
,
$r_2$
. According to the second law of thermodynamics (Greiner et al. Reference Greiner, Neise and Stöcker1995), the Gibbs free energy of the system gradually decreases until the minimum is reached. Moreover, the overall bubble volume during the ripening process remains unchanged. This is a typical minimisation problem under fixed constraints, which could be solved by the method of Lagrange multipliers (Hoffmann & Bradley Reference Hoffmann and Bradley2003). Since the theoretical derivations in the above sections show that the ripening mechanisms for 2-D and 3-D cases are very similar, in the following work, we only consider the 2-D results and analyses. The objective function
$f$
and the constraint
$g$
are given by
\begin{align} \begin{array}{*{20}{c}} {\min \left \{ {f({\theta _1},{\theta _2})} \right \}\!,\quad f({\theta _1},{\theta _2}) = \dfrac {{2{r_1}{\theta _1}}}{{\sin {\theta _1}}} + \dfrac {{2{r_2}{\theta _2}}}{{\sin {\theta _2}}}},\\[12pt]{\text{s.t.} \quad g({\theta _1},{\theta _2}) = \dfrac {{\left ( {{\theta _1} - \sin {\theta _1}\cos {\theta _1}} \right )r_1^2}}{{\sin ^2 \theta _1}} + \dfrac {{\left ( {{\theta _2} - \sin {\theta _2}\cos {\theta _2}} \right )r_2^2}}{{\sin ^2 \theta _2}} - {C_0} = 0}, \end{array} \end{align}
where
$C_0$
is the initial overall volume of the dual bubbles. By substituting
$f$
and
$g$
in (2.33) into the Lagrangian function
$\mathcal{L}(\theta _1, \theta _2, \lambda ) = f(\theta _1, \theta _2) + \lambda \times g(\theta _1, \theta _2)$
, and calculating its derivative against the contact angles
$\theta _1$
,
$\theta _2$
and the Lagrange multiplier
$\lambda$
, we obtain the equations
\begin{align} \begin{cases} \lambda r_1^2\big[(1 - \cos 2{\theta _1})\sin ^2{\theta _1} - \sin 2{\theta _1}({\theta _1} - 0.5\sin 2{\theta _1})\big] \\ \quad + 2{r_1}\sin ^2{\theta _1}(\sin {\theta _1} - {\theta _1}\cos {\theta _1}) = 0,\\[6pt]\lambda r_2^2\big[(1 - \cos 2{\theta _2})\sin ^2{\theta _2} - \sin 2{\theta _2}({\theta _2} - 0.5\sin 2{\theta _2})\big] \\ \quad + 2{r_2}{\sin ^2}{\theta _2}(\sin {\theta _2} - {\theta _2}\cos {\theta _2}) = 0,\\[6pt]r_1^2({\theta _1} - 0.5\sin 2{\theta _1}){\sin ^2}{\theta _2} + r_2^2({\theta _2} - 0.5\sin 2{\theta _2}){\sin ^2}{\theta _1} \\ \quad - {C_0}{\sin ^2}{\theta _1}{\sin ^2}{\theta _2} =0. \end{cases} \end{align}
Unfortunately, due to the complexity and nonlinearity of (2.34), the direct obtaining of the analytical solutions for
$\theta _1$
and
$\theta _2$
is currently infeasible. Despite this, an inspiring piece of news is that we can approximately obtain the solution using the numerical method of sequential quadratic programming (Nocedal & Wright Reference Nocedal and Wright2006). In the next section, by comparison with LBM in 2-D cases, we will show that the curvature radii (
$R_i = r_i / \sin \theta _i$
) obtained from the numerical solution of (2.34) strictly achieve egalitarianism, as we analysed from the aspect of chemical potential at the beginning of this subsection.
3. Numerical modelling by two-phase LBM
In this section, we apply the two-phase LBM (Qin et al. Reference Qin, Fei, Zhao, Kang, Derome and Carmeliet2023) to investigate the dual-bubble ripening process. The model is first validated with Ostwald ripening of dual bubbles in the free space and on the homogeneous substrate with CCA, as given in Appendix B. Afterwards, we model the ripening of the dual-bubble system on the heterogeneous substrate with CCR under the six initial conditions as shown in figure 1. The initial set-ups of bubble contact angles
$\theta$
, contact radii
$r$
and curvature radii
$R$
under the six conditions are also given in table 1. Figure 5 illustrates the evolution of normalised growth rate
$k_R/k$
against contact angle
$\theta$
for the two bubbles under the six conditions, which generally agrees well with the theoretical equation (2.10) as we proposed. The small deviations between the simulations and theoretical predictions are mainly due to two reasons. First, the initial numerical oscillations make the bubble deviate from the spherical shape to some extent, as shown at approximately
$\theta = 112^\circ$
in figure 5. Second, the contact angle scheme to realise contact angle hysteresis is not sufficiently accurate at low bubble contact angle (Ding & Spelt Reference Ding and Spelt2008), as shown at
$\theta \in (50^\circ ,60^\circ )$
in figure 5.
Comparison of bubble contact angle (
$\theta$
), contact radii (
$r$
) and curvature radii (
$R$
) between LBM simulation and theory after ripening on the heterogeneous substrate with CCR under diverse initial conditions.
Evolution of normalised growth rate
$k_R/k$
against contact angle
$\theta$
during dual-bubble ripening on the heterogeneous substrate with CCR under six different conditions.
Comparison of bubble profiles between LBM simulation and theory after ripening on the heterogeneous substrate with CCR under diverse conditions.
After ripening completion, the final bubble profiles under the six conditions are shown in figure 6. The initial bubble profiles are indicated with green solid curves, while the final bubble profiles obtained by LBM simulations and theoretical predictions are indicated with blue solid and red dashed curves, respectively. Figure 6 shows that excellent agreement is observed between the LBM simulations and theoretical predictions. As we analysed in § 2.4, the theoretical predictions of the final profiles can be obtained by solving the Lagrangian function of (2.34) to reach the minimum free energy of the system. In the current LBM simulations, we also recorded the normalised change of the total free interfacial area (
$({{S - {S_0}}})/{{{S_0}}}$
), which is an indicator of the system Gibbs free energy evolution during the ripening process, as shown in figure 7. We can see under all the six different conditions considered,
$({{S - {S_0}}})/{{{S_0}}}$
first gradually decreases and then reaches a minimum value, and remains approximately unchanged afterwards. Figure 7 shows that the free energy evolution governed by (2.34) is well reproduced by the LBM simulations.
The change of total interfacial area of the dual-bubble system during ripening with CCR under different conditions.
Finally, we investigate the curvature radii
$R_B$
of the two bubbles after the ripening process. In figure 8, the initial and final
$R_B$
values, as well as their evolution processes, are illustrated for the six conditions considered in figure 6. We clearly observe that the bubbles initially experienced O-R in cases 1 and 2, R-R in cases 3 and 4, and C-R in cases 5 and 6. Despite different ripening processes, the curvature radii
$R_B$
of the dual bubbles finally reached egalitarianism, i.e. the final
$R_1$
and
$R_2$
situated on the linear line of
$R_1=R_2$
, as shown in figure 8. The exact values of the final curvature radii
$R_B$
, as well as the contact angles
$\theta$
of the six cases, are given in table 1, where the maximum error is only 2.58 %, while the average error is as low as 0.60 %, which mutually validated both the accuracies of LBM and the theoretical predictions. To further validate the proposed theory, we conduct 13 more simulations with different bubble contact radii (
$r_b$
) and contact angles (
$\theta$
). As shown in figure 8, after experiencing various ripening processes (one, two or three ripening processes), the final curvature radii of the two bubbles also reached egalitarianism at
$R_1=R_2$
. We further emphasise that the various ripening processes analysed in § 2.3 may transit from one to another. For instance, in case 8, the bubbles first experienced the O-R process from point A to point B, then experienced C-R from B to C, and finally experienced R-R from C to D to reach egalitarianism.
The evolutions of curvature radii
$R_1$
and
$R_2$
during dual-bubble ripening with CCR undergoing (a) one process, (b) two processes, (c) three processes.
4. Conclusions
We investigated the physical processes of dual bubbles ripening on a solid substrate under constant contact radius (CCR) mode. Based on the kinetic theory, we first derived a theoretical equation for the bubble curvature radius. Then we identified three distinct phenomena under six bubble conditions: classical Ostwald ripening, reversed ripening and consistent ripening. Following the principle of minimum free energy, all three ripening processes are found to achieve the egalitarianism of bubble curvature radii. Furthermore, under different two-bubble configurations, the three ripening processes can dynamically transition from one to another during bubble evolution. The two-phase LBM simulations validated the theoretical predictions and analyses very well. This study provides strong theoretical and numerical foundations for exploring complex bubble ripening phenomena on solid substrates. Crucially, our study addresses the theoretical gap concerning bubble ripening under the combined effects of interface control and CCR on a substrate, offering insights and predictions to guide future experimental investigations in this less-explored regime.
Future work directions could include extending this analysis to 3-D systems to better approximate real-world bubble geometries. Another important avenue for future research is to conduct a more in-depth theoretical analysis of diffusion-controlled ripening of dual bubbles under the CCR mode, and its validation through LBM simulations. Although some numerical studies have already been conducted on attached bubbles under the diffusion-controlled regime (Zhu et al. Reference Zhu, Zhang and Lohse2018), a systematic exploration of bubbles ripening under CCR in the diffusion-controlled regime, particularly at the theoretical level as conducted in this work, is still required. Such studies would complement the present work and contribute to a comprehensive understanding of the impact of the CCR constraint under different ripening mechanisms.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10852.
Acknowledgements
The authors acknowledge Professor C. Sun at Tsinghua University for valuable discussions. They acknowledge the use of Google’s large language model for assistance in the preparation of this paper. They take full responsibility for the accuracy and originality of the final content.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos 11872315, 12302337, 12388101), the Guangdong Basic and Applied Basic Research Foundation (grant no. 2022A1515011201), Postdoctoral Special Project (grant no. K5002), Shaanxi Province Science Foundation for Youths (grant no. K1004), the 111 Project of China (grant no. B17037), and the Fundamental Research Funds for the Central Universities (grant no. G2022KY05105).
Declaration of interests
The authors report no conflict of interest.
Author contributions
F.Q.: theoretical derivation, numerical modelling, data analysis, original manuscript, revision, funding. Q.X.: numerical modelling, data analysis, original manuscript, revision. X.S.: data analysis, revision, supervision. X.C.: conception, data analysis, original manuscript, revision, supervision, funding. H.H.: revision, supervision.
Appendix A. Ostwald ripening of dual bubbles in free space
In reaction-controlled ripening in 2-D free space, the kinetic equation governing the ripening process reads as
\begin{align} \frac {{{\textrm {d}}{V_B}}}{{{\textrm {d}}t}} = {S_B} \times k\big ( {{C_P} - {C_{{R_B}}}} \big ), \end{align}
where
$V_B=\pi R_B^2$
and
$S_B=2 \pi R_B$
are bubble volume and interfacial area, respectively. Here,
$R_B$
denotes the radius of the bubble,
$C_P$
is the actual solute concentration in the solution near the bubble, and
$C_{R_B}$
is the solute concentration on the surface of the bubble. The parameter
$k$
denotes the rate of dissolution or phase transition. Similar to the derivations in § 2.1, we can obtain the governing equation for the dual bubbles ripening in the free space as
\begin{align} \dot {R_1}=k\lambda \left (\frac {1}{R_c}-\frac {1}{R_1}\right )\!,\quad \dot {R_2}=k\lambda \left (\frac {1}{R_c}-\frac {1}{R_2}\right )\!,\quad R_c=\frac {R_1+R_2}{2}. \end{align}
Equation (A2) strictly follows the Lifshitz–Slyozov–Wagner theory (Voorhees Reference Voorhees1985; Bray Reference Bray1994).
Results of dual bubbles ripening in the free space by LBM. (a) The evolution of bubble radii against simulation time. The circles represent the original data, while the curves represent the spline interpolation. (b) Linear relation between the bubble radius change rate and the driving force, i.e. the difference of the reciprocals of critical radius
$R_c$
and current bubble radius
$R$
.
To validate the theoretical solution as given in (A2), we conduct a lattice Boltzmann simulation (Qin et al. Reference Qin, Fei, Zhao, Kang, Derome and Carmeliet2023) to investigate the ripening process of dual bubbles in the free space. The computational domain is set as
$1000 \times 500$
$\textit{lattice}^2$
. The initial radii of the two bubbles are set as
$R_1=57$
and
$R_2=60$
lattices, respectively. In the initial stage of the simulation, as shown in figure 9(a), the bubble radii illustrate some oscillations. But they gradually decayed after iterations
$t=0.5\times 10^4$
. Afterwards, the larger bubble with
$R_2=60$
gradually grows, while the smaller bubble with
$R_1=57$
gradually diminishes. Until iteration
$t=2.2\times 10^4$
, the smaller bubble completely dissolves into the solution, and the larger bubble reaches the maximum radius approximately
$R_{2,max}=82$
lattices.
Quantitatively, we plotted the bubble radii change rate against the driving force, i.e. the difference of the reciprocals of critical radius
$R_c$
and current bubble radius
$R$
. As shown in figure 9(b), a linear dependence between
${{{\rm d}R}}/{{{\rm d}t}}$
and
${1}/{{{R_c}}} - {1}/{{{R}}}$
is obtained for both bubbles. The result agrees well with the theoretical solution of (A2). Moreover, we also obtained the linear coefficient
$k\lambda = 0.52$
in (A2) for the current simulation. For a given system, the coefficients
$k$
and
$\lambda$
are fixed. Therefore, they can be further used in the paper for the comparison of ripening in the free space and on the solid substrate.
Appendix B. Ripening of dual bubbles on the homogeneous substrate with CCA
Due to different surface wettabilities, bubbles or droplets may exhibit different contact angles when they sit on a solid surface. When the surface is smooth without chemical treatment, the contact angle is normally constant. In this appendix, we study the dual bubbles ripening on the homogeneous substrate with CCA. Similar to the derivations in Appendix A, here we derive the reaction-controlled ripening under CCA for a two-bubble system. Under this condition, the 2-D bubbles are regarded as deficient circles. By analogising to (A1), we can write the kinetic equation as
\begin{align} \frac {\mathrm{d} V_{\textit{cap}}}{\mathrm{d} t}=S_{\textit{cap}} \times k\lambda \left (\frac {1}{R_c}-\frac {1}{R_B}\right )\!, \end{align}
where the volume and interfacial area of the deficient circle are
\begin{align} V_{\textit{cap}}=\theta R_B{ }^2-\sin \theta \cos \theta\, R_B{ }^2,\quad S_{\textit{cap}}=2 \theta R_B. \end{align}
Similarly, we can obtain the growth rate of the bubble as
\begin{align} \dot {R_B}=\frac {\theta k \lambda }{\theta -\sin \theta \cos \theta } \left (\frac {1}{R_c}-\frac {1}{R_B}\right )=k_{\theta } \lambda \left (\frac {1}{R_c}-\frac {1}{R_B}\right )\!. \end{align}
By comparing the following (B4) with (A2) above, we can see the bubble growth rate
$k_{\theta }$
under the CCA model is a function of
$k$
in the free space and the contact angle
$\theta$
, i.e.
\begin{align} k_{\theta }=\frac {\theta }{\theta -\sin \theta \cos \theta }\times k. \end{align}
We note that since the contact angle
$\theta$
is a fixed value in CCA mode, the bubble growth rate
$k_{\theta }$
is also a fixed value related to
$\theta$
. The variation of normalised growth rate
$k_{\theta }/k$
against the contact angle
$\theta$
is shown in figure 10, where we observe a monotonic decreasing trend.
Comparison of the normalised growth rate
$k_{\theta }/k$
against contact angle
$\theta$
for bubbles ripening on the homogeneous substrate with CCA.
Similar to Appendix A, we also conducted LBM simulations to investigate the dual bubbles ripening on the homogeneous substrate with CCA. The simulation set-ups are the same, and we only vary the bubble contact angle. Five cases for
$\theta$
, i.e.
$33^\circ$
,
$60^\circ$
,
$90^\circ$
to
$121^\circ$
and
$160^\circ$
, are conducted. The obtained normalised bubble growth rates
$k_{\theta }/k$
are shown in figure 10. From the comparison between theoretical prediction (B4) and the current simulation results, we can see that excellent agreement is achieved.
