2 Formulation of the problem

Consider a system of two superposed layers of immiscible liquids with different physical properties (see figure 1). The bottom layer rests on a solid substrate; the top layer is in contact with the adjacent gas phase. The temperature of the solid substrate is $T_{s}(x,y)$ , and the temperature of the gas is $T_{g}$ . Herein all variables referring to the bottom layer are marked by subscript 1, and all variables referring to the top layer are marked by subscript 2. The coordinates of the interfaces in a quiescent state are $z=H_{m}^{0}$ , with $m=1,2$ . The deformable interfaces are described by the equations $z=H_{1}(x,y,t)$ (liquid–liquid interface) and $z=H_{2}(x,y,t)$ (liquid–gas interface). The $m$ th fluid has density $\unicode[STIX]{x1D70C}_{m}$ , kinematic viscosity $\unicode[STIX]{x1D708}_{m}$ , dynamic viscosity $\unicode[STIX]{x1D702}_{m}=\unicode[STIX]{x1D70C}_{m}\unicode[STIX]{x1D708}_{m}$ , thermal diffusivity $\unicode[STIX]{x1D712}_{m}$ and heat conductivity $\unicode[STIX]{x1D705}_{m}$ . The surface tension coefficients on the lower and upper interfaces, $\unicode[STIX]{x1D70E}_{1}$ and $\unicode[STIX]{x1D70E}_{2}$ , are linear functions of temperature $T$ , i.e. $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{1}^{0}-\unicode[STIX]{x1D6FC}_{1}T$ and $\unicode[STIX]{x1D70E}_{2}=\unicode[STIX]{x1D70E}_{2}^{0}-\unicode[STIX]{x1D6FC}_{2}T$ . The problem is considered in the framework of the Boussinesq approximation, i.e. the variation of parameters is neglected (Gershuni & Zhukhovitsky Reference Gershuni and Zhukhovitsky1976).

Figure 1.

Geometric configuration of the region and coordinate axes.

Recall that there are two kinds of Marangoni instabilities in liquid layers: a short-wave instability, with critical wavenumber $k_{c}\neq 0$ , and a long-wave deformational instability, with critical wavenumber $k_{c}=0$ (see Smith Reference Smith1966; Takashima Reference Takashima1981; Golovin, Nepomnyashchy & Pismen Reference Golovin, Nepomnyashchy and Pismen1994). For the latter kind of instability, the instability threshold is proportional to the Galileo number $Ga=gH_{1}^{0}/(\unicode[STIX]{x1D708}_{1}\unicode[STIX]{x1D712}_{1})$ (Nepomnyashchy, Simanovskii & Legros Reference Nepomnyashchy, Simanovskii and Legros2012), where $H_{1}^{0}$ is the equilibrium thickness of the lower layer; hence it prevails under microgravity conditions and for thin layers. Below, we disregard the effect of gravity, which is negligible for sufficiently thin layers. In the absence of gravity, the long-wave instability has a zero threshold; hence it is always more dangerous than the short-wave one.

We assume that the characteristic spatial scale of the substrate temperature inhomogeneity is large with respect to the layers’ thicknesses. The temperature inhomogeneity generates a thermocapillary flow and surface deformations which have the same characteristic spatial scale. That circumstance allows one to apply the long-wave approach, formerly used for studying nonlinear waves created by a long-wave Marangoni instability (Danov et al. Reference Danov, Alleborn, Raszillier and Durst1998a ,Reference Danov, Paunov, Alleborn, Raszillier and Durst b ,Reference Danov, Paunov, Stoyanov, Alleborn, Raszillier and Durst c ; Pototsky et al. Reference Pototsky, Bestehorn, Merkt and Thiele2005; Nepomnyashchy & Simanovskii Reference Nepomnyashchy and Simanovskii2007, Reference Nepomnyashchy and Simanovskii2012).

