2.1 Apparatus and experimental procedure

All experiments were conducted using a custom-built dynamic fluid-film interferometer (DFI). The assembly of the DFI is described in Frostad et al. (Reference Frostad, Tammaro, Santollani, Bochner de Araujo and Fuller2016). We give a brief description of the apparatus and modifications made to the equipment. Figure 1 shows a schematic of the DFI. A chamber 2 cm wide by 2.3 cm long by 1.8 cm deep is mounted on a platform that is connected to a motorized actuator (Newport TRA12PPD). Two of the chamber sidewalls are made of glass for the visualization of the needle, the bubble and the position of the liquid–air interface. The bubble is supported by a 16-gauge capillary with an inner diameter of 1.2 mm. The capillary is connected to a $100~\unicode[STIX]{x03BC}\text{l}$ syringe (Hamilton 1710CX) that is controlled by a syringe pump (Harvard Apparatus Pump 11 Elite HA1100W). A control valve is added between the capillary and the syringe pump to reduce the dead volume and to minimize dead volume fluctuations due to air compressibility.

Figure 1.

Schematic of the dynamic fluid-film interferometer (DFI).

At the start of a set of experiments, the chamber is mounted onto the platform with the top of the capillary positioned at the centre of the chamber. An amount of 4 ml of silicone oil is then added to the chamber to submerge the capillary. The top of the chamber is then covered with a glass slide to prevent evaporation and the light source (CCS Inc. LAV-80SW2) is turned on. The system is allowed to equilibrate for an hour. During this time, the chamber and the liquid are warmed due to the top light and the temperature of the chamber and the liquid equilibrates to $25\,^{\circ }\text{C}$. A side convection glass wall (5 cm tall) is mounted to the motor platform. The convection wall is added to reduce ambient disturbances on the sides; however, there are still gaps to allow evaporation to take place.

After the temperature has stabilized, a bubble is generated by pumping out $4~\unicode[STIX]{x03BC}\text{l}$ of air. The bubble radius typically ranges from 0.8 to 1 mm and due to the small radii, the initial bubble shapes are near-spherical. Once the bubble is formed, the control valve is closed. Using videos captured by the side camera (ThorLabs DCU223C), the motor stage is positioned to achieve a small initial clearance between the top of the bubble and the free liquid–air interface. The initial clearance in the experiments ranges between 5 % and 15 % of the bubble radius. In post-processing, the initial clearance and the bubble radius are determined by analysing the captured video from the side camera. Variations in the initial positioning of the bubble are accounted for by conducting simulations with average values of the bubble radius and the initial clearance.

Once in position, the system is allowed to equilibrate for another minute. At the start of an experiment, the top glass slide is removed and the motor stage is set to move downward at a speed of $0.075~\text{mm}~\text{s}^{-1}$ for three seconds, after which the stage is held fixed. During this step, the two interfaces move towards each other and at the final motor position, the top of the bubble penetrates the initially flat free liquid–air interface. The top camera (Imaging Development Systems UI-3080CP) starts recording when the motor starts to move. When the liquid film sandwiched between the air–liquid interfaces is less than $5~\unicode[STIX]{x03BC}\text{m}$ in thickness, the interference pattern can be observed through the top camera and its thickness analysed. Image analysis procedures and example interference patterns can be found in Frostad et al. (Reference Frostad, Tammaro, Santollani, Bochner de Araujo and Fuller2016).

2.2 Silicone oil samples

Silicone oils (polydimethylsiloxane, ShinEtsu DM-Fluid) at three different kinematic viscosities were used in this study: 1, 20 and 5000 cSt. The high-viscosity silicone oil was used for model validation and the other two oils were blended to study the solutal–thermal Marangoni flows in a binary mixture of 0.1 vol% of the 1 cSt silicone oil with 99.9 vol% of the 20 cSt silicone oil.

Table 1 shows the material properties of the silicone oils reported by the manufacturer. Kinematic viscosity ($\unicode[STIX]{x1D708}$), density ($\unicode[STIX]{x1D70C}$), surface tension ($\unicode[STIX]{x1D70F}^{\ast }$), thermal conductivity ($\unicode[STIX]{x1D706}$) and heat capacity ($C_{p}$) are used for non-dimensionalization and scaling. The refractive index ($n$) is used for analysing the interference patterns. For the binary mixture interference pattern, the refractive index of 20 cSt silicone oil is used because the mixture is predominantly made up of this material. The binary diffusivity ($D$) of 1 cSt silicone oil in 20 cSt silicone oil is estimated to be $0.0012~\text{mm}^{2}~\text{s}^{-1}$ (Walls, Meiburg & Fuller Reference Walls, Meiburg and Fuller2018). The heat of vaporization ($\unicode[STIX]{x0394}h_{ev}$) of 1 cSt silicone oil is $36.9~\text{kJ}~\text{mol}^{-1}$ (Chickos & Acree Jr Reference Chickos and Acree2003). Thermal dependence of surface tension for the 20 cSt silicone oil is estimated to be $\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}^{\ast }/\unicode[STIX]{x2202}T=-0.08~\text{mN}~\text{m}^{-1}~\text{K}^{-1}$ (Lechner, Wohlfarth & Wohlfarth Reference Lechner, Wohlfarth and Wohlfarth2015).

