2 Mathematical model

2.1 Governing equations

Consider a two-layer system of viscous thermally conductive immiscible fluids (figure 1). Here $G_{1}$ is the layer occupied by the gas, $G_{2}$ is the liquid layer, $f(t,x)$ is the interface between them. Below and above, the system is bounded by solid walls $f_{2}(t,x)$ and $f_{0}(t,x)$ , respectively; the side boundaries are solid planes $x=0$ and $x=L$ . The gravity vector $\boldsymbol{g}=(0,-g)$ is parallel to the $y$ axis in the selected coordinate system and is directed downward. One or more heaters are located on the substrate $f_{2}(t,x)$ . When heating, they excite convective motion in the layers and deformation of the interface. In the general case the $j$ th thermal element occupies non-empty area $Q_{j}$ that is part of the lower boundary defined by the equation $y=f_{2}(t,x)$ and $q_{j}(t)$ is the heater’s temperature. Areas $Q_{j}$ and values $q_{j}$ are determined according to the requirements of specific problems. In the case of steady heating $q_{j}(t)\equiv \text{const}$ .

Each layer $G_{i}$ ( $i=1,2$ ) is characterized by its thickness $\,h_{i}$ , density $\,\unicode[STIX]{x1D70C}_{i}$ and kinematic viscosity $\,\unicode[STIX]{x1D708}_{i}$ . Here, the viscosity $\,\unicode[STIX]{x1D708}_{i}$ in each layer is considered to be a constant value, but the change in density $\,\unicode[STIX]{x1D70C}_{i}$ due to the temperature heterogeneity in the layers is taken into account. We choose $x_{0}=h_{1}$ and $v_{0}=\unicode[STIX]{x1D708}_{1}/h_{1}$ as characteristic scales of length and velocity, respectively, and assume them as well as the pressure scale $p_{0}=\unicode[STIX]{x1D70C}_{1}v_{0}^{2}$ to be identical for both layers of the system. Here, $h_{1}$ is the layer thickness of $G_{1}$ at rest.

The fluid flow and heat transfer in each of the layers $G_{i}$ is described by a system of Navier–Stokes equations in the Oberbeck–Boussinesq approximation, as well as the heat conductivity equation expressed in the variables $\,\unicode[STIX]{x1D713}$ , $\,\unicode[STIX]{x1D714}$ , $\,\unicode[STIX]{x1D703}\,$ (streamfunction, vorticity and temperature),

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x2202}t}\,+\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\bigg(\unicode[STIX]{x1D714}_{i}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}y}\bigg)-\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\bigg(\unicode[STIX]{x1D714}_{i}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}x}\bigg)=\,\frac{1}{Re_{i}}\,\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{i}+\,\frac{Gr_{i}}{Re_{i}^{2}}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{i}}{\unicode[STIX]{x2202}x}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{i}+\unicode[STIX]{x1D714}_{i}=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{i}}{\unicode[STIX]{x2202}t}\,+\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\bigg(\unicode[STIX]{x1D703}_{i}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}y}\bigg)-\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\bigg(\unicode[STIX]{x1D703}_{i}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}x}\bigg)=\frac{1}{Pr_{i}Re_{i}}\,\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{i}. & \displaystyle\end{eqnarray}$$

The required functions are associated with the original unknown characteristics through the following relations:

(2.4a-d ) $$\begin{eqnarray}u_{i}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}y},\quad v_{i}=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}x},\quad \unicode[STIX]{x1D714}_{i}=\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}y},\quad \unicode[STIX]{x1D703}_{i}=\frac{T_{i}-T_{0}}{\unicode[STIX]{x1D6FF}T},\end{eqnarray}$$

where $T_{0}$ is the characteristic temperature, $\unicode[STIX]{x1D6FF}T$ is the characteristic temperature difference for the whole system. Dimensionless Reynolds, Grashoff and Prandtl numbers in the $i$ th medium are introduced in the usual way,

(2.5a-c ) $$\begin{eqnarray}Re_{i}=\,\frac{x_{0}v_{0}}{\unicode[STIX]{x1D708}_{i}},\quad Gr_{i}=\,\frac{\unicode[STIX]{x1D6FD}_{i}\unicode[STIX]{x1D6FF}Tgx_{0}^{3}}{\unicode[STIX]{x1D708}_{i}^{2}},\quad Pr_{i}=\,\frac{\unicode[STIX]{x1D708}_{i}}{\unicode[STIX]{x1D712}_{i}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{i}$ , $\unicode[STIX]{x1D712}_{i}$ are the thermal expansion and thermal diffusivity coefficients of the $i$ th fluid, respectively.

The liquid layer $G_{2}$ is characterized by its surface tension coefficient $\,\unicode[STIX]{x1D70E}(T)$ at the boundary of this liquid with air or its own steam. It is assumed that the surface tension coefficient of the liquid is a linear function of temperature decreasing with the growth of  $T$ ,

(2.6a-d ) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(T)=\unicode[STIX]{x1D70E}_{0}(1-\unicode[STIX]{x1D70E}_{T}(T-T_{0})),\quad \unicode[STIX]{x1D70E}_{0}=\unicode[STIX]{x1D70E}(T_{0}),\quad \unicode[STIX]{x1D70E}_{T}=-\frac{1}{\unicode[STIX]{x1D70E}_{0}}\frac{\text{d}\unicode[STIX]{x1D70E}}{\text{d}T}\mid _{T=T_{0}},\quad \unicode[STIX]{x1D70E}_{0},\unicode[STIX]{x1D70E}_{T}>0.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D70E}_{T}$ is the temperature coefficient characterizing the rate of change of the surface tension at temperature deviations with respect to the characteristic value $T_{0}$ .

