2.1. Governing equations

We suppose that the liquid is an incompressible fluid, whose dynamic viscosity $\mu $ and thermal diffusivity $\chi $ are constants. Both a Newtonian fluid and a viscoelastic fluid are considered in the following. The latter is modelled by an Oldroyd-B fluid, which can be derived from the kinetic theory for a dilute suspension of polymer molecules in a Newtonian solvent (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987).

The dimensionless numbers can be given as follows:

(2.2a)\begin{gather}Ma = \frac{{b\gamma {d^2}}}{{\chi \mu }},\quad Pr = \frac{\mu }{{\rho \chi }},\quad R = \frac{{Ma}}{{Pr}} = \frac{{\rho Ud}}{\mu },\quad S = \frac{{\rho d{{\sigma ^{\prime}}_0}}}{{{\mu ^2}}},\end{gather}

(2.2b)\begin{gather}Ca = \frac{{Ma}}{{Pr{\mkern 1mu} {\mkern 1mu} S}},\quad G = \frac{{\rho g{d^2}}}{{\mu U}},\quad Bo = \frac{{{\rho _0}g\kappa {d^2}}}{\gamma }.\end{gather}

Here, Ma is the Marangoni number, which measures the thermocapillary effect; b is the temperature gradient on the plane; Pr is the Prandtl number; R is the Reynolds number; S is the non-dimensional surface-tension number; $U = b\gamma d/\mu $ is the velocity scale; and Ca is the capillary number, which measures the magnitude of the surface deformation (Smith & Davis Reference Smith and Davis1983a). When S is large enough, $Ca \to 0$ and the surface deformation can be neglected, and CI rather than SWI is preferred. However, when $O(Ca) \ge 0.001$, the surface wave can be the preferred mode in some cases (Smith & Davis Reference Smith and Davis1983b). Also above, G is the Galileo number, which measures the gravity effect; and Bo is the dynamic Bond number, which measures the buoyancy effect. Suppose the fluid density varies linearly with the temperature,

(2.3)\begin{equation}\rho = {\rho _0}[1 - \kappa (T - {T_0})],\end{equation}

where κ is the thermal expansion coefficient. We can see that Bo is independent of the surface tension gradient. Although the fluid is assumed to be incompressible, the variation of density with temperature in the presence of gravity can lead to a buoyancy effect, which may change the flow significantly. So the Boussinesq approximation is used.

The dimensionless governing equations can be obtained as follows, which are the continuity equation, the momentum equation and the energy equation, respectively (Yan et al. Reference Yan, Hu and Chen2018); the scales of non-dimensionalization are presented in table 1:

(2.4)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0,\end{gather}

(2.5)\begin{gather}R\left( {\frac{{\partial \boldsymbol{u}}}{{\partial t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}} \right) ={-} \boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau } + (Bo{\mkern 1mu} {\mkern 1mu} T - G)\boldsymbol{\nabla }z,\end{gather}

(2.6)\begin{gather}\frac{{\partial {\kern 1pt} T}}{{\partial {\kern 1pt} t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }T = \frac{1}{{Ma}}{\nabla ^2}T.\end{gather}

Table 1.

The scales of non-dimensionalization.

Here, u = (u, v, w), p, T and $\boldsymbol{\tau }$ stand for the velocity, pressure, temperature and stress, respectively. The buoyancy effect is taken into account in the momentum equation (Hu et al. Reference Hu, He, Chen and Liu2018). When $G \gg L/d$, we have $\rho gd \gg (\mu U/{d^2})L$, which means that the variation of pressure along the x direction is much smaller than the pressure caused by gravity. Thus, the thickness expression $d = 2{\hat{\kappa }^{ - 1}}\textrm{sin}(\varphi /2)$ derived by Brochard (Reference Brochard1989) is still valid for droplet migration.

For example, the values of physical and dimensionless parameters for liquid silicon (Smith & Davis Reference Smith and Davis1983b) are stated in table 2. Here, we set the values of $\varphi $, b and L. It can be seen that the dimensionless parameters Pr, R, S, Ca, Ma and G are within the ranges we will consider in this paper. In addition, $L \gg {\hat{\kappa }^{ - 1}},d$ and $G \gg L/d$, while the variation of surface tension along the x direction is small: $bL\gamma \ll {\sigma ^{\prime}_0}$.

Table 2.

The values of physical and dimensionless parameters for liquid silicon.

The dimensionless constitutive equation for a Newtonian fluid is

(2.7)\begin{equation}\boldsymbol{\tau } = \boldsymbol{S},\end{equation}

where ${\textstyle{1 \over 2}}\boldsymbol{S} = {\textstyle{1 \over 2}}[\boldsymbol{\nabla }\boldsymbol{u} + {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}}]$ is the strain-rate tensor. For an Oldroyd-B fluid, one has