The evaporation rate of the 1 cSt silicone oil is determined by recording the mass loss of silicone oil over time on a scale with a configuration that closely matches that in the experiments. The net volume flux per unit area for the 1 cSt silicone oil is determined to be $0.075~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$. It is assumed that in a binary mixture, the evaporation rate of the mixture is linearly dependent on the mole fraction of evaporative species. Because $\unicode[STIX]{x1D719}$ represents the volume fraction of the 1 cSt silicone oil in the mixture, a unit conversion between the mole-fraction-based evaporation rate and the volume-fraction-based evaporation was conducted. By accounting for the molecular weight ratio and the density ratio of the two silicone oils, the volume-fraction-based evaporation rate for the 1 cSt silicone oil is determined to be $E^{\ast }=0.546~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$. The 20 cSt silicone oil evaporates on a time scale that is much longer than the duration of an experiment and it is assumed that the 20 cSt silicone oil is non-evaporative.

Table 1.

Selected material properties of the silicone oils used. Values for kinematic viscosity $\unicode[STIX]{x1D708}$, density $\unicode[STIX]{x1D70C}$, surface tension $\unicode[STIX]{x1D70F}^{\ast }$, thermal conductivity $\unicode[STIX]{x1D706}$, heat capacity $C_{p}$, refractive index $n$ and molecular weight $MW$ are provided by the manufacturer (at $25\,^{\circ }\text{C}$).

3.1 Lubrication scaling and governing equations

The small initial clearance atop the bubble allows the lubrication approximation to hold for the problem of interest ($\unicode[STIX]{x1D716}\ll 1$). The lubrication length scales in the vertical and radial directions are $b$ and $\sqrt{ab}$, respectively. The velocity in the vertical direction ($v_{z}^{\ast }$) is scaled by the motor speed $U$ and according to the mass conservation equation, the velocity in the radial direction ($v_{r}^{\ast }$) is scaled by $U/\sqrt{\unicode[STIX]{x1D716}}$. Time is non-dimensionalized by $b/U$. The pressure scale is $(1/\unicode[STIX]{x1D716}^{2})(\unicode[STIX]{x1D707}U/a)$. The temperature deviation from the reference temperature is non-dimensionalized by $(\unicode[STIX]{x1D70C}_{e}/\unicode[STIX]{x1D70C}_{ne})(\unicode[STIX]{x0394}h_{ev}/C_{p})$.

Based on these scalings, the dimensionless governing equation for the total film thickness $h(t,r)$ is

(3.2)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(rh\left<v\right>\right)=-{\mathcal{E}}{\mathcal{v}}\unicode[STIX]{x1D719},\end{eqnarray}$$

where ${\mathcal{E}}{\mathcal{v}}=E^{\ast }/U$ is the evaporation number. The leading-order dimensionless solute species conservation equation and the energy conservation equation are

(3.3)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\left<v\right>\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}-\frac{1}{{\mathcal{P}}{\mathcal{e}}}\frac{1}{h}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}-\frac{1}{{\mathcal{P}}{\mathcal{e}}}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}\right)=-\frac{{\mathcal{E}}{\mathcal{v}}}{h}\unicode[STIX]{x1D719}(1-\unicode[STIX]{x1D719}) & & \displaystyle\end{eqnarray}$$

and

(3.4)$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}t}+\left<v\right>\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}-\frac{1}{{\mathcal{P}}{\mathcal{e}}_{T}}\frac{1}{h}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}-\frac{1}{{\mathcal{P}}{\mathcal{e}}_{T}}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\right)=-\frac{{\mathcal{E}}{\mathcal{v}}}{h}\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}(t,r)$ is the volume fraction of the evaporative species and $\unicode[STIX]{x1D6E9}(t,r)$ is the non-dimensionalized temperature deviation from $T_{ref}$. The mass and thermal Péclet numbers reflect the relative rates of convection to mass and thermal diffusion. They are defined as ${\mathcal{P}}{\mathcal{e}}=aU/D$ and ${\mathcal{P}}{\mathcal{e}}_{T}=aU/\unicode[STIX]{x1D705}$.

The interface deformations are introduced as

(3.5)$$\begin{eqnarray}\displaystyle \bar{h}_{1}=h_{1}+{\mathcal{E}}{\mathcal{v}}\unicode[STIX]{x1D719}_{0}t & & \displaystyle\end{eqnarray}$$

and

(3.6)$$\begin{eqnarray}\displaystyle \bar{h}_{2}=h_{2}+1+\frac{r^{2}}{2}-Ht,\quad \text{where }H=\left\{\begin{array}{@{}l@{}}1\quad (t<t_{stop})\quad \\ 0\quad (t\geqslant t_{stop}),\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

and the modified pressure $\bar{\mathscr{P}}=P+({\mathcal{B}}{\mathcal{o}}/{\mathcal{C}}{\mathcal{a}})z$, where the capillary number is ${\mathcal{C}}{\mathcal{a}}=(1/\unicode[STIX]{x1D716}^{2})(\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D70F}_{ne})$ and the Bond number is ${\mathcal{B}}{\mathcal{o}}=\unicode[STIX]{x1D70C}gab/\unicode[STIX]{x1D70F}_{ne}$. The normal stress balances at $O(1)$ and $O(\unicode[STIX]{x1D716})$ on the two deformed surfaces $\bar{h}_{1}$ and $\bar{h}_{2}$ are