The natural vertical length scale is $H_{1}^{0}$ . The choice of the horizontal length scale, $L^{\ast }$ , is arbitrary (see Nepomnyashchy & Simanovskii Reference Nepomnyashchy and Simanovskii2012). We choose

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}^{\ast }=\frac{\unicode[STIX]{x1D702}_{1}(L^{\ast })^{4}}{\unicode[STIX]{x1D70E}_{1}^{0}(H_{1}^{0})^{3}} & & \displaystyle\end{eqnarray}$$

as a time scale and

(2.2) $$\begin{eqnarray}\displaystyle p^{\ast }=\frac{\unicode[STIX]{x1D70E}_{1}^{0}H_{1}^{0}}{(L^{\ast })^{2}} & & \displaystyle\end{eqnarray}$$

as a pressure scale.

Figure 2.

Isolines of $h_{1}(X,Y)-1$ , for $M=-2$ , $\unicode[STIX]{x1D6FF}_{X}=0.15$ and $\unicode[STIX]{x1D6FF}_{Y}=0.25$ .

Figure 3.

Isolines of $h_{1}(X,Y)-1$ , for $M=-2$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.25$ .

Figure 4.

General diagram of regimes on the plane $(\unicode[STIX]{x1D6FF}_{X},\unicode[STIX]{x1D6FF}_{Y})$ for $M=-2$ : open squares, stationary pattern; black squares, oscillatory flow at small temperature modulation; open circles, oscillatory flow at large temperature modulation; open triangles, rupture.

Figure 5.

Phase trajectory on the plane $(s_{10},s_{01})$ , for $M=-2$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.1$ .

Figure 6.

Phase trajectory in the plane $(s_{10},s_{01})$ , for $M=-2$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.1$ .

Figure 7.

Snapshots of the isolines of $h_{1}(X,Y,\unicode[STIX]{x1D70F})-1$ , for $M=-2$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.1$ : (a)  $\unicode[STIX]{x1D70F}=2410\,097$ ; (b)  $\unicode[STIX]{x1D70F}=2410\,200$ ; (c)  $\unicode[STIX]{x1D70F}=2410\,385$ ; (d)  $\unicode[STIX]{x1D70F}=2410\,500$ .

Figure 8.

Decay of Fourier component $c_{10}(\unicode[STIX]{x1D70F})$ , $c_{01}(\unicode[STIX]{x1D70F})$ .

Figure 9.

Phase trajectory in the plane $(s_{10},s_{01})$ , for $M=-2$ , $\unicode[STIX]{x1D6FF}_{X}=0.25$ and $\unicode[STIX]{x1D6FF}_{Y}=0.4$ .

Figure 10.

Snapshots of the isolines of $h_{1}(X,Y,\unicode[STIX]{x1D70F})-1$ , for $M=-2$ , $\unicode[STIX]{x1D6FF}_{X}=0.25$ and $\unicode[STIX]{x1D6FF}_{Y}=0.4$ : (a)  $\unicode[STIX]{x1D70F}=2010\,099$ ; (b)  $\unicode[STIX]{x1D70F}=2010\,185$ ; (c)  $\unicode[STIX]{x1D70F}=2010\,195$ ; (d)  $\unicode[STIX]{x1D70F}=2010\,220$ ; (e)  $\unicode[STIX]{x1D70F}=2010\,366$ ; (f)  $\unicode[STIX]{x1D70F}=2010\,507$ .

Figure 11.

Oscillations of $s_{10}(\unicode[STIX]{x1D70F})$ , $s_{01}(\unicode[STIX]{x1D70F})$ , for $M=-6$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.1$ .

Figure 12.

Phase trajectory on the plane $(s_{10},s_{01})$ , for $M=-6$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.1$ .

Figure 13.

Dependence $h_{max,2}(\unicode[STIX]{x1D70F})$ , for $M=-6$ and $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0.1$ .

Figure 14.

General diagram of regimes on the plane $(\unicode[STIX]{x1D6FF}_{X},\unicode[STIX]{x1D6FF}_{Y})$ , for $M=-6$ : stars, periodic oscillatory flow; black circles, quasi-periodic oscillatory flow; open diamonds, periodic two-cell flow; black diamonds, quasi-periodic two-cell flow; asterisks, periodic three-cell flow.