2.2 Interface. Functional dependencies

Consider the interface $y=f(t,x)$ between the layers $G_{1}$ and $G_{2}$ and define the unit tangent and normal vectors on this interface

(2.7a,b ) $$\begin{eqnarray}\boldsymbol{s}=\left\{\frac{1}{\sqrt{1+f_{x}^{2}}},\frac{f_{x}}{\sqrt{1+f_{x}^{2}}}\right\},\quad \boldsymbol{n}=\left\{-\frac{f_{x}}{\sqrt{1+f_{x}^{2}}},\frac{1}{\sqrt{1+f_{x}^{2}}}\right\}.\end{eqnarray}$$

Note that for a liquid $G_{2}$ , where $\boldsymbol{n}$ is the unit vector of the outer normal at the interface between the layers, while for a fluid $G_{1}$ for the same surface $y=f(t,x)$ this normal is directed into the gas layer.

2.2.1 Functional dependencies for temperature

In this paper, we do not consider the evaporation (condensation) processes at the internal interface, and we also neglect the energy consumption associated with its deformation (since they are very insignificant for the considered layer thicknesses). Therefore, the temperature conditions at the interface of the layers $f(t,x)$ can be written as continuity conditions of the temperature and heat fluxes

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{1}=\unicode[STIX]{x1D703}_{2},\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{2}}{\unicode[STIX]{x2202}n}-\bar{\unicode[STIX]{x1D705}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{1}}{\unicode[STIX]{x2202}n}=0,\end{eqnarray}$$

where $\bar{\unicode[STIX]{x1D705}}=\unicode[STIX]{x1D705}_{1}/\unicode[STIX]{x1D705}_{2}$ is the ratio of the fluid thermal conductivities.

2.2.2 Functional dependencies for the streamfunction

From the equality of the tangential velocities of the fluids $G_{1}$ and $G_{2}$ at their common interface $f(t,x)$ and volume preserving conditions for each medium, it follows that

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{1}=\unicode[STIX]{x1D713}_{2};\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}n}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}n}=0.\end{eqnarray}$$

2.2.3 Functional dependencies for vorticity

To derive these conditions the kinematic and dynamic conditions at the interface were used (Napolitano Reference Napolitano1979; Andreev et al. Reference Andreev, Gaponenko, Goncharova and Pukhnachov2012). The dynamic condition is based on the continuity retention of normal and tangential stresses on a free surface, which is relevant if the medium itself retains continuity. If evaporation (condensation) processes are neglected, the dynamic condition at the interface $f(t,x)$ can be presented as

(2.10) $$\begin{eqnarray}(\unicode[STIX]{x1D63F}_{2}-\unicode[STIX]{x1D63F}_{1})\boldsymbol{n}=\frac{\unicode[STIX]{x1D70E}(T)}{R}\boldsymbol{n}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}(T),\end{eqnarray}$$

where $\unicode[STIX]{x1D63F}_{1}$ , $\unicode[STIX]{x1D63F}_{2}$ are the stress tensors in the layers $G_{1}$ and $G_{2}$ , respectively, $\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}(T)$ is the surface tension gradient along the interface, $R$ is the interface curvature radius $f(t,x)$ . Thus, we have

(2.11) $$\begin{eqnarray}\frac{1}{R}=\frac{f_{xx}}{\sqrt{(1+f_{x}^{2})^{3}}}.\end{eqnarray}$$

Two consequences can be obtained from the dynamic condition. Carrying out scalar multiplication of (2.10) by the vector $\boldsymbol{s}$ we obtain

(2.12) $$\begin{eqnarray}\boldsymbol{s}(\unicode[STIX]{x1D63F}_{2}-\unicode[STIX]{x1D63F}_{1})\boldsymbol{n}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}(T)}{\unicode[STIX]{x2202}s}\end{eqnarray}$$

or

(2.13) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D708}_{2}\left[\frac{1-f_{x}^{2}}{2}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}y^{2}}-\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}x^{2}}\right)-2f_{x}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}x}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad \displaystyle -\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D708}_{1}\left[\frac{1-f_{x}^{2}}{2}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}y^{2}}-\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}x^{2}}\right)-2f_{x}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}_{1}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}x}\right]=\frac{1+f_{x}^{2}}{2}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}(T)}{\unicode[STIX]{x2202}s}.\end{eqnarray}$$

