Long-wave-instability-induced pattern formation in an evaporating sessile or pendent liquid layer

We investigate the nonlinear dynamics and stability of an evaporating liquid layer subject to vapor recoil, capillarity, thermocapillarity, ambient cooling, viscosity, and negative or positive gravity combined with buoyancy effects in the lubrication approximation. Using linear theory, we identify the mechanisms of finite-time rupture, independent of thermocapillarity and direction of gravity, and predict the effective growth rate of an interfacial perturbation which reveals competition among the mechanisms. A stability diagram is predicted for the onset of long-wave (LW) evaporative convection. In the two-dimensional simulation, we observe well-defined capillary ridges on both sides of the valley under positive gravity and main and secondary droplets under negative gravity, while a ridge can be trapped in a large-scale drained region in both cases. Neglecting the other non-Boussinesq effects, buoyancy does not have a significant influence on interfacial evolution and rupture time but makes contributions to the evaporation-driven convection and heat transfer. The average Nusselt number is found to increase with a stronger buoyancy effect. The flow field and interface profile jointly manifest the LW Marangoni-Rayleigh-Bénard convection under positive gravity and the LW Marangoni convection under negative gravity. In the three-dimensional simulation of moderate evaporation with a random perturbation, the rupture patterns are characterized by irregular ridge networks with distinct height scales for positive and negative gravity. A variety of interfacial and internal dynamics are displayed, depending on evaporation conditions, gravity, Marangoni effect, and ambient cooling. Reasonable agreement is found between the present results and the reported experiments and simulations. The concept of dissipative compacton also sheds light on the properties of interfacial fractalization.

Figure 1

(a) Schematic of an unstable evaporating liquid layer subject to gravity combined with buoyancy, covering a heated horizontal substrate with surface wavelength $\tilde{\lambda }\gg h_0$. Here $\mathbf{g}=-g{\mathbf{e}}_z,{\mathbf{e}}_z$ is the unit vector in the $z$ direction, an arrow with a solid line shows $g>0$, and an arrow with a dashed line shows $g<0$. The smaller dot in the layer denotes a lighter fluid particle. (b) Cases considered and study outline.

Figure 2

Solutions obtained with $(1+1)\mathrm{D}$ and $(2+1)\mathrm{D}$ versions of Eq. (1) in Ref. [54] using the parameters (in the notation therein) $\mathrm{Bi}=1,M=35.1,P=7.02,G=\frac{1}{3},S=100$, and $k_m=0.0677$. (a) Interface evolution at an interval of $100{\tau }_{\text{visc}}$ up to $t_r$ with the IC $H(x,0)=1+0.1cos\left(\frac{1}{4}{k}_{m}x\right)$ and a mesh of 601 points. (b) Surface and contour plots at $t_r$ on the domain of $l_x=l_y=185.62$ with the IC $H(x,y,0)=1+{\varepsilon }_1[cos\left(\frac{1}{2}{k}_{m}x\right)+sin\left(\frac{1}{2}{k}_{m}x\right)]cos\left(\frac{1}{2}{k}_{m}y\right)$ and an $81\times 81$ mesh.

Figure 3

Interfacial and convective stabilities of evaporating layers. (a) Plot of $s_{\text{eff}}$ as a function of $t$ and $k$ for hanging layers. The dispersion curves present an increase in $s_{\text{eff},m}$ and broadening in the band of unstable modes with $t$. Here $\mathcal{G}=-0.023,\mathcal{M}=0.29$, and $K=8.7\times {10}^{-5}$ (upper curves) are obtained from Table 2 and $\mathcal{D}=1$ is for illustration of the vapor recoil effect, feasible for a relatively small ambient pressure. The lower curves with $K=0.02$ are appropriate for a certain experimental situation (e.g., smaller $a)$. (b) Stability diagram of convective instabilities in sessile layers: ${BA}_i$, neutral curves without evaporation $\frac{1}{48}{M}_c+\frac{1}{320}{R}_c=1$ [62, 63]; ${Ob}_i$, neutral curves expressed in Eq. (34b) with a vapor recoil of $\mathcal{D}=-0.03$ (subscripts $i=0,1,2$ corresponding to $\mathcal{G}=0,-0.1,-0.5$, respectively); ${Om}_{j=1,2,3},\text{Ma}=-\frac{1}{3}\mathcal{M}\text{Ga}$ (with $\mathcal{M}=-0.01,-1,-5$, respectively); $A_0{D}_{0}Ox$ and $A_{1,2}{D}_{1,2}O$, stability; ${yBD}_{1,2}{b}_{1,2}$, Marangoni-dominated instability; ${BOD}_{i=0,1,2}$, vapor-recoil-dominated instability; and $b_{1,2}{D}_{1,2}{A}_{1,2}x$, buoyancy-dominated instability. Logarithmic scale is used for the abscissa for large values of Ga (e.g., ground conditions). All the curves are straight on a linear scale.