(3.7)$$\begin{eqnarray}\displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\bar{h}_{1}}{\unicode[STIX]{x2202}r}\right)+{\mathcal{C}}{\mathcal{a}}\bar{\mathscr{P}}+2\unicode[STIX]{x1D716}{\mathcal{C}}{\mathcal{a}}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\left<v\right>\right)={\mathcal{B}}{\mathcal{o}}\bar{h}_{1} & & \displaystyle\end{eqnarray}$$

and

(3.8)$$\begin{eqnarray}\displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\bar{h}_{2}}{\unicode[STIX]{x2202}r}\right)-{\mathcal{C}}{\mathcal{a}}\bar{\mathscr{P}}-2\unicode[STIX]{x1D716}{\mathcal{C}}{\mathcal{a}}\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\left<v\right>\right)=0. & & \displaystyle\end{eqnarray}$$

The gravity term is included at the top surface such that in the far field gravity ensures $\bar{h}_{1}$ approaches zero. The lower surface is affected by the pinning of the needle and its far-field shape is matched to that of a pinned bubble in static fluid (see appendix A for details).

Finally to close the set of equations, we introduce the dimensionless surface tension, which is scaled by $\unicode[STIX]{x1D70F}_{ne}$:

(3.9)$$\begin{eqnarray}\unicode[STIX]{x1D70F}=1-\unicode[STIX]{x1D716}{\mathcal{C}}{\mathcal{a}}({\mathcal{M}}{\mathcal{a}}\unicode[STIX]{x1D719}+{\mathcal{M}}{\mathcal{a}}_{T}\unicode[STIX]{x1D6E9}),\end{eqnarray}$$

where the solutal and thermal Marangoni numbers are ${\mathcal{M}}{\mathcal{a}}=-(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D707}U)(\unicode[STIX]{x1D70F}_{e}-\unicode[STIX]{x1D70F}_{ne})|_{T=T_{ref}}$ and ${\mathcal{M}}{\mathcal{a}}_{T}=-(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D707}U)\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ne}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}|_{T=T_{ref}}$. The negative signs in the definitions of the two Marangoni numbers are used to make the values of the appropriate material constants positive. For silicone oils, the more evaporative species has a lower surface tension and, as temperature increases, the surface tension of the silicone oil mixture decreases. The combined tangential stress balance at $O(1)$ and $O(\unicode[STIX]{x1D716})$ is

(3.10)$$\begin{eqnarray}\displaystyle \frac{1}{2}\frac{\unicode[STIX]{x2202}\bar{\mathscr{P}}}{\unicode[STIX]{x2202}r} & = & \displaystyle -\frac{{\mathcal{M}}{\mathcal{a}}}{h}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}-\frac{{\mathcal{M}}{\mathcal{a}}_{T}}{h}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\left<v\right>)\right)+\frac{\unicode[STIX]{x1D716}}{h}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\left<v\right>\right)+\frac{\unicode[STIX]{x2202}\left<v\right>}{\unicode[STIX]{x2202}r}\right).\end{eqnarray}$$

In this problem, we need to consider the combined effect of the $O(1)$ and the $O(\unicode[STIX]{x1D716})$ terms to account for the solutions obtained both in the presence and in the absence of the Marangoni flow. Specifically, when there are no Marangoni effects and evaporation present, the pressure scales as $1/\unicode[STIX]{x1D716}(\unicode[STIX]{x1D707}U/a)$. In the normal stress balance, the pressure for a clean bubble is balanced by the viscous stress associated with the velocity field. Thus, effectively variables $\bar{h}_{1},\bar{h}_{2},\bar{\mathscr{P}}$ and $\left<v\right>$ are composites of their $O(1)$ and $O(\unicode[STIX]{x1D716})$ values, i.e. $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}^{(0)}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D713}^{(1)}$, where $\unicode[STIX]{x1D713}$ stands for the four variables of interest. Variables $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6E9}$ are solved to $O(1)$. More details can be found in appendix A.

The initial and boundary conditions associated with (3.2), (3.3), (3.4), (3.7), (3.8) and (3.10) are

(3.11a-f)$$\begin{eqnarray}\text{at }t=0:\quad \bar{h}_{1}=0,\bar{h}_{2}=0,\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0},\unicode[STIX]{x1D6E9}=0,\bar{\mathscr{P}}=0,\left<v\right>=0;\end{eqnarray}$$

(3.12a-f)$$\begin{eqnarray}\text{at }r=0:\quad \frac{\unicode[STIX]{x2202}\bar{h}_{1}}{\unicode[STIX]{x2202}r}=0,\frac{\unicode[STIX]{x2202}\bar{h}_{2}}{\unicode[STIX]{x2202}r}=0,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}=0,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}=0,\frac{\unicode[STIX]{x2202}\bar{\mathscr{P}}}{\unicode[STIX]{x2202}r}=0,\left<v\right>=0;\end{eqnarray}$$

