4.1. Symmetry breaking

In this section, we describe symmetry breaking associated with the dynamics of an evaporating binary liquid film over a spherical glass dome. Expanding on our previous study (Rodríguez-Hakim et al. Reference Rodríguez-Hakim, Barakat, Shi, Shaqfeh and Fuller2019), the present work is focused on binary blends of 1 cSt and 5 cSt silicone oils at a higher volume fraction of the non-volatile, high surface tension species (5 cSt, $\phi _0 \geqslant 0.0035$). Similar to the thin film dynamics that we previously described as being in the ‘capillary regime’ (1 cSt and 5 cSt blend, $0.0002 \lesssim \phi _0 \leqslant 0.003$), the film dynamics for higher concentrations of the non-volatile species demonstrates initial drainage due to the squeezing motion of the glass dome. A concentration gradient is simultaneously developed as the volatile species evaporates. This concentration gradient drives a Marangoni flow toward the apex of the glass dome, forming a solute-rich mound, bounded by the rim, $R(t,\theta )$, where the film thickness is at its minimum. Because the solute-enriched mound region is characterized by a higher surface tension than that in the bulk surrounding region, the liquid is continuously drawn into the mound. In the previously described capillary regime, the mound volume continues to grow in an axisymmetric manner throughout this process. However, at higher initial non-volatile solute concentration, $\phi _0$, while the mound volume is growing, the film thickness at the rim of the mound grows unevenly, leading to symmetry breaking and non-axisymmetric liquid discharge from the mound into the bulk.

Figure 2 contains interferograms and their corresponding film thickness maps, before, during and after a mound discharge event for a 1.00–5.00 cSt ($\phi _0 = 0.005$) blend. Figure 3 provides the parameter matched simulation film thickness and solute concentration surfaces at the corresponding time points. At the start of visible symmetry breaking, the film thickness and solute concentration fields remain near axisymmetric (figures 2d and 3a,d). As the disturbances grow, one side of the film thins while its opposite side thickens (figures 2e and 3b), leading to a visible loss of axisymmetry. As the asymmetry becomes even more apparent, the apex of the mound shifts away from $r = 0$ and toward the direction of film thickening. Simultaneously, the solute-rich region is convected toward the region of film thickening (figure 3e). The spatial shift in the high surface tension region leads to an azimuthal Marangoni flow that pulls more liquid toward the location of the radially shifted mound. Eventually, the region of increased film thickness grows to such an extent that the previously well-defined rim becomes a shoulder and there is no clear distinction between the mound region and the thicker, bulk region in the direction of film thickening (figures 2f and 3c). This process completes one cycle of mound discharge.

Figure 2.

Interference patterns (a–c, $t = 165$, $t = 218$, $t = 236$) and the corresponding experimental dimensionless film thickness profiles (d–f) during a mound discharge event for a binary mixture of 1 cSt silicone oil with 0.5 vol% of 5 cSt silicone oil. The diameter of the interferograms corresponds to 1 mm and the horizontal bars in the film thickness profiles indicate a dimensionless radial length scale of 1.

Figure 3.

Film thickness and concentration profiles during a mound discharge event of a simulation conducted at $\phi _0 = 0.005$ and $Pe = 200$. The other simulation parameters are $Ca = 2.426 \times 10^{-6}$, $Bo = 0.0427$, $Ma/\phi _0 = 106.4$, $Ev = 0.0026$, $t_{stop} = 29.57$ and $Ha = 1.438\times 10^{-5}$. For the rest of the paper, these parameters are kept constant. (a–c) Film thickness surface plots at time points 356, 369 and 381. The horizontal bar indicates a dimensionless radial length scale of 1. (d–f) Contour plots of the non-evaporative species concentration (in filled colour) overlaid with the location of the rim (white lines). The directions of mound discharge (‘Ds’), film thinning (‘Th’) and the locations between the two points (‘$\textrm {Ds} + {\rm \pi}/2$’ and ‘$\textrm {Th} + {\rm \pi}/2$’) are labelled with white dots. The arrows plotted over the four positions indicate the general flow direction at that point. The length of the arrows gives a general indication of the growth in flow magnitude. The arrow lengths are not drawn to scale.

In the experiments, ambient disturbance events such as ambient air disturbances and vibrations can induce a mound discharge event; in the numerical simulation shown in figure 3, the accumulated numerical error leads to the eventual mound discharge event. Despite the difference in the source of the disturbance, the simulation and the experiment agree quantitatively in the form of the break in symmetry, suggesting that the most dangerous mode follows the form of the observed disturbance in the simulations.