Figure 4

Interactions of thermocapillarity and gravity with ${\mathrm{\Sigma }}_0=100$ and $K=0.003$. (a) Plot of $k_c$ as a function of $\overline{H}$ represented by four groups of curves with $\text{Ma}=2000,1000,100,0$ and $\text{Bi}=1$. The solid lines show $\left|G\right|=10$ and $\left|\text{Gr}\right|=1$, the lower (upper) family for $G>0\left(G<0\right)$; the dashed lines show $G=0$; and the dotted line shows the nonvolatile weightless case with $K\to \infty ,E=G=0,\text{Ma}=1,{\mathrm{\Sigma }}_0=1/3$, and $\text{Pr}=2$ to recover the dotted line in Fig. 2 of Ref. [41]. (b) Dispersion curves in the case of RTI for $G=-10,\text{Bi}=\text{Gr}=0$, and $\overline{H}=1,3/4$, and $1/2$. The other parameters are taken from Table 2 in Appendix pp1.

Figure 5

Three pairs of dispersion curves with $\overline{H}\left(0\right)=1$ and $G=10$ for $E=0.03,0.06,0.09$ with $K=0.003$ or for $K=0.002,0.003,0.004$; the other parameters are taken from Table 2. Influences of (a) and (c) ambient cooling and (b) and (d) buoyancy are shown.

Figure 6

Stability diagram with $K=0.06,\text{Ra}=0$, and $E=0.011$ for $G>0$: (a) $\overline{H}=1$ and ${\mathrm{\Sigma }}_0=0$ with different $\text{Bi}$, (b) as layer thins with $\text{Bi}=1$ and ${\mathrm{\Sigma }}_0=0\left(T_e=53.636\right)$, (c) $\overline{H}=1$ and $\text{Bi}=5$ with different ${\mathrm{\Sigma }}_0$ and $k$, and (d) as layer thins with $\text{Bi}=5,k=6$, and ${\mathrm{\Sigma }}_0=1/3\left(T_e=64.545\right)$.

Figure 7

Rupture patterns of case I evaporating layers subject to RTI for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\mathcal{D}=1,\mathcal{G}=-0.2$, and different $\text{Bi}$ with a mesh of 201 points. (a) $\mathcal{M}=1$. (b) $\mathcal{M}=5$. (c) Lines represent numerical solutions of Eq. (25) and circles show the corresponding DCs with Eq. (44). The inset shows a close-up of the window of (a) illustrating the transition to (b) (see the text in Sec. 6b). Here a large value of Bi could be suitable with these considerations: (i) It could be achieved in an experiment where a liquid layer of low $k_{\text{th}}$ (e.g., HFE-7300) is surrounded by a thin gas layer with a high $k_{\text{th},g}$ (e.g., helium, similar to the experiment of VanHook et al. [8]), (ii) evaporation can give an explanation through a thermal dissipation mechanism along interface as if the gas had a large $k_{\text{th},g}$ [12], and (iii) a large Bi could be taken in the parameter study [6].

Figure 8

Streamlines and isotherms of the case I evaporating layer subject to RTI at three representative moments, calculated for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\mathcal{D}=1,\mathcal{M}=5,\mathcal{G}=-0.2$, and $\text{Bi}=1$ with a mesh of 201 points. Contours are separated by intervals of $\mathrm{\Delta }\mathrm{\Psi }$ and $\mathrm{\Delta }\mathrm{\Theta }=0.076$. (a) Initial stage of valley formation $\left(\mathrm{\Delta }\mathrm{\Psi }=0.016\right)$ at $t=4.2$. (b) Emergence of a secondary droplet $\left(\mathrm{\Delta }\mathrm{\Psi }=0.064\right)$ at $t=5.6$. (c) Formation of drained regions just before rupture $\left(\mathrm{\Delta }\mathrm{\Psi }=0.12\right)$ at $t=6.08$. The inset shows an enlarged intervening film. The initial disturbance is shown by the black dashed line. The red dots highlight intersections where streamlines $\mathrm{\Psi }=0$ cross the free surface, although some of them seem to be tangential to the interface owing to the influences of negative gravity on interfacial velocity and its slope under weak mass loss.