Using the kinematic condition for each of the fluids on their common interface, condition (2.13) can be reduced to

(2.14a ) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{2}-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D714}_{1}=F_{1}(t,x),\end{eqnarray}$$

where

(2.14b ) $$\begin{eqnarray}F_{1}=Mn\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}s}+2(1-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}})\left(\frac{\unicode[STIX]{x2202}v_{n}}{\unicode[STIX]{x2202}s}+\frac{v_{s}}{R}\right).\end{eqnarray}$$

Here, $v_{s}$ , $v_{n}$ are the tangent and normal components of the velocity for the points lying at the interface $f(x,t)$ , $Mn=(\unicode[STIX]{x1D70E}_{T}\unicode[STIX]{x1D6FF}T)/(\unicode[STIX]{x1D70C}_{2}v_{0}\unicode[STIX]{x1D708}_{2})$ is the Marangoni number, $\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{1}/\unicode[STIX]{x1D70C}_{2}$ , $\bar{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D708}_{1}/\unicode[STIX]{x1D708}_{2}$ are the ratios of densities and viscosities, respectively. Carrying out scalar multiplication of (2.10) by the vector $\boldsymbol{n}$ , we obtain

(2.15) $$\begin{eqnarray}\boldsymbol{n}(\unicode[STIX]{x1D63F}_{2}-\unicode[STIX]{x1D63F}_{1})\boldsymbol{n}=\frac{\unicode[STIX]{x1D70E}(T)}{R}\end{eqnarray}$$

or transforming the expression to a more convenient form, we obtain

(2.16) $$\begin{eqnarray}-(p_{2}-p_{1})=-2\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D708}_{2}\left[\bigg(\frac{\unicode[STIX]{x2202}v_{n}}{\unicode[STIX]{x2202}n}\bigg)_{2}-\frac{f_{x}v_{s}}{R}\right]+2\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D708}_{1}\left[\bigg(\frac{\unicode[STIX]{x2202}v_{n}}{\unicode[STIX]{x2202}n}\bigg)_{1}-\frac{f_{x}v_{s}}{R}\right]+\frac{\unicode[STIX]{x1D70E}(T)}{R}.\end{eqnarray}$$

Next, we consider the motion equations in each of the $G_{i}$ layers in the Oberbeck–Boussinesq approximation written in the ‘velocity–pressure’ variables,

(2.17) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{v}_{i}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}_{i}\unicode[STIX]{x1D735})\boldsymbol{v}_{i}\,=-\frac{1}{\unicode[STIX]{x1D70C}_{i}}\,\unicode[STIX]{x1D735}p_{i}+\unicode[STIX]{x0394}\boldsymbol{v}_{i}-\boldsymbol{g}\unicode[STIX]{x1D6FD}_{i}\unicode[STIX]{x1D703}_{i}.\end{eqnarray}$$

A solution should be found so that $v_{s}$ at the interface $f(t,x)$ meets conditions (2.9) (the equality of the tangential velocities of the points lying at the interface). This is a requirement of the physical model of the concerned problem. Here, the dynamic condition has its own requirement, namely, the pressure jump at this interface must satisfy balance relations (2.10).

While carrying out scalar multiplication of vector equation (2.17) for $i=1$ by the vector $-\unicode[STIX]{x1D70C}_{1}\boldsymbol{s}(t,x)$ , and the equation for $i=2$ by the vector $\unicode[STIX]{x1D70C}_{2}\boldsymbol{s}(t,x)$ , and summing them up, we obtain the following relation:

(2.18) $$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})\,\frac{\unicode[STIX]{x2202}v_{s}}{\unicode[STIX]{x2202}t}\,+(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})v_{s}\,\frac{\unicode[STIX]{x2202}v_{s}}{\unicode[STIX]{x2202}s}=-\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(p_{2}-p_{1})-\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}(\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D708}_{2}\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D708}_{1}\unicode[STIX]{x1D714}_{1})\nonumber\\ \displaystyle & & \displaystyle \qquad -\,(\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D6FD}_{2}\unicode[STIX]{x1D703}_{2}-\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D6FD}_{1}\unicode[STIX]{x1D703}_{1})(\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{s})+(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})\,\frac{f_{x}v_{n}^{2}}{R}+v_{n}(\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D714}_{2}-\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D714}_{1}),\end{eqnarray}$$

in which $p_{i}$ is the deviation of the actual pressure from the hydrostatic one. Excluding $\unicode[STIX]{x2202}(p_{2}-p_{1})/\unicode[STIX]{x2202}s$ by means of (2.16) and given the equality of the temperatures at the interface of the layers, after reducing to the dimensionless form, we obtain

(2.19a ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{2}}{\unicode[STIX]{x2202}n}-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{1}}{\unicode[STIX]{x2202}n}=F_{2}(t,x),\end{eqnarray}$$