Based on the observations above, we conducted an analysis of the evolution of the rim thickness and the apex location to further characterize the mound discharge event and compare the simulation against the experiment. The rim film thickness ($h_R$) corresponding to the panels shown in figures 2 and 3 is extracted and plotted in figures 4(a) and 4(d). When plotted against the azimuthal angle, the rim film thickness takes the approximate form of a cosine function. Fourier transform of $h_R$ confirms this observation (figure 4b,e): as the mound discharges, the amplitude of the signal with an azimuthal wavenumber of one grows ($m = 1$, vertical dashed line). Thus, for this system, apparently the disturbance that is amplified follows the form of $\tilde {\psi }(t,r,\theta ) = \cos (\theta ) \hat {\psi }(t,r)$, where $\psi = h_1, \phi$. During the discharge, the apex of the mound moves away from $r = 0$ at an exponential rate (figure 4c,f). In the simulation, during $0 < t < 310$, the flat line in the radial position of the apex signifies the smallest radial grid point away from $r = 0$ (for details on meshing, see the numerical methods section).

Figure 4.

Comparison of selected quantities associated with experimental data (a–c) presented in figure 2 and the parameter matched simulation (d–f) results presented in figure 3. (a,d) Rim film thickness before, during and after mound discharge. (b,e) Fourier transform of the rim film thickness presented in (a,d). The vertical dashed line is plotted over the azimuthal wavenumber $m = 1$ to guide the eye. (c,f) Radial position of the apex of the mound as a function of time.

4.2. Flow characteristics during symmetry breaking

In this section, we examine flow characteristics during the symmetry breaking event by comparing the results of two-dimensional and one-dimensional (1-D) (i.e. axisymmetric) simulations conducted with the same set of parameters (i.e. those presented in figure 3). The most visible indication of a mound discharge event is the deviation in the mound volume ($V_m$ defined below) in the 2-D simulation as compared to the 1-D simulation. In the 1-D simulation, the assumption of axisymmetry restricts the solution of the governing equations to be only a function of $t$ and $r$. Therefore, the solution does not have a second spatial dimension into which momentum can be transferred. In the 2-D simulation, the liquid can move in both the radial and the azimuthal directions. Thus in these two cases, there is a distinct difference in the evolution of the ‘mound volume’, defined as

(4.1)\begin{equation} V_m(t) = \int_0^{2{\rm \pi}} \int_0^{R(t,\theta)} h(t,r,\theta) r \,\textrm{d}r \,\textrm{d}\theta, \end{equation}

where the rim position, $R(t,\theta )$, is the radial location of the minimum film thickness in a given azimuthal direction. For the 1-D simulation, this quantity simply becomes

(4.2)\begin{equation} V_{m,1D}(t) = 2{\rm \pi} \int_0^{R(t)} h(t,r) r \,\textrm{d}r. \end{equation}

Figure 5(a) shows the comparison of the mound volume evolution of the 2-D and 1-D simulations. In both sets of simulations, the mound volume first undergoes a transient growth, after which the mound enters a monotonic growth region beyond $t = 76$. From $t = 76$ until approximately $t = 356$, the mound volume grows at the same rate for both 1-D and 2-D simulations, and there is apparently no break in symmetry. At $t = 368$, the disturbance in the 2-D simulation has grown to a large enough amplitude to create an asymmetric mound discharge event, thus leading to a sharp drop in the mound volume. In the 1-D simulation, the mound volume continues to grow.

Figure 5.

Comparison of the time evolutions of selected quantities associated with the simulation shown in figure 3 and the corresponding 1-D axisymmetric simulation. (a) Mound volume evolution. (b) Apex and rim film thickness. In the 1-D simulation, the apex is located at $r = 0$ and the rim is located at the position of minimum film thickness. For the 2-D simulation, the film thicknesses at the apex and at the rim on the opposing sides along the direction of the mound discharge are plotted; ‘Ds’ and ‘Th’ stand for ‘discharging side’ and ‘thinning side’ respectively. (c) Radial velocity at the rim locations. (d) Pressure (‘$P$’) and Marangoni (‘$M$’) contributions to the radial velocity at the rim locations. (e) Azimuthal velocities at the rim locations that are ${\rm \pi}$/2 offset from the discharging and thinning sides. (f) Pressure and Marangoni contributions to the azimuthal velocity at the rim locations presented in (e).