and

(3.13a-f)$$\begin{eqnarray}\displaystyle & & \displaystyle \text{as }r\rightarrow \infty :\quad \bar{h}_{1}\rightarrow 0,\frac{\unicode[STIX]{x2202}\bar{h}_{2}}{\unicode[STIX]{x2202}r}={\displaystyle \frac{\bar{h}_{2}}{r\left(\log [r]+1-{\displaystyle \frac{1}{\cos [\unicode[STIX]{x1D703}_{b}]}}-\log [2\tan [\unicode[STIX]{x1D703}_{b}/2]]\right)}},\nonumber\\ \displaystyle & & \displaystyle \quad \unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{0},\quad \unicode[STIX]{x1D6E9}\rightarrow 0,\quad \bar{\mathscr{P}}\rightarrow 0\quad \text{and}\quad \frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\left<v\right>\right)\rightarrow 0.\end{eqnarray}$$

The far-field boundary condition on the lower surface $h_{2}$ is derived by matching the shape function in the lubrication region to the bulk region static shape of a pinned bubble (Yiantsios & Davis Reference Yiantsios and Davis1990). Details regarding the derivation of the governing equations and the associated boundary conditions can be found in appendix A.

In the set of six equations, there are ten non-dimensional parameters, summarized in table 2. The capillary number is the ratio of viscous to capillary forces; it reflects the deformability of the interface. The Bond number is the square of the gravity to radial characteristic length scales; it reflects the strength of gravity. The Péclet number compares convection to mass diffusion and the thermal Péclet number compares convection to thermal diffusivity. The solutal Marangoni number is the ratio of surface tension difference to the viscous forces; it reflects the maximum scale of surface tension gradient that can be generated with a given liquid pair. The thermal Marangoni number is the ratio of the thermal dependence of the non-evaporative species to the viscous forces. Because only 1 and 20 cSt silicone oils were used in the binary system, the physical properties are fixed. When probing the parameter space, the experimentally controlled values of $a$, $b$, $d$ and $\unicode[STIX]{x1D719}_{0}$ are varied. All the parameters are simultaneously varied when $b$ or $a$ is varied, since $\unicode[STIX]{x1D716}=b/a$.

Table 2.

Non-dimensional parameters used. In the definition of $t_{stop}$, the quantity $d$ stands for the dimensional distance that the motor moves. The last column of the table lists the value range or the order of magnitude of the parameters that are reported in this study for a liquid mixture of 1 cSt (dilute) and 20 cSt silicone oils.

3.2 Numerical methods

The six governing equations are numerically solved using the finite difference method. To enhance numerical stability, equation (3.3) is converted into a natural log form and the variable $\log [\unicode[STIX]{x1D719}/\unicode[STIX]{x1D719}_{0}]$ is used instead of $\unicode[STIX]{x1D719}$. The time evolution follows an iterative Crank–Nicholson scheme with adaptive time stepping. At each time step, the linearized set of equations is solved using the biconjugate gradient-stabilized method via an external sparse matrix package (Guennebaud et al. Reference Guennebaud, Jacob and Nuentsa-Wakam2010). Spatially, all unknowns are solved on a smoothly stretched map that clusters 50 % of the grid points over $r\in [0,1]$ and the remaining points over $r\in (1,R_{max})$. Verification simulations were completed and found to match the analytical solutions in the small ${\mathcal{C}}{\mathcal{a}}$ limit following the procedures described in Rodríguez-Hakim et al. (Reference Rodríguez-Hakim, Barakat, Shi, Shaqfeh and Fuller2019).

Numerical validations are done by comparing experimental and simulation apex film drainage rates of a 5000 cSt silicone oil in the absence of evaporation and Marangoni flows (figure 3). The higher-viscosity oil was used to produce slower dynamics than that of a 20 cSt silicone oil. In the absence of Marangoni driving forces, the apex film thins at an exponential rate as dictated by the gravitational–viscous time scale of $\unicode[STIX]{x1D708}/(a\cdot g)$ reported by Debrégeas, De Gennes & Brochard-Wyart (Reference Debrégeas, De Gennes and Brochard-Wyart1998). In the small set of validation experiments done by the authors, only exponential drainage rates were observed. It is possible that in another experimental parameter space there may be algebraic drainage rates as predicted by Howell (Reference Howell1999). Simulations and experiments show good agreement when the film thickness is above 100 nm. When the film becomes thinner than 100 nm, molecular-level interactions, such as van der Waals interactions, become non-negligible. In the experimental data, this phenomenon is reflected in a change in the slope of the apex drainage data. However, in the simulations, there is no distinct change in the slope because the model does not account for these molecular-level interactions.

Figure 3.

Experimental (red circles) and numerical (black solid lines) apex film thickness evolution for three experiments of 5000 cSt silicone oil in the absence of evaporation. The motor speeds are at $0.01$ and $0.05~\text{mm}~\text{s}^{-1}$. In the simulations, it is assumed that there is no evaporative mass loss or Marangoni flow, i.e. $\unicode[STIX]{x1D719}_{0}={\mathcal{E}}{\mathcal{v}}={\mathcal{M}}{\mathcal{a}}={\mathcal{M}}{\mathcal{a}}_{T}=0$.