where

(2.19b ) $$\begin{eqnarray}\displaystyle F_{2}(t,x) & = & \displaystyle -2\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}v_{n}}{\unicode[STIX]{x2202}n}\right)_{2}-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{\unicode[STIX]{x2202}v_{n}}{\unicode[STIX]{x2202}n}\right)_{1}\right]+2\left[\left(1-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\right)\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{f_{x}v_{s}}{R}\right)\right]\nonumber\\ \displaystyle & & \displaystyle +\,Ca^{-1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left[\frac{1}{R}\left(1-Mn\,Ca\unicode[STIX]{x1D703}\right)\right]+\left[\left(\frac{Gr_{2}}{Re_{2}}\,-\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\,\frac{Gr_{1}}{Re_{1}}\right)\unicode[STIX]{x1D703}+\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}\,\frac{G}{Re_{1}}\right]\frac{f_{x}}{\sqrt{1+f_{x}^{2}}}\nonumber\\ \displaystyle & & \displaystyle \displaystyle +\,Re_{2}\left[\!(\bar{\unicode[STIX]{x1D70C}}-1)\frac{\unicode[STIX]{x2202}v_{s}}{\unicode[STIX]{x2202}t}+(\bar{\unicode[STIX]{x1D70C}}-1)v_{s}\,\frac{\unicode[STIX]{x2202}v_{s}}{\unicode[STIX]{x2202}s}+(1-\bar{\unicode[STIX]{x1D70C}})\,\frac{f_{x}v_{n}^{2}}{R}+v_{n}(\unicode[STIX]{x1D714}_{2}-\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D714}_{1})\!\right],\end{eqnarray}$$

(2.19c,d ) $$\begin{eqnarray}Ca=\frac{\unicode[STIX]{x1D70C}_{2}v_{0}\unicode[STIX]{x1D708}_{2}}{\unicode[STIX]{x1D70E}_{0}},\quad G=\frac{gx_{0}}{v_{0}^{2}}.\end{eqnarray}$$

Here, $Ca$ is the capillary number, $G$ is the Galileo number. Note that conditions (2.14a ) and (2.19a ) include the parameter $\bar{\unicode[STIX]{x1D70C}}\bar{\unicode[STIX]{x1D708}}$ , which plays an important role in the formation of convective flows in the layers (Goncharova et al. Reference Goncharova, Kabov and Pukhnachov2012), especially in the cases where there is no thermal load, and convective flows are caused by other factors (Kushnir et al. Reference Kushnir, Segal, Ullmann and Brauner2014).

The position of the interface $f(t,x)$ is determined from the kinematic condition on the discontinuity surface, for example, for the lower layer. In the case under consideration this condition has the form

(2.20) $$\begin{eqnarray}f_{t}+\sqrt{1+f_{x}^{2}}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}}{\unicode[STIX]{x2202}s}=0.\end{eqnarray}$$

2.3 Boundary conditions at the external boundaries of the system

Boundary conditions of the first kind for the temperature function are set at the lateral and upper boundaries,

(2.21a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{i\,}|_{\,x=0}=0,\quad \unicode[STIX]{x1D703}_{i\,}|_{\,x=L}=0,\quad \unicode[STIX]{x1D703}_{1\,}|_{\,y=f_{0}(x)}=0.\end{eqnarray}$$

Boundary conditions for the temperature at the lower boundary, written in general form taking into account the presence of local heaters, are as follows:

(2.22a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{2\,}|_{\,y=f_{2}(x),\,x\notin Q_{j}^{s}}=0,\quad \unicode[STIX]{x1D703}_{2\,}|_{\,y=f_{2}(x),\,x\in Q_{j}^{s}}=q_{j}^{s}(t),\end{eqnarray}$$

where $Q_{j}^{s}$ is the area occupied by the $j$ th heater, $q_{j}^{s}$ is the temperature of the $j$ th heater. The exact form of conditions (2.22) specifying the areas $Q_{j}^{s}$ and values $q_{j}^{s}$ is given in the corresponding sections along with the results of the numerical calculations.

Boundary conditions for the current function on solid surfaces $x=0$ , $x=L$ , $y=f_{0}(x)$ , $y=f_{2}(x)$ correspond to the no-slip condition for viscous fluids

(2.23) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{i}=0,\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}}{\unicode[STIX]{x2202}n}=0.\end{eqnarray}$$