A second set of measures that indicate symmetry breaking are the film thicknesses at both the apex and the rim. In a 2-D simulation, the apex is defined as the location of the maximum film thickness within the mound region; whereas in a 1-D simulation, the apex is always located at $r = 0$ (figure 5b). Similar to the mound volume evolution, the apex film thickness in the 2-D simulation decreases at around $t = 368$ (i.e. at the onset of a discharge event), while its 1-D counterpart grows indefinitely due to the symmetry constraint. For the 2-D simulation, the locations of the maximum and minimum rim film thickness can be examined. These opposing rim locations are of interest because they are positioned along the axis of mound discharge, and thus where film thickness changes are most rapid and pronounced. The direction in which the film thickens during discharge is referred to as the discharging side (‘Ds’). This is the side where the apex and rim will eventually merge, thus signalling the end of a mound discharge event. The opposing side is referred to as the thinning side (‘Th’), where the film thickness decreases during discharge.

To examine the total liquid flow into and out of the mound region, we look at the depth averaged radial and azimuthal volumetric flux per unit circumference evaluated at the rim

(4.3a,b)\begin{align} \left\langle {v_r} \right\rangle_R = \left.\left( {-\frac{h^2}{3}\frac{\partial {P}}{\partial {r}} + \frac{Ma}{2}h\frac{\partial {\phi}}{\partial {r}}} \right) \right|_R,\quad \text{and}\quad \left\langle {v_\theta} \right\rangle_R = \left.\left( { -\frac{h^2}{3}\frac{1}{r}\frac{\partial P}{\partial \theta} + \frac{Ma}{2} h \frac{1}{r}\frac{\partial \phi}{\partial \theta}} \right) \right|_R. \end{align}

Both velocities have two contributing terms: the pressure contribution and the Marangoni contribution (figure 5c–f). For the pressure contribution of the radial velocity, the dominant term comes from capillary forces. Based on the definition of the coordinate system, a positive value in $\left \langle {v_r} \right \rangle _R$ signifies liquid flow away from $r = 0$ towards the bulk. During the transient mound growth ($t_{stop} < t < 33$), the film thickness grows monotonically in all azimuthal directions. The volume of the mound is then reduced quickly as liquid flows out axisymmetrically due to the dominant pressure contributions (figure 5(d), $33 < t < 76$). During this period, a concentration gradient is established such that the Marangoni term overwhelms the pressure term, leading to an inward radial flux toward $r = 0$ ($\left \langle {v_r} \right \rangle _R < 0$, $76 < t < 178$). The film thickens as a result of Marangoni regeneration, thereby increasing the magnitude of both contributions and reducing the magnitude of the total inward flow (figure 5(c), $178 < t < 356$). Soon after, the difference in rim film thickness between the discharging and thinning side increases. On the discharging side, the film thickens and the pressure contribution to the radial flow eventually overwhelms the Marangoni flow, leading to the eventual outward discharge flow (figure 5(c) when the dash-dotted line goes from being negative to positive). On the thinning side, the film thickness is further reduced and both contributions to the radial flow decrease in magnitude; however, the Marangoni contribution decreases at a slower rate, leading to a stronger inward flow (in figure 5(c) the solid line grows more negative). In addition to the change in flow magnitude on the two opposing sides, as the mound is shifted in the direction of discharge, the solute rich mound gets convected. The off-centre mound creates secondary Marangoni flows in the azimuthal direction (thus a non-zero $\left \langle {v_\theta } \right \rangle _R$) that draws more liquid into the mound and speeds its radial discharge (figure 3d–f).

Finally, we look at the azimuthal flow $\left \langle {v_\theta } \right \rangle _R$ at selected rim locations. Figure 5(e,f) shows the magnitude of the azimuthal flow at rim locations that are ${\rm \pi} /2$ away from the discharging and thinning locations. These locations were chosen because they have the maximum azimuthal flow amplitude; whereas, at the discharging and thinning locations, the magnitude of the azimuthal flow is small due to the mirror symmetry along the axis of mound discharge. The azimuthal flow has negligible magnitude until significant symmetry breaking occurs ($0 < t < 330$). When the mound discharges along the axis of symmetry breaking, it convects the high surface tension liquid away from the $z$-axis, thus rendering the concentration/surface tension profile non-axisymmetric. The disruption in the concentration profile leads to an increase in the magnitude of the azimuthal Marangoni flow, which is slightly damped by the pressure contributions (figure 5f). Overall, at ${\rm \pi} /2$ away from the discharging and thinning rim locations, the azimuthal flows drive fluid toward the discharging side, thus further thickening the film at that location (figure 5(f), and see the arrows in figure 3).

4.3. Effects of solute concentration and diffusivity on the onset of symmetry breaking

As one can deduce, of primary interest are the effects of $\phi _0$ and $Pe$ on the onset of symmetry breaking. The experimental set-up is subjected to ambient disturbances throughout an experiment. In simulations, we can inject a known form of disturbance at a specific time point to determine its effect on the subsequent mound discharge. Because the sources and amplitudes of the disturbance in the experiments and the simulations are different, we choose to compare the trends in how the experiments and simulations respond to changes in one input parameter. To this end, we introduce the following form of disturbance at $t = 127$ for all the simulations shown in this section:

(4.4a,b)\begin{equation} \tilde{h}_{1,ini} = C_h \cos(\theta) \frac{r}{2R^3 + r^3}, \quad\text{and}\quad \tilde{\phi}_{ini} = C_\phi \cos(\theta) \frac{r}{4R^5 + r^5}, \end{equation}

where $C$ is a normalization coefficient. The analytical form of the injected disturbance was found by fitting a smooth surface to the difference between the 2-D solution during mound discharge and the corresponding solution in the 1-D simulation. Noise-injection time of $t = 127$ is chosen such that all simulations studied are in the Marangoni regeneration regime and have yet to break axisymmetric flow. The rim location, $R$, is found for each simulation at the time of the noise injection. The normalization coefficients, $C_h$ and $C_\phi$, are chosen such that the injected disturbance is sufficiently large to create a system response that is not overshadowed by the intrinsic numerical difference between the 2-D and the axisymmetric simulations, while small enough such that there is agreement with the linearized disturbance evolution at short times. We will devote more discussion to the linear stability analysis in the next section. (See Appendix A for the specific values of the normalization coefficients and a comparison of the radial profiles of the disturbance quantities.)

We study the binary mixtures at four initial concentrations ($\phi _0 = 0.0035, 0.005, 0.01, 0.02$). These concentrations were chosen such that the onset of the discharge time in the experiments is individually distinguishable. Figure 6 shows experimental and simulation apex film thickness and mound volume evolution up to the end of the first mound discharge event. Marangoni driving forces increase in magnitude as the initial concentration increases. As a result, Marangoni regeneration happens at an earlier time and mound discharge also happens more quickly. Figure 7 shows the experimental and simulation mound discharge time plotted as a function of initial concentration. In the simulations, three Péclet numbers were examined to account for the uncertainty in diffusivity. The diffusivity that corresponds to $Pe = 200$ is extrapolated from four binary diffusivities of 1 cSt silicone oil at 1000, 10 000, 60 000 and 100 000 cSt, reported by Walls et al. (Reference Walls, Meiburg and Fuller2018). Clearly the extrapolation has a high level of uncertainty associated with it. Thus we examine three Péclet numbers to account for this uncertainty.

Figure 6.

Time evolution of experimental (a) and simulation (b) apex film thickness at four concentrations. The error bars in the experimental data represent the standard deviation taken from fifteen experiments conducted for each concentration. The simulations were conducted at $Pe = 200$, $Ca = 2.426 \times 10^{-6}$, $Bo = 0.0427$, $Ma/\phi _0 = 106.4$, $Ev = 0.0026$, $t_{stop} = 29.57$ and $Ha = 1.438\times 10^{-5}$. For each simulation, the same disturbance was introduced at $t = 127$.

Figure 7.

Effects of initial concentration on experimental (a) and simulation (b) mound discharge time. For each simulation, the same disturbance was introduced at $t = 127$. Except for $\phi _0$ and $Pe$, all other simulation parameters are the same as the ones used in figure 3.

For all three Péclet numbers, as the concentration increases, the time for the mound to discharge decreases, in agreement with the experimental observations. For the simulation conducted at $\phi _0 = 0.0035, Pe = 100$, the disturbance ultimately decayed, indicating a stability boundary in the parameter space of $\phi _0$ and $Pe$. This strongly suggests that diffusion plays a crucial role in stabilizing the system. In the experiments, the film thickness profile remained axisymmetric for $\phi _0 \leqslant 0.003$ and symmetry breaking is observed for $\phi _0 \geqslant 0.0035$. This observation of the stability cross-over in the experiments further validates the predictions of our theoretical model. Note that no attempt has been made to examine the complete stability boundary in $\phi _0$ and $Pe$ space, and this is a fruitful direction for future studies.

In the range of concentrations and Péclet numbers studied, it is unclear if there is a critical film thickness or a critical mound volume at which the system will discharge. A closer look at the disturbance evolution is necessary to elucidate the mechanism behind this instability. In this context, we therefore turn to linear stability analysis.

4.4. Linear stability analysis