Figure 15.

Phase trajectory on the plane $(s_{10},s_{01})$ , for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.15$ .

Figure 16.

Snapshots of the isolines of $h_{2}(X,Y,\unicode[STIX]{x1D70F})-h$ , for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.25$ : (a)  $\unicode[STIX]{x1D70F}=2100\,000$ ; (b)  $\unicode[STIX]{x1D70F}=2100\,117$ ; (c)  $\unicode[STIX]{x1D70F}=2100\,185$ ; (d)  $\unicode[STIX]{x1D70F}=2100\,255$ ; (e)  $\unicode[STIX]{x1D70F}=2100\,364$ ; (f)  $\unicode[STIX]{x1D70F}=2100\,399$ .

Figure 17.

Snapshots of the isolines of $h_{1}(X,Y,\unicode[STIX]{x1D70F})-1$ , for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.25$ : (a)  $\unicode[STIX]{x1D70F}=2100\,000$ ; (b)  $\unicode[STIX]{x1D70F}=2100\,117$ ; (c)  $\unicode[STIX]{x1D70F}=2100\,185$ ; (d)  $\unicode[STIX]{x1D70F}=2100\,255$ ; (e)  $\unicode[STIX]{x1D70F}=2100\,364$ ; (f)  $\unicode[STIX]{x1D70F}=2100\,399$ .

Figure 18.

Oscillations of $s_{20}(\unicode[STIX]{x1D70F})$ (solid line) and $s_{01}(\unicode[STIX]{x1D70F})$ (dashed line), for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.25$ .

Figure 19.

Oscillations of $s_{20}(\unicode[STIX]{x1D70F})$ (solid line) and $c_{20}(\unicode[STIX]{x1D70F})$ (dashed line), for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.05$ and $\unicode[STIX]{x1D6FF}_{Y}=0.35$ .

Figure 20.

Snapshots of the isolines of $h_{2}(X,Y,\unicode[STIX]{x1D70F})-h$ , $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.01$ and $\unicode[STIX]{x1D6FF}_{Y}=0.4$ : (a)  $\unicode[STIX]{x1D70F}=2020\,940$ ; (b)  $\unicode[STIX]{x1D70F}=2020\,960$ ; (c)  $\unicode[STIX]{x1D70F}=2020\,968$ ; (d)  $\unicode[STIX]{x1D70F}=2020\,971$ ; (e)  $\unicode[STIX]{x1D70F}=2020\,976$ ; (f)  $\unicode[STIX]{x1D70F}=2020\,989$ .

Figure 21.

Snapshots of the isolines of $h_{1}(X,Y,\unicode[STIX]{x1D70F})-1$ , for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.01$ and $\unicode[STIX]{x1D6FF}_{Y}=0.4$ : (a)  $\unicode[STIX]{x1D70F}=2020\,940$ ; (b)  $\unicode[STIX]{x1D70F}=2020\,960$ ; (c)  $\unicode[STIX]{x1D70F}=2020\,968$ ; (d)  $\unicode[STIX]{x1D70F}=2020\,971$ ; (e)  $\unicode[STIX]{x1D70F}=2020\,976$ ; (f)  $\unicode[STIX]{x1D70F}=2020\,989$ .

Figure 22.

Shapes of (a) the upper surface $z=h_{2}(X,Y,\unicode[STIX]{x1D70F})$ and (b) the interface $z=h_{1}(X,Y,\unicode[STIX]{x1D70F})$ , for $M=-6$ , $\unicode[STIX]{x1D6FF}_{X}=0.01$ and $\unicode[STIX]{x1D6FF}_{Y}=0.4$ .

The non-dimensional parameters of the problem are the Biot number,

(2.3) $$\begin{eqnarray}\displaystyle Bi=\frac{qH_{1}^{0}}{\unicode[STIX]{x1D705}_{2}}, & & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{1}/\unicode[STIX]{x1D702}_{2}$ , $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{1}/\unicode[STIX]{x1D705}_{2}$ , $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{2}^{0}/\unicode[STIX]{x1D70E}_{1}^{0}$ and $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{2}/\unicode[STIX]{x1D6FC}_{1}$ . Also, we define the local modified Marangoni number,