Here, we have two conditions for the function $\unicode[STIX]{x1D713}$ and no conditions for vorticity $\unicode[STIX]{x1D714}$ . To find the value of $\unicode[STIX]{x1D714}$ on a solid wall we use Tom’s condition (Roache Reference Roache1976). It should be emphasized that in this model the closing relations (i.e. boundary conditions) are not only derived from the conservation laws, but are themselves written in the form of conservation laws in accordance with the assumptions and constraints which form the basis of the model. In conclusion, we note that this model can be adapted to solve a wide range of problems, such as calculations for media with variable viscosity (Ovcharova Reference Ovcharova2006), for conditions that take into account the evaporation processes from one layer to another across the interface (Kuznetsov et al. Reference Kuznetsov, Bartashevich and Kabov2012; Bekezhanova & Goncharova Reference Bekezhanova and Goncharova2016) or changes in the internal energy of the interface associated with the Marangoni forces (Bekezhanova & Kabov Reference Bekezhanova and Kabov2016), etc. The numerical method for solving problems similar to the present task is described in detail in Ovcharova (Reference Ovcharova2017).

In the present work we performed a series of calculations for one and two heaters arranged at the lower boundary of the liquid–gas system and working in different modes and with various heating intensities. Since the configuration of the entire thermal stabilization system is very complicated, at the initial stage, the issues of heater arrangement and their operation in a separate working section of the fluid path are studied. In the onboard equipment there are both thermal elements operating in the stationary mode, and periodically switched on/off elements; this fact explains and determines the choice of the regimes under study. In the commutated operation regime, some heaters can be turned off completely and then be switched on, or they can be switched to operating regimes with large heat fluxes by activating some system units (antenna correction, telemetry data transmission, etc.); then they can be switched to the initial regime with constant temperature (or magnitude of the heat flux). Thus, the configurations in question are the model ones.

As mentioned above, the scaling of the whole system should be unified, i.e. the characteristic values of length $x_{0}$ , velocity $v_{0}$ and pressure $p_{0}$ should be the same for each layer of the system. Given the chosen method of non-dimensionalization, where the first fluid layer of the system is taken as a reference, the Reynolds numbers $Re_{1}=1$ , and $Re_{2}$ , corresponding to the first and second medium, respectively, are determined unambiguously. Nitrogen (gas) and ethanol (liquid) are selected as working media. The thicknesses of the layers $G_{1}$ and $G_{2}$ are assumed to be the same, $h_{1}=h_{2}=5\times 10^{-3}$  m. The thermophysical properties of the system are given below in the form of $\{(1),(2)\}$ , where the value $(i)$ corresponds to the $i$ th fluid occupying the region $G_{i}$ , while the value of temperature coefficient of the surface tension is specified only for ethanol: $\unicode[STIX]{x1D70C}=\{1.25$ , $0.79\times 10^{3}\}~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D708}=\{0.15\times 10^{-4}$ , $0.15\times 10^{-5}\}~\text{m}^{2}~\text{s}^{-1}$ , $\unicode[STIX]{x1D6FD}=\{3.67\times 10^{-3}$ , $0.108\times 10^{-2}\}~\text{K}^{-1}$ , $\unicode[STIX]{x1D712}=\{0.3\times 10^{-4}$ , $0.89\times 10^{-7}\}~\text{m}^{2}~\text{s}^{-1}$ , $\unicode[STIX]{x1D705}=\{0.02717$ , $0.1672\}~\text{W}~(\text{m}^{-1}~\text{K}^{-1})$ , $\unicode[STIX]{x1D70E}_{T}=0.8\times 10^{-4}~\text{N}~(\text{m}^{-1}~\text{K}^{-1})$ . The characteristic temperature drop is taken as $\unicode[STIX]{x1D6FF}T=10$  K, the length of the test section is $L=0.2$  m. In both cases, the size of the heaters is the same $6h_{1}=0.03$  m and they are symmetrically located relative to the middle of the test section. Calculations are carried out with the following values of the defining dimensionless parameters corresponding to the reduced gravity conditions $g=g_{0}\times 10^{-1}$ ( $g_{0}=9.81~\text{m}~\text{s}^{-2}$ ): $Ca=2\times 10^{-3}$ , $Mn=112$ , $G=545$ , $Gr_{1}=200$ , $Gr_{2}=5751$ , $Pr_{1}=0.5$ , $Pr_{2}=16.8$ , $Re_{2}=10$ . In both cases, at the initial state, both media are at rest, and the interface between them is flat.

Figure 2. Temporal alteration of the flow structure in the system with the heater with constant temperature on the substrate, $Q_{1}^{s}=\{(x,y):x\in [8.5;11.5],\,y=0\}$ , $q_{1}^{s}=0.1$ . Flow patterns at $t=2$  s (a), $t=8.5$  s (b), $t=26.5$  s (c), $t=50$  s (d).

Figure 3. Evolution of the thermal field (a,c,e) and the interface location (b,d,f) in the system with the heater with constant temperature on the substrate at $t=8.5$  s (a,b), $t=26.5$  s (c,d), $t=50$  s (e,f), $Q_{1}^{s}=\{(x,y):x\in [8.5;11.5],\,y=0\}$ , $q_{1}^{s}=0.1$ .

3 Steady heating from below by a single heater