In this section, we examine the system response to an asymptotically small perturbation, by looking at the linearized disturbance evolution. Based on experimental observations, the initially axisymmetric film evolves until it transitions to a non-axisymmetric state and the Marangoni-regenerated mound discharges along a given axis. The constant presence of evaporation indicates that the system continuously evolves in time. Thus, the base state in this case is the axisymmetric, 1-D solution that evolves in time and is a function of only the radial coordinate. The complexity of the base state suggests that we formulate a linearized initial value problem to study the transient response of the system to small perturbations (Schmid Reference Schmid2007), without making the assumption of exponential disturbance growth in time (the usual eigenvalue problem). The linearized disturbance equations are derived by taking the difference in the governing equations for the 2-D and the 1-D systems, followed by neglecting all nonlinear terms in the disturbance quantities, $\tilde {\psi } \equiv \psi _{2D} - \psi _{1D}$, where $\psi$ represents $h_1$ (position of the top air/oil interface), $\phi$, $\left \langle {v_r} \right \rangle$, $\left \langle {v_\theta } \right \rangle$ and $P$. We seek solutions in the form of $\tilde {\psi } = \exp (\textrm {i}m\theta )\hat {\psi }(t,r)$, i.e. the disturbance quantities in the azimuthal direction follow a sinusoidal form. After some algebra, the disturbance quantities are governed by the following set of linearized disturbance equations:

(4.5)\begin{align} &\frac{\partial {\hat{h}_1}}{\partial {t}} + \frac{1}{r}\frac{\partial}{\partial r}\left(rh_{1D}\left\langle {\hat{v}_r} \right\rangle\right) + \frac{1}{r}\frac{\partial}{\partial r}\left(r\left\langle {v_r} \right\rangle_{1D}\hat{h}_1\right) + h_{1D} \frac{\textrm{i}m}{r} \left\langle {\hat{v}_\theta} \right\rangle - Ev\phi_0 \hat{\phi} = 0, \end{align}

(4.6)\begin{align} &\frac{\partial {\hat{\phi}}}{\partial {t}} +\left\langle {v_r} \right\rangle_{1D} \frac{\partial {\hat{\phi}}}{\partial {r}} +\frac{\partial {\phi_{1D}}}{\partial {r}} \left\langle {\hat{v}_r} \right\rangle -\frac{1}{Pe}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial {\hat{\phi}} }{\partial r}\right) +\frac{1}{Pe}\frac{m^2}{r^2}\hat{\phi} +\cdots \nonumber\\ &\quad-\frac{1}{Pe} \frac{1}{h_{1D}} \frac{\partial {\phi_{1D}}}{\partial {r}} \frac{\partial {\hat{h}_1}}{\partial {r}} - \frac{1}{Pe} \frac{1}{h_{1D}} \frac{\partial {h_{1D}}}{\partial {r}} \frac{\partial {\hat{\phi}}}{\partial {r}} +\frac{1}{Pe} \frac{1}{h_{1D}^2} \frac{\partial {h_{1D}}}{\partial {r}} \frac{\partial {\phi_{1D}}}{\partial {r}} \hat{h}_1 \nonumber\\ &\qquad = Ev \frac{1 - 2\phi_0\phi_{1D}}{h_{1D}}\hat{\phi} - Ev \frac{\phi_{1D}(1 - \phi_0\phi_{1D})}{h_{1D}^2} \hat{h}_1, \end{align}

where

(4.7)\begin{gather} \left\langle {\hat{v}_r} \right\rangle ={-}\frac{h_{1D}^2}{3} \frac{\partial {\hat{P}}}{\partial {r}} -\frac{2}{3} h_{1D} \frac{\partial {P_{1D}}}{\partial {r}} \hat{h}_1 +\frac{Ma}{2}\frac{\partial {\phi_{1D}}}{\partial {r}} \hat{h}_1 +\frac{Ma}{2}h_{1D}\frac{\partial {\hat{\phi}}}{\partial {r}}, \end{gather}

(4.8)\begin{gather}\left\langle {\hat{v}_\theta} \right\rangle ={-}\frac{h_{1D}^2}{3} \frac{\textrm{i}m}{r} \hat{P} + \frac{Ma}{2} h_{1D} \frac{\textrm{i}m}{r} \hat{\phi},\quad \text{and} \end{gather}

(4.9)\begin{gather}\hat{P} ={-}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial {\hat{h}_1} }{\partial r}\right) + \frac{m^2}{r^2} \hat{h}_1 + Bo \hat{h}_1 + \frac{3Ha}{h_{1D}^4} \hat{h}_1. \end{gather}

In the numerical implementation, the above equations are rearranged such that only $\hat {h}_1$ and $\hat {\phi }$ are solved with an initial condition that captures the form, as described above (4.4a,b), of what is apparently the ‘most dangerous mode’ – i.e. a smooth approximation to the difference between the fully 2-D simulations during mound discharge and the corresponding 1-D axisymmetric state. By setting $m$ to an integer value and applying the following initial and boundary conditions, we complete the initial value problem