(2.4) $$\begin{eqnarray}\displaystyle M(x,y)=\frac{\unicode[STIX]{x1D6FC}_{1}(T_{s}(x,y)-T_{g})}{\unicode[STIX]{x1D70E}_{1}^{0}}\left(\frac{L^{\ast }}{H_{1}^{0}}\right)^{2}. & & \displaystyle\end{eqnarray}$$

In the framework of the long-wave approach based on the appropriate rescaling of variables, one can derive a closed system of strongly nonlinear equations for non-dimensional interfacial deformations $h_{j}=H_{j}/H_{1}^{0}$ , with $j=1,2$ (Nepomnyashchy & Simanovskii Reference Nepomnyashchy and Simanovskii2016),

(2.5a,b ) $$\begin{eqnarray}\displaystyle h_{1\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{q}_{1}=0,\quad h_{2\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{q}_{2}=0. & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D70F}$ is the rescaled non-dimensional time variable, and $\unicode[STIX]{x1D735}_{\bot }$ denotes the differentiation with respect to the rescaled non-dimensional horizontal coordinate $\boldsymbol{X}=(X,Y)$ . Each flux $\boldsymbol{q}_{j}$ , with $j=1,2$ , is a sum of two terms,

(2.6) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{j}=\boldsymbol{q}_{j}^{\unicode[STIX]{x1D70E}}+\boldsymbol{q}_{j}^{T},\quad j=1,2. & & \displaystyle\end{eqnarray}$$

The first term, $\boldsymbol{q}_{j}^{\unicode[STIX]{x1D70E}}$ , is the flux generated by the gradients of capillary pressures,

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{j}^{\unicode[STIX]{x1D70E}}=\mathop{\sum }_{k=1}^{2}f_{jk}\unicode[STIX]{x1D735}_{\bot }p_{k},\quad j=1,2 & & \displaystyle\end{eqnarray}$$

(Pototsky et al. Reference Pototsky, Bestehorn, Merkt and Thiele2004; Bandyopadhyay, Gulabani & Sharma Reference Bandyopadhyay, Gulabani and Sharma2005; Fisher & Golovin Reference Fisher and Golovin2005), where

(2.8a,b ) $$\begin{eqnarray}\displaystyle f_{11}=-{\textstyle \frac{1}{3}}h_{1}^{3},\quad f_{12}=-{\textstyle \frac{1}{2}}h_{1}^{2}(h_{2}-h_{1}), & & \displaystyle\end{eqnarray}$$

(2.9a,b ) $$\begin{eqnarray}\displaystyle f_{21}={\textstyle \frac{1}{6}}h_{1}^{3}-{\textstyle \frac{1}{2}}h_{1}^{2}h_{2},\quad f_{22}=(h_{2}-h_{1})[h_{1}^{2}({\textstyle \frac{1}{2}}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D702})+h_{1}h_{2}(-1+{\textstyle \frac{2}{3}}\unicode[STIX]{x1D702})-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D702}h_{2}^{2}], & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle p_{1}=-\unicode[STIX]{x1D6FB}_{\bot }^{2}h_{1}-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}_{\bot }^{2}h_{2}, & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle p_{2}=-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}_{\bot }^{2}h_{2}. & \displaystyle\end{eqnarray}$$

The second term, $\boldsymbol{q}_{j}^{T}$ , corresponds to the flux generated by the thermocapillary effect,

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}_{1}^{T}=-{\textstyle \frac{1}{2}}h_{1}^{2}\unicode[STIX]{x1D735}_{\bot }\{M[1+d(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D705}-Bi\,h_{1})]\}, & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}_{2}^{T}=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D705}h_{2}^{2}\unicode[STIX]{x1D735}_{\bot }(Md)+{\textstyle \frac{1}{2}}(2h_{2}-h_{1})h_{1}\unicode[STIX]{x1D735}_{\bot }\{M[-1+Bi\,h_{1}d-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D705}(1-\unicode[STIX]{x1D702})d]\}, & \displaystyle\end{eqnarray}$$