4.1 Oscillatory film thickness

Figure 4 shows data collected from four experiments of binary silicone oil mixtures (1 cSt 0.1 vol%–20 cSt 99.9 vol%). Note that the experimentally measured film profile can deviate from the axisymmetric state, so the apex film thickness is defined as the maximum film thickness in the dimple region. During the initial second in the experiments, the motor moves the lower surface towards the upper interface and the film quickly thins due to the squeezing motion. The fast dynamics and the relatively large thickness of the film make it difficult to determine the apex film thickness through the interference patterns. For $1~\text{s}<t^{\ast }<t_{stop}^{\ast }=3~\text{s}$, the film thickens (see dashed lines in figure 4a). After the motor stops, the film is allowed to evolve under the effects of evaporation. Thereafter, the apex film thickness fluctuates between 0.5 and $2~\unicode[STIX]{x03BC}\text{m}$ for many minutes. Apex film thickness oscillations are observed and the amplitudes of those oscillations vary between 0.2 and $1~\unicode[STIX]{x03BC}\text{m}$. The time it takes for one cycle of film thickness oscillation to occur in the experiment varies between 1 and 4 s. In the first 20 s of the experiments, the film stays nearly axisymmetric (figure 4b (i–v)); at longer times, there is significant asymmetry in the film profile (figure 4b (vi)).

Figure 4.

Experimental film thickness evolution profiles. (a) Time evolution of the apex film thickness from four binary silicone oil experiments (1 cSt 0.1 vol% and 20 cSt 99.9 vol%). The open symbols mark the experimental apex film thickness and the dashed lines are there to guide the eye. The vertical dash-dot line represents the motor stop time ($t_{stop}^{\ast }$) for all four experiments shown. The six filled symbols in experiment 1 correspond to (b) the six interference patterns at various time points. The diameter of each circular image corresponds to 0.5 mm. Panels (i)–(v) show the growth and decay for one oscillation cycle. Panel (vi) shows symmetry breaking at 20 s for that experiment; note the foot that developed in the two o’clock position. The reference colour map for the interference patterns is provided on the left-hand side of the plot.

Figure 5(a) shows the time evolution of apex film thickness, $h(t,r=0)$, from three numerical simulations. The numerical simulations were performed at three initial concentrations: $\unicode[STIX]{x1D719}_{0}=0.001,0.003$ and $0.005$. The three conditions were chosen to examine the effect of initial concentration on the film thickness oscillation. The initial clearance and bubble radius were chosen to match the corresponding average values in the experiments. All material properties were matched to those of the materials used in the experiments. The motor stop time in the experiments was 3 s, while in the simulations the motor stop time was 1.5 s. During the motor movement period, the apex film thickness in the simulation grows apparently at three times the rate observed in the experiments. It was necessary to shorten the motor stop time in the simulations to avoid further overshoot. This overshoot of the film growth rate is a result of the global deformation of the bubble during approach, which is not included in the theory. In developing the lubrication theory, the far-field deformation of the lower surface is matched to the static shape of a pinned bubble. Whereas in the experiments, as the bubble moves closer to the top interface, it is flattened globally: evidence that the theory is under-accounting for the lower surface deformation. After the motor stops, the film in the simulations quickly thins, followed by oscillations in the film thickness. For all three simulations, the average apex film thickness is in the range of film thickness as observed across the four experiments. At $\unicode[STIX]{x1D719}_{0}=0.001$, the apex film thickness shows clear periodic oscillations, with a single dominant frequency of around 1 Hz (figure 5c). The simulation stops when the minimum film thickness becomes nanoscopic. At $\unicode[STIX]{x1D719}_{0}=0.003$, the film oscillation is first damped where the amplitude of the oscillation decreases from $1~\unicode[STIX]{x03BC}\text{m}$ to a small value. At longer times, the film oscillation then becomes more pronounced, with an oscillation amplitude growing again to $1~\unicode[STIX]{x03BC}\text{m}$. The dominant frequency for the apex film thickness oscillation is around 1.2 Hz (figure 5d). At $\unicode[STIX]{x1D719}_{0}=0.005$, the simulated film thickness oscillates only once after $t_{stop}^{\ast }$ and the subsequent oscillations are damped. No long-time oscillation was observed in this case.

Figure 5.

Apex film thickness evolution profiles and their frequency spectra. (a) Time evolution of the apex film thickness from numerical simulations conducted at initial volume fractions of $\unicode[STIX]{x1D719}_{0}=0.001,0.003$ and $0.005$. The vertical dash-dot line represents the motor stop time ($t_{stop}^{\ast }$). Other non-dimensional parameters corresponding to the three simulations are $\unicode[STIX]{x1D716}=0.067,{\mathcal{C}}{\mathcal{a}}=0.0156,{\mathcal{B}}{\mathcal{o}}=0.0678,{\mathcal{M}}{\mathcal{a}}=173,{\mathcal{M}}{\mathcal{a}}_{T}=255,{\mathcal{P}}{\mathcal{e}}=94,{\mathcal{P}}{\mathcal{e}}_{T}=1.15,{\mathcal{E}}{\mathcal{v}}=0.00728,\unicode[STIX]{x1D703}_{b}=130^{\circ }$. (b–d) Frequency spectra for selected apex film thickness evolution profiles in figure 4(a) and panel (a). The amplitude of the signal is normalized by the amplitude of the maximum signal.