(2.8)\begin{equation}\left( {1 + \lambda \frac{{\delta {\kern 1pt} }}{{\delta {\kern 1pt} t}}} \right)\boldsymbol{\tau } = \left( {1 + \lambda \tilde{\beta }\frac{{\delta {\kern 1pt} }}{{\delta {\kern 1pt} t}}} \right)\boldsymbol{S},\end{equation}

where $\delta {\kern 1pt} /\delta t$ is the upper convected derivative with the expression

(2.9)\begin{equation}\frac{{\delta \boldsymbol{\tau }}}{{\delta t}} = \frac{{\partial \boldsymbol{\tau }}}{{\partial t}} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau } - {(\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T}}\boldsymbol{\cdot }\boldsymbol{\tau } - \boldsymbol{\tau }\boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{u}).\end{equation}

Here, $\lambda = (\mu /\tilde{G})(U/d)$ is the Weissenberg number, $\tilde{G}$ is the elastic modulus and $\tilde{\beta }$ is the ratio of solvent viscosity to total viscosity. The elasticity number $\varepsilon = \lambda /R$ is used in the following, as it depends only on the properties of the fluid and the geometry. When $\lambda = 0$ or $\tilde{\beta } = 1$, the Oldroyd-B fluid recovers the Newtonian fluid.

We choose the reference frame travelling with the droplet, while the plane moves at the speed of $\zeta $ in the negative x direction. Thus, the flux in the x direction is zero: $\int_0^1 {u\,\textrm{d}z} = 0$. The boundary conditions on the rigid plane are

(2.10a,b)\begin{equation}u{|_{z = 0}} ={-} \zeta ,\quad T{|_{z = 0}} ={-} x,\end{equation}

where $\zeta $ is the dimensionless migration velocity and a temperature gradient is imposed on the plane. We suppose that a deformable gas–liquid interface is located at $z = 1 + \xi (x,y,t)$, where $\xi (x,y,t)$ is the small displacement of the interface from its undisturbed position $z = 1$. The boundary conditions at the interface are presented as follows, which are the kinematic boundary condition, the tangential and normal components of the stress balance (Pérez-García & Carneiro Reference Pérez-García and Carneiro1991) and the continuity of heat flux, respectively:

(2.11a)\begin{gather}{\partial _t}\xi + {\boldsymbol{u}_ \bot }\boldsymbol{\cdot }\boldsymbol{\nabla }\xi = w,\end{gather}

(2.11b)\begin{gather}{\boldsymbol{t}_j}\boldsymbol{\cdot }\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{n} ={-} \boldsymbol{\nabla }T\boldsymbol{\cdot }{\boldsymbol{t}_j},\end{gather}

(2.11c)\begin{gather}- p + \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{n} ={-} C{a^{ - 1}}({\boldsymbol{\nabla }_s}\boldsymbol{\cdot }\boldsymbol{n}) - G\xi ,\end{gather}

(2.11d)\begin{gather}- \boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{n} = 0.\end{gather}

In (2.11a), ${\boldsymbol{u}_ \bot } = (u,v)$ is the two-dimensional velocity in the x–y plane; (2.11b) stands for the shear stress caused by the thermocapillary effect; and ${\boldsymbol{t}_j}$ and $\boldsymbol{n}$ are the unit tangent and unit normal vectors to the free surface, respectively. Their linearized expressions in the perturbed state are (Patne et al. Reference Patne, Agnon and Oron2021a)

(2.12a–c)\begin{equation}\boldsymbol{n} ={-} {\partial _x}\xi {\boldsymbol{e}_x} + {\boldsymbol{e}_z} - {\partial _y}\xi {\boldsymbol{e}_y},\quad {\boldsymbol{t}_{1}} = {\boldsymbol{e}_x} + {\partial _x}\xi {\boldsymbol{e}_z},\quad {\boldsymbol{t}_{2}} = {\partial _y}\xi {\boldsymbol{e}_z} + {\boldsymbol{e}_y}.\end{equation}