where

(2.14) $$\begin{eqnarray}\displaystyle d(\boldsymbol{X},\unicode[STIX]{x1D70F})=[\unicode[STIX]{x1D705}+Bi(1-\unicode[STIX]{x1D705})h_{1}(\boldsymbol{X},\unicode[STIX]{x1D70F})+Bi\,\unicode[STIX]{x1D705}h_{2}(\boldsymbol{X},\unicode[STIX]{x1D70F})]^{-1}. & & \displaystyle\end{eqnarray}$$

3 Marangoni waves in the absence of modulation

The stability of the mechanical equilibrium and the nonlinear flow regimes in the absence of the substrate temperature modulation ( $M(X,Y)=M=\text{const.}$ ) have been studied by Nepomnyashchy & Simanovskii (Reference Nepomnyashchy and Simanovskii2012). It has been shown that, in addition to a monotonic instability leading to the rupture of the film, the two-layer film is subject to an oscillatory instability that generates self-sustained Marangoni waves. The oscillatory instability is developed by heating from below ( $M>0$ ) if $Bi<Bi_{c}$ , and by heating from above ( $M<0$ ) if $Bi>Bi_{c}$ , where

(3.1) $$\begin{eqnarray}\displaystyle Bi_{c}=\frac{1+\unicode[STIX]{x1D6FC}(1+\unicode[STIX]{x1D702}\unicode[STIX]{x1D705}a^{2}+2\unicode[STIX]{x1D705}a)}{a}. & & \displaystyle\end{eqnarray}$$

The oscillatory instability threshold is determined by the formula

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle M_{o}=M_{2}k^{2}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle M_{2}=\frac{2(\unicode[STIX]{x1D705}+Bi+Bi\,\unicode[STIX]{x1D705}a)^{2}}{Bi\,\unicode[STIX]{x1D705}a(Bi_{c}-Bi)}\left[\frac{1}{3}+\unicode[STIX]{x1D70E}\left[a(a+1)+\frac{1}{3}(1+\unicode[STIX]{x1D702}a^{3})\right]\right], & \displaystyle\end{eqnarray}$$

where $k$ is the wavenumber of the disturbance.

If the problem is solved in a finite region $L\times L$ with periodic boundary conditions, then the set of admissible wavevectors is discrete:

(3.4a-c ) $$\begin{eqnarray}\displaystyle k_{x}=mk_{0},\quad k_{y}=nk_{0},\quad k_{0}=2\unicode[STIX]{x03C0}/L, & & \displaystyle\end{eqnarray}$$

where $m$ and $n$ are integers. Therefore,

(3.5) $$\begin{eqnarray}\displaystyle k^{2}=(m^{2}+n^{2})k_{0}^{2}. & & \displaystyle\end{eqnarray}$$

Because of the conservation of liquid volumes, the Fourier components of $h_{j}(X,Y)$ with $m=n=0$ do not change in time. For an unstable mode, the minimum value of the wavenumber compatible with periodic boundary conditions is $k_{0}$ . The corresponding modes ( $m=\pm 1$ , $n=0$ ) and ( $m=0$ , $n=\pm 1$ ) describe waves with the isolines of $h_{j}(X,Y)$ parallel to axes $X$ and  $Y$ .