The oscillation frequencies in the experiments and the simulations are clearly different. In the experiments, the oscillations do not show a single dominant frequency (figure 5b shows an example frequency spectrum). For the simulations with sustained oscillations, there is clear periodicity in the oscillations (figure 5c,d). The difference may be caused by sources of noise in the experiments. For example, the disturbances in the evaporation rate can impose temporally and spatially dependent fluctuations to the film thickness. In addition, imperfect mixing in the binary silicone oil mixtures can lead to spatial concentration inhomogeneities that can shift the oscillatory states between the stable oscillatory state and the damped state. Despite the differences, the model does provide a mechanism for the spontaneous oscillations and allows for an estimate regarding the relative contributions from the solutal and the thermal Marangoni flows. Furthermore, we are presenting a free, nonlinear oscillator with physical relevance, and we will discuss how various experimental parameters may affect the oscillation frequency.

4.2 Oscillation mechanism

We introduce the modified surface tension ($\bar{\unicode[STIX]{x1D70F}}$) to explain the mechanism behind the spontaneous oscillations:

(4.1)$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70F}}\equiv \unicode[STIX]{x1D70F}-\unicode[STIX]{x1D70F}(t=0)=-\unicode[STIX]{x1D716}{\mathcal{C}}{\mathcal{a}}({\mathcal{M}}{\mathcal{a}}(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{0})+{\mathcal{M}}{\mathcal{a}}_{T}\unicode[STIX]{x1D6E9}).\end{eqnarray}$$

Throughout the process, the species with the lower surface tension evaporates across the top interface and $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6E9}$ decrease due to the coupled mass and heat transfer process. Thus, surface tension in the thin film will increase as a result of evaporation, i.e. $\bar{\unicode[STIX]{x1D70F}}>0$. The rate at which $\bar{\unicode[STIX]{x1D70F}}$ increases will vary across the film based on the relative rates of change in $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6E9}$, which are coupled to film thickness and curvature.

Figure 6.

Time evolution of interface positions (a–c), film thickness (d–f), surface tension (g–i), and difference between capillary and Marangoni pressure gradients (j–l) for a simulation that corresponds to $\unicode[STIX]{x1D716}=0.067$, ${\mathcal{C}}{\mathcal{a}}=0.0156$, ${\mathcal{B}}{\mathcal{o}}=0.0678$, ${\mathcal{M}}{\mathcal{a}}=173$, ${\mathcal{M}}{\mathcal{a}}_{T}=255$, ${\mathcal{P}}{\mathcal{e}}=94$, ${\mathcal{P}}{\mathcal{e}}_{T}=1.15$, ${\mathcal{E}}{\mathcal{v}}=0.00728$, $\unicode[STIX]{x1D719}_{0}=0.001$ and $t_{stop}=1.5$ (corresponding to the $\unicode[STIX]{x1D719}_{0}=0.001$ case in figure 5a). The legend represents non-dimensional time and the grey arrows indicate increasing time. Panel (a,d,g,j) shows the dimple formation while the needle is moving upward. Panel (b,e,h,k) shows the first capillary discharge after the needle has stopped moving. Panel (c,f,i,l) shows the first film regeneration associated with the Marangoni flows.

Figure 6 shows simulation results for the first cycle of film thickness oscillation, the associated changes in $\bar{\unicode[STIX]{x1D70F}}$ and the pressure gradients, for the $\unicode[STIX]{x1D719}_{0}=0.001$ case presented in figure 5(a). Before $t_{stop}$, the bubble moves towards the initially flat top interface, squeezing liquid away from the centreline into the bulk region, thereby creating a film whose thickness varies in the radial direction. The evaporative species volume fraction depletes faster in regions where the film is thinner; consequently, a surface tension gradient is established.

During $0<t<1$, the minimum film thickness is at $r=0$, where the surface tension is largest. The established surface tension gradient draws liquid from the bulk towards the centreline, generating a Marangoni flow ($1<t<1.61$). Due to the inflow, the film thickness becomes non-monotonic in the radial direction, forming the so-called spontaneous dimple. The formation of the dimple leads to a change in the capillary pressure and the capillary pressure gradient:

(4.2)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{\mathscr{P}}_{{\mathcal{C}}{\mathcal{a}}}}{\unicode[STIX]{x2202}r}=-\frac{1}{4}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}\right)\right).\end{eqnarray}$$

Here $(\unicode[STIX]{x2202}\bar{\mathscr{P}}_{{\mathcal{C}}{\mathcal{a}}}/\unicode[STIX]{x2202}r)$ is the dominant contributing term to the dynamic pressure gradient and its value is predominantly positive due to the gradient in the curvature. Balancing the capillary pressure gradient is the depth-averaged surface tension gradient (equation (3.10)):

(4.3)$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D716}}\frac{1}{h}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}r}=-\frac{{\mathcal{C}}{\mathcal{a}}}{h}\left({\mathcal{M}}{\mathcal{a}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}+{\mathcal{M}}{\mathcal{a}}_{T}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\right).\end{eqnarray}$$

