1. Introduction
The study of transient gaseous bulk nanobubbles is of interest across various fields, such as water electrolysis, catalysis involving dissolved gases, water treatment, agriculture, medicine, and the flotation of fine particles (Zhu et al. Reference Zhu, An, Alheshibri, Liu, Terpstra, Liu and Craig2016; Wu et al. Reference Wu, Lyu, Yue, Tonoli, Verderio, Ma and Pan2019; Batchelor et al. Reference Batchelor, Armistead, Ingram, Peyman, McLaughlan, Coletta and Evans2022; Li & Zhang Reference Li and Zhang2022; Zhang et al. Reference Zhang, Zhu, Wood and Lohse2024; Yadav, Nirmalkar & Ohl Reference Yadav, Nirmalkar and Ohl2024). For instance, in the flotation process, nanobubbles offer a higher surface area to volume ratio and longer contact times for attachment to particles compared to macroscale bubbles, thereby enhancing particle–bubble collision and attachment efficiency (Li & Zhang Reference Li and Zhang2022). In medicine therapy, drugs or therapeutic agents can be encapsulated within the shell of coated nanobubbles and then be released in targeted positions (Batchelor et al. Reference Batchelor, Armistead, Ingram, Peyman, McLaughlan, Coletta and Evans2022). Therefore, achieving active control over the movements of bulk nanobubbles is essential for improving the efficiency of applications involving these bubbles.
For free surface flows, the gradient of surface tension can often lead to considerable shear stress at the free surface, significantly altering flow behaviours. A concentration gradient of solutes in a solvent can cause a surface tension gradient, and the resulting flow is known as solutal Marangoni flow. For instance, a free surface can be contaminated by surface-active agents such as surfactants. When the transport of surfactants on the free surface is dominated by advection, indicated by a relatively large Marangoni number or Péclet number, the distribution of surfactants tends to be non-uniform. This non-uniformity results in a surface tension gradient and a shear stress on the free surface (Manikantan & Squires Reference Manikantan and Squires2020). This Marangoni stress often acts against the advection flow, making surfactants crucial in stabilizing emulsions (Parker, Claesson & Attard Reference Parker, Claesson and Attard1994; Constante-Amores et al. Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui, Juric and Matar2021), damping waves on ocean surfaces (Erinin et al. Reference Erinin, Liu, Liu, Mostert, Deike and Duncan2023; Lohse Reference Lohse2023), translating drops and bubbles (Chen & Stebe Reference Chen and Stebe1997; Lohse & Zhang Reference Lohse and Zhang2020), and deteriorating superhydrophobic surfaces (Rodriguez-Broadbent & Crowdy Reference Rodriguez-Broadbent and Crowdy2023). Notably, the effects of surfactants on the free surface, as described by the Boussinesq–Scriven model (Scriven Reference Scriven1960), not only alter the surface tension but also introduce a new property, i.e. surface viscosity (Zhang & Ding Reference Zhang and Ding2023). Apart from surfactants, ions in water can also increase or decrease the surface tension according to the empirical Hofmeister series (Boström et al. Reference Boström, Williams and Ninham2001). Recent studies have shown that ion-related solutal Marangoni effects significantly impact the coalescence of hydrogen bubbles and thereby the current density during water electrolysis (Park et al. Reference Park, Liu, Demirkır, van der Heijden, Lohse, Krug and Koper2023; Liu et al. Reference Liu, Manica, Liu, Xu, Klaseboer and Yang2023).
In non-isothermal systems, a surface tension gradient and the resulting Marangoni flow can arise from temperature variations on the free surface, as the surface tension of a liquid can change significantly with temperature. Thermal Marangoni effects are common in processes such as water electrolysis (Park et al. Reference Park, Liu, Demirkır, van der Heijden, Lohse, Krug and Koper2023; Bashkatov et al. Reference Bashkatov, Babich, Hossain, Yang, Mutschke and Eckert2023), drop evaporation (Ristenpart et al. Reference Ristenpart, Kim, Domingues, Wan and Stone2007; Shiri et al. Reference Shiri, Sinha, Baumgartner and Cira2021; Kant et al. Reference Kant, Souzy, Kim, Van Der Meer and Lohse2024), vapour bubble detachment during boiling (Christopher & Wang Reference Christopher and Wang2001), the stability of liquid films (Kalliadasis, Kiyashko & Demekhin Reference Kalliadasis, Kiyashko and Demekhin2003), the motion of drops or bubbles in liquids (Young, Goldstein & Block Reference Young, Goldstein and Block1959), and plasmonic bubbles (Namura et al. Reference Namura, Nakajima, Kimura and Suzuki2015; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021). The surface tension of most liquids decreases with increasing temperature. However, there are so-called ‘self-rewetting’ fluids whose surface tension increases with increasing temperature. During boiling, where vapour bubbles are produced on the heating solid, these ‘self-rewetting’ fluids flow spontaneously to hotter regions on the solid due to the Marangoni force, removing dry spots on the solid and thus enhancing heat transfer (Abe Reference Abe2006; Hu et al. Reference Hu, Liu, Li and Wang2014). In many non-isothermal and multicomponent liquid systems, thermal Marangoni and solutal Marangoni effects coexist. For example, during water electrolysis, Ohmic heating in the electrolyte and the electrode generates a hot spot near the electrode. This results in a thermal Marangoni force that keeps bubbles attached to the electrode (Bashkatov et al. Reference Bashkatov, Babich, Hossain, Yang, Mutschke and Eckert2023). In the meantime, solutal Marangoni forces due to ions may generate forces whose directions depend on the type of ions (Park et al. Reference Park, Liu, Demirkır, van der Heijden, Lohse, Krug and Koper2023). Another example is the periodic bouncing of plasmonic bubbles produced in a binary liquid, which is due to the competition between solutal and thermal Marangoni forces (Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021).
The study of macroscopic bubble motion by thermal Marangoni forces was pioneered by Young et al. (Reference Young, Goldstein and Block1959) (referred to as YGB59 theory hereafter). Through experiments, they showed that small bubbles could be held stationary or even driven downwards against gravity by a strong negative temperature gradient. This is due to the competition between the Marangoni force and the buoyancy force. They also analytically derived the bubble velocity, under the assumption of small Reynolds numbers, small Marangoni numbers, and a spherical bubble shape. Later, this analytical solution was examined extensively by other experiments (Thompson & DeWitt Reference Thompson and DeWitt1979; Hardy Reference Hardy1979; Thompson, DeWitt & Labus Reference Thompson, DeWitt and Labus1980; Morick & Woermann Reference Morick and Woermann1993; Treuner et al. Reference Treuner, Galindo, Gerbeth, Langbein and Rath1996). Under terrestrial conditions, the Marangoni effect is often masked by the buoyancy force exerted on the bubble due to gravity. Special interest is thus given to experiments of Marangoni-driven bubbles in micro-gravity environments (Thompson & DeWitt Reference Thompson and DeWitt1979; Thompson et al. Reference Thompson, DeWitt and Labus1980; Treuner et al. Reference Treuner, Galindo, Gerbeth, Langbein and Rath1996) where the buoyancy force is considerably deSimulated three cases with different settings.creased. However, as the YGB59 work did not present the expressions of Marangoni forces and drag forces exerted on the bubble, approximate expressions were used in literature (Morick & Woermann Reference Morick and Woermann1993; Lubetkin Reference Lubetkin2003; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021). The YGB59 solution was also generalized to large Reynolds numbers (Balasubramaniam & Chai Reference Balasubramaniam and Chai1987) and large Marangoni numbers (Balasubramaniam & Subramanian Reference Balasubramaniam and Subramanian2000). With the development of computational methods and power, numerical simulations of the Navier–Stokes equations and the heat advection–diffusion equation were performed (Haj-Hariri, Shi & Borhan Reference Haj-Hariri, Shi and Borhan1997; Liu & Zhang Reference Liu and Zhang2015; Abu-Al-Saud, Popinet & Tchelepi Reference Abu-Al-Saud, Popinet and Tchelepi2018; Meulenbroek, Vreman & Deen Reference Meulenbroek, Vreman and Deen2021) to study thermocapillary motions of bubbles under shape deformations in terms of large Reynolds numbers.
Bulk nanobubbles, characterized by radii smaller than a micrometre, exhibit behaviour similar to bubbles in reduced-gravity environments, as they are largely unaffected by gravity. Consequently, Marangoni forces may play a dominant role in controlling their motion. While there are extensive studies on the behaviour of macroscopic bubbles driven by temperature gradients, research on nanobubbles influenced by such gradients remains relatively sparse. Nanoscale heat transfer and flow can often behave in a manner that is different from that at the macroscale. For example, interfacial thermal resistance plays a crucial role in nanoscale heat transfer problems. Many works have shown the existence of a temperature jump across the liquid–solid surface (Barrat & Chiaruttini Reference Barrat and Chiaruttini2003; Alosious et al. Reference Alosious, Kannam, Sathian and Todd2019; Hadjiconstantinou & Swisher Reference Hadjiconstantinou and Swisher2022), the liquid–liquid interface, and the liquid–gas interface (Patel, Garde & Keblinski Reference Patel, Garde and Keblinski2005; Plascencia, Bird & Liang Reference Plascencia, Bird and Liang2022; Dockar, Gibelli & Borg Reference Dockar, Gibelli and Borg2023). The slip length has also been shown to play a major role in microflows and nanoflows (Lauga, Brenner & Stone Reference Lauga, Brenner and Stone2007; Bocquet & Charlaix Reference Bocquet and Charlaix2010). Compared to the well-studied liquid–solid slip, recently there has been more interest in studying the slip length at the liquid–liquid interface (Poesio, Damone & Matar Reference Poesio, Damone and Matar2017; Telari, Tinti & Giacomello Reference Telari, Tinti and Giacomello2022; Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023). However, it remains unclear whether the YGB59 theory (Young et al. Reference Young, Goldstein and Block1959), which assumes temperature and velocity continuity at the liquid–gas interface, works at the nanoscale.
Due to the small spatial and temporal scales of nanobubbles, and due to the unavoidable contaminants, studying clean nanobubble dynamics in experiments remains challenging (Lohse & Zhang Reference Lohse and Zhang2015). As a result, molecular dynamics simulations have become a valuable tool for investigating nanobubbles, including nanobubble evolution in water electrolysis (Maheshwari et al. Reference Maheshwari, van der Hoef, Zhang and Lohse2016; Zhang & Lohse Reference Zhang and Lohse2023; Zhang et al. Reference Zhang, Zhu, Wood and Lohse2024), nanobubble cavitation (Vedadi et al. Reference Vedadi, Choubey, Nomura, Kalia, Nakano, Vashishta and Van Duin2010; Shekhar et al. Reference Shekhar, Nomura, Kalia, Nakano and Vashishta2013; Menzl et al. Reference Menzl, Gonzalez, Geiger, Caupin, Abascal, Valeriani and Dellago2016; Gao, Wu & Wang Reference Gao, Wu and Wang2021), vapour nanobubble nucleation (Gallo et al. Reference Gallo, Magaletti, Cocco and Casciola2020; Tinti et al. Reference Tinti, Giacomello, Meloni and Casciola2023) and the coalescence of bulk nanobubbles (Bird, Smith & Liang Reference Bird, Smith and Liang2021). Note that clean bulk nanobubbles in saturated or undersaturated environments tend to dissolve within microseconds, as predicted by the Epstein–Plesset theory (Epstein & Plesset Reference Epstein and Plesset1950; Lohse & Zhang Reference Lohse and Zhang2015). However, compared to the time scales of other nanobubble hydrodynamic phenomena, such as cavitation and coalescence, the time scale of dissolution is relatively longer. Therefore, on very short time scales, the effect of nanobubble dissolution, which causes a time-dependent bubble radius, may be neglected for simplicity, as demonstrated in previous studies.
To unveil the mechanism of nanobubble motions driven by thermal Marangoni flows, here the movement of a gas nanobubble between two solid plates with a temperature gradient is investigated by well-controlled molecular dynamics simulations. The previous thermal Marangoni flow theory for macroscale gas bubbles is extended by considering potentially enhancing nanoscopic interfacial effects such as the thermal resistance and the slip length across the liquid–gas interface. Expressions of the Marangoni force and the drag force are derived. By balancing the Marangoni force and the drag force, the theoretical velocity of bubble migration is determined. Finally, the simulation results are compared with the theory, finding good agreement.
This paper is organized as follows. In § 2, the procedure of molecular simulations is described. Section 3 shows the bubble trajectories from simulations. This section also includes a systematic measurement of the transport properties of water, gas, and their interfaces. Section 4 derives a model for nanobubble displacement. The comparison between simulation results and theories is then made. We conclude our paper in § 5 with a summary and an outlook.
2. Molecular simulations of nanobubble motion in a temperature gradient
(a) A snapshot of a gas nanobubble in water between two solid plates in MD simulations. The top plate has a lower temperature
$T_l$
, and the bottom plate has a higher temperature
$T_h$
. (b) Sketch of a gas nanobubble in water between two solid plates. Here,
$\gamma$
is the local surface tension of the bubble surface, and the bubble portion with a lower (higher) temperature has a higher (lower) surface tension, denoted as
$\gamma +$
(
$\gamma -$
). This creates a Marangoni flow around the bubble, resulting in a Marangoni force
$F_{Ma}$
on the bubble.
Molecular dynamics (MD) simulations are used to simulate the motion of bulk gas nanobubbles in liquid between two solid plates whose temperatures are different. We adopt the open-source code LAMMPS (Plimpton Reference Plimpton1995). As shown in figure 1(a), the molecular system consists of water molecules (represented in orange), gas atoms (represented in green), atoms of the bottom solid (represented in red) with a higher temperature
$T_h$
, and atoms of the top solid (represented in blue) with a lower temperature
$T_l$
.
The mW water potential is adopted to model water. The mW water model is a monatomic water model proposed by Molinero & Moore (Reference Molinero and Moore2009). It uses the Stillinger–Weber potential:
\begin{align} E&=\sum \limits _{i}{\sum \limits _{j\gt i}{{{\phi }_{2}}\left ( {{r}_{ij}} \right )}}+\sum \limits _{i}{\sum \limits _{j\ne i}{\sum \limits _{k\gt j}{{{\phi }_{3}}\left ( {{r}_{ij}},{{r}_{ik}},{{\theta }_{ijk}} \right )}}}, \nonumber \\ {{\phi }_{2}}\left ( r \right )&=A\varepsilon \left [ B{{\left ( \frac {\sigma }{r} \right )}^{p}}-{{\left ( \frac {\sigma }{r} \right )}^{q}} \right ]\exp \left ( \frac {\sigma }{r-a\sigma } \right ), \nonumber \\ {{\phi }_{3}}\left ( r,s,\theta \right )&=\lambda \varepsilon {{\left ( \cos \theta -\cos {{\theta }_{0}} \right )}^{2}}\exp \left ( \frac {\gamma \sigma }{r-a\sigma } \right )\exp \left ( \frac {\gamma \sigma }{s-a\sigma } \right ). \end{align}
Here,
$A=7.049556277$
,
$B=0.6022245584$
,
$\varepsilon =6.189$
kcal mol
$^{-1}$
,
$\sigma =2.3925 {{\unicode{x00C5}}}$
,
$a=1.80$
,
$\lambda =23.15$
,
$\gamma =1.20$
,
$\cos \theta _0= -1/3$
,
$p=4$
and
$q=0$
. Compared to other water models such as TIP4P/2005 and SPC/E, the mW water model has a good accuracy of surface tension and a low computational cost (Molinero & Moore Reference Molinero and Moore2009). For example, a large time step can be chosen for the mW water model (Molinero & Moore Reference Molinero and Moore2009). In the current simulations, it is 8 fs.
Except for water itself, the intermolecular potentials
$U$
between
$i$
-type atoms and
$j$
-type atoms are simulated with the standard Lennard-Jones 12–6 potential:
\begin{equation} U({{r}_{ij}}) = \begin{cases} 4\varepsilon _{ij} \left [ {{\left ( \dfrac {\sigma _{ij} }{{{r}_{ij}}} \right )}^{12}}-{{\left ( \dfrac {\sigma _{ij} }{{{r}_{ij}}} \right )}^{6}} \right ] & \text {if} \,\,\,{{r}_{ij}}\leqslant {{r}_{c,ij}},\\ 0 & \text {if}\,\,\, {{r}_{ij}}\gt {{r}_{c,ij}}.\\ \end{cases} \end{equation}
Here,
$r_{ij}$
,
$\varepsilon _{ij}$
,
$\sigma _{ij}$
and
${r}_{c,ij}$
are the pairwise distance, energy parameter, length parameter and cut-off distance, respectively. The complete list of parameters among water (W), gas (G) and solid (S) is given in table 1. The gas modelled by the standard 12–6 Lennard-Jones potential represents nitrogen, and it has density
$\rho _{\infty }=1.147$
kg m−3 at 1 atm and 300 K. The molar mass of the gas is 28 g mol−1. As noted earlier, clean bulk nanobubbles in saturated or undersaturated environments tend to dissolve in microseconds, as predicted by the Epstein–Plesset theory (Epstein & Plesset Reference Epstein and Plesset1950; Lohse & Zhang Reference Lohse and Zhang2015). We indeed find that in our molecular simulations, a nanobubble dissolves over time, leading to a time-dependent bubble radius. This obviously complicates the theoretical modelling and goes beyond the scope of current work to validate the YGB59 theory (where a constant bubble radius is assumed) to the nanoscale. To circumvent this problem, one can tune the energy parameter, distance parameter and cut-off distance between gas and water to obtain a small gas solubility. In current simulations, the cut-off distance between water and gas is set to be
$2^{1/6}\sigma _{GW}=3.4$
nm so that only the repulsive force between gas and water is effective, leading to an extremely small gas solubility. The same simulation strategy was used when simulating nanobubble cavitation (Shekhar et al. Reference Shekhar, Nomura, Kalia, Nakano and Vashishta2013).
Interaction parameters among water (W), gas (G) and solid (S).
The box has lateral size
$L_x=19.2$
nm and
$L_y=19.2$
nm; see figure 2(b). The height
$L_z$
of this box is fixed during the bubble movement, but it varies for different cases with different bubble radii, as discussed below and shown in table 2. The height of each solid plate at the lower and upper ends of the domain is
$L_p=0.96$
nm. Both the top and bottom plates are made of five layers of atoms with the same density as water for simplicity. Periodic boundary conditions are applied in the lateral directions. On the temperature conditions, the bottom hot plate is held at 350 K or 325 K, and the top cold plate is held at 300 K using the Langevin thermostat with damping factor 100 fs. A constant spring force 10 kcal mol−1
${{\unicode{x00C5}}}$
−2 is applied to each atom in the plates to tether them to their initial position. The positions and velocities of the fluid atoms (liquid and gas) are updated based on Newton’s second law without extra thermostating.
Simulated three cases with different settings.
(a) Snapshots of the bubble movement in the MD simulations at three different times for case 1. (b) For the same simulation, the three-dimensional trace of the bubble motion is shown.
The tank is initially filled with 467 200 water atoms and is isothermal with temperature 300 K. Then
$N$
water atoms next to each other centred at
$h_0=28.8$
nm are switched into gas atoms, which forms a nanobubble. Notably, the nanobubble is initially tethered to its initial position through a spring force so that it can only expand or shrink without Brownian motions. During this process, the top plate can move vertically to maintain the ambient liquid pressure at 1 atm in response to the appearance of the nanobubble. Two different situations with
$N=1196$
and
$N=3316$
are simulated. The corresponding radii of the obtained nanobubbles at equilibrium are measured to be 3.1 nm and 5.0 nm, consistent with the prediction of the real gas law (Zhang & Lohse Reference Zhang and Lohse2023). At equilibrium, the position of the top plate barely moves, and
$L_z=38.4$
nm for
$R=3.1$
nm, and
$L_z=39$
nm for
$R=5.0$
nm, are measured and listed in table 2. After obtaining the initial configuration, the Langevin thermostat is applied to the two plates. The formed nanobubbles are still tethered to their initial positions for another 100 ns, during which a steady temperature gradient is produced between the two solid plates. After 100 ns, the nanobubbles are released so that they can move freely in the liquid. Due to the thermal motions of nanobubbles, for each case, several realizations (see table 2) have to be run to improve the statistics.
3. Results from the MD simulations
3.1. Nanobubble displacement
Figure 2(a) shows snapshots of the moving bubble (projected at the
$x{-}z$
plane) obtained from the molecular simulations at three different times for case 1. Just as for the case of macroscopic bubbles driven by a temperature difference, the nanobubble moves to the hot plate due to the thermal Marangoni force. Figure 2(b) plots the three-dimensional trajectory of the moving bubble. It can be seen that though on average the nanobubble moves downwards towards the hot plate, it experiences considerable thermal motions in all three directions during the movement. Thus, as shown in figures 3(a–c), a number of realizations have been performed to get the mean bubble trajectory, from ensemble averaging. For example, for case 1, five realizations represented by grey lines in Figure 3(a) have been conducted. The red line represents the mean value of the five realizations. The error bars stand for the standard deviations of the five samples. For case 2, the performed numbers of independent realizations are 9, as smaller bubbles are more vulnerable to thermal fluctuations, evidenced by the larger standard deviations of the samples, as shown in figure 3(b). For case 3 in figure 3(c), five realizations were performed. The obtained mean bubble distance as a function of time for the three cases is shown in figure 3(d) for comparison. The dependence of bubble distances on time is linear when the bubbles are away from the plates. This means that the Marangoni force that drives the bubble motion is balanced by the drag force. Thus thermocapillarity works down to the nanoscale. By a linear fitting of the data from 10 ns to 20 ns, the corresponding bubble velocity
$U_{MD}$
is obtained, as shown in table 3. To quantitatively compare the simulation results with the YGB59 theory (Young et al. Reference Young, Goldstein and Block1959), one has to systematically measure the transport properties of the simulated liquid and gas.
Comparisons among the bubble velocity results of MD simulations (
$U_{MD}$
), the theoretical velocity ignoring the gas thermal conductivity and the viscosity (
$U_{Macro}$
), and the revisited theoretical velocity
$U_{Nano-NS}$
(
$U_{Nano-FS}$
) considering nanoscale effects and without (with) the slip length. Here,
$\Delta T_{eff}$
is the effective temperature difference after considering the solid–liquid thermal resistance.
Trajectories of the moving bubbles in the molecular simulations for the three different cases. Grey lines represent different realizations in each case. The solid coloured lines are the mean trajectories of all the realizations, with the error bars representing the standard deviations: (a) case 1, bubble radius
$R=5.0$
nm and
$\Delta T=50$
K; (b) case 2, bubble radius
$R=3.1$
nm and
$\Delta T=50$
K; (c) case 3, bubble radius
$R=5.0$
nm and
$\Delta T=25$
K; (d) mean bubble trajectory for the mentioned three cases, for comparisons.
3.2. Transport properties of the simulated water and gas
Though the mW water model (Molinero & Moore Reference Molinero and Moore2009) adopted here is a standard model, to our knowledge, its transport properties such as surface tension and viscosity as a function of temperature are not available in the literature. A systematic measurement of these mentioned properties is thus performed. The surface tension of water is measured by the standard method using the mechanical definition of surface tension at the atomic level (see Appendix A). The value of surface tension at different temperatures ranging from
$300$
K to
$350$
K is shown in figure 4(a). In this temperature range, water surface tension decreases linearly with temperature, and
${\rm d} \gamma /{\rm d}T=-0.082$
mN m−1 K−1 is obtained by a linear fitting as shown in figure 4(a). The shear viscosity of water is measured by the Green–Kubo method as also explained in Appendix A. We obtain that the dynamical viscosity
$\mu ^o$
(the superscript
$^o$
is for the liquid, and the superscript
$^i$
is for the inner fluid) of the mW water model at 300 K is approximately 0.35 mPa s, which underpredicts the experimental value compared to other water models such as TIP4P/2005 and SPC/E (González & Abascal Reference González and Abascal2010). The dynamic viscosity of water also exhibits a temperature dependence like surface tension. However, the dependence is weak in the temperature range
$300$
–
$350$
K. Therefore we take the dynamic viscosity of water to be the mean dynamic viscosity
$\mu ^o=0.27$
mPa s in this temperature range, for simplicity. In principle, the thermal conductivity of water can also be measured by the Green–Kubo method, but this method is known to converge poorly for thermal conductivity (Sirk, Moore & Brown Reference Sirk, Moore and Brown2013). Thus here the thermal conductivity of water
$k^o$
is measured by the non-equilibrium heat flux method (Wirnsberger, Frenkel & Dellago Reference Wirnsberger, Frenkel and Dellago2015), as shown in Appendix A. Within the simulated temperature range (300–500 K), the thermal conductivity of water hardly varies, thus the mean value
$k^o=0.34$
W m−1 K−1 is used. The mW water model underpredicts the experimental value of water thermal conductivity, while other water models TIP4P/2005 and SPC/E overpredict the experimental value (Sirk et al. Reference Sirk, Moore and Brown2013).
(a) Surface tension of the mW water model as a function of temperature. Symbols are MD results. The dashed lines show the linear fitting. (b) Viscosity (left-hand axis) and thermal conductivity (right-hand axis) for the mW water model. (c) Viscosity (left-hand axis) and thermal conductivity (right-hand axis) of gas as a function of gas density. (d) Liquid–gas slip (left-hand axis) and thermal resistance (right-hand axis) as a function of gas density.
Normally for macroscopic gas bubbles, the heat conduction and the flow inside the bubble are ignored because the heat conductivity and the viscosity of the gas are both small compared to those for water outside. For a dilute gas, based on the kinetic theory, its thermal conductivity and viscosity depend only weakly on temperature, and are independent of gas pressure. Under these conditions, one may ignore the gas and simply take the property of water to calculate the bubble velocity. Based on the YGB59 theory (Young et al. Reference Young, Goldstein and Block1959), i.e. (4.28), the theoretical bubble velocities for the three cases are then calculated to be 0.91 m s−1, 0.56 m s−1 and 0.45 m s−1 (see
$U_{Macro}$
in table 3), which, however, considerably overpredict the simulated bubble velocity
$U_{MD}$
.
There are multiple factors that may cause the mismatch between the MD results and the YGB59 theory. For nanobubbles whose pressure can be as high as dozens of MPa, the predictions of the gas viscosity and thermal conductivity, based on the kinetic theory, may be wrong. Indeed, we find that the thermal conductivity and viscosity of gas strongly depend on the gas density as shown in figure 4(c), measured by the method discussed below. Another reason may be that the temperature and velocity continuity condition at the bubble surface breaks down.
3.3. Thermal resistance and slip length across the liquid–gas interface
(a) Snapshot of the liquid–gas system in MD simulations to measure the liquid–gas thermal resistance. (b) The temperature profile for gas density
$0.0036\ {{{\unicode{x00C5}}}}^3$
. (c) The temperature profile for gas density
$0.0063\ {{{\unicode{x00C5}}}}^3$
. (d) The temperature profile for gas density
$0.0096\ {{{\unicode{x00C5}}}}^3$
. The symbols are MD results, and the red lines are linear fittings.
Many researchers have shown the existence of a temperature jump across a liquid–solid interface (Barrat & Chiaruttini Reference Barrat and Chiaruttini2003; Alosious et al. Reference Alosious, Kannam, Sathian and Todd2019; Hadjiconstantinou & Swisher Reference Hadjiconstantinou and Swisher2022), a liquid–liquid interface (Patel et al. Reference Patel, Garde and Keblinski2005) and a liquid–gas interface (Patel et al. Reference Patel, Garde and Keblinski2005; Plascencia et al. Reference Plascencia, Bird and Liang2022; Dockar et al. Reference Dockar, Gibelli and Borg2023). A velocity jump may also exist at a liquid–liquid interface (Poesio et al. Reference Poesio, Damone and Matar2017; Telari et al. Reference Telari, Tinti and Giacomello2022; Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023).
To check whether there is thermal resistance between the water–gas interface and the water–solid interface, we use the setting as shown in figure 5(a), where both gas (in green) and water (in dark yellow) are confined between two plates held at different temperatures. While the number of water atoms is fixed in the box, the number of gas atoms is varied to have different gas densities. After equilibrium, the temperature profile is obtained. For a low gas number density
$\rho =0.0036\ {{{\unicode{x00C5}}}}^3$
, as shown in figure 5(b), we observe noticeable temperature jumps at the liquid–solid, liquid–gas and gas–solid interfaces. As the gas–solid interface is not relevant for the Marangoni-driven bubble motion, only the thermal resistance for the liquid–solid interface and the liquid–gas interface are calculated. To obtain the temperature jump, we interpolate the linear fitting of the temperature profile into the interfaces: a temperature jump 3.4 K is obtained for the liquid–solid interface (at
$z=0$
nm), and 10.0 K is obtained for the liquid–gas interface (at
$z=6.2$
nm). By its definition (4.7), the thermal resistance of the liquid–solid interface is calculated to be
$1.30\times 10^{-8}$
m2K W−1 , and for the liquid–gas interface, it is
$4.16\times 10^{-8}$
m2 K W−1. The order of obtained thermal resistance is consistent with results reported in the literature (Plascencia et al. Reference Plascencia, Bird and Liang2022). Using the ratio of the slopes of temperature profiles in the liquid and in the gas, the thermal conductivity of the gas of density
$\rho =0.0036\ {{{\unicode{x00C5}}}}^3$
is 0.028 W m−1 K−1.
(a) Snapshot of the Couette flow to measure the gas–liquid slip. (b) The velocity profile for gas density
$0.0036{\unicode{x00C5}}^3$
. (c) The velocity profile for gas density
$0.0063{\unicode{x00C5}}^3$
. (d) The velocity profile for gas density
$0.0096{\unicode{x00C5}}^3$
. The symbols are MD results, and the red lines are linear fittings.
For the gas density the same as case 1 (
$\rho =0.0063 {{\unicode{x00C5}}}^3$
) and case 2 (
$\rho =0.0096{{{\unicode{x00C5}}}}^3$
), the temperature profiles in figures 5(c) and 5(d) show smaller temperature jumps
$7.7$
K and
$5.3$
K at the liquid–gas interface. The thermal resistances of the liquid–gas are calculated to be
$2.0\times 10^{-8}$
m2 K W−1 and
$1.1\times 10^{-8}$
m2K W−1. The thermal resistance is thus reduced by increasing the gas density, which is similar to the theoretical prediction of the solid–gas thermal resistance (Chen et al. Reference Chen, Xu, Zhou and Li2022) and molecular simulations of liquid–vapour thermal resistance (Holyst & Litniewski Reference Holyst and Litniewski2008). The thermal conductivity of the gas with a higher density is, however, increased. For
$\rho =0.0063{{{\unicode{x00C5}}}}^3$
, it is 0.0466 W m−1 K−1 , and for
$\rho =0.0096{{{\unicode{x00C5}}}}^3$
, it is 0.0593 W m−1 K−1. The relation between gas thermal conductivity or liquid–gas thermal resistance and gas density is shown in figures 4(c) and 4(d). One can also use the heat flux method (Wirnsberger et al. Reference Wirnsberger, Frenkel and Dellago2015) to measure the gas thermal conductivity as a function of gas density, and we have checked that both methods produce the same results. We also note that the temperature jumps across the liquid–solid interface indeed increase with an increasing heat flux, as shown by figures 5(b–d), and the mean liquid–solid thermal resistance of the three cases is
$1.22\times 10^{-8}$
m2 K W−1, which is used later to estimate the effective temperature difference that drives nanobubble motions (see
$T_{eff}$
in table 3).
In terms of the slip length between liquid and gas, we use a similar configuration, as shown in figure 6(a). Here, the left-hand plate can move in the
$x$
direction with velocity
$5\times 10^{-4}$
m s-1 , while the right-hand plate is fixed. The velocity profiles of this Couette flow are shown in figures 6(b–d) for three different gas densities. One can see that there are indeed velocity jumps at the liquid–gas interface. The velocity jump decreases with increasing gas density. By its definition (B1) shown in Appendix B, the slip length for the gas (liquid) phase
$b^i$
(
$b^o$
) is measured to be 1.1 nm (1.3 nm), 0.61 nm (0.3 nm), and 0.52 nm (0.15 nm), respectively for the three different gas densities. So the slip lengths
$b^i$
and
$b^o$
decrease with increasing gas densities, as shown in figure 4(d). The gas viscosity can be obtained by the ratio of the slopes of the red fitting lines based on the continuity of the shear stress. The obtained gas viscosities are 0.0229 mPa s, 0.0270 mPa s and 0.0326 mPa s, respectively, for the three different gas densities, thus the gas viscosities increase with increasing gas densities. This relation is also plotted in figure 4(c).
After all, we discovered that the thermal conductivity and viscosity of gas nanobubbles have been enhanced a lot by the large gas density on the nanoscale. It was also found that temperature jumps and velocity jumps coexist at the bubble surface, which invalidates the assumption of temperature and velocity continuities at the bubble surface adopted by the classic YGB59 theory (Young et al. Reference Young, Goldstein and Block1959). Thus the YGB59 theory has to be generalized. Notably, though we have obtained thermal resistance and slip length from our simulations whose values are of the same order as reported values in the literature, we have used an ideal potential to describe the water–gas intermolecular interactions compared to other more realistic systems (Plascencia et al. Reference Plascencia, Bird and Liang2022; Dockar et al. Reference Dockar, Gibelli and Borg2023; Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023). However, the existence of thermal resistance and slip at the liquid–gas interface is generally true.
4. Extension of the YGB59 theory towards the nanoscale
Though Young et al. (Reference Young, Goldstein and Block1959) (YGB59 theory) derived the correct expression of bubble velocity, their work did not include the expression of Marangoni forces and drag forces. This has led to some approximate expressions in the literature (Morick & Woermann Reference Morick and Woermann1993; Lubetkin Reference Lubetkin2003; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021), which are valid only at the macroscale. Here, the YGB59 equation is revisited, and more general expressions for Marangoni forces and drag forces are obtained, which show important nanoscale effects, such as those of the thermal resistance and the slip at the interface.
The governing dimensionless numbers for both fluids (liquid and gas) are the Reynolds number
$Re=\rho u_0 R/\mu$
and the Marangoni number
$Ma=({{\rm d}\gamma }/{{\rm d}T})({{\rm d}T}/{{\rm d}z})({{R}^{2}}/ ( \mu \alpha ))$
. Here,
$\rho$
,
$\mu$
,
$u_0$
,
$R$
,
$\gamma$
and
$\alpha$
are the density, dynamic viscosity, velocity, radius of the nanobubble, surface tension and thermal diffusivity, respectively. Assuming that
$Re$
and
$Ma$
are much smaller than unity for both fluids, which is valid at the small scale (small
$R$
) of nanobubbles, the momentum equation and continuity equation for both fluids are
\begin{equation} \mu\, {\boldsymbol {\nabla }^{2}} {\boldsymbol u}=\nabla p, \end{equation}
\begin{equation} \boldsymbol {\nabla } \cdot {\boldsymbol u}=0, \end{equation}
respectively. In the momentum equation, the unsteady term
$\rho\, \partial u / \partial t$
is dropped. This is justified because the ratio of the time scale required for the water velocity to adjust to the changes in the bubble position (
$\rho R^2/ \mu$
) to the time scale of bubble movement (
$R/u_0$
), i.e.
$Re$
, is small. In fact, the largest Reynolds number in our simulations is evaluated to be 0.01 for case 1. Therefore, the water velocity field may be assumed to be quasi-static. The heat transfer equation for both fluids is
\begin{equation} {{\nabla }^{2}}T=0. \end{equation}
The general solution to the temperature field for both fluids in spherical coordinates is
\begin{equation} T={{A}_{1}}+\frac {{{A}_{2}}}{r}+{{A}_{3}}r\cos \theta +\frac {{{A}_{4}}}{{{r}^{2}}}\cos \theta . \end{equation}
Here,
$A_i$
(
$i=1,2,3,4$
) is the coefficient to be determined by the boundary conditions. At the far field, the temperature distribution
$T^o$
for the outer fluid (i.e. water here – the superscript
$^o$
is for the liquid, and the superscript
$^i$
is for the inner fluid, i.e. gas hereafter) is
\begin{equation} T^o=T_0+\Gamma r\cos {\theta }, \end{equation}
where
$T_0$
is the temperature at the centre of the bubble, and
$\Gamma$
is the temperature gradient. At the fluid–fluid interface, the continuity of the heat flux requires
\begin{equation} k^o\frac {\partial T^o}{\partial r}= k^i\frac {\partial T^i}{\partial r}, \end{equation}
where we have assumed a spherical bubble with the coordinate system at its centre. For the temperature condition across the interface, motivated by the observation from our simulations, we adopt the temperature jump condition across the interface instead of the continuity of temperature, so that one can have
\begin{equation} T^o-T^i=k^o\frac {\partial T^o}{\partial r}G, \end{equation}
where
$G$
is the thermal resistance.
Using the above temperature boundary conditions, the temperature distribution can be obtained for both fluids:
\begin{align} & {{T}^{o}}={{T}_{0}}+{\Gamma }r\cos \theta +\frac { ( {{k}^{o}}-{{k}^{i}} )R+{{k}^{o}}{{k}^{i}}G}{( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {{\Gamma }{{R}^{3}}\cos \theta }{{{r}^{2}}}, \quad r\geqslant R,\nonumber \\ & {{T}^{i}}={{T}_{0}}+\frac {3{{k}^{o}}R}{ ( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o}}{{k}^{i}}G}{\Gamma }r\cos \theta, \quad r \lt R. \end{align}
In terms of the velocity field, the usual way to solve (4.1) and (4.2) is by introducing the Stokes stream function
$\psi (r,\theta )$
(Bird, Stewai & Lightfoot Reference Bird, Stewai and Lightfoot2006)
\begin{align} & {{u}_{r}}=-\frac {1}{{{r}^{2}}\sin \theta }\frac {\partial \psi }{\partial \theta },\nonumber \\ & {{u}_{\theta }}=\frac {1}{r\sin \theta }\frac {\partial \psi }{\partial r}. \end{align}
Equation (4.1) is then rewritten in terms of the Stokes stream function
$\psi$
as
\begin{equation} {{E}^{4}}\psi =0, \end{equation}
where
${{E}^{2}}=({\partial }/{\partial {{r}^{2}}})+({\sin \theta }/{{{r}^{2}}})({\partial }/{\partial \theta }) ( ({1}/{\sin \theta })({\partial }/{\partial \theta }) )$
. The general solution to (4.10) for both fluids is
\begin{align} & {{\psi }^{o}}=\left ( {{C}_{1}}{{r}^{-1}}+{{C}_{2}}r+{{C}_{3}}{{r}^{2}}+{{C}_{4}}{{r}^{4}} \right ){{\sin }^{2}}\theta, \quad r\geqslant R, \end{align}
\begin{align} & {{\psi }^{i}}=\left ( {{D}_{1}}{{r}^{-1}}+{{D}_{2}}r+{{D}_{3}}{{r}^{2}}+{{D}_{4}}{{r}^{4}} \right ){{\sin }^{2}}\theta,\quad r\lt R. \end{align}
The resulting expressions for the velocities are
\begin{align} & u_{r}^{o}=-2\cos \theta \left ( \frac {{{C}_{1}}}{{{r}^{3}}}+\frac {{{C}_{2}}}{r}+{{C}_{3}}+{{C}_{4}}{{r}^{2}} \right ), \end{align}
\begin{align} & u_{\theta }^{o}=\sin \theta \left ( -\frac {{C}_{1}}{{r}^{3}}+\frac {{{C}_{2}}}{r}+2{{C}_{3}}+4{{C}_{4}}{{r}^{2}} \right ), \end{align}
\begin{align} & u_{r}^{i}=-2\cos \theta \left ( \frac {{{D}_{1}}}{{{r}^{3}}}+\frac {{{D}_{2}}}{r}+{{D}_{3}}+{{D}_{4}}{{r}^{2}} \right ), \end{align}
\begin{align} & u_{\theta }^{i}=\sin \theta \left ( -\frac {{D}_{1}}{{r}^{3}}+\frac {{{D}_{2}}}{r}+2{{D}_{3}}+4{{D}_{4}}{{r}^{2}} \right ) . \end{align}
These factors
$C$
and
$D$
are determined by the velocity boundary conditions. We assume that the bubble has velocity
$U\equiv -u_\infty$
(
$u_\infty \gt 0$
). Taking the reference frame with the bubble, the far-field liquid velocity will be
$u_r=u_\infty \cos {\theta }$
and
$u_\theta =-u_\infty \sin {\theta }$
. Thus the far-field stream function will be
\begin{equation} \psi ^{o}=-\tfrac {1}{2}u_\infty r^2 \sin ^2 \theta, \quad r\rightarrow \infty, \end{equation}
so
$C_4=0$
and
$C_3=- ({1}/{2})u_\infty$
using (4.11). For the flow inside the nanobubble, the velocity at
$r=0$
has to be finite, so that
$D_1=0$
and
$D_2=0$
using (4.12). Furthermore, there are no masses across the bubble surface, and no bubble shape changes, so that
$u_{r}^{o}=0$
and
$u_{r}^{i}=0$
at
$r=R$
. This leads to
$C_1= (u_{\infty }R-2C_2 )R^2/2$
and
$D_3=-R^2D_4$
.
For boundary conditions at the liquid–gas interface, the shear stress continuity is obeyed (Leal Reference Leal2007):
\begin{equation} \tau _{r\theta }^{o}-\tau _{r\theta }^{i}=-{{\nabla }_{s}}\gamma. \end{equation}
Here, one can obtain
$-{{\nabla }_{s}}\gamma =-({1}/{R})({\partial \gamma }/{\partial \theta })=- ({1}/{R}) ({\partial \gamma }/{\partial T}) ({\partial T}/{\partial \theta })= ({\partial \gamma }/{\partial T})$${(3{{k}^{o}}R+3k^ok^iG)}/{ (( 2{{k}^{o}}+{{k}^{i}} )R+2{{k}^{o})}{{k}^{i}}G}{\Gamma }\sin \theta$
, knowing the temperature distribution by (4.8). Later, we define
$A=({\partial \gamma }/{\partial T}) {(3{{k}^{o}}R+3k^ok^iG)}/{ (( 2{{k}^{o}}+{{k}^{i}} )R+ 2{{k}^{o}}{{k}^{i}}G}{\Gamma )}$
such that we can simply write
$-{{\nabla }_{s}}\gamma =A \sin \theta$
.
A number of recent works have shown that there exists slip at a fluid–fluid interface (Poesio et al. Reference Poesio, Damone and Matar2017; Telari et al. Reference Telari, Tinti and Giacomello2022; Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023). As seen from the simulation of the two-phase Couette flow in figure 6, this is indeed the case for the liquid–gas interface in our system. Thus it is natural to include slip in the theoretical model for the bubble velocity. However, how to model the fluid–fluid slip condition is not a well-explored topic. Even for the liquid–solid slip, the frequently used Navier’s slip condition (Zhang Reference Zhang2024), which assumes that the slip velocity is proportional to shear stress at the wall, is empirical (Leal Reference Leal2007). Thus at first, we still use the classic no-slip condition to derive the theoretical bubble velocity. Later, we then try to use empirical fluid–fluid slip models to see any improvements in the predictions.
The no-slip condition at the fluid–fluid interface is
\begin{equation} u_{\theta }^o=u_{\theta }^i. \end{equation}
Using the shear stress condition (4.18) and the no-slip condition (4.19), one can obtain
\begin{equation} {{D}_{4}}=\frac {4{{C}_{2}}-3{{u}_{\infty }}R}{4{{R}^{3}}}, \end{equation}
\begin{equation} {{C}_{2}}=\frac {6{{\mu }^{o}}{{u}_{\infty }}{R}^{2}+9{{\mu }^{i}}{{u}_{\infty }}{{R}^{2}}-2A {{R}^{3}}}{12{{\mu }^{o}}R+12{{\mu }^{i}}R}. \end{equation}
Now the Marangoni force exerted on the bubble can be evaluated by integrating the shear stress:
\begin{align} {{F}_{Ma}}&=-2\pi {{R}^{2}}{{\int _{0}^{\pi }{\tau _{r\theta } }}}\sin^{2}\theta\, {\rm d}\theta =-8\pi {{\mu }^{o}}\left ({{u}_{\infty }R-2{{C}_{2}}} \right ) \nonumber \\ & =-\frac {4\pi }{3}{{\mu }^{o}}\frac {2A{{R}^{2}}-3{{\mu }^{i}}{{u}_{\infty }}{{R}}}{{{\mu }^{o}}+{{\mu }^{i}}}. \end{align}
Similarly, by integrating the normal stress
$\sigma _{rr}=-p+2\mu ({\partial u_r}/{ \partial r})$
, the pressure drag force experienced by the bubble is
\begin{align} {{F}_{P}} &=2\pi {{R}^{2}}\int _{0}^{\pi }{\sigma _{rr} \sin \theta }\cos \theta\, {\rm d}\theta =8\pi \mu ^o \left ( {{u}_{\infty }}R-{{C}_{2}} \right ) \nonumber \\ & =2\pi \mu ^o \frac {2AR^2+6{{\mu }^{o}}{{u}_{\infty }} R+3{{\mu }^{i}}{{u}_{\infty }}{{R}}}{3{{\mu }^{o}}+3{{\mu }^{i}}}. \end{align}
One can see that the eventual expression of the pressure drag depends on
$C_2$
, which again depends on the boundary condition at the liquid–gas interface. For a clean bubble with a negligible gas viscosity,
$\tau _{r\theta }=0$
so that
$C_2=u_{\infty }R/2$
and
$F_P=4\pi \mu ^o u_{\infty }R$
, which recovers the classic drag on a conventional bubble. For our case,
$C_2$
is different. By letting
$F_{Ma}+F_{P}=0$
, one finds that
$C_2=0$
. Using (C2), the velocity of the bubble (
$U\equiv -u_\infty$
) is obtained as
\begin{equation} U=-\frac {2R}{2{{\mu }^{o}}+3{{\mu }^{i}}}\,\frac {{{k}^{o}}R+k^ok^iG}{\left ( 2{{k}^{o}}+{{k}^{i}} \right )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}
If there is no temperature jump (
$ G=0 $
), then (4.24) reduces to the result of YGB59 theory (Young et al. Reference Young, Goldstein and Block1959):
\begin{equation} U=-\frac {2R}{2{{\mu }^{o}}+3{{\mu }^{i}}}\,\frac {{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}
In this case, the Marangoni force and drag force are
\begin{equation} {{F}_{Ma}}=-{4\pi }{{\mu }^{o}}\frac {{{R}^{2}}\dfrac {\partial \gamma }{\partial T}\dfrac {2{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}{{\Gamma }}-{{\mu }^{i}}{{u}_{\infty }}R}{{{\mu }^{o}}+{{\mu }^{i}}}, \end{equation}
\begin{equation} {{F}_{P}}=2\pi {{\mu }^{o}}\frac {{{R}^{2}}\dfrac {\partial \gamma }{\partial T}\dfrac {2{{k}^{o}}}{2{{k}^{o}}+{{k}^{i}}}{{\Gamma }}+2{{\mu }^{o}}{{u}_{\infty }}R+{{\mu }^{i}}{{u}_{\infty }}R}{{{\mu }^{o}}+{{\mu }^{i}}}. \end{equation}
For macroscale gas bubbles whose viscosity and thermal conductivity are usually ignored, the velocity becomes
\begin{equation} U=-\frac {R}{2{{\mu }^{o}}}\frac {\partial \gamma }{\partial T}{{\Gamma }}, \end{equation}
and the forces are
${{F}_{Ma}}=-4\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$
and
${{F}_{P}}=4\pi {{\mu }^{o}}{{u}_{\infty }}R+2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$
. Thus the drag force is increased by the presence of the surface tension gradient. This extra term
$2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$
was not included in previous works about bubble motions in temperature gradients on the macroscale (Morick & Woermann Reference Morick and Woermann1993; Lubetkin Reference Lubetkin2003; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021). Reassuringly, if one investigates the total force
${{F}_{Ma}}+{{F}_{P}}$
on the bubble, then it is
$4\pi {{\mu }^{o}}{{u}_{\infty }}R-2\pi ({\partial \gamma }/{\partial T}){{\Gamma }}{{R}^{2}}$
, which is the same for the previous study and the current study.
With all the measured transport properties of liquid, gas, and their interfaces, as summarized in table 3, now the theoretical velocity of nanobubbles without slip
$U_{Nano-NS}$
(see (4.24)) can be evaluated. As shown in table 3,
$U_{Nano-NS}$
values for the three cases agree reasonably well with the simulation results
$U_{MD}$
. This highlights the consideration of the enhanced nanoscopic effects such as gas thermal conductivity, gas viscosity and liquid–gas thermal resistance for accurate descriptions of the nanobubble motion in a temperature gradient.
However, will the consideration of the slip improve the prediction even more? To answer this question, we have to use ‘empirical’ slip conditions for the fluid–fluid interface to include slip in the theoretical model of bubble velocity. In analogy to Navier’s slip condition for a liquid–solid interface, the slip condition for a fluid–fluid interface is proposed to be (Hilaire et al. Reference Hilaire, Siboulet, Charton and Dufrêche2023)
\begin{equation} {{u}_{1}}-{u}_{2}=\frac {{{b}_{1}}}{{{\mu }_{1}}}\tau _{1}=\frac {{{b}_{2}}}{{{\mu }_{2}}}\tau _{2 }, \end{equation}
where
${u}_{1}$
,
$b_1$
and
$\tau _1$
(
${u}_{2}$
,
$b_2$
and
$\tau _2$
) are the interfacial velocity, slip length and interfacial shear stress of fluid 1 (fluid 2). For a clean surface where
$\tau _1 = \tau _2$
, (4.29) means
$b_{1}/b_{2}=\mu _{1}/\mu _{2}$
. However, in our case,
$\tau _1 \neq \tau _2$
due to the extra Marangoni stress at the interface. Thus to account for this inconsistency, an extended slip model for the fluid–fluid interface may be (see Appendix B)
\begin{equation} {{u_\theta }^{o}}-{{u_\theta }^{i}}=\frac {{{b}^{o}}}{{{\mu }^{o}}}\tau _{r\theta }^{o}+\frac {{{b}^{i}}}{{{\mu }^{i}}}\tau _{r\theta }^{i}. \end{equation}
Note that the definition of slip length in (4.30) is different from (4.29).
Using this boundary condition, the obtained bubble velocity with finite slip
$U_{Nano-FS}$
is (see Appendix C)
\begin{equation} U=-\frac {2R}{2{{\mu }^{0}}+3{{\mu }^{i}}\dfrac {1+2{{b}^{o}}/R}{1+3{{b}^{i}}/R}}\,\frac {{{k}^{o}}R+k^ok^iG}{\left ( 2{{k}^{o}}+{{k}^{i}} \right )R+2{{k}^{o}}{{k}^{i}}G}\,\frac {\partial \gamma }{\partial T}\,{\Gamma }. \end{equation}
One can see that the slip from the inner phase
$b^i$
may speed up the bubble movement since the inner circulation is reduced by the slip. In contrast, the slip from the outer phase may slow down the bubble since the also reduced outer flow is the reason for bubble movement. However, using the slip lengths
$b^i$
and
$b^o$
measured from independent MD simulations (see table 3), the consideration of slip does not change the theoretical prediction too much in our case (see
$U_{Nano-FS}$
in table 3) because
$b/R$
is small. But a large liquid–liquid slip may also exist in other systems, such as at the interface between different types of polymer melts (Himmel & Wagner Reference Himmel and Wagner2013), so the consideration of slip may then become necessary.
Notably, the small difference between
$U_{Nano}$
and
$U_{MD}$
may be attributed to several reasons. First, the water viscosity shows a weak temperature dependence in the temperature range from 300 K to 350 K. However, we have used the mean viscosity in this temperature range for simplicity. Second, the Brownian motions of bubbles are strong in MD simulations, and only a limited number of realizations have been performed due to the high computational costs. Third, the solid plates may affect the flow field, though we have used the data of bubble displacements away from the plates.
5. Conclusions and outlooks
In summary, motivated by the fact that macroscopic bubbles often originate from nucleated nanobubbles, the nanobubble motion in a temperature gradient is investigated by molecular dynamics simulations and approximate analytical calculations. Unlike macroscopic bubbles, whose thermal conductivity and viscosity are normally ignored compared to those of water, the simulation results show that the gas thermal conductivity and viscosity of nanobubbles are enhanced to a great extent due to the high internal pressure and the high gas density, and have to be considered. The thermal resistance and slip length across the interface are also found to exist on the liquid–gas interface. The thermocapillary theory by Young et al. (Reference Young, Goldstein and Block1959) for macroscale bubbles is thus extended by considering those effects. For the liquid–gas slip condition, an empirical model is used. Expressions for the Marangoni force, the drag force and the bubble migration velocity are derived. The theoretical bubble velocity evaluated using the measured transport properties of liquid, gas, and their interfaces matches well with the results of molecular simulations. This work highlights the need to consider nanoscopic effects for accurate predictions of nanobubble motions driven by thermal Marangoni flows. However, the small slip length in the current liquid–gas system does not have considerable effects on nanobubble motions.
Our findings have important implications for controlling nanobubble motions during water electrolysis and boiling where current models have not taken the nanoscopic effects into consideration. For example, for the nanobubble detachment, one has to estimate the thermal Marangoni force that keeps nanobubbles attached to the solid (Park et al. Reference Park, Liu, Demirkır, van der Heijden, Lohse, Krug and Koper2023; Bashkatov et al. Reference Bashkatov, Babich, Hossain, Yang, Mutschke and Eckert2023). Our results above have proven that nanoscopic effects may change the thermal Marangoni force significantly at the nanoscale. Another example is the shear flow over the substrate, which is often used to detach bubbles (Darband, Aliofkhazraei & Shanmugam Reference Darband, Aliofkhazraei and Shanmugam2019). The enhanced gas viscosity at the nanoscale will certainly change the drag force on surface nanobubbles. Also, for the coalescence of nanobubbles, the enhanced gas viscosity may have some effects on the film drainage between two nanobubbles and then the bridge growth dynamics (Eggers, Sprittles & Snoeijer Reference Eggers, Sprittles and Snoeijer2024). In the future, it is also important to see the dynamics of growing or dissolving nanobubbles under a temperature gradient. Finally, the current slip model is empirical. More research is required to clarify the liquid–liquid slip condition in the future as well.
Acknowledgements
We are grateful for the computational resources provided by Dutch National supercomputer Snellius and EuroHPC supercomputer Discoverer.
Funding
We acknowledge the financial support from the Advanced Research Center Chemical Building Blocks Consortium (ARC CBBC) under project ARC CBBC 2021.038.UT.15.
Declaration of interests
The authors report no conflict of interest.
Appendix A: Measurements of transport properties
(a) Temporal correlations of the stress when measuring the mW water model with temperature 275.6 K. (b) Integrating the temporal correlation. (c) Temperature distribution in the heat flux method.
The surface tension
$\gamma$
is obtained by its mechanic definition. For a planar liquid–vapour interface perpendicular to the
$z$
axis of the molecular simulation,
$\gamma$
is given by
\begin{equation} \gamma =\int _{-\infty }^{\infty }\left (S_{N}-S_{T}\right )\,{\rm d}z. \end{equation}
Here,
$S$
is the atomic definition of the pressure tensor in molecular simulations,
$S_N$
is the normal component of
$S$
, and
$S_T$
is the tangential component. The dynamic viscosity is obtained from the Green–Kubo expression
\begin{equation} \mu =\frac {V}{{{k}_{B}}T}\int _{0}^{\infty }{\left \langle {{S }_{zx}}(t)\,{{S }_{zx}}(0) \right \rangle }\,{\rm d}t, \end{equation}
where
$V$
is the volume of the liquid, and
$S_{zx}$
is the off-diagonal stress. To calculate the viscosity numerically in MD simulations, a cubic box of molecules undergoing thermal motions is simulated with periodic conditions in all three directions. The temporal correlation of the shear stress, which decays to zero with time, is then obtained from the MD simulation (see figure 7
a) and integrated to obtain the viscosity (see figure 7
b).
The heat flux method is used to calculate the thermal conductivity of water, as described by Wirnsberger et al. (Reference Wirnsberger, Frenkel and Dellago2015). A periodic cubic box with length 4.8 nm filled with 3703 water atoms is simulated. An energy current
$F$
$F=1.38\times 10^{-7}$
J s−1 is subtracted and added into the specific region of the system to sustain a temperature gradient as shown in figure 7(c). The thermal conductivity can then be calculated as
$k=F/(2A\,{\rm d}T/{\rm d}z)$
, where
$A$
is surface area, and the factor
$1/2$
is due to the periodicity condition (Wirnsberger et al. Reference Wirnsberger, Frenkel and Dellago2015).
Appendix B: Fluid–fluid slip condition
For a clean surface, the stress continuity of
$\tau _1 = \tau _2$
means
$b_{1}/b_{2}=\mu _{1}/\mu _{2}$
in (4.29). However, in this way, the defined
$b_1$
and
$b_2$
may not be the intrinsic slip of the fluid–fluid interface, since
$\tau _1 \neq \tau _2$
if there is other shear stress such as the Marangoni stress at the interface. Here, an extension of (4.29) is proposed. Unlike (4.29), where the definition of slip velocity is the same for both fluids (i.e.
$u_1-u_2$
), it may be different for two fluids, in fact. For example, consider a simple two-phase Couette flow shown in figure 8. The interfacial velocity is
$u_1=u_2=u_0$
when there is no slip at the fluid–fluid interface, as shown in figure 8(a). When there is a slip, as shown in figure 8(b),
$u_1\gt u_0$
and
$u_2\lt u_0$
. So the slip velocity for fluid 1 is
$u_1-u_0$
, while it is
$u_2-u_0$
for fluid 2. Thus the slip condition may be
\begin{align} &u_1-u_0=\frac {b_1}{\mu _1}\tau _1, \nonumber \\ &u_2-u_0=-\frac {b_2}{\mu _2}\tau _2. \end{align}
The difference between
$u_1$
and
$u_2$
is thus
\begin{equation} u_1-u_2=\frac {b_1}{\mu _1}\tau _1+\frac {b_2}{\mu _2}\tau _2. \end{equation}
For a clean surface,
$\tau _1=\tau _2$
, (B2) becomes
\begin{equation} u_1-u_2=\left (\frac {b_1}{\mu _1}+\frac {b_2}{\mu _2}\right )\tau _1=\left (\frac {b_1}{\mu _1}+\frac {b_2}{\mu _2}\right )\tau _2, \end{equation}
which is effectively the same as (4.29). Note that the value of the slip length depends on how exactly it is defined. The factor in front of the stress will be the same for (4.29) and (B2) when
$\tau _1=\tau _2$
. The extended model may be general as there is no requirement of
$\tau _1=\tau _2$
and
$b_1/\mu _1 =b_2 /\mu _2$
.
To extract the slip length from our independent MD simulations shown in figure 6(a), we need to know
$u_0$
. It is simply given by
${{u}_{0}}= {U}/({1+{{\mu }_{2}}{{h}_{1}}/ ( {{\mu }_{1}}{{h}_{2}} ))}$
, where
$h_1=6.2$
nm and
$h_2=3.4$
nm are the heights of the two fluid domains. For the case of a gas density
$0.0096{{{\unicode{x00C5}}}}^3$
shown in figure 6(d),
$\mu _2/\mu _1=0.12$
so that
$u_0=4.1\times 10^{-4}$
m s−1. By the definition of slip,
$b^o=0.15$
nm and
$b^i=0.52$
nm are obtained.
Two-phase Couette flow between two plates: (a) no-slip, (b) finite slip.
Appendix C: Bubble velocity with finite slip
Using the shear stress condition (4.18) and the finite slip condition (4.30), one can obtain
\begin{equation} {{D}_{4}}=\frac {{{C}_{2}}\left ( \dfrac {2}{R}+\dfrac {6{{b}^{o}}}{{{R}^{2}}} \right )-\left ( \dfrac {3}{2}+\dfrac {3{{b}^{o}}}{R} \right ){{u}_{\infty }}}{6{{b}^{i}}R+2{{R}^{2}}}, \end{equation}
\begin{equation} {{C}_{2}}=\frac {6{{\mu }^{0}}{{u}_{\infty }}\left ( 3{{b}^{i}}R+{{R}^{2}} \right )+3{{\mu }^{i}}{{u}_{\infty }}\left ( 3{{R}^{2}}+6{{b}^{o}}R \right )-2A\left ( R+3{{b}^{i}} \right ){{R}^{2}}}{12{{\mu }^{0}}\left ( R+3{{b}^{i}} \right )+12{{\mu }^{i}}\left ( R+3{{b}^{o}} \right )}. \end{equation}
Now the Marangoni force exerted on the bubble can be evaluated by integrating the shear stress:
\begin{align} {{F}_{Ma}}&=-2\pi {{R}^{2}}{{\int _{0}^{\pi }{\tau _{r\theta } }}}\sin^{2}\theta\, {\rm d}\theta \nonumber =-8\pi {{\mu }^{o}}\left ({{u}_{\infty }R-2{{C}_{2}}} \right ) \nonumber \\ & =-\frac {4\pi }{3}{{\mu }^{o}}\frac {2A\left ( R+3{{b}^{i}} \right ){{R}^{2}}-3{{\mu }^{i}}{{u}_{\infty }}{{R}^{2}}}{{{\mu }^{0}}\left ( R+3{{b}^{i}} \right )+{{\mu }^{i}}\left ( R+3{{b}^{o}} \right )}. \end{align}
The pressure drag force experienced by the bubble is
\begin{align} {{F}_{P}} &=2\pi {{R}^{2}}\int _{0}^{\pi }{\sigma _{rr} \sin \theta }\cos \theta\, {\rm d}\theta \nonumber =8\pi \mu ^o \left ( {{u}_{\infty }}R-{{C}_{2}} \right ) \nonumber \\ & =2\pi {{\mu }^{o}}\frac {2A\left ( R+3{{b}^{i}} \right ){{R}^{2}}+6{{\mu }^{o}}{{u}_{\infty }}\left ( {{R}^{2}}+3{{b}^{i}}R \right )+3{{\mu }^{i}}{{u}_{\infty }}\left ( {{R}^{2}}+6{{b}^{o}}R \right )}{3{{\mu }^{0}}\left ( R+3{{b}^{i}} \right )+3{{\mu }^{i}}\left ( R+3{{b}^{o}} \right )}. \end{align}
The condition
$F_{Ma}+F_{P}=0$
leads to
$C_2=0$
and thus the expression of bubble velocity with finite slip (4.31).
Open access