We investigate the behaviour of the system characteristics under the conditions where one heater is located in the centre of the substrate. When the heater with constant temperature $q_{1}^{s}=0.1$ is switched on, we first observe the formation of convective structures generated by the action of mass forces. Under the effect of buoyancy forces, the hotter liquid rises up (figure 2 a). A dual-vortex flow structure is observed in each layer. Heat from the thermal source is transferred into the fluid bulk by convection (see the temperature distribution in figure 3) since the Prandtl number for the chosen pair of working media is quite large. This directly follows from the energy equation (2.3). Diffusive transfer in any direction is too slight. The formation of the vortices in the upper gas layer, which are weaker in comparison with the swirls in the lower liquid layer, is conditioned by continuity of the liquid–gas system provided by the equality of tangential velocities of the media at the interface. Once heat reaches the interface, the surface forces come into action. Under the thermocapillary effect, the liquid spreads along the interface, transferring heat from hot spots to the areas with lower temperature. The action of Marangoni forces leads to alteration of the temperature field in the liquid. As a result, three indicative temperature peaks appear (figure 3 a). In other words, the temperature gradient is observed both at the interface and in the volume of the liquid. It results in the development of four-vortex motion in each layer (figure 2 b). Thus, such a transition from the two-vortex flow structure to the quadruple-vortex pattern observed in both fluids is induced by the action of both buoyancy forces and the thermocapillary effect. The surface forces additionally cause the deformation of the interface. Note that at the initial stage, the thermocapillary deflection of the interface in the zone of the heater has two local dimples (concave meniscus) located above the areas of ascending flows in each pair of vortices with opposite circulation. Between these dimples there is an inflection zone with the positive curvature (the surface is convex outward figure 3 b).

Over time, the temperatures in the internal part of the layer and at the interface above the heater gradually equalize (figure 3 e), and the temperature gradient becomes close to zero. It excludes the action of both mass forces and surface ones in the area above the heater. Therefore, the pair of ‘internal’ swirls, caused by the action of these forces, weakens (figure 2 c) and then vanishes completely (figure 2 d), since the reasons for their existence disappear. Strong heating in the area of the heater action attenuates the surface tension of the interface and it ‘hangs down’. Such an alteration of the velocity field and character of the interface deformation, induced by local thermal exposure, can be explained by analysing equation (2.1) and the right-hand side (2.19b ) of the boundary condition for vorticity (2.19a ). The summand $MnCa\unicode[STIX]{x1D703}$ in (2.19b ) characterizes the degree of deformability of the interface with varying temperature. With the growth of $\unicode[STIX]{x1D703}$ on the surface the value $MnCa\unicode[STIX]{x1D703}$ is increased. This means that the surface tension is decreased since the term is included on the right-hand side with the negative sign. At the same time, the equalization of the temperature in the area of the heater action and subsequent decrease of the temperature gradient result in the last term on the right-hand side of momentum transfer equation (2.1) becoming small, whereas the convective terms remain large. In such a way the reverse transition to the double-vortex flow in each layer occurs with time. When the ‘internal’ pair of vortices vanishes, then the ‘external’ ones get close to each other. Beyond the domain of the local thermal load a temperature gradient is observed, and therefore, provides additional pressure. Here, the real pressure is higher than the hydrostatic one. But the pressure within the area above the heater is lower than the hydrostatic one. The interface flexure is indicative of this. Otherwise, if the real pressure is equal to the hydrostatic one, then there will be no deflection. Thus, the vortices move from the domain of higher pressure into the region with the lower one. The process continues until the equilibrium between the force of buoyancy, the Marangoni effect and the temperature on the entire outer contour of the layer system is reached. The upper gas layer completely repeats the movement of the liquid layer due to the requirement on equality of the tangential velocities. In the course of time, the flow pattern is not changed. The dual-vortex flow structure is retained in each layer, which is supported by thermocapillary spreading of the liquid along the interface from the hot pole to the region with the lower temperature (figures 2 d and 3 e). Here, the steady thermocapillary deflection with the negative curvature is formed in the area of the heater (the surface is convex inward, figure 3 f). Such a quadruple-vortex flow in the liquid–gas system represents an asymptotic steady regime.

Thus, successive transitions from the double-vortex mode to the four-vortex one and back in each layer (we call it the 2–4–2 transition) of the steady locally heated bilayered system are observed when activating the heater. The transient regimes are caused by the joint action of convective and thermocapillary mechanisms. Later, the flow stabilization occurs where the steady thermocapillary deflection in the area of heating is formed and the quadruple-vortex flow in the liquid–gas system is reached. It should be noted that such a scenario of system behaviour is typical also for systems with several heaters having equal constant temperature that are well spaced. Here, if there are $n$ heaters on the substrate we observe $2n$ – $4n$ – $2n$ transition in each layer of the system. The issue of mutual location of the heaters requires a particular study.

4 Non-uniform heating from below by two heaters