(4.10)\begin{gather} \hat{h}_1(t_{ini},r) = \frac{r}{2R^3 + r^3}, \end{gather}

(4.11)\begin{gather}\hat{\phi}(t_{ini},r) = \frac{r}{4R^5 + r^5}, \quad\text{and} \end{gather}

(4.12)\begin{gather}\text{as } r\rightarrow\infty, \quad\hat{h}_1 \rightarrow 0, \quad\frac{\partial {\hat{h}_1}}{\partial {r}} \rightarrow 0, \quad\frac{\partial {^2\hat{h}_1}}{\partial {r^2}} \rightarrow 0, \quad\frac{\partial {^3\hat{h}_1}}{\partial {r^3}} \rightarrow 0,\quad \hat{\phi} \rightarrow 0, \quad\frac{\partial {\hat{\phi}}}{\partial {r}} \rightarrow 0. \end{gather}

In the linear stability analysis, we are particularly interested in the ‘gain’ of the linearized disturbance quantities, specifically,

(4.13)\begin{equation} G(t) \equiv \frac{\displaystyle\int_0^\infty \left( {\hat{h}_1^2 + \hat{\phi}^2} \right) r \,\textrm{d}r}{ \displaystyle\int_0^\infty\left( {\hat{h}_{1,ini}^2 + \hat{\phi}_{ini}^2} \right) r \,\textrm{d}r}. \end{equation}

Correspondingly, the gain associated with the disturbances injected into the 2-D simulation is defined as

(4.14)\begin{equation} G_{2D}(t) = \frac{\displaystyle\int_0^{2{\rm \pi}} \displaystyle\int_0^\infty \left( {(h_{1,2D} - h_{1,1D})^2 + (\phi_{2D} - \phi_{1D})^2} \right) r \,\textrm{d}r \,\textrm{d}\theta} {\int_0^{2{\rm \pi}} \int_0^\infty \left( {\tilde{h}_{1,ini}^2 + \tilde{\phi}_{ini}^2} \right) r \,\textrm{d}r \,\textrm{d}\theta}. \end{equation}

Note that the gains, as defined in (4.13) and (4.14), are suitably scaled inner products of the fluctuations in $h_1$ and $\phi$. When the gain of the disturbances decays to a value well below one, the system is stable, whereas a value above one indicates disturbance growth. By comparing the gain associated with the linearized disturbances, $G(t)$ and $G_{2D}$ , we can also establish the time range in which the linear stability analysis is applicable. To such an end, we examine two simulations in the $\phi _0 - Pe$ parameter space: $\phi _0 = 0.0035, Pe = 100$ and $\phi _0 = 0.005, Pe = 200$. In the full 2-D simulations $\phi _0 = 0.0035, Pe = 100$ is the only parameter set with a long-time stable response, whereas the other parameter set gives an unstable response.

Figure 8 shows the evolution of the gain of the disturbances for these two parameter sets. In the 2-D simulations, the disturbance defined by (4.4a,b) is introduced at $t = 127$ as a one-time noise injection. The 2-D simulations are evolved and we track the deviation from the corresponding axisymmetric, 1-D simulation without noise injection. We examine the most dangerous azimuthal mode $m = 1$. Higher-order modes are found to be more stable (see Appendix A).

Figure 8.

(a,b) Contain comparisons of mound volume evolution between the corresponding noise-injected 2-D simulations and the axisymmetric 1-D simulations. (c,d) Contain comparisons of the gain of the linearized disturbance (circles) to the gain associated with the injected disturbances in the 2-D simulations (lines). (a,c) Correspond to the simulations conducted at $\phi _0 = 0.0035$ and $Pe = 100$. (b,d) Correspond to the simulations conducted at $\phi _0 = 0.005$ and $Pe = 200$. The noise is injected at $t = 127$.

The evolution in mound volume gives a direct visual clue to the onset of mound discharge. For the stable case ($\phi _0 = 0.0035$, $Pe = 100$, figure 8a), up to $t = 600$, there is no significant difference between the mound volume in the 2-D simulation with noise injection and that in the axisymmetric simulation, signifying that the injected noise decays over time. For the unstable case ($\phi _0 = 0.005$, $Pe = 200$, figure 8b), the mound discharged by $t = 200$, shortly after the noise injection, signifying that the injected noise brings about the onset of symmetry breaking.