In (2.11c), ${\boldsymbol{\nabla }_s} = (\boldsymbol{I} - \boldsymbol{nn})\boldsymbol{\cdot }\boldsymbol{\nabla }$ is the surface gradient operator; and finally (2.11d) stands for zero heat flux.

We restrict our attention to flow that is not near the rim of the droplet. Thus, the basic flow is assumed to be parallel, while the temperature is linear in x plus a distribution in z as follows:

(2.13a,b)\begin{equation}\boldsymbol{u} = ({U_0}(z),0,0),\quad {T_0}(x,z) ={-} x + {T_b}(z).\end{equation}

The velocity and temperature distributions for the base state can be determined as follows (Hu et al. Reference Hu, Yan and Chen2019):

(2.14a)\begin{gather}{U_0}(z) ={-} \zeta \left( {1 - 3z + \frac{3}{2}{z^2}} \right) + \left( { - \frac{1}{2}z + \frac{3}{4}{z^2}} \right) + \frac{{Bo}}{4}\left( { - \frac{1}{2}z + \frac{5}{4}{z^2} - \frac{2}{3}{z^3}} \right),\end{gather}

(2.14b)\begin{gather}{T_b}(z) = {z^2}Ma\left[ {\frac{\zeta }{8}{{( - 2 + z)}^2} + \frac{z}{4}\left( {\frac{1}{3} - \frac{1}{4}z} \right) + \frac{{Bo}}{{24}}\left( {\frac{1}{2}z - \frac{5}{8}{z^2} + \frac{1}{5}{z^3}} \right)} \right].\end{gather}

It is found that the average of the velocity gradient and the surface temperature increase with Bo and $\zeta $, while the upper surface is hotter than the bottom.

When Bo = 0, the velocity distribution (2.14a) agrees with equation (8) of Pratap et al. (Reference Pratap, Moumen and Subramanian2008). The above solutions are the same for the Newtonian fluid and the Oldroyd-B fluid. However, due to the elasticity, there is a normal stress component for the latter. The total stress tensor accounting for both Newtonian and elastic contributions can be given as follows:

(2.15)\begin{equation}{\boldsymbol{\tau }_0} = \left[ {\begin{array}{*{20}{c}} {{{\bar{\tau }}_{11}}}&0&{{{\bar{\tau }}_{13}}}\\ 0&0&0\\ {{{\bar{\tau }}_{13}}}&0&0 \end{array}} \right],\end{equation}

where ${\bar{\tau }_{13}} = \textrm{D}{U_0}(z)$ and ${\bar{\tau }_{11}} = 2\lambda (1 - \tilde{\beta }){[\textrm{D}{U_0}(z)]^2}$ with $\textrm{D}{U_0}(z) = (\textrm{d}/\textrm{d}z){U_0}(z)$.

The migration velocity can be derived as follows:

(2.16)\begin{equation}\zeta = \frac{{2\textrm{cos}\phi + 1}}{6} + \frac{1}{{24}}Bo,\end{equation}

where the details can be found in Hu et al. (Reference Hu, Yan and Chen2019) and Dai et al. (Reference Dai, Khonsari, Shen, Huang and Wang2016).

2.2. Perturbed state

Now we apply an infinitesimal three-dimensional disturbance in the normal mode form to the base state:

(2.17a)\begin{gather}(\boldsymbol{u},T,P,\boldsymbol{\tau },z) = ({\boldsymbol{u}_0},{T_0},{P_0},{\boldsymbol{\tau }_0},\textrm{1}) + \boldsymbol{\delta },\end{gather}

(2.17b)\begin{gather}\boldsymbol{\delta } = (\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over v} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over w} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over T} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over P} ,\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \tau } },\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \xi } )\textrm{exp[}\sigma t + \textrm{i(}\alpha x + \beta y\textrm{)]},\end{gather}