Computations have been carried out by Nepomnyashchy & Simanovskii (Reference Nepomnyashchy and Simanovskii2012) for the system of Fluorinert FC70 (liquid 1) and silicone oil 10 (liquid 2). The Marangoni convection in that system was studied formerly in microgravity experiments (see e.g. Géoris et al. Reference Géoris, Hennenberg, Lebon and Legros1999) and by means of a theoretical analysis (Liu et al. Reference Liu, Zhou, Nguyen Thi and Hu2004; Nepomnyashchy et al. Reference Nepomnyashchy, Simanovskii and Legros2012). It is characterized by the following set of parameters: $\unicode[STIX]{x1D702}=3.04$ , $\unicode[STIX]{x1D705}=0.522$ , $\unicode[STIX]{x1D6FC}=2$ , $\unicode[STIX]{x1D70C}=0.482$ , $\unicode[STIX]{x1D70E}=2.6$ . For the system under consideration $\unicode[STIX]{x1D70C}<1$ , therefore, there is no Rayleigh–Taylor instability. The computations have been performed for $h=2.5$ , $Bi=10$ and $L=240$ . Equations (2.5)–(2.13) have been discretized by central differences for spatial derivatives on a grid of $80\times 80$ and solved using an explicit scheme. The Fourier components are

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle c_{mn}(\unicode[STIX]{x1D70F})=\frac{2}{L^{2}}\int _{0}^{L}\!\int _{0}^{L}h_{1}(X,Y,\unicode[STIX]{x1D70F})\cos \left[\frac{2\unicode[STIX]{x03C0}}{L}(mX+nY)\right]\text{d}X\,\text{d}Y, & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle s_{mn}(\unicode[STIX]{x1D70F})=\frac{2}{L^{2}}\int _{0}^{L}\!\int _{0}^{L}h_{1}(X,Y,\unicode[STIX]{x1D70F})\sin \left[\frac{2\unicode[STIX]{x03C0}}{L}(mX+nY)\right]\text{d}X\,\text{d}Y, & \displaystyle\end{eqnarray}$$

where $m$ and $n$ are integers. Also, variables

(3.8) $$\begin{eqnarray}\displaystyle r_{mn}(\unicode[STIX]{x1D70F})=\sqrt{c_{mn}^{2}(\unicode[STIX]{x1D70F})+s_{mn}^{2}(\unicode[STIX]{x1D70F})}, & & \displaystyle\end{eqnarray}$$

which characterize the amplitudes of the corresponding complex Fourier harmonics, and quantities

(3.9) $$\begin{eqnarray}\displaystyle h_{max,j}(\unicode[STIX]{x1D70F})=\max _{X,Y}h_{j}(X,Y,\unicode[STIX]{x1D70F}),\quad j=1,2, & & \displaystyle\end{eqnarray}$$

which describe the deformations of the surfaces, have been used.

For $Bi=10$ , the mechanical equilibrium state is unstable with respect to oscillatory disturbances by heating from above. In the region $L\times L$ with periodic boundary conditions, the linear stability theory (see (3.2) and (3.3)) predicts the instability threshold $M_{\ast }=M_{o}(k_{0})\approx -1.02$ . Below the threshold, the system tends to the mechanical equilibrium state. Above the threshold, the nonlinear development of the instability creates a wave pattern known as alternating rolls. That motion is a superposition of two periodic standing waves with wavevectors directed along the axes $X$ and $Y$ . The oscillations of waves with orthogonal wavevectors have a mutual phase shift corresponding to a quarter of a period, $T/4$ . The mean values of the Fourier components $s_{10}(\unicode[STIX]{x1D70F})$ , $c_{10}(\unicode[STIX]{x1D70F})$ , $s_{01}(\unicode[STIX]{x1D70F})$ and $c_{01}(\unicode[STIX]{x1D70F})$ are equal to zero; the amplitudes of both standing waves are equal,

(3.10) $$\begin{eqnarray}\displaystyle \max _{\unicode[STIX]{x1D70F}}r_{10}(\unicode[STIX]{x1D70F})=\max _{\unicode[STIX]{x1D70F}}r_{01}(\unicode[STIX]{x1D70F}). & & \displaystyle\end{eqnarray}$$

The quantities $h_{max,j}(\unicode[STIX]{x1D70F})$ , $j=1,2$ , oscillate with period  $T/4$ .

4 Nonlinear flow regimes for moderate $|M|$

In this section, we discuss the nonlinear regimes observed in the case of a spatially periodic temperature modulation. The most interesting dynamics is expected when the problem has a definite symmetry, which leads to the enhancement of the dimension of the central manifold near the instability point. In the present paper, we investigate the case $L_{x}=L_{y}=L$ when the computational region has a symmetry $X\rightarrow Y$ and $Y\rightarrow X$ . Correspondingly, we apply a special disturbance of the kind