From the radial distribution of the surface tension plot, it is evident that the surface tension gradient is predominantly negative. The difference between the two gradients will give a clear indication as to which is the key mechanism that is driving the flow (figure 6j–l). As the lower surface is forced towards the top surface, the surface tension gradient dominates over the capillary pressure gradient, indicated by the positive radial distribution of the pressure gradient difference.

Eventually the capillary pressure gradient becomes large enough to overcome the Marangoni forces and to dispel liquid from the centreline region into the bulk ($1.61<t<1.86$, figure 6k, $t=1.70$). Because the film near the centreline is thin and the bubble small, gravitational drainage contributes very little as compared to the capillary discharge. After the discharge, the film near the centreline becomes flat and the capillary pressure gradient is again overcome by the surface tension gradient and the difference between the gradients again become positive. Furthermore, the discharge process also re-establishes the surface tension gradient because the centreline region again has the thinnest film and thus the fastest change in $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D6E9}$ – creating a high surface tension again near $r=0$.

Finally, with the new surface tension gradient, the combined solutal and thermal Marangoni flow will again draw liquid from the bulk (where surface tension is lower than that in the centreline region) into the centreline, thereby reducing the surface tension gradient, increasing the capillary pressure gradient and providing the conditions for another fluid discharge and subsequent inflow.

4.3 Relative contributions of solutal and thermal Marangoni flows

In this section, we examine the relative contributions of the solutal and thermal Marangoni flows by comparing the relative contributions of the concentration gradient and the temperature gradient to the overall surface tension gradient:

(4.4)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}r}=-\unicode[STIX]{x1D716}{\mathcal{C}}{\mathcal{a}}\left({\mathcal{M}}{\mathcal{a}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}+{\mathcal{M}}{\mathcal{a}}_{T}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\right).\end{eqnarray}$$

Shown in figure 7 are radial profiles of surface tension gradients at $t=t_{stop}$ for three different initial concentrations that correspond to the simulations presented in figure 5(a). The motion of the motor rising is the initial transient behaviour that starts the subsequent oscillations. The time at which the motor stops is the starting point of the subsequent long-term behaviour, and thus we can gain insights about the system by examining the dynamic variables at $t_{stop}$.

Figure 7.

Radial profiles of surface tension gradient and its components at $t=t_{stop}$ for the three simulations shown in figure 5(a).

At $\unicode[STIX]{x1D719}_{0}=0.001$, the maximum values for solutal and thermal gradients are similar in magnitude, signifying a near-equal contribution from the two sources. The two gradients have different radial shapes: the concentration gradient has a broad peak whereas the thermal gradient grows and decays more sharply. Thus, in this case, the presence of the thermal Marangoni flow significantly modifies the curvature of the film. At $\unicode[STIX]{x1D719}_{0}=0.003$, the maximum magnitude of the concentration gradient is greater than that of the thermal gradient. The maxima of the two gradients occur near $r=0.4$ and the thermal gradient at that point contributes to one-third of the overall surface tension gradient. At $\unicode[STIX]{x1D719}_{0}=0.005$, solutal Marangoni flow contributes roughly three-quarters of the overall Marangoni flow at the maxima. In the latter two cases, the presence of the thermal Marangoni flow serves to enhance the overall Marangoni flow, but it does not significantly affect the curvature of the film, as the radial profiles of the gradients are similar in shape. Comparing results across the three concentrations, the thermal Marangoni flow is not affected by the change in the initial concentration, whereas the solutal Marangoni flow scales with the initial concentration. As a result, the overall surface tension gradient increases as the initial concentration increases and the oscillator becomes more damped as the initial concentration increases.

4.4 Damped oscillators

We examine the nature of these oscillations by looking at the phase portraits of the dimple volume, whose motion describes the composite effect of the Marangoni flows, capillarity, gravity and evaporation. The dimple volume is defined as

(4.5)$$\begin{eqnarray}V=2\unicode[STIX]{x03C0}\int _{0}^{R}hr\,\text{d}r,\end{eqnarray}$$

where $R$ is the radial position of the minimum film thickness. For the three simulations presented in figure 5(a), we plot the dimple volume against its rate of change in time (figure 8). The colour represents non-dimensional time. At $\unicode[STIX]{x1D719}_{0}=0.001$, as the motor moves upward, a dimple forms and by $t=t_{stop}$ the dimple volume reaches a maximum of 0.003, after which the dimple continues to oscillate at a near-constant amplitude and frequency – evidenced by the constant trajectory on the phase portrait. At $\unicode[STIX]{x1D719}_{0}=0.003$, the oscillator first enters an attractor region (around $t=5$) then traverses to a repeller region ($t>10$) where the oscillation amplitude grows and then enters a stable orbit. At $\unicode[STIX]{x1D719}_{0}=0.005$, the dimple grows and only oscillates a few times before entering a damped, dynamic steady state where draining due to capillarity and gravity is balancing the regenerative Marangoni flow. Across the three simulations, the maximum dimple volume increases with concentration, because the dimple-generating Marangoni flow strengthens with initial concentration (figure 7).