where the subscript 0 stands for the basic flow while the variables without subscript 0 refer to the perturbation. The complex eigenvalue $\sigma = {\sigma _r} + \textrm{i}{\sigma _i}$ consists of the growth rate ${\sigma _r}$ and the frequency ${\sigma _i}$. The parameters $\alpha $ and $\beta $ are the wavenumbers in the x and y directions, respectively. The total wavenumber, wave speed and wave propagation angle are defined as $k = \sqrt {{\alpha ^2} + {\beta ^2}} $, $c = |{{\sigma_i}} |/k$ and $\theta = \textrm{ta}{\textrm{n}^{ - 1}}(\beta /\alpha )$, respectively. Owing to the symmetry, we only need to consider the case $\theta \in \textrm{[}0\mathrm{^\circ },180\mathrm{^\circ }\textrm{)}$. The perturbation equations have been presented in Hu et al. (Reference Hu, He, Chen and Liu2018) and the boundary conditions are described in the Appendix.

The Chebyshev collocation method (Schmid & Henningson Reference Schmid and Henningson2001) is used to solve the general eigenvalue problem in the form of $\boldsymbol{\mathsf{W}}\boldsymbol{g} = \sigma \boldsymbol{\mathsf{Z}}\boldsymbol{g}$, where W and Z are two matrices, and $\boldsymbol{g}$ is the eigenvector. We set N_c Gauss–Lobatto points for the governing equations in the flow region $z = [1 - \cos (j{\rm \pi} /({N_c} + 1))]/2$, $j = 1- {N_c}$, and two points for the boundary conditions at z = 0,1 (2.25). The results are sufficiently accurate when ${N_c} = 80-100$.

In order to validate our code, we compute some parameters of the two-dimensional SWI in the thermocapillary liquid layer in table 3 and the critical parameters of CI in the thermocapillary–buoyancy convection in table 4. The code for the constitutive equations of an Oldroyd-B fluid is validated by computing the critical parameters of plane Poiseuille flow for the Oldroyd-B fluid in table 5. Most of perturbation equations in the liquid layer and the channel flow are the same as those in the droplet migration, so we only need to change the basic flow and some boundary conditions in the code. Comparisons are made with the results of Smith & Davis (Reference Smith and Davis1983b), Riley & Neitzel (Reference Riley and Neitzel1998), Chan & Chen (Reference Chan and Chen2010) and Sureshkumar & Beris (Reference Sureshkumar and Beris1995). Our results are comparable to those in the previous works.

Table 3.

Some parameters of the two-dimensional SWI in thermocapillary liquid layers: (a) critical parameters of the linear flow; and (b) neutral Reynolds numbers of the return flow. Here, ‘SD’ and ‘PW’ stand for the results of Smith & Davis (Reference Smith and Davis1983b) and the present work, respectively.

Table 4.

The critical parameters of CI in the thermocapillary–buoyancy convection at $Pr = 13.9$ and $Bo = 0.142$. Here, ‘RN’, ‘CC’ and ‘PW’ stand for the results of Riley & Neitzel (Reference Riley and Neitzel1998), Chan & Chen (Reference Chan and Chen2010) and the present work, respectively.

Table 5.

The critical parameters of plane Poiseuille flow for the Oldroyd-B fluid ($\tilde{\beta } = 0.5$, $\beta = 0$). Our results are exactly the same as those in Sureshkumar & Beris (Reference Sureshkumar and Beris1995).

3.1. Newtonian fluid

3.1.1. The effect of the Galileo number

First, we study the effect of the Galileo number G on the instability. The variation of Ma_c with Pr at $\zeta = 0.49$, Ca = 0.001 and Bo = 0 is displayed in figure 2. It can be seen that Ma_c always increases with Pr. SWI is preferred at Pr < 19.1 for G = 10 and Pr < 2.05 for G = 100. The preferred mode is the downstream streamwise wave ($\theta = 0\mathrm{^\circ }$, curves a and b in figure 2), and Ma_c increases with G. When Pr becomes larger, the preferred mode changes to the oblique wave of CI (curve c in figure 2), which is independent of G. When G = 1000, we can only find CI in the flow.

Figure 2.

The variation of Ma_c with Pr at $\zeta = 0.49$, Ca = 0.001 and Bo = 0. The curves correspond to streamwise waves of SWI (a and b) and the oblique wave of CI (c).

In figure 3, we show the wavenumber, propagation angle and wave speed of the preferred modes. It can be seen that, for SWI, the wavenumbers are less than 0.6. The wave speeds increase significantly with Pr. For CI, the wavenumbers are of the order of 2, the oblique wave travels upstream $(90\mathrm{^\circ } < \theta < 180\mathrm{^\circ })$ and the wave speed is much smaller than those of SWI. We will confine ourselves to the case of G = 100 in the following study.