Figure 9

Evolution of the case III evaporating layer subject to RTI for $K=1,\mathcal{E}=0.1,\mathcal{D}=0.1,\mathcal{M}=5,\mathcal{G}=-0.1$, and $\text{Bi}=1$. The side length of the periodic domain is $4\pi$. A uniform $151\times 151$ mesh is employed. The number of semidiscrete ODEs to be integrated at each step is $22500$. Successive snapshots of the interface contour with the minimum and maximum elevations $(H_{\text{min}}$ and $H_{\text{max}})$ are shown at (a) $t=1.5$ (0.9039 and 0.9885), (b) $t=4$ (0.4312 and 1.3231), and (c) $t_r=5.781(5.9019\times {10}^{-5}$ and 2.5040). Red arrows denote the orientation of the breakup of annular dry patches. (d) Surface plot at $t_r$. Also shown is the evolution of representative profiles for (e) $x=2.1$ and 7.2 and (f) $x=10$.

Figure 10

Rupture patterns of the case I evaporating layers lying on substrates for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\mathcal{D}=-1,\mathcal{G}=-0.2$, and different $\text{Bi}$ with a mesh of 201 points. (a) $\mathcal{M}=-1$. (b) $\mathcal{M}=-5$. (c) $\mathcal{M}=-5$, with lines representing numerical solution of Eq. (25) and circles showing the DC with Eq. (44). The inset shows a close-up of the window illustrating the next rupture (see the text in Sec. 6b).

Figure 11

Streamlines and isotherms of the case I evaporating layer for $G>0$ at three representative moments, calculated for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\mathcal{D}=-1,\mathcal{M}=-5,\mathcal{G}=-0.2$, and $\text{Bi}=1$ with a mesh of 201 points. Contours are separated by intervals of $\mathrm{\Delta }\mathrm{\Psi }$ and $\mathrm{\Delta }\mathrm{\Theta }$. (a) Initial stage of valley formation $(\mathrm{\Delta }\mathrm{\Psi }=0.1$ and $\mathrm{\Delta }\mathrm{\Theta }=0.064)$ at $t=10$. (b) Emergence of capillary ridges $(\mathrm{\Delta }\mathrm{\Psi }=0.12$ and $\mathrm{\Delta }\mathrm{\Theta }=0.062)$ at $t=11.35$. (c) Formation of the drained region just before rupture $(\mathrm{\Delta }\mathrm{\Psi }=0.096$ and $\mathrm{\Delta }\mathrm{\Theta }=0.062)$ at $t=11.42$. At the same level liquid is warmer under a convex surface than that under a concave one. The smaller dot denotes lighter fluid particle, illustrating the destabilizing effect of buoyancy. The hexagon wavelength [24] ${\lambda }_6=\frac{3}{2}{l}_6\approx 3.9$ with the side length $l_6=\frac{2}{\sqrt{3}}\times 2.27\approx 2.62$. The insets show the enlarged middle ridges.

Figure 12

Evolution of the case III evaporating layer lying on a substrate for $K=1,\mathcal{E}=0.1,\mathcal{D}=-0.1,\mathcal{M}=-5,\mathcal{G}=-0.1$, and $\text{Bi}=1$ with a $171\times 171$ mesh. The side length of the periodic domain is $6\pi$. The amount of coupled ODEs to be solved at each step is $28900$. Successive snapshots of the interface contour with the minimum and maximum elevations $(H_{\text{min}}$ and $H_{\text{max}})$ are shown at (a) $t=1$ (0.9559 and 0.9766), (b) $t=7$ (0.7330 and 0.7586), (c) $t=12.5$ (0.4688 and 0.5350), and (d) $t_r=18.741(5.5066\times {10}^{-5}$ and 0.2526). (e) Surface plot at $t_r$. (f) Evolution of representative profiles for $x=1.7$ and 16.2.

Figure 13

Streamlines and pressure contours of case I at the emergence of a valley, calculated for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\left|\mathcal{D}\right|=1,\left|\mathcal{M}\right|=5,\mathcal{G}=-0.2$, and $\text{Bi}=1$ with a mesh of 201 points. The $\widehat{P}$ contours are dashed with an interval of $\mathrm{\Delta }\widehat{P}$. The parameters on the left are $G>0,\mathrm{\Delta }\mathrm{\Psi }=0.1$, and $\mathrm{\Delta }\widehat{P}=0.2$; on the right $G<0,\mathrm{\Delta }\mathrm{\Psi }=0.016$, and $\mathrm{\Delta }\widehat{P}=0.25$. A smaller dot along a streamline denotes a lighter fluid particle due to buoyancy. The distinct differences when streamlines cross the interfaces can be understood from a mechanical point of view, as in the insets, which are qualitative sketches illustrating the composition of forces on the interfacial particles.