(4.1) $$\begin{eqnarray}\displaystyle M(\boldsymbol{X})=M\left(1+\unicode[STIX]{x1D6FF}_{X}\sin \frac{2\unicode[STIX]{x03C0}X}{L}+\unicode[STIX]{x1D6FF}_{Y}\sin \frac{2\unicode[STIX]{x03C0}Y}{L}\right). & & \displaystyle\end{eqnarray}$$

Note that function (4.1) is symmetric with respect to axes $X=L/4$ , $X=3L/4$ , $Y=L/4$ and $Y=3L/4$ .

Let us note that the problem is characterized by a large number of non-dimensional parameters, and the investigation of their influence is beyond the scope of the present paper. The computations are carried out for a typical system of liquids appropriate for experiments, namely, the system of liquids Fluorinert FC70/silicone oil 10, and the values of parameters ( $h=2.5$ , $Bi=10$ , $L=240$ ) formerly used by Nepomnyashchy & Simanovskii (Reference Nepomnyashchy and Simanovskii2012). The same numerical method is applied.

We are interested in the case of heating from above ( $M<0$ ), where there is no film rupture in the absence of a temperature modulation. It is sufficient to take $\unicode[STIX]{x1D6FF}_{X}>0$ and $\unicode[STIX]{x1D6FF}_{Y}>0$ , because the change of the sign of $\unicode[STIX]{x1D6FF}_{X}$ or $\unicode[STIX]{x1D6FF}_{Y}$ is obtained by a translation $X\rightarrow X+L/2$ or $Y\rightarrow Y+L/2$ , correspondingly. Because of the symmetry of the problem with respect to the transformation $X\rightarrow Y$ and $Y\rightarrow X$ , it is sufficient to consider the case $\unicode[STIX]{x1D6FF}_{Y}\geqslant \unicode[STIX]{x1D6FF}_{X}$ .

4.3 Oscillatory flows for larger values of $\unicode[STIX]{x1D6FF}_{X}$ and $\unicode[STIX]{x1D6FF}_{Y}$

With the growth of $\unicode[STIX]{x1D6FF}_{Y}$ , another oscillatory instability of stationary flows takes place. The Fourier harmonics $s_{10}$ and $s_{01}$ oscillate mostly in phase (see figure 9), while $c_{10}=c_{01}=0$ . The characteristic snapshots of height isolines are shown in figure 10. The corresponding region in the plane $(\unicode[STIX]{x1D6FF}_{X},\unicode[STIX]{x1D6FF}_{Y})$ is shown in figure 4 by open circles. In the limit $\unicode[STIX]{x1D6FF}_{X}\rightarrow 0$ , this type of oscillation tends to a wave travelling in the $X$ -direction.

For even larger values of $\unicode[STIX]{x1D6FF}_{X}$ and $\unicode[STIX]{x1D6FF}_{Y}$ , the thickness of the bottom layer becomes equal to zero in the course of the development of instability. Note that a self-consistent consideration of a liquid layer rupture cannot be carried out in the framework of model (2.5)–(2.14): the contributions of intermolecular forces have to be included (Fisher & Golovin Reference Fisher and Golovin2005; Pototsky et al. Reference Pototsky, Bestehorn, Merkt and Thiele2005; Nepomnyashchy & Simanovskii Reference Nepomnyashchy and Simanovskii2006). Therefore, we stopped the simulations when

(4.2) $$\begin{eqnarray}\displaystyle h_{min,1}(\unicode[STIX]{x1D70F})=\min _{X,Y}h_{1}(X,Y,\unicode[STIX]{x1D70F}) & & \displaystyle\end{eqnarray}$$

approached zero ( $h_{min,1}(\unicode[STIX]{x1D70F})<10^{-19}$ ). The region where the model used predicts film rupture is shown in figure 4 by triangles. In our computations, we never observed the rupture of the upper layer preceding the rupture of the lower layer. Note that in the case of interfacial heating, both kinds of ruptures were observed (Nepomnyashchy & Simanovskii Reference Nepomnyashchy and Simanovskii2017).