Figure 3.

The (a) wavenumber and (b) wave speed of the modes in figure 2.

3.1.2. The effect of the surface-tension number

In this section, we study the effect of the surface-tension number S. The variation of Ma_c with S is displayed in figure 4. The dashed portions of curves a and e are where the condition (A3) is not satisfied. When Pr = 0.01, the preferred mode includes the streamwise waves (curves a and c in figure 4) and the oblique waves (curves b and d in figure 4). Here, curve b belongs to the SWI whose propagation angle has $\theta \in (0\mathrm{^\circ },40\mathrm{^\circ })$. When S is large enough, the preferred mode changes from SWI to CI. When Pr = 1, the SWI is more prominent at S < 3.39 × 10⁵ (curve e in figure 4). When Pr = 100, only the CI is found. There are two kinds of streamwise waves, which travel upstream ($\theta = 180\mathrm{^\circ }$, curve a) and downstream ($\theta = 0\mathrm{^\circ }$, curves c and e), respectively.

Figure 4.

The variation of Ma_c with S at $\zeta = 0.49$ and Bo = 0. The curves correspond to streamwise waves (a, c and e) and oblique waves (b, d,f and g). SWI includes curves a, b, c and e, while CI includes curves d,f and g. The dashed portions of curves a, e and g are where the condition (A3) is not satisfied.

The wavenumber and wave speed of the preferred modes are displayed in figure 5. In the following, we will confine ourselves to the case Ca = R/S = 0.001, so that the condition (A3) can always be satisfied.

Figure 5.

The (a) wavenumber and (b) wave speed of the modes in figure 4.

3.1.3. The effects of the migration velocity and the Bond number

The effects of the migration velocity $\zeta $ and the Bond number Bo are discussed in this section.

(A) $Pr=0.01$

The variations of Ma_c with $\zeta $ and Bo at Pr = 0.01 are displayed in figure 6. It can be seen from (2.16) that the region of $\zeta $ depends on Bo. The Ma_c value always decreases with Bo. When Bo = 0, the variation of Ma_c with $\zeta $ is non-monotonic. The preferred modes include the streamwise waves of SWI (curves a and c in figure 6) and the oblique wave of CI (curve b in figure 6). Here, curves a and b correspond to the wave travelling upstream $(\theta > 90\mathrm{^\circ })$ while curve c corresponds to the wave travelling downstream $(\theta < 90\mathrm{^\circ })$. When Bo = 5 and 10, Ma_c always decreases with $\zeta $. The preferred modes are the downstream streamwise waves of SWI (curves d and e in figure 6).

Figure 6.

The variation of Ma_c with $\zeta $ and Bo at Pr = 0.01. The curves correspond to streamwise waves of SWI (a, c, d and e) and the oblique wave of CI (b).

The wavenumber and wave speed of the preferred modes in figure 6 are displayed in figure 7. It can be seen that the wavenumbers of SWI (curves a, c, d and e in figure 7) are much smaller than those of CI (curve b in figure 7), while the reverse is true for the wave speed.

Figure 7.

The (a) wavenumber and (b) wave speed of the modes in figure 6.

The wavenumbers of curves a, d and e in figure 7 tend to zero, which means that they belong to the long-wave instability. For these modes, Ma_N increases with k. However, the wavenumber of SWI cannot be arbitrarily small. Owing to (A3), $k \gg Ca \times \mathrm{2}{\rm \pi} \approx 0.006$. In addition, the largest wavelength is the width L, which leads to $k \ge \mathrm{2}{\rm \pi} d/L$. Therefore, the actual Ma_c will be larger than the value shown in figure 6. When Bo = 0 and $\zeta = 0$, the Ma_c of CI is 19.38. Thus, SWI is preferred when $k < 0.094$. If this value does not satisfy the previous requirements, CI will be preferred.

Here, the perturbation flow fields of the preferred modes at Pr = 0.01 are plotted in figure 8. It can be seen that there are hot spots on the surface. In figure 8(a) (respectively, figure 8b,c), the vertical up-flow on the surface at the phase $\phi = 3.5$ make $\xi $ increase, so the wave travels upstream (respectively, downstream), while the hot spot is near the wave trough (respectively, crest). There is no roll in the flow field, which is different from the case of CI. The modes in figures 8(a) and 8(b) belong to the long-wave instability.