Figure 14

The $H_{\text{min}}$ evolution of the case I evaporating layers with different thermocapillarity for $\lambda =2\sqrt{2}\pi ,K=0.1,\mathcal{E}=0.05,\left|\mathcal{D}\right|=1,\mathcal{G}=0,\text{Bi}=0$ or 0.5, and (a) and (c) $G<0$ and (b) and (d) $G>0$.

Figure 15

Influences of buoyancy on the case III evaporating layers lying on a substrate with $K=1,\mathcal{E}=0.1,\mathcal{D}=-0.1,\mathcal{M}=-5$, and $\text{Bi}=1$. (a) Variation of $\langle \text{Nu}\rangle$ at $t_r$ with $\mathcal{G}$. Squares are for $l=8\pi$ with a $181\times 181$ mesh, circles are for $l=6\pi$ with a $171\times 171$ mesh, and the markers are labeled with values of $t_r$. (b) Line integral convolution plot for the case with $\mathcal{G}=-0.2$ and $l=8\pi$ at $t_r$.

Figure 16

Results for $G<0,K=0.1$, and $\mathcal{D}=1$ with a mesh of 201 points. Case I evolution with $\mathcal{E}=0.05$ is shown for (a) $\mathcal{M}=1$ and $\mathcal{G}=-0.2$ or $-0.3$ at $t=1,2,3,4,5,6,6.5,6.8,7,t_r$ and (b) $\mathcal{M}=5$ and $\mathcal{G}=-0.2$ at $t=1,2,3,4,4.6,5,5.2,5.3,5.4,t_r$. (c) Rupture patterns of cases I and II with $\mathcal{G}=0$.

Figure 17

Results for $G>0,K=0.1$, and $\mathcal{D}=-1$ with a mesh of 201 points. Case I evolution with $\mathcal{E}=0.05$ is shown for (a) $\mathcal{M}=-1$ and $\mathcal{G}=-0.2$ or $-0.4$, corresponding to $0\le t\le 11$ with $\mathrm{\Delta }t=1$ and instants of 11.4 and $t_r$, and (b) $\mathcal{M}=-5$ and $\mathcal{G}=-0.2$ at $t=1,2,3,4,5,6,7,8,9,9.4,9.7,9.8,9.92,t_r$. (c) Rupture patterns of cases I and II with $\mathcal{G}=0$.

Figure 18

Hierarchical DCs indicate the self-affine structure of Eq. (44). (a)–(c) Evolution of $H(x,t)$ illustrates droplet formation for $K=1,\left|\mathcal{M}\right|=5,\left|\mathcal{D}\right|=0.1$, and $\text{Bi}=0.1$ using the IC ((27a) with ${\varepsilon }_0=0.1$. Lines represent numerical results of Eq. (42) and circles show the corresponding DCs with Eq. (44). (b) and (c) show close-ups of the imminent DCs in (a). (d) Plot of $L_n$ and ${\mathcal{H}}_m$ versus $2l_n$.

Figure 19

Generalized basic-state behavior for (a) $\overline{H}$, (b) $\overline{J}$, and (c) ${\overline{\mathrm{\Theta }}}_w-{\overline{\mathrm{\Theta }}}_I$. In (b) and (c) time is scaled with $T_e$. In (b) the horizontal dotted line corresponding to $K^{-1}$ denotes the value of ${\overline{J}}_{1,2}$ at $T=T_e$.

Figure 20

Potential energy for $K=0.1,\left|\mathcal{G}\right|=0.2,\mathrm{Bi}=0.5,\left|\mathcal{D}\right|=0.2$ (dashed line) or 1 (solid line), and different values of $\mathcal{M}$ showing global and local minima: (a) $V\left(H\right)$ for condensing layers hanging from a cooled substrate $\left(\mathcal{D},\mathcal{G}>0\right)$ with $I_0=-1$ and (b) $-V\left(H\right)$ for evaporating layers lying on a heated substrate with $I_0=-3$. The horizontal dotted lines indicate the potential wells for several curves.