One way to control the effective operation of the thermal stabilization system is the commutated mode of work for the heaters, which ensures long-term continuous operation of the equipment. Such conditions allow for increasing of the heater life time. We have investigated the characteristics of the two-layer system in the following successive thermal load conditions: (i) at the initial time the temperature of both heaters is the same $q_{A}^{s}=q_{B}^{s}=0.1$ , and then, the temperature of the heater $B$ increases abruptly according to the following law $q_{B}^{s}(10)=0.25$ , $q_{B}^{s}(20)=0.5$ , $q_{B}^{s}(30)=0.75$ , $q_{B}^{s}(40)=1$ , while the heater $A$ produces heat at a constant initial temperature; (ii) the temperature of the heater $A$ increases abruptly, while the temperature of the second heater $B$ decreases according to the following law: $q_{A}^{s}(50)=0.25$ , $q_{B}^{s}(50)=0.75$ , $q_{A}^{s}(60)=0.5$ , $q_{B}^{s}(60)=0.5$ , $q_{A}^{s}(70)=0.75$ , $q_{B}^{s}(70)=0.25$ , $q_{A}^{s}(90)=1$ , $q_{B}^{s}(90)=0.1$ . Reducing the temperature of heater $B$ simulates the cooling of the thermal element when it is turned off.

Figure 4. Temporal alteration of the flow structure in the system with switchable heating by two heaters on the substrate, $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5],\,y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],\,y=0\}$ . Flow patterns at $t=23.5$  s (a), $t=32$  s (b), $t=44$  s (c).

Figure 5. Evolution of the thermal field (a,c,e,g) and the interface location (b,d,f,h) in the system with combined heating by two heaters on the substrate at $t=23.5$  s (a,b), $t=32$  s (c,d), $t=33.5$  s (e,f), $t=44$  s (g,h), $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5]$ , $y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],\,y=0\}$ .

With each increase in the temperature of the heater, there is a certain time delay during which the system responds to the change in the thermal load. In this case, we observe the regimes at which the interface undergoes significant oscillations with the change in the flexure form and deformation amplitude not only in the area above the heater, but also along the entire length of the tested area (figure 5). In such regimes, small vortices are formed near the hot spot in each layer, which are split and drift to the side walls (figure 4). During the oscillations of the interface, certain surface shapes with two or more menisci appear periodically (figure 5 b,d,f). The formation of these modes with the interface oscillations accompanied by the alteration of the deflection form and appearance of travelling stretched eddies is a characteristic feature of the switchable heating regime. The interface oscillations and vortex drift with the change of their form occur only if the heater temperature $q^{s}$ increases over time.

Analysis of the hydrodynamic field structure shows that a boundary layer arises near the interface. The boundary layer flattens a core where vortex formations appear. It leads to changes in the form of the vortices, they become ‘flat’. The boundary layer is particularly evident within and close by the heating zone, and it vanishes with distance from region of local thermal exposure. In the heating domain the boundary layer is subjected to thermal attack, it is distorted along with the interface and deforms the flux core. When there is a steady thermal load the interface deformation also occurs, but it is unidirectional and the interface does not undergo any oscillations. The surface tension becomes weaker under each next temperature jump. At the same time the buoyancy force pushes the interface up while the Marangoni forces, acting with some delay, pull the interface in different directions. In such a way the conditions emerge to destabilize the interface, boundary layer and core of the flux as a whole. The vortices break away and drift due to all these factors. As soon as the heater temperature reaches a maximal constant value (in our case $q_{s}=1$ ) the breakaway and drift of the swirls cease, and the motion is stabilized; it is supported by the balance of all forces. However, vortices, having already changed their shape and being flattened, no longer take another shape. Constant tangent velocity and the interface curvature are set at the interface. They dictate the form of the vortex formations.

The presence of the weak heater $A$ breaks the flow symmetry, leading to the deformation of the vortex structures, and formation of short-term asymmetric menisci above the more powerful heater with a certain time delay after each abrupt change in temperature $q_{B}^{s}$ at the stage (i) (figure 5 d,f). The lag time corresponds to the time in which heat from the thermal element $B$ reaches the interface. The convective transfer predominates over the thermocapillary one in this time period. If the liquid layer thickness is decreased, then the lag time is reduced. Over time, the vortex structures arising in the zone of the weak heater $A$ are merged by the flow, which is induced by the second heater $B$ (figure 4 c). If the temperature of the heater $B$ does not increase further, the resulting mode is stabilized, and a stable thermocapillary deflection is formed in the zone of this heater (with a higher temperature). Note, that the stretched form of the vortices caused by the oscillations of the interface will be kept even after stabilization of the free surface form.