Figures 8(c) and 8(d) compare the disturbance gain obtained from the 2-D simulation with noise injection and from linearized disturbance equations. For the stable case, there is good agreement in the gain from the noise injection at $t = 127$ until $t = 250$ (figure 8c). The amplitude of the gain first grows, reaching its peak value at $t = 150$. It then decays to its initial value of one at $t = 200$. Thereafter, the disturbance further decays, showing that the system is linearly stable for the parameter set of $\phi _0 = 0.0035$ and $Pe = 100$. In the unstable case, the disturbance gain grows immediately after noise injection (figure 8d). The linearized disturbance growth follows that of the 2-D simulation until $t = 190$, at which point the mound in the 2-D simulation starts to discharge. From $t = 190$ onward, the amplitude of the gain of the linearized disturbances grows exponentially, while the gain in the 2-D simulation is affected by secondary, nonlinear flows associated with symmetry breaking.

The linear stability analysis conducted for these two parameter sets has shown that, at $t = 127$, one is linearly stable while the other one is linearly unstable. This analysis will now be carried out at different time points, i.e. the initial value problem will be solved at different starting points during Marangoni regeneration. For the linearly stable case ($\phi _0 = 0.0035$ and $Pe = 100$, figure 9a), the system shows a linearly stable response at $t_{ini} = 127$, 152, 178, 203, 229 and $254$. The growth rate and magnitude of the gain differ for each $t_{ini}$, but they all eventually decay, showing that this parameter set is indeed linearly stable. For the other parameter set, the disturbances immediately grow when the linearized disturbance equations are solved (figure 9b), confirming that the system is unstable. From these two figures we conclude that the linear stability analysis can accurately depict the stability for the two parameter sets. We now present the factors that affect the evolution of the gain of the linearized disturbance quantities and reveal the mechanism behind this instability.

Figure 9.

Evolution of the gain of the linearized disturbance quantities. (a) Corresponds to simulations conducted at $\phi _0 = 0.0035$ and $Pe = 100$. (b) Corresponds to simulations conducted at $\phi _0 = 0.005$ and $Pe = 200$. Each simulation is a solution to the linearized initial value problem, starting at $t_{ini}$. The horizontal dashed lines are reference lines with values of one. Values above the dashed line indicate disturbance growth, when compared to the initial condition.

4.5. Symmetry breaking mechanism

To elucidate the mechanism behind the observed instability, we analyse the factors that affect the evolution of the gain of the linearized disturbance quantities. By algebraically manipulating the linearized disturbance equations, we obtain a scalar equation that relates the time rate of change in the gain, $G$, to the four main physical factors that affect the disturbances, namely Marangoni flow ($Ma$), pressure-driven flow ($P$), diffusion ($Pe$) and evaporation ($Ev$), i.e.

(4.15)\begin{equation} \frac{1}{2}\frac{\textrm{d}G}{\textrm{d}t} = \frac{\textrm{d}G_{Ma}}{\textrm{d}t} + \frac{\textrm{d}G_{P}}{\textrm{d}t} + \frac{\textrm{d}G_{Pe}}{\textrm{d}t} + \frac{\textrm{d}G_{Ev}}{\textrm{d}t}. \end{equation}

A positive value contributes to disturbance growth, while a negative value indicates a stabilizing influence on the disturbances. Note that, while we have denoted the terms ${\textrm {d}G_{P}}/{\textrm {d}t}$ as ‘pressure-driven flow’, we found in all cases that the pressure was almost entirely created by capillary forces. Thus, we could have equally called this term ‘capillary-driven flow’. The detailed expressions for each of these terms can be found in Appendix A. Herein, we focus on the evolution of each of these four terms and the key contributing terms to ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ and ${\textrm {d}G_{P}}/{\textrm {d}t}$

(4.16)\begin{gather} A_{i} \equiv{-}\frac{\displaystyle\int_0^\infty \hat{h}_1 \frac{1}{r}\frac{\partial}{\partial r}\left(rh_{1D}\left\langle {\hat{v}_r} \right\rangle_{i}\right) r \,\textrm{d}r} {\displaystyle\int_0^\infty \left( {\hat{h}_{1,ini}^2 + \hat{\phi}_{ini}^2} \right) r \,\textrm{d}r} ,\quad \text{where}\ i = Ma, P \end{gather}

(4.17)\begin{gather}B_{Ma} \equiv \frac{\displaystyle\int_0^\infty \hat{h}_1 \frac{Ma}{2}\frac{m^2}{r^2} h_{1D}^2 \hat{\phi} r \,\textrm{d}r} {\displaystyle\int_0^\infty \left( {\hat{h}_{1,ini}^2 + \hat{\phi}_{ini}^2} \right) r \,\textrm{d}r}, \end{gather}

(4.18)\begin{gather}B_{P} \equiv{-}\frac{\displaystyle\int_0^\infty \hat{h}_1 \frac{m^2}{r^2}\frac{h_{1D}^3}{3}\hat{P} \,\textrm{d}r} {\displaystyle\int_0^\infty \left( {\hat{h}_{1,ini}^2 + \hat{\phi}_{ini}^2} \right) r \,\textrm{d}r}. \end{gather}