Figure 8.

The perturbation flow field and surface deformation at Pr = 0.01, Ca = 0.001 and G = 100: (a) Bo = 0, $\zeta = 0$; (b) Bo = 10, $\zeta = 0.5$; and (c) Bo = 0, $\zeta = 0.49$. The perturbation phase is defined as $\phi = \alpha x + \beta y$.

(B) $Pr=1$

Figure 9 shows Ma_c at Pr = 1. The variation of Ma_c with Bo is non-monotonic; and Ma_c decreases with $\zeta $ except for the case at Bo = 0, $\zeta < - 0.025$. The preferred modes include the streamwise waves of SWI (curves a, c, e and f in figure 9) and the oblique waves of CI (curves b and d in figure 9). Here, curves a, b and d travel upstream, while curves c, e and f travel downstream.

Figure 9.

The variation of Ma_c with $\zeta $ and Bo at Pr = 1. The curves correspond to streamwise waves of SWI (a, c, e and f) and oblique waves of CI (b and d).

The wavenumber and wave speed of the preferred modes in figure 9 are displayed in figure 10. Here, curves a, e and f in figure 10 belong to the long-wave instability. Once again, the wavenumbers (wave speeds) of SWI are much smaller (larger) than those of the CI.

Figure 10.

The (a) wavenumber and (b) wave speed of the modes in figure 9.

(C) $Pr=100$

The variations of Ma_c with $\zeta $ and Bo at Pr = 100 are displayed in figure 11. The preferred modes include the streamwise waves of SWI (curves a, c and d in figure 11) and the oblique wave of CI (curve b in figure 11). Here, curves a and b travel upstream, while curves c and d travel downstream. SWI is dominant at Bo > 0.

Figure 11.

The variation of Ma_c with $\zeta $and Bo at Pr = 100. The curves correspond to streamwise waves of SWI (a, c and d) and the oblique wave of CI (b).

The wavenumber and wave speed of the preferred modes in figure 11 are displayed in figure 12. Here, curves a and d in figure 12 belong to the long-wave instability. For SWI, the wave speed increases with $\zeta $ when Bo > 0. The relationship between SWI and CI is similar to the cases at Pr = 0.01 and 1.

Figure 12.

The (a) wavenumber and (b) wave speed of the modes in figure 11.

When Bo > 0, there is a new type of perturbation field, which is displayed in figure 13. We can find two hot spots in the vertical direction (z = 0.5 and 0.9). There are two rolls in one cycle $(\phi \in [0,2{\rm \pi}])$, which almost coincide with the isothermals.

Figure 13.

The perturbation flow field and surface deformation at Pr = 100, $\zeta = 0.4$, Bo = 5, Ca = 0.001 and G = 100. The preferred mode is the downstream streamwise wave. Here, T_s is the perturbation temperature on the surface.

3.1.4. Energy analysis

The energy mechanism of SWI is studied in this section. The rate of change for the perturbation energy can be derived as follows (Wanschura et al. Reference Wanschura, Shevtsova, Kuhlmann and Rath1995; Hu, Peng & Zhu Reference Hu, Peng and Zhu2012):

(3.2)\begin{align} \dfrac{{\partial {E_{kin}}}}{{\partial t}} & ={-} \dfrac{1}{R}\int {p\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,{\textrm{d}^2}r} - \dfrac{1}{{2R}}\int {\textrm{(}\boldsymbol{\tau }\ {\bf :}\ \boldsymbol{S}\textrm{)}\,{\textrm{d}^3}r} + \dfrac{1}{R}\int {\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{n}\,{{\rm d}^2}r} \nonumber\\ & \quad - \int {\boldsymbol{u}\boldsymbol{\cdot }[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol{u}_0}]\,{\textrm{d}^3}r} + \int {\left( {\dfrac{{Bo}}{R}T{\kern 1pt} {\boldsymbol{e}_3}\boldsymbol{\cdot }\boldsymbol{u}} \right)\,{\textrm{d}^3}r} \nonumber\\ & = P - N + M + I + B. \end{align}