5 Nonlinear flow regimes for large $|M|$

Far from the instability threshold, the sequence of nonlinear flow regimes is quite different. At $M=-6$ , we did not find a region where the flow is stationary. At sufficiently small $\unicode[STIX]{x1D6FF}_{X}$ and $\unicode[STIX]{x1D6FF}_{Y}$ , we obtain a time-periodic oscillatory flow with a phase shift between oscillations of $s_{10}(\unicode[STIX]{x1D70F})$ and $s_{01}(\unicode[STIX]{x1D70F})$ close to a quarter of a period, $T/4$ (see figure 11), which is quite similar to alternating rolls observed as $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}=0$ . However, now the phase trajectory in the plane $(s_{10},s_{01})$ is not symmetric even in the case $\unicode[STIX]{x1D6FF}_{X}=\unicode[STIX]{x1D6FF}_{Y}$ (see figure 12). While for ‘usual’ alternating rolls the oscillation period of quantities $h_{max,j}(\unicode[STIX]{x1D70F})$ is equal to $T/4$ , that period is now $T$ : the maxima of $h_{max,j}(\unicode[STIX]{x1D70F})$ corresponding to maxima of $|s_{10}|(\unicode[STIX]{x1D70F})$ and $|s_{01}|(\unicode[STIX]{x1D70F})$ are not equal to each other (see figure 13). Note that $c_{10}(\unicode[STIX]{x1D70F})$ and $c_{01}(\unicode[STIX]{x1D70F})$ decay; therefore, the structure keeps the symmetry with respect to two vertical axes and two horizontal axes, corresponding to the symmetry of the temperature modulation (4.1). The flow regime described above is shown in figure 14 by stars.

In a certain region of parameters $\unicode[STIX]{x1D6FF}_{X}$ and $\unicode[STIX]{x1D6FF}_{Y}$ , indicated by black circles in figure 14, quasi-periodic oscillations are observed (see figure 15). The Fourier component $s_{10}(\unicode[STIX]{x1D70F})$ oscillates with the mean frequency $\unicode[STIX]{x1D714}_{10}$ slightly larger than the mean frequency $\unicode[STIX]{x1D714}_{01}$ of $s_{01}(\unicode[STIX]{x1D70F})$ . Their nonlinear interaction creates modulation of oscillations with the characteristic time scale $2\unicode[STIX]{x03C0}/(\unicode[STIX]{x1D714}_{10}-\unicode[STIX]{x1D714}_{01})$ , which is much larger than $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{10}$ and $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}_{01}$ .

For larger values of $\unicode[STIX]{x1D6FF}_{Y}$ , a ‘two-cell’ time-periodic regime with $|s_{20}|>|s_{10}|$ is developed (see figures 16 and 17). The symmetry with respect to horizontal axes is retained, but the symmetry with respect to vertical axes is broken. The oscillations of $c_{20}(\unicode[STIX]{x1D70F})$ and $s_{20}(\unicode[STIX]{x1D70F})$ are proportional, that is, a standing wave. Owing to the nonlinearity of the equations, oscillations of $c_{20}(\unicode[STIX]{x1D70F})$ and $s_{20}(\unicode[STIX]{x1D70F})$ generate oscillations of $s_{01}$ with a doubled frequency (see figure 18). That regime is indicated in figure 14 by open diamonds.

With a further growth of $\unicode[STIX]{x1D6FF}_{Y}$ , the standing waves become quasi-periodic (see figure 19). The region of that regime is shown in figure 14 by black diamonds.

Even more short-scale time-periodic standing waves are observed for larger values of $\unicode[STIX]{x1D6FF}_{Y}$ (see figures 20–22). That ‘three-cell’ regime is quite similar to the periodic regime described above, but now the proportional Fourier components $s_{30}(\unicode[STIX]{x1D70F})$ and $c_{30}(\unicode[STIX]{x1D70F})$ prevail. Its existence region is indicated by asterisks in figure 14.