Here, $A_{Ma}$ and $A_{P}$ are associated with the radial advection of fluid mass, and $B_{Ma}$ and $B_{P}$ are associated with flows in the azimuthal directions.

Figure 10 presents the time evolution of these quantities, for the disturbance quantities defined by (4.10) and (4.11). The linearized disturbance equations are solved starting at $t_{ini} = 127$, with an azimuthal wavenumber of $m = 1$, i.e. the same disturbance introduced in the linearized stability analysis presented in figure 8. The solution to the linearized disturbance equations is then integrated to yield the terms in (4.15). The analysis is conducted for the two cases presented in the previous section: the first one corresponds to the linearly stable case ($\phi _0 = 0.0035, Pe = 100$, figure 10a–c), and the second one corresponds to the unstable case ($\phi _0 = 0.005, Pe = 200$, figure 10d–f).

Figure 10.

Evolution of the time rate of change in the gain ($G$) for two parameter sets: $\phi _0 = 0.0035, Pe = 100$ (a–c) and $\phi _0 = 0.005, Pe = 200$ (d–f). The linearized disturbance equations are solved starting at $t = 127$, with $m = 1$ (same noise introduced in figure 8); ${\textrm {d}G_{Ma}}/{\textrm {d}t}$, ${\textrm {d}G_{P}}/{\textrm {d}t}$, ${\textrm {d}G_{Pe}}/{\textrm {d}t}$, and ${\textrm {d}G_{Ev}}/{\textrm {d}t}$ represent the total contributions of Marangoni flow, capillarity, diffusion and evaporation to the time rate change of the gain. In panels (a,d) and (b,e), the lines labelled ‘$A_{Ma}$’, ‘$B_{Ma}$’ and ‘$A_{P}$’, ‘$B_{P}$’, respectively, represent the top two contributing terms for Marangoni flow and for capillarity. Their analytical expressions are defined in (4.16), (4.17), and (4.18).

In the first case, ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ is positive and contributes to the growth of the gain while ${\textrm {d}G_{P}}/{\textrm {d}t}$ is negative and serves to stabilize the system. From $t_{ini} = 127$ to $t = 140$, the growth in ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ is dominated by $A_{Ma}$ and the decay in ${\textrm {d}G_{P}}/{\textrm {d}t}$ is driven by $A_{P}$. As the disturbances evolve, the terms related to the azimuthal flow grow in magnitude (see the peaks of $B_{Ma}$ and $B_{P}$ after $t = 140$ in figure 10a,b). These observations suggest that, initially, the radial advective Marangoni flow acts to destabilize while the radial pressure-driven flow counteracts and stabilizes. These interactions lead to a growth in the azimuthal flows, which further accentuates the destabilizing nature of Marangoni flows and the stabilizing effects of the pressure-driven flow. The total sum of these two effects are overall destabilizing (red dashed line in figure 10c), which is offset by the stabilizing effect of diffusion primarily at long times ($t = 160\text {--}240$), leading to an eventual decay in the gain. Evaporation is destabilizing for short times, until $t = 175$, and thereafter is nearly negligible. Thus, diffusion is the primary stabilizing factor, as we had anticipated from the results shown in figure 7.

In the unstable case, the evolution in ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ and ${\textrm {d}G_{P}}/{\textrm {d}t}$ is similar: Marangoni terms act to destabilize, while the pressure/capillary terms act to stabilize. The initial growths of ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ and ${\textrm {d}G_{P}}/{\textrm {d}t}$ are driven by the advective terms, and the azimuthal contribution later fuels the growth of each term. Similar to the previous case, the overall effect of the Marangoni and pressure terms are destabilizing (figure 10f). In contrast to the previous case, the contributions from the Marangoni and pressure terms dominate over the stabilizing effects of evaporation and diffusion, leading to the eventual exponential growth in the gain.

Through the gain-component analysis, we have elucidated the mechanism behind the instability. In the two cases studied, ${\textrm {d}G_{Ma}}/{\textrm {d}t}$ acts to increase the gain (figure 10a,d), while ${\textrm {d}G_{P}}/{\textrm {d}t}$ serves to stabilize the system (figure 10b,e). For the linearly stable case, diffusion has a significant stabilizing effect that contributes to the eventual decay in the gain (figure 10c). For the unstable parameter set, the contributions from diffusion and evaporation are negligible, and the destabilizing effect of Marangoni contributions overwhelms the stabilizing effects of the pressure/capillary contribution, leading to the exponential growth of the disturbance quantities (figure 10f).