Here, P is the work done by the pressure on the surface, N is the work done by stress in the bulk of the layer, M is the work done by stress on the surface, I is the energy from the basic flow, and B is the work done by the buoyancy (Hu et al. Reference Hu, He, Chen and Liu2018); and $\int {f\,{\textrm{d}^2}r} $ and $\int {f\,{\textrm{d}^3}r} $ stand for the surface and volume integrals, respectively. For the Newtonian fluid, $\boldsymbol{\tau } = \boldsymbol{S}$, so N > 0, which stands for the viscous dissipation. The perturbation is normalized as $\int {{{|\boldsymbol{u} |}^2}\,{\textrm{d}^3}r} = 1$.

As the stress on the deformed surface has horizontal and vertical components, M can be decomposed into two terms as follows:

(3.3)\begin{equation}M = \frac{1}{R}\int {(u{\tau _{13}} + v{\tau _{23}})\,{\textrm{d}^2}r} + \frac{1}{R}\int {w{\tau _{33}}\,{\textrm{d}^2}r} = {M_h} + {M_v}.\end{equation}

In addition, we can see from (A2b) and (A2c) that shear stress can be induced by the Marangoni effect and the surface deformation. Thus, M_h is decomposed into the following two terms:

(3.4)\begin{gather}{M_h} = {M_{h,T}} + {M_{h,\xi }},\end{gather}

(3.5)\begin{gather}{M_{h,T}} = \frac{1}{R}\int {\left[ {u\left( { - \frac{{\partial T}}{{\partial x}}} \right) + v\left( { - \frac{{\partial T}}{{\partial x}}} \right)} \right]\,{\textrm{d}^2}r} ,\end{gather}

(3.6)\begin{gather}{M_{h,\xi }} = \frac{1}{R}\int {\left[ {u( - {\partial_z}{T_0} + {{\bar{\tau }}_{11}})\frac{{\partial \xi }}{{\partial x}} - u\frac{{\partial {{\bar{\tau }}_{13}}}}{{\partial z}}\xi + v( - {\partial_z}{T_0})\frac{{\partial \xi }}{{\partial y}}} \right]\,{\textrm{d}^2}r} .\end{gather}

Table 6 shows the terms of perturbation energy growth for SWI. It can be seen that M is often the main energy source, I is important in some cases at small Pr, while B is negligible; P always leads to dissipation. The maximum of $|{P/N} |$ is 12 %, while $|{P/N} |< 3\,\%$ in most cases, so P is not essential for the energy mechanism. Also, M_h is much larger than M_v. For the long-wave mode (cases 1 and 5), M_h is caused by the surface deformation, while for the mode with finite wavelength (cases 2, 3 and 4), it is caused by the Marangoni effect.

Table 6.

The terms of the perturbation energy growth for SWI: case 1, $Pr = 0.01,\, Bo = 0,\, \zeta = 0$; case 2, $Pr = 0.01,\, Bo = 0,\, \zeta = 0.49$; case 3, $Pr = 1,\, Bo = 0,\, \zeta = 0.35$; case 4, $Pr = 100,\, Bo = 5,\, \zeta = 0.4$; and case 5, $Pr = 100,\, Bo = 10,\, \zeta = 0.6$.

3.2. Oldroyd-B fluid

We study the effect of elasticity in this section. The variation of Ma_c with $\varepsilon $ at $\zeta = 0.49$, Bo = 0 and $\tilde{\beta } = 0.1$ is displayed in figure 14. The choice of $\zeta = 0.49$ is close to the complete wetting case $\zeta = 0.\textrm{5}$ (equivalently $\varphi = \textrm{0}\mathrm{^\circ }$). The downstream streamwise wave of SWI is preferred at small $\varepsilon $ (curves a, d and f in figure 14). We can see that Ma_c increases slightly with $\varepsilon $ at Bo = 0, while the reverse is true for Bo > 0. When $\varepsilon $ is large enough, the preferred mode changes to the waves of CI (curves c, e and g in figure 14), while Ma_c decreases with $\varepsilon $. For SWI, the perturbation flow field and energy mechanism of the Oldroyd-B fluid are similar to those of the Newtonian fluid.

Figure 14.

The variation of Ma_c with $\varepsilon $ at $\zeta = 0.49$, Bo = 0 and $\tilde{\beta } = 0.1$. The curves correspond to downstream streamwise waves (a, b, c, d and f) and upstream oblique waves (e and g). SWI includes curves a, b, d and f, while CI includes curves c, e and g.