Figure 8.

Phase portraits of the dimple volume ($V$) for the three simulations shown in figure 5(a). The colour represents non-dimensional time, as indicated by the colour bars.

4.5 Oscillation frequency

In the experiments, it is relatively easy to vary $\unicode[STIX]{x1D719}_{0}$ and $t_{stop}$. In this section, we look at how these two parameters affect the behaviour of the oscillator. The nonlinear nature of the governing equations makes it prohibitive to derive specific frequency scalings with respect to the parameter space. In this section, we focus on how the experimentally adjustable parameters affect the behaviour of the oscillator through the aid of numerical simulations (figure 9). With the previous set of simulations, it is clear that as $\unicode[STIX]{x1D719}_{0}$ increases, there is an increase in the Marangoni driving force and the system is damped more quickly. The simulations presented in this section all have sustained oscillations, which means that parameters such as the bubble radius and the initial apex clearance have been tuned to avoid cases where the system exhibits a dynamic steady state due to Marangoni damping.

When $\unicode[STIX]{x1D719}_{0}$ is the only parameter that is varied, the amplitude of the dimple volume oscillation is nearly constant (figure 9a). At lower concentrations ($\unicode[STIX]{x1D719}_{0}=0.003$–$0.005$), the dominant frequency increases with concentration and the waveform of the dimple volume does not change significantly. The waveform changes shape when $\unicode[STIX]{x1D719}_{0}$ increases from 0.005 to 0.006. This change may be caused by a shift in the balance between the two contributing Marangoni flows for this specific set of simulation parameters. Figure 9(c) shows the period associated with the dominant frequency when each of the four simulations has reached steady oscillation, i.e. when the effect of the initial motor movement has worn off. As $\unicode[STIX]{x1D719}_{0}$ increases there is no clear trend in the period of the dimple volume oscillation. The period for the dimple to regenerate and to discharge can be further split into the Marangoni regeneration segment (from valley to peak, $1/f_{{\mathcal{M}}{\mathcal{a}}}$) and the capillary discharge segment (from peak to valley, $1/f_{{\mathcal{C}}{\mathcal{a}}}$). As the initial concentration increases, the time it takes for the dimple to regenerate decreases, because the overall Marangoni flow strength scales with $\unicode[STIX]{x1D719}_{0}$. On the other hand, $1/f_{{\mathcal{C}}{\mathcal{a}}}$ shows no clear relation to $\unicode[STIX]{x1D719}_{0}$.

Figure 9.

(a,b) Time evolution of the dimple volumes (left-hand column) and their frequency spectra (right-hand column) at different $\unicode[STIX]{x1D719}_{0}$ and $t_{stop}$. The amplitude of the Fourier signals are normalized by the amplitude of the maximum signal. The frequency is non-dimensionalized by $U/b$. (c,d) Corresponding steady oscillation time scales associated with the dominant frequency ($f_{D}$) and the time scales associated with Marangoni regeneration ($1/f_{{\mathcal{M}}{\mathcal{a}}}$) and capillary discharge ($1/f_{{\mathcal{C}}{\mathcal{a}}}$). (a) Variation in $\unicode[STIX]{x1D719}_{0}$. All other parameters are kept constant at $\unicode[STIX]{x1D716}=0.15$, ${\mathcal{C}}{\mathcal{a}}=0.00307$, ${\mathcal{B}}{\mathcal{o}}=0.0434$, ${\mathcal{M}}{\mathcal{a}}=389$, ${\mathcal{M}}{\mathcal{a}}_{T}=573$, ${\mathcal{P}}{\mathcal{e}}=50$, ${\mathcal{P}}{\mathcal{e}}_{T}=0.612$ and ${\mathcal{E}}{\mathcal{v}}=0.00728$. (b) Variation in $t_{stop}$. All other parameters are kept constant at $\unicode[STIX]{x1D716}=0.067$, ${\mathcal{C}}{\mathcal{a}}=0.0156$, ${\mathcal{B}}{\mathcal{o}}=0.0678$, ${\mathcal{M}}{\mathcal{a}}=173$, ${\mathcal{M}}{\mathcal{a}}_{T}=255$, ${\mathcal{P}}{\mathcal{e}}=93.8$, ${\mathcal{P}}{\mathcal{e}}_{T}=1.15$ and ${\mathcal{E}}{\mathcal{v}}=0.00728$. (c) Time scales for steady oscillations shown in (a). (d) Time scales for steady oscillations shown in (b).

When $t_{stop}$ is the only increasing parameter, the amplitude of the oscillations increases because the duration of the motor movement is longer (figure 9b). When the motor is moving, Marangoni flows dominate and capillarity is not strong enough to discharge any liquid. This behaviour is reflected in the increase in apex film thickness during $1<t<t_{stop}$ and it is seen in both the simulations and the experiments. Thus, a longer time period for the initial inward flow will lead to a larger initial dimple and an increase in oscillation amplitude. The dominant frequency decreases with increasing $t_{stop}$ value, because a larger dimple volume takes longer to discharge. This trend is also reflected in the regeneration and discharge time scales (figure 9d). As the motor stop time increases, both processes take longer to complete, leading to the overall increase in the dominant period.