Further, we describe the characteristics of the mode arising at stage (ii). At the first change in the temperature of the heaters (at the time point $t=50$  s), weak oscillations (with a small amplitude) of the interface are most noticeable in the area above the heater $B$ . Within the zone of the heater $A$ we observe the flow intensification, resulting in the growth of the thermal gradient, and increase in the size of the drifting vortices (compare figures 4(c) and 6 a). The next temperature switching (at the time point $t=60$  s) leads to the transition from the mode with splitting vortices to the triple-vortex regime (figure 6 b) and then, to the quadruple-vortex flow (figure 6 c) in each layer. In such transient modes, interface shapes with three or two concave menisci and inflection zones (figure 7) can be observed, while the oscillations of the interface are preserved. At two subsequent changes in the temperature of the heaters, the modes are still observed in the system. They are accompanied by fluctuations in the interface (figure 9 b,d), and by the formation of splitting drifting vortices, now within the area of heater $A$ , which has a higher temperature than heater $B$ , as well as transitions from quadruple-vortex flows to triple-vortex ones, and finally to dual-vortex flows in each layer (figure 8), similar to those regimes which are typical for the mode arising at stage (i).

Figure 6. Temporal alteration of the flow structure in the system with two switched on heaters on the substrate, $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5],\,y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],\,y=0\}$ . Flow patterns at $t=59$  s (a), $t=64.5$  s (b), $t=66.5$  s (c).

Figure 7. Evolution of the thermal field (a,c) and the interface location (b,d) in the system with two switched on heaters on the substrate at $t=64.5$  s (a,b), $t=66.5$  s (c,d), $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5],\,y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],\,y=0\}$ .

Figure 8. Temporal alteration of the flow structure in the system with two switched on heaters on the substrate, $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5],\,y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],\,y=0\}$ . Flow patterns at $t=72$  s (a), $t=78.5$  s (b), $t=94$  s (c).

Figure 9. Evolution of the thermal field (a,c,e) and the interface location (b,d,f) in the system with two switched on heaters on the substrate at $t=72$  s (a,b), $t=78.5$  s (c,d), $t=94$  s (e,f), $Q_{A}^{s}=\{(x,y):x\in [5.5;8.5],\,y=0\}$ , $Q_{B}^{s}=\{(x,y):x\in [11.5;14.5],y=0\}$ .

Additional thermal elements can be used to prevent the appearance of regimes that are close to the critical one (with the critical decrease of the liquid layer thickness or thermal hysteresis where the parameters of an electronic device do not revert to the initial values after its heating and subsequent cooling). Switching to these additional heaters allows one to ensure continuity of the operation of onboard equipment without losing effective power. The lifetime of thermal elements can be prolonged by placing a duplicate heater and periodic switching of the electric circuit from one heater to another while the system relaxes (compare figures 5 h and 9 f, where the interface forms are presented before the commutation to the second heater (figure 5 h) and after that (figure 9 f)).

It should be noted that oscillations of the interface with the increase in the deformation amplitude and short-term changes of the thermocapillary deflection form are observed only under conditions of commutated heating with an abrupt increase of the temperature in switching. If the heater temperature is increased continuously, then an inertial behaviour of the free boundary will be observed and no oscillations of the interface, hence regimes with splitting and drifting to the periphery vortices will appear. The interface will relax monotonically under an abrupt decreasing of the heater temperature.

5 Concluding remarks

The problem of thermocapillary convection in a two-layer system undergoing local heating has been considered. The use of the Oberbeck–Boussinesq equations to describe the motion and heat transfer, and the formulation of conditions at the interface in the form of conservation laws, allow us to comprehend the role and degree of influence of particular mechanisms and their combinations on the characteristics of convective regimes arising in the system under consideration.

Within the framework of the proposed approach, we described the formation of complex cellular motions of different symmetries caused by the combined effect of convective and thermocapillary mechanisms, competing structures arising under the thermal load from the heaters arranged on the bottom boundary of the system and various scenarios of the phase interface behaviour. Based on numerical simulation, we have shown the existence of successive transient regimes of the 2–4–2 type in the system with steady heating by one heater, which are characterized by the appearance of multi-vortex structures, depending on the features of reciprocity of the basic mechanisms. The resulting effect of joint action of these mechanisms in the case of the steady (uniform) heating leads to the formation of stable regimes. Under the conditions of non-uniform thermal load the system responds to the local temperature change with a certain time delay. Abrupt increase of the heater temperature results in oscillations of the free surface accompanied by a change in the deformation amplitude and its form, and the appearance of small eddies travelling from the hot spot to the lateral walls of the container.

The results for the case of one heater with constant temperature obtained in the framework of the complete formulation of the thermocapillary convection problem are in good agreement with the results obtained in the model of the thin layer approximation. The pattern of the asymptotic double-vortex regime forming in the liquid (lower) layer of a bilayered system subjected to steady local heating by one heater and the form of the thermocapillary deflection which are predicted by our model are identical to the topological structure of the flow and configuration of the free surface obtained in Marchuk (Reference Marchuk2009).

The results obtained allow one to better understand the formation and interaction mechanisms of vortex, thermocapillary and convective structures in systems with internal interfaces. The results of the present study can aid in designing thermal control systems, involving local heaters with the specified characteristics of heat removal and operation modes, and optimization of fluidic section configurations.