2. Formulation of the problem

An inextensible thin sheet of infinite length, thickness $ h$ , density $ \rho _{\textit{sh}}$ , and bending modulus $ B$ , lies on top of a Newtonian fluid layer of depth $ d$ and is subjected to uniaxial compression by a lateral force per unit length $ F_x$ . The bottom of the fluid, assumed to be a rigid surface, and the sheet are maintained at constant temperatures $ T_{{h}}$ and $ T_{{c}}$ , respectively, where $ T_{{h}} \gt T_{{c}}$ . The fluid has kinematic viscosity $\nu$ , thermal diffusivity $\kappa$ and temperature-dependent density

(2.1) \begin{equation} \rho _{\ell }=\rho _{\ell }^0\left [1-\alpha (T-T_{{c}})\right ]\! , \end{equation}

where $\alpha$ is the coefficient of volume expansion, $T$ is the temperature of the fluid, and $\rho _{\ell }^0$ is the density of the fluid at temperature $T_{{c}}$ . The system is subjected to gravitational acceleration $g$ in the negative $y$ -direction, where we place a Cartesian coordinate system at the undeformed configuration of the sheet, as in figure 1.

Initially, the fluid is at rest, and the sheet is flat but laterally compressed by the force $ F_x$ . The temperature difference $ \Delta T \equiv T_{{h}} - T_{{c}}$ is then increased until instability occurs in the form of natural convection. In this analysis, we examine the critical conditions at which instability arises, and the characteristics of the fastest-growing unstable mode.

Figure 1. The system set-up consists of a fluid layer of depth $ d$ , covered by a thin elastic sheet that is subjected to lateral compression. The bottom surface of the fluid and the sheet are maintained at constant temperatures $ T_{{h}}$ and $ T_{{c}}$ , respectively. The sheet remains flat, and the fluid is at rest before instability sets in (not shown in the figure). However, after instability has set in, the fluid flow may induce buckling of the upper sheet.

Hereafter, we use the standard normalisation of the Rayleigh–Bénard set-up: we normalise all lengths to the depth of the fluid layer $d$ , and time to the thermal diffusive time scale ${d}^2/\kappa$ . In addition, we define the temperature relative to $T_{{c}}$ , and normalise it by the temperature difference $\Delta T$ , i.e. $T\rightarrow (T-T_{{c}})/\Delta T$ . Taking into consideration the physical parameters of the system, this gives rise to four independent dimensionless numbers:

(2.2) \begin{equation} P_r=\frac {\nu }{\kappa } , \quad R_g=\frac {gd^3}{\kappa \nu }, \quad \varLambda =\frac {B}{\kappa ^2 \rho _{\textit{sh}}h}, \quad \lambda =\frac {\rho _{\textit{sh}}h}{\rho _{\ell }^0d}, \end{equation}

where $ P_r$ is the Prandtl number, $ R_g$ represents the ratio between gravity and thermal and viscous diffusivity (sometimes referred to as the Galilei number; Or & Kelly Reference Or and Kelly2002), $ \varLambda$ denotes the ratio of the thermal diffusivity time scale to the inertial time scale of the sheet, and $ \lambda$ is the sheet-to-fluid mass ratio. In addition, there are two control parameters:

(2.3) \begin{equation} R_a=\frac {gd^3 \alpha\, \Delta T}{\kappa \nu }, \quad f_x=\frac {F_x {d}^2}{B}, \end{equation}

where $R_a$ is the Rayleigh number, which is proportional to the temperature difference $\Delta T$ , and $f_x$ is the normalised lateral force.

The state of the fluid is determined by four fields, two of which are the components of the fluid’s velocity vector $\boldsymbol { u}(x,y,t)=(u_x,u_y)$ ; the other two are the pressure $p(x,y,t)$ and temperature $T(x,y,t)$ fields. The temporal evolution of these fields is determined by equations of motion that describe the balance of mass, momentum and energy. In the Boussinesq approximation, these equations read (Chandrasekhar Reference Chandrasekhar1961)

(2.4a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol { u}&=0, \end{align}

(2.4b) \begin{align} \frac {1}{P_r}\frac {{\textrm{D}}\boldsymbol { u}}{{\textrm{D}}t} &=\boldsymbol{\nabla }\boldsymbol{\cdot }\sigma + (R_a T-R_g)\hat{\boldsymbol{y}}, \\[-12pt] \nonumber \end{align}

(2.4c) \begin{align} \frac {{\textrm{D}}T}{{\textrm{D}}t} &={\nabla} ^2 T, \\[12pt] \nonumber \end{align}

where ${\textrm{D}}/{\textrm{D}}t=\partial /\partial t+\boldsymbol { u}\boldsymbol{\cdot }\boldsymbol{\nabla}$ is the convective derivative, $\boldsymbol{\nabla}$ is the two-dimensional gradient operator, and $\sigma _{ij}=-p\delta _{ij}+ ({\partial u_i}/{\partial x_j})+({\partial u_j}/{\partial x_i})$ is the stress tensor of the fluid, where $\delta _{ij}$ is the isotropic tensor.

At $y=-1$ , the fluid is at rest and is held at a constant temperature. This gives the boundary conditions

(2.5) \begin{equation} u_x(x,-1,t)=0, \quad u_y(x,-1,t)=0, \quad T(x,-1,t)=1. \end{equation}

At the upper layer, the fluid is in contact with the sheet. To describe this contact, we first define the function $y_{\textit{sh}}(x,t)$ , which represents the height of the sheet relative to the flat configuration. Additionally, we restrict the formulation to shallow deflections of the sheet (Landau & Lifshitz Reference Landau and Lifshitz1986) and neglect deformations in the $x$ -direction. (For an inextensible sheet, deformations in the $x$ -direction are proportional to $y_{\textit{sh}}^2$ and are later negligible in the linear stability analysis.) Since the sheet and the fluid have the same velocity at their interface, and the sheet is maintained at a constant temperature, we prescribe the following boundary conditions:

(2.6) \begin{equation} u_x(x,y_{\textit{sh}},t)=0, \quad u_y(x,y_{\textit{sh}},t)=\frac {\partial y_{\textit{sh}}}{\partial t}, \quad T(x,y_{\textit{sh}},t)=0. \end{equation}

To close the system of equations, we must also balance the sheet’s inertia in the $y$ -direction with the three distinct forces acting along this direction – namely, the bending force, the lateral compression and the vertical component of the hydrodynamic force. After normalisation, this equation takes the form (Brau et al. Reference Brau, Damman, Diamant and Witten2013; Diamant Reference Diamant2021; Goncharuk, Feldman & Oshri Reference Goncharuk, Feldman and Oshri2023; Oshri, Goncharuk & Feldman Reference Oshri, Goncharuk and Feldman2024)

(2.7) \begin{equation} \frac {\partial ^2y_{\textit{sh}}}{\partial t^2}+\varLambda \frac {\partial ^4 y_{\textit{sh}}}{\partial x^4}+ \varLambda f_x \frac {\partial ^2y_{\textit{sh}} }{\partial x^2}+\frac {P_r}{\lambda }(\sigma _{yy})_{y=y_{\textit{sh}}}=0. \end{equation}

This completes the formulation of the problem. In summary, given the four dimensionless parameters $ P_r$ , $ R_g$ , $ \varLambda$ and $ \lambda$ in (2.2), and setting values for the control parameters, i.e. $ R_a$ and $ f_x$ , we can, in principle, solve (2.4)–( 2.7) by assuming periodic boundary conditions in the $ x$ -direction. However, the following analysis is focused not on finding general solutions to these equations, but rather on obtaining the critical conditions at which instability occurs. To determine these conditions, we need to add one more step to this formulation, and linearise the equations around their static base solution. This linearisation is described in the next section.

3. Linear stability analysis

Before the onset of instability, the sheet is flat and the fluid is at rest. To maintain this stationary state, while satisfying the balance of momentum and energy (2.4b ) and (2.4c ), a quadratic pressure field and a linear temperature field are held constant along the $y$ -direction. After the instability sets in, these static conditions are slightly perturbed. Therefore, to describe the instability, we introduce a small perturbation of size $\epsilon \ll 1$ that grows exponentially over time with a growth rate $\sigma$ , on top of the static solution:

(3.1a) \begin{align} \boldsymbol { u}(x,y,t)&=\epsilon \int _{-\infty }^{\infty }{\textrm{e}}^{\sigma t +{\textrm{i}}kx}\,\boldsymbol { U}(k,y)\,{\textrm{d}}k, \end{align}

(3.1b) \begin{align} p(x,y,t)&=-R_g y-\frac {R_a}{2}y^2+\epsilon \int _{-\infty }^{\infty }{\textrm{e}}^{\sigma t+{\textrm{i}}kx}\,P(k,y)\,{\textrm{d}}k, \end{align}

(3.1c) \begin{align} T(x,y,t)&=-y+\epsilon \int _{-\infty }^{\infty }{\textrm{e}}^{\sigma t+{\textrm{i}}kx}\,\varTheta (k,y)\,{\textrm{d}}k, \end{align}

(3.1d) \begin{align} y_{\textit{sh}}(x,t)&=\epsilon \int _{-\infty }^{\infty }{\textrm{e}}^{\sigma t+{\textrm{i}}kx}\,Y_{\textit{sh}}(k)\,{\textrm{d}}k, \end{align}

where $\boldsymbol { U}(k,y)$ , $P(k,y)$ and $\varTheta (k,y)$ represent, respectively, perturbations of the fluid’s velocity, pressure and temperature fields with the Fourier transform applied in the $x$ -direction. Correspondingly, $Y_{\textit{sh}}(k)$ is the perturbation of the sheet’s height function with wavenumber $k$ around the flat solution. Substituting the above perturbed fields in the nonlinear equations, and expanding to linear order in $\epsilon$ , yields, after some simplifications, the following linear equations:

(3.2a) \begin{align} \left (\frac {{\textrm{d}}^2}{{\textrm{d}}y^2}-k^2\right )\left (\frac {{\textrm{d}}^2}{{\textrm{d}}y^2}-\frac {\sigma }{P_r}-k^2\right )U_y&=R_a k^2 \varTheta , \end{align}

(3.2b) \begin{align} \left (\frac {{\textrm{d}}^2}{{\textrm{d}}y^2}-k^2-\sigma \right )\varTheta &=-U_y, \end{align}

(3.2c) \begin{align} \left (\sigma ^2+\varLambda k^4-\varLambda f_x k^2+\frac {P_r R_g}{\lambda }\right )Y_{\textit{sh}}&=\frac {P_r}{\lambda k^2}\left (\frac {{\textrm{d}}^3U_y}{{\textrm{d}}y^3}\right )_{y=0}\!, \end{align}

where we have eliminated $P(k,y)$ and $U_x(k,y)$ in favour of $U_y(k,y)$ ; see Appendix A for details. These equations are supplemented by the following boundary conditions:

(3.3a) \begin{align} U_y(-1) &=0, \quad \frac {{\textrm{d}}U_y}{{\textrm{d}}y}(-1)=0, \quad \varTheta (-1)=0, \end{align}

(3.3b) \begin{align} U_y(0) &=\sigma Y_{\textit{sh}}, \quad \frac {{\textrm{d}}U_y}{{\textrm{d}}y}(0)=0, \quad \varTheta (0)=Y_{\textit{sh}}. \end{align}

Note that the boundary condition $\varTheta (0) = Y_{\textit{sh}}$ arises from the fact that the temperature is held constant on the sheet’s surface, rather than at $y=0$ .

Equations (3.2) and (3.3) form a closed system of linear and homogeneous equations, and thereby complete the linear stability formulation. For a given set of parameters (2.2), these equations have a non-trivial solution if their corresponding determinant equals zero. The critical conditions for the onset of instability are obtained, for example, by fixing $f_x$ , setting the real part of $\sigma$ to zero, and then finding the smallest $R_a$ at which the determinant vanishes to zero. This typically results in two algebraic equations that must be solved simultaneously for the wavenumber $k$ and the imaginary part of the growth rate. Nonetheless, in all cases considered below, we find – through the numerical solution of (3.2) and (3.3) – that the imaginary part of $\sigma$ equals zero when instability occurs (see Appendix A.1 for details on the numerical method). This finding simplifies the analysis considerably, and we use the simplified analysis in the analytical derivations that follow.

Before we continue to explore the solutions of (3.2) and (3.3), it is instructive to recap two familiar cases of instability from this formulation: one driven purely by thermal effects, as in the classical Rayleigh–Bénard set-up, and the other driven purely by elastic effects. In the first case, the sheet is assumed to be very stiff, and the fluid is essentially enclosed between two rigid plates that are held at different temperatures. In the second case, the temperature difference is zero, and the sheet is laterally compressed until the wrinkling instability occurs. We recap these two limiting scenarios in the next subsection.

3.1. Recap of previous studies

This subsection is divided into two parts. In the first part, we recap the Rayleigh–Bénard instability that is driven by buoyancy and thermal effects, and in the second part, we recap the wrinkling instability that emerges in a laterally compressed sheet lying on a fluid.

3.1.1. The thermal instability

To recap the scenario of thermal instability, we must first eliminate the elastic response of the sheet from the model. To this end, we first set the lateral force $ f_x$ to zero, and increase the bending stiffness ( $ \varLambda k^4 \gg 1$ ), so as to inhibit elastic deformations. Under these conditions, it readily follows from (3.2c ) that the height function falls to zero and can be neglected from the formulation.

In this limit, the marginal stability problem ( $\sigma =0$ ) depends on a single parameter, the Rayleigh number $R_a$ , and recaps the Rayleigh–Bénard set-up of two rigid plates at different temperatures that enclose a fluid layer. Although the solution to this problem is provided in many textbooks (see e.g. Chandrasekhar Reference Chandrasekhar1961; Drazin Reference Drazin2002; Goluskin Reference Goluskin2016), it is instructive at this point to introduce an alternative approximated approach for the solution, as we later use this approach in the analysis of the more general case. We start the derivation by guessing the following function for the perturbed temperature:

(3.4) \begin{equation} \varTheta (k,y)=A \sin (\unicode{x03C0} y), \end{equation}

where $A$ is as yet an unknown constant. While (3.4) satisfies the boundary conditions on the two rigid plates, it provides only an approximate solution to the governing equations. It should be understood as the leading term in a Fourier sine expansion, which, as we now show, captures the solution with sufficient accuracy. Next, we substitute this expansion in (3.2a ) and solve for $U_y(k,y)$ with the boundary conditions given in (3.3). (The solution for the vertical velocity is $U_y(k,y)= {\textit{AR}_a k^2} {[(k+\sinh k)\sin (\unicode{x03C0} y)-\unicode{x03C0} (1+y)\sinh (k y)-\unicode{x03C0} y\sinh (k(1+y) ) ]}/{ (k^2+\unicode{x03C0} ^2 )^2(k+\sinh k)}$ ). Finally, we substitute $U_y(k,y)$ and $\varTheta (k,y)$ in (3.2b ), multiply the equation by $\sin (\unicode{x03C0} y)$ , and integrate over the domain, i.e. solve the weak form of (3.2b ). This procedure yields a linear and homogeneous equation for the amplitude $A$ , which has a non-trivial solution when its coefficient equals zero. This gives the following relation between $R_a$ and $k$ :

(3.5) \begin{equation} R_a=\frac {(k^2+\unicode{x03C0} ^2)^5\left (k+\sinh k\right )}{-8\unicode{x03C0} ^2k^3\left (1+\cosh k\right )+k^2 (k^2+\unicode{x03C0} ^2)^2\left (k+\sinh k\right )}. \end{equation}

The instability occurs at the wavenumber at which $R_a$ is minimised. This minimisation gives

(3.6) \begin{equation} k\equiv k_{\textit{RB}}\approx 3.114 \quad \text{and} \quad R_a\equiv R_a^{\textit{RB}}\approx 1715, \end{equation}

where we use subscript or superscript ‘RB’ to denote quantities related to the Rayleigh–Bénard set-up. Notably, this approximated solution is very close to the numerical solution of the equations, which yields $k\simeq 3.117$ and $R_a\simeq 1707$ (Chandrasekhar Reference Chandrasekhar1961).

In figure 2(a), we plot the marginal stability line on the $ (k, R_a)$ plane, obtained from the numerical solution of (3.2) and (3.3), and compare it with the approximated analytical prediction (3.5). As may be seen, the two solutions are in excellent agreement. To obtain the marginal stability numerically, we first fix the wavenumber $ k$ , set the Rayleigh number to zero, and calculate the growth rate $ \sigma$ . Then we increase the Rayleigh number by some increment, and repeat the calculation until the real part of the growth rate becomes positive. When this first occurs, we identify the corresponding value as the critical Rayleigh number.

At the instability, the flow field exhibits a periodic pattern of vortices, where a unit cell of this pattern consists of two counter-rotating vortices. The vertical length of the cell coincides with the separation between the plates, and the horizontal length is given by $ \lambda _{\textit{RB}} = 2\unicode{x03C0} /k_{\textit{RB}} \simeq 2.01$ ; see figure 2(b). When dimensions are restored in this analysis, it becomes clear that the emerging pattern at the instability depends solely on the distance $d$ between the two rigid plates. However, this length scale does not play any role in the elastic instability, as demonstrated in the next section.

Figure 2. Solution of the linear stability analysis in the Rayleigh–Bénard set-up. (a) The marginal stability diagram on the $(k, R_a)$ plane. The solid black line corresponds to (3.5), and the open blue circles correspond to the numerical solution of (3.2) and (3.3). Instability first occurs when the Rayleigh number reaches $ R_a^{\textit{RB}}$ . At this value, the emerging pattern has a wavenumber $ k_{\textit{RB}}$ . (b) The flow field in the marginal stability solution $(k_{\textit{RB}}, R_a^{\textit{RB}})$ consists of a periodic structure of vortices, whose typical length scales with the separation between the plates $ d$ . Arrows represent the streamlines, and colours represent the magnitude of the perturbed velocity. (c) The perturbed temperature at the instability $ \varTheta (k_{\textit{RB}}, y)$ . The solid line corresponds to the profile given by (3.4) with the normalisation $ A = 1$ , and the open circles correspond to the numerical data.

3.1.2. The wrinkling instability

To derive the wrinkling instability from our model, we must first eliminate all thermal effects. Since early in this formulation we normalised the temperature by $ \Delta T$ , we now reintroduce the temperature dimensions into the model and set $ \Delta T = 0$ . This procedure removes $ \varTheta (k, y)$ from the formulation, and leaves us to solve (3.2a ) and (3.2c ), along with the boundary conditions in (3.3).

When $ \sigma = 0$ , these equations allow a non-trivial solution when the velocity field vanishes, and $ f_x = k^2 + ({P_r R_g}/{\varLambda \lambda k^2})$ . The instability occurs at the wavenumber at which $ f_x$ is minimised, giving

(3.7) \begin{equation} k\equiv k_{\textit{sh}}=\left (\frac {P_r R_g}{\varLambda \lambda }\right )^{1/4} \quad \text{and} \quad f_x^{\textit{sh}}=2k_{\textit{sh}}^2, \end{equation}

where we use subscript or superscript ‘sh’ to denote quantities related to the elastic instability.

Equation (3.7) may initially appear to be counterintuitive, because the elastic instability depends neither on thermal effects nor on fluid dynamics. However, thermal conductivity and viscosity appear in the emergent wavelength. This outcome is simply a result of our choice of normalisation. When reintroducing the dimensions into the model, we find that the wavenumber $ \tilde {k}_{\textit{sh}} = (\rho _{\ell }^0 g / B)^{1/4}$ and the force $ F_x^{\textit{sh}} = 2(\rho _{\ell }^0 g B)^{1/2}$ depend only on the sheet’s elasticity ( $ B$ ) and the specific weight of the fluid ( $ \rho _{\ell }^0 g$ ). (Here, we use a tilde to distinguish the dimensional wavenumber from its dimensionless counterpart.)

When $f_x\lt f_x^{\textit{sh}}$ , the sheet remains flat, and the system is stable against out-of-plane deformations. However, at $f_x^{\textit{sh}}$ , the flat solution becomes unstable, and buckling occurs. Since this analysis is concerned with the instability of the system around the flat configuration, we must limit the lateral force to the range $0\leqslant f_x\lt f_x^{\textit{sh}}$ .

In contrast to the thermal instability, the dominant length scale in the elastic scenario is entirely independent of the depth of the fluid. Both length scales, $ (B/\rho _{\ell }^0 g)^{1/4}$ and $ d$ , can influence the emergent pattern only when the instability is driven by thermal and elastic effects. Exploring these effects is the main focus of this work, as discussed further in the next section.

4. Investigation of the thermoelastic instability

We now address the general scenario, in which convection is driven by both thermal and elastic effects. Given the relatively large number of parameters defining this problem, it would be useful, at this point, to focus on a specific range of parameters that are more experimentally accessible. For typical properties of the sheet, we refer, for example, to the experiments in Pocivavsek et al. (Reference Pocivavsek, Dellsy, Kern, Johnson, Lin, Lee and Cerda2008), and assume that the sheet is made of a polyester-like material, with Young’s modulus of the order of ${\sim}1\,\text{GPa}$ , density , Poisson’s ratio $ \sim 0.3$ , and thickness in the range $h \in [0.1, 10]\ \unicode{x03BC} \textrm{m}$ . Additionally, we take the fluid substrate to be water, with density , thermal diffusivity $\kappa \sim 10^{-7}\, {\textrm{m}}^2\,{s}^{-1}$ and kinematic viscosity $\nu \sim 10^{-6}\, {\textrm{m}}^2\,{s}^{-1}$ . The distance between the bottom of the fluid and the sheet is assumed to range from the sub-millimetre to the centimetre scale, $d \in [0.5, 100] \,\text{mm}$ (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000; Schatz & Neitzel Reference Schatz and Neitzel2001). From these considerations, we can fix the Prandtl number at $P_r \sim 7$ , as is typical for water, and treat the mass ratio $\lambda$ as a small number, $\lambda \in [10^{-5},10^{-2}]$ . Additionally, $R_g \in [10^3, 10^7]$ and $\varLambda \in [10^5, 10^8]$ are typically large but can vary significantly due to changes in the depth $d$ and the thickness of the sheet. This gives an elastic wavelength that ranges from approximately a millimetre to a centimetre, and can be either much larger or even smaller than the thermal wavelength.

Consider, for example, the marginal stability diagram in figure 3, where the parameters are chosen such that $ k_{\textit{sh}} \simeq 15.4 \gt k_{\textit{RB}} \simeq 3.11$ . When $ f_x/f_x^{\textit{sh}} \ll 1$ , the system exhibits behaviour very similar to that of two rigid plates, i.e. the instability occurs only slightly below $R_a^{\textit{RB}}$ , and the wavenumber is $ k_{\textit{RB}}$ . However, if the lateral force is further increased, then a new minimum appears in the diagram ( $ f_x/f_x^{\textit{sh}} = 0.82$ in figure 3), with a wavenumber that coincides approximately with that of the sheet. Initially, this minimum has a larger Rayleigh number than $ R_a^{\textit{RB}}$ , therefore the system is still governed by thermal instability. But as the lateral force increases further, this minimum pre-empts the Rayleigh–Bénard solution, and the instability first appears at the elastic wavelength ( $ f_x/f_x^{\textit{sh}} = 0.995$ in figure 3). This example essentially demonstrates that adding an elastic sheet on top of the fluid layer can, under some conditions, have a dramatic effect in the pattern that emerges at the onset of instability.

Figure 3. Marginal stability diagram for a fluid-supported sheet. Solid lines correspond to the analytical approach, while open blue circles represent the numerical solution of (3.2) and (3.3). The following parameters are assumed: $R_g = 4500$ , $\lambda = 5 \times 10^{-5}$ , $P_r = 7$ and $\varLambda = 11 \times 10^3$ . From (3.7), the elastic wavenumber is $k_{\textit{sh}} \simeq 15.46$ . For very small values of the lateral force, instability occurs slightly below $R_{a}^{\textit{RB}}$ at the wavenumber $k_{\textit{RB}}\simeq 3.11$ . However, as the lateral force increases, a new minimum appears with wavenumber $k_{\textit{sh}}$ , which eventually pre-empts the Rayleigh–Bénard solution. In this example, the exchange of stabilities occurs just before $ f_x / f_x^{\textit{sh}} = 0.995$ . The dashed grey lines correspond to the analytical approximation in the case $k_{\textit{sh}}\gg k_{\textit{RB}}$ , (4.3).

When elastic effects dominate the instability, the flow pattern still consists of a periodic structure of vortices (figure 4 a), but their typical length no longer scales with the distance $ d$ ; instead, it scales with the wavelength associated with the sheet. At the onset of instability, wrinkles with wavelength approximately $ 2\unicode{x03C0} /k_{\textit{sh}}$ appear on the surface of the sheet. Below each wrinkle, two counter-rotating vortices emerge in the flow, with a hill or a valley synchronised with the upstream or downstream flow of the vortices, respectively. The strength of these vortices decays to zero away from the sheet over length scale $ 1/k_{\textit{sh}}$ . This is evident from the centre of the vortices, which are much closer to the surface of the sheet than to the bottom wall. Another indication of this localisation may be seen in the eigenfunction $\varTheta (k_{\textit{sh}},y)$ plotted in figure 4(b). Clearly, the perturbed temperature is maximised close to the sheet, and decays to zero over a length proportional to the elastic length scale.

In contrast, when thermal effects dominate the system, i.e. when the global minimiser of the stability line corresponds to the Rayleigh–Bénard length scale $ k_{\textit{RB}}$ , the flow resembles that shown in figure 2(b). The only difference is that the upper plate forms small periodic undulations with wavelength $ 2\unicode{x03C0} /k_{\textit{RB}}$ (see figure 4 c). This additional flexibility of the system slightly reduces the critical Rayleigh number below $ R_a^{\textit{RB}}$ , but does not have any other prominent effect. In this case, the eigenfunction $ \varTheta (k_{\textit{RB}},y)$ exhibits a pattern similar to that seen in the Rayleigh–Bénard instability (compare figures 4(d) and 2(c)), with some discrepancies near $ y = 0$ due to the different boundary conditions in the two problems.

Figure 4. Eigenfunctions of the emerging pattern at the onset of instability. The parameters used are the same as those in figure 3. All eigenfunctions are obtained from the numerical solution of (3.2) and (3.3), and for visualisation they are normalised such that $Y_{\textit{sh}}=0.02$ . (a) The flow field at the critical Rayleigh number ( $R_a \simeq 965$ ) for $f_x/f_x^{\textit{sh}}=0.995$ (see figure 3). Here and in (c), arrows represent the streamlines, and colours represent the magnitude of the perturbed velocity. The flow forms a periodic structure of vortices, with two counter-rotating vortices making up each unit cell. The upward and downward streams of the vortices are synchronised with the maxima and minima points of the sheet’s wrinkle pattern. The magnitude of the fluid velocity gradually falls to zero away from the sheet. (b) The perturbed temperature under the conditions considered in (a) exhibits a localised profile that decays to zero away from the sheet. (c) The flow pattern at the critical Rayleigh number ( $R_a \simeq 1655$ ) when $f_x/f_x^{\textit{sh}}=0$ . The length scale of the emergent flow and wrinkle pattern scale with the distance $d$ . (d) The perturbed temperature under the same conditions exhibits a profile similar to that of the Rayleigh–Bénard set-up (see figure 4 c). In (b) and (d), solid lines correspond to (4.1), and open blue circles correspond to the numerical solution.

These results suggest that the perturbed temperature profile can be estimated by

(4.1) \begin{equation} \text{thermo-elastic instability} \quad \varTheta (k,y)=A \sin (\unicode{x03C0} y)+\frac {Y_{\textit{sh}}}{\sinh k}\sinh [k(y+1) ] , \end{equation}

where $A$ is as yet an unknown amplitude. As one can easily verify, this profile satisfies the boundary conditions at the bottom surface and on the sheet (3.3). Additionally, when the amplitude $ A$ approaches zero, the temperature exhibits a localised pattern, as expected when elastic effects dominate. However, if the first term in (4.1) is relatively larger than the second, then the temperature profile resembles the Rayleigh–Bénard solution (3.4), with small corrections due to the different boundary conditions. Therefore, in contrast to the thermal instability analysis (3.4), where the amplitude $A$ is, to leading order, independent of $k$ , here the leading-order approximation must depend on the wavenumber.

To obtain an analytic approximation of the solution, we follow the procedure described in § 3.1.1. First, we substitute the temperature profile (4.1) into (3.2a ), and solve for the vertical velocity $ U_y(k, y)$ , assuming $ \sigma = 0$ . (We use Mathematica for the symbolic analysis; see Wolfram-Research (2024).) Then we substitute this solution along with $ \varTheta (k, y)$ into (3.2b ), and satisfy the weak form of the equation by projecting it onto $ \sin (\unicode{x03C0} y)$ . This yields the amplitude $ A$ as a function of the wavenumber $ k$ and the Rayleigh number, as shown in figure 5. Indeed, when $ k \gg k_{\textit{RB}}$ , the amplitude $ A$ approaches zero, and the temperature becomes localised. However, near $ k_{\textit{RB}}$ , the amplitude is minimised, and the temperature profile closely resembles the Rayleigh–Bénard solution. Finally, we solve (3.2c ) using $ U_y(k, y)$ and $ A$ , which we have determined immediately above. This equation provides a relationship between the Rayleigh number and the system’s parameters.

The solid lines in the marginal stability diagram, figure 3, and the temperature profiles, figures 4(b) and 4(d), are based on this analytic approach. While the former marginal stability diagram shows very good agreement with the numerical data, some quantitative discrepancies appear in the temperature profile, primarily visible in figure 4(d). We attribute these differences to the limitations of the ansatz approach.

Since the analytical expressions for the Rayleigh number and the amplitude $A$ are cumbersome, we do not present them explicitly here. However, this analytic approach yields tractable predictions within certain limits, such as when the elastic wavenumber is much larger than or comparable to the thermal wavenumber. We explore these limits in the following subsections.

Figure 5. The amplitude $A/Y_{\textit{sh}}$ as a function of the wavenumber $k$ . When $k\gg 1$ , the amplitude decays to zero, and the temperature profile (4.1) exhibits a localised profile. However, close to the Rayleigh–Bénard wavenumber, $A/Y_{\textit{sh}}$ is minimised, and the temperature restores the solution given by figure 2(c).

4.1. The limit $k_{\textit{sh}}\gg k_{\textit{RB}}$

When $k_{\textit{sh}}/k_{\textit{RB}}\gg 1$ , i.e. $ (\rho _{\ell }^0 gd^4/B )^{1/4}\gg 1$ , the perturbed temperature is localised and differs from zero only close to the region of the sheet (figure 4 b). We can therefore neglect the amplitude $A$ in (4.1), and estimate the temperature solely by

(4.2) \begin{equation} \varTheta (k,y)=Y_{\textit{sh}}\frac {\sinh [k(y+1) ]}{\sinh k}. \end{equation}

The marginal stability curve is derived from this approximated profile by substituting (4.2) into (3.2a ) and solving for the perturbed velocity $ U_y(k, y)$ using the boundary conditions given in (3.3). Then the solution of (3.2c ) yields $R_a$ . In the limit $k\gg 1$ , this solution reads

(4.3) \begin{equation} R_a\simeq \frac {4R_g k}{3}\left (1-2\frac {k^2}{k_{\textit{sh}}^2}\frac {f_x}{f_x^{\textit{sh}}}+\frac {k^4}{k_{\textit{sh}}^4}\right )\!. \end{equation}

In figure 3, we plot this approximation (dashed grey line) alongside the numerical solution, and show that the approximation and the numerical solution agree well near the minimum arising from the elastic instability.

Therefore, we can exploit the above approximation to estimate the wavenumber and the critical Rayleigh number at which this minimum appears. To this end, we minimise (4.3) with respect to $k$ and obtain

(4.4) \begin{equation} \text{elastic minima} \quad k_{\textit{min}}^{\pm }/k_{\textit{sh}}=\left [\frac {3f_x}{5f_x^{\textit{sh}}}\pm \frac {3}{5}\left (\left (\frac {f_x}{f_x^{\textit{sh}}}\right )^2-\frac {5}{9}\right )^{1/2}\right ]^{1/2}. \end{equation}

Below $f_x/f_x^{\textit{sh}}\lt \sqrt {5}/3\simeq 0.74$ , the minimum does not exist, and the marginal stability line is minimised only at $k_{\textit{RB}}$ . However, above this critical value, two wavenumbers emerge: one increases ( $+$ ) as $ f_x/f_x^{\textit{sh}}$ increases, while the other decreases ( $-$ ) as $ f_x/f_x^{\textit{sh}}$ increases; see figure 6(a). Upon substitution of $ k_{\textit{min}}^{\pm }$ back into (4.3), it can be verified that the positive solution ( $+$ ) corresponds to a value of $ R_a$ that is lower than the solution with the minus sign. Therefore, we continue the analysis with the solution corresponding to the plus sign.

Figure 6. The wavenumber at the elastic minimum and an example of the stability diagram. (a) The wavenumbers $k_{\textit{min}}^{\pm }$ for the elastic minimum as a function of $f_x / f_x^{\textit{sh}}$ , as given by (4.4). No solution exists for $f_x / f_x^{\textit{sh}} \lt \sqrt {5} / 3$ . Above this value, two wavenumbers emerge: $k_{\textit{min}}^{+}$ (solid line) and $k_{\textit{min}}^{-}$ (dashed line). The former solution with the plus sign corresponds to a smaller value of $R_a$ . The dotted line represents the approximation given by (4.5). (b) The state diagram on the $ (f_x / f_x^{\textit{sh}}, R_a)$ plane, where the same parameters as those in figure 3 are used. The light blue region corresponds to stable states, where the fluid remains at rest and the sheet remains flat. The light yellow region indicates instability. For $ f_x \lt f_x^{{cr}}$ , where $ f_x^{{cr}} / f_x^{\textit{sh}} \simeq 0.99$ from (4.6), the transition between these regions is given approximately by $ R_a^{\textit{RB}}$ . When the temperature difference increases, the system progresses along a vertical line in this region, and the instability emerges with the wavenumber $ k_{\textit{RB}}$ . For $ f_x \gt f_x^{{cr}}$ , the critical Rayleigh number is given by (4.5) and depends on $ f_x$ . In this region, when the temperature difference exceeds the critical value of $R_a$ , the instability emerges with the wavenumber $ k_{\textit{sh}}$ .

Note that this approximated approach predicts a critical value of $ f_x / f_x^{\textit{sh}}$ that is slightly lower compared to the numerical solution in our example; see figure 3, where a minimum point appears when $ f_x / f_x^{\textit{sh}} \gtrsim 0.8$ . This discrepancy arises because $k_{\textit{sh}}$ is not significantly larger than $k_{\textit{RB}}$ in this example. Nonetheless, even in this case, the critical value $ f_x / f_x^{\textit{sh}} \simeq 0.74$ provides a lower bound of the lateral compression, only above which another minimum appears in the marginal stability curve.

This new minimum becomes the global minimiser only when $f_x/f_x^{\textit{sh}}$ gets closer to 1. In this region, we can estimate the wavenumber and the Rayleigh number at the minimum point by expanding (4.3) and (4.4) around $1-f_x/f_x^{\textit{sh}}\ll 1$ . This gives

(4.5) \begin{equation} \text{elastic minima} \ (f_x/f_x^{\textit{sh}}\rightarrow 1) \,\,\, k_{\textit{min}}^{+}/k_{\textit{sh}}\simeq 1-\frac {3}{4}\left (1\!-\!\frac {f_x}{f_x^{\textit{sh}}}\right )\!, \,\,\, R_a\simeq \frac {8}{3}R_g k_{\textit{sh}}\left (1\!-\!\frac {f_x}{f_x^{\textit{sh}}}\right )\!. \end{equation}

This approximation for the wavenumber is plotted in figure 6(a), and is shown to fit well with the more general solution when $f_x/f_x^{\textit{sh}}\gtrsim 0.9$ , which is sufficient for our current exchange of stability analysis. When $ f_x / f_x^{\textit{sh}}$ approaches 1, the wavenumber converges to $ k_{\textit{sh}}$ , and $ R_a$ approaches zero. Thus we recover the wrinkling instability. However, for any $ f_x / f_x^{\textit{sh}}$ slightly less than 1, the Rayleigh number takes on a finite value, which can be quite large, because $ R_g$ and $ k_{\textit{sh}}$ are considered as large numbers. The instability first occurs with $k_{\textit{sh}}$ when this finite $ R_a$ becomes smaller than the Rayleigh–Bénard minimum.

The Rayleigh–Bénard minimum occurs very close to $k_{\textit{RB}}$ , with $R_a \simeq R_a^{\textit{RB}}$ . Deviations from these values are primarily attributed to variations in the parameter $R_g$ . This is because in the marginal stability state ( $\sigma =0$ ) and under the limit $k \ll k_{\textit{sh}}$ , (3.2c ) simplifies to $R_g k^2 Y_{\textit{sh}} \simeq ({\textrm{d}}^3U_y/{\textrm{d}}y^3 )_{y=0}$ , and all parameters associated with the sheet vanish from the formulation. Within the range of our parameter space where $R_g \gg 1$ , it can be shown that deviations from the Rayleigh–Bénard values scale as $1/R_g$ and have only minor effects on the following analysis. Consequently, we choose to neglect them.

In summary, to approximate the critical force at which the elastic minimum becomes smaller than the Rayleigh–Bénard value, we equate (4.5) to $R_a^{\textit{RB}}$ , and obtain

(4.6) \begin{equation} \frac {f_x^{\textit{cr}}}{f_x^{\textit{sh}}}\simeq 1-\frac {3R_a^{\textit{RB}}}{8R_g k_{\textit{sh}}}. \end{equation}

Since both $ R_g$ and $ k_{\textit{sh}}$ are assumed to be large in this analysis, the elastic minimum becomes the global minimiser only when the lateral force approaches very close to its value in the static wrinkling instability.

Figure 6(b) illustrates the stability regions on the $ (f_x / f_x^{\textit{sh}}, R_a)$ plane, according to the parameters used in the example of figure 3. For $ f_x \lt f_x^{{cr}}$ , the marginal stability line is independent of the lateral force and occurs approximately when $ R_a$ exceeds $ R_a^{\textit{RB}}$ . In this region, the emergent wavenumber in the vicinity of the instability is $ k_{\textit{RB}}$ . However, for $f_x^{{cr}}\lt f_x \lt f_x^{\textit{sh}}$ , the transition occurs at an $ R_a$ value that is given by (4.5), with the wavenumber $k_{\textit{sh}}$ .

Figure 7. Investigation of marginal stability when $ k_{\textit{sh}} \simeq k_{\textit{RB}}$ . (a) The marginal stability curves on the $(k, R_a)$ plane for three different values of the lateral compression. A logarithmic scale is used on the vertical axis. The parameters $ R_g = 4500$ , $ \lambda = 5 \times 10^{-3}$ , $ P_r = 7$ and $ \varLambda = 60 \times 10^3$ , chosen are such that $ k_{\textit{sh}} \simeq 3.2$ . Open blue circles represent the numerical solution of (3.2) and (3.3) with $ \sigma = 0$ , and the corresponding lines (solid, dotted and dashed) show the analytical approximations. The triangles indicate the approximated minima of each curve, as given by (4.8). The insets show the perturbed temperature profiles at two different minima, where the open blue circles correspond to the numerical solution, and the solid lines represent (4.7). (b) Regions of stability on the $(f_x / f_x^{\textit{sh}}, R_a)$ plane, where the light blue region corresponds to stable states, and the light yellow region corresponds to unstable states. The marginal stability line is given by (4.8), with the same parameters as those used in (a).

4.2. The limit $k_{\textit{sh}}\simeq k_{\textit{RB}}$

By modifying the system’s parameters, for example, by using softer and thicker sheets, we can reduce the elastic wavenumber to the order of $k_{\textit{RB}}$ , i.e. $ (\rho _{\ell }^0 gd^4/B )^{1/4}\simeq 3.11$ . An example is shown in figure 7(a), where the sheet’s parameters are selected such that $ k_{\textit{sh}} \simeq 3.2$ , which is very close to the Rayleigh–Bénard wavenumber $ k_{\textit{RB}} \simeq 3.114$ . Marginal stability curves for several values of the lateral compression are plotted in the figure. When the lateral force approaches zero, the solution converges approximately to the Rayleigh–Bénard stability line. However, as compression increases, instability occurs at smaller values of $ R_a$ . In the limit $ f_x / f_x^{\textit{sh}} \rightarrow 1$ , the minimum value of $ R_a$ falls to zero, and the elastic instability is recovered.

Despite this significant change in the minimum of the Rayleigh number, the critical wavenumber remains nearly unchanged during the gradual variation in $ f_x$ . This invariance of $ k_{\textit{min}} \simeq k_{\textit{RB}}$ at the onset of instability allows us to derive some approximate analytical predictions, as we now show.

When the elastic and thermal wavenumbers coincide, it is no longer possible to estimate the perturbed temperature using only one term in (4.1). This is because the perturbed temperature is expected to recover the Rayleigh–Bénard solution when $f_x/f_x^{\textit{sh}} \simeq 0$ , and to converge to the localised solution when $f_x/f_x^{\textit{sh}} \simeq 1$ . Therefore, both terms must be included to describe the solution accurately. Keeping these two terms, and substituting $k=k_{\textit{RB}}$ in the perturbed temperature profile, we obtain

(4.7) \begin{equation} \varTheta (k_{\textit{RB}},y)/Y_{\textit{sh}}\simeq \frac {\tilde {A} R_a}{R_a-R_a^{\textit{RB}}}\sin (\unicode{x03C0} y)+\frac {\sinh [k_{\textit{RB}}(y+1)]}{\sinh k_{\textit{RB}}}, \end{equation}

where $\tilde {A}\simeq 0.29$ is a numerical factor obtained from the substitution of $k_{\textit{RB}}$ in the expression for $A$ . In addition, $ R_a$ corresponds to the Rayleigh number at $ k_{\textit{RB}}$ , which is yet to be determined as a function of the system’s parameters. When $ f_x / f_x^{\textit{sh}} \rightarrow 0$ , the Rayleigh number at the minimum point approaches $ R_a^{\textit{RB}}$ , and the first term in (4.7) dominates the second. Conversely, when $ f_x / f_x^{\textit{sh}} \rightarrow 1$ , the Rayleigh number approaches zero, and the second term becomes dominant. The insets in figure 7(a) depict these two scenarios at the minima associated with the two extreme values of the lateral force. The numerical eigenfunctions shown in these insets agree relatively well with the prediction given by (4.7).

Furthermore, the value of $ R_a$ at the minimum point of the stability line can be estimated by substituting $ k = k_{\textit{RB}}$ into the approximated analytical analysis. This gives

(4.8a) \begin{align} \text{for}\ k=k_{\textit{RB}}\simeq k_{\textit{sh}}, & \quad R_a=g_1\left (1-\sqrt {1-g_2/g_1^2}\right )\!, \end{align}

(4.8b) \begin{align} g_1=c_1+c_2\left (1-\frac {f_x}{f_x^{\textit{sh}}}\right )R_g, & \quad g_2=c_3\left (1-\frac {f_x}{f_x^{\textit{sh}}}\right )R_g, \end{align}

where $ c_1 \simeq 1362.7$ , $ c_2 \simeq 7.2$ and $ c_3 \simeq 24\,643.5$ are numerical factors derived by substituting $ k = k_{\textit{RB}}$ into the general solution. The triangles plotted near the minima of the stability curves in figure 7(a) are based on (4.8) and illustrate the agreement between this approximated solution and the numerical values.

Finally, figure 7(b) summarises the system’s stability regions on the $(f_x / f_x^{\textit{sh}}, R_a)$ plane for the case $ k_{\textit{sh}} \simeq k_{\textit{RB}}$ . For a given lateral compression, the system is stable below (4.8) and unstable above it. As $ f_x / f_x^{\textit{sh}} \rightarrow 0$ , the critical Rayleigh number converges approximately to the value predicted by the Rayleigh–Bénard solution. However, instability occurs at significantly smaller Rayleigh numbers in the opposite limit, $ f_x / f_x^{\textit{sh}} \rightarrow 1$ . In all cases, the emergent wavenumber at the onset of instability is approximately $ k_{\textit{RB}}$ , and the flow pattern is similar in shape to that shown in figure 4(c).

4.3. The limit $k_{\textit{sh}}\ll k_{\textit{RB}}$

Since the wavenumber $ k_{\textit{sh}}$ is inversely proportional to the thickness of the sheet, thicker sheets may, in principle, exhibit a much smaller wavenumber than $ k_{\textit{RB}}$ , i.e. $ (\rho _{\ell }^0 gd^4/B )^{1/4}\ll 1$ . One such example is shown in figure 8(a), where the system parameters are chosen such that $ k_{\textit{sh}} \simeq 0.28$ . As in previous cases, the marginal stability line always exhibits a local minimum at approximately $ k_{\textit{RB}}$ . Instability occurs at this minimum only if the lateral force decreases below a certain threshold. Above this threshold, instability develops with a wavenumber $ k_{\textit{sh}}$ .

The elastic minimum can be further quantified using our analytical approach. Away from the Rayleigh–Bénard minimum, the perturbed temperature is primarily described by the second term in (4.1). In this region, where $ k \ll 1$ , the perturbed temperature can be estimated solely by

(4.9) \begin{equation} \varTheta (k,y)/Y_{\textit{sh}}= \frac {\sinh [k(y+1)]}{\sinh k}\simeq y+1. \end{equation}

The inset in figure 8(b) shows that this approximated profile matches well with the numerical solution at the elastic minimum. To obtain the marginal stability line, we substitute (4.9) in (3.2a ), and solve for the perturbed velocity $U_y(k,y)$ . Then we substitute this solution in (3.2c ), and solve for $R_a$ . To fourth order in $k$ we obtain

(4.10) \begin{equation} R_a=\frac {20R_g}{7}\left (1-2\frac {k^2}{k_{\textit{sh}}^2}\frac {f_x}{f_x^{\textit{sh}}}+\frac {k^4}{k_{\textit{sh}}^4}\right )\!. \end{equation}

This profile is plotted in figure 8(b) (dashed grey lines) and shows good agreement with the numerical solution near the elastic minimum. Equation (4.10) implies that in addition to the Rayleigh–Bénard minimum, the marginal stability line is also minimised at

(4.11) \begin{equation} k_{\textit{min}}=\left (\frac {f_x}{f_x^{\textit{sh}}}\right )^{1/2}k_{\textit{sh}}, \quad R_a=\frac {20R_g}{7}\left [1-\left (\frac {f_x}{f_x^{\textit{sh}}}\right )^2\right ]\! . \end{equation}

When $ f_x/f_x^{\textit{sh}} \to 0$ , the wavenumber $ k_{\textit{min}}$ falls to zero, and the minimum occurs at $ R_a \simeq 20R_g/7$ . Although this solution arises from our analysis of the elastic minimum, it is not caused by bending effects, since $ \varLambda k^4 \to 0$ when the wavenumber is zero. Instead, it results solely from gravitational waves that modify the surface of the sheet. Similar phenomena, although in systems involving Marangoni effects, are discussed by Scriven & Sternling (Reference Scriven and Sternling1964), Smith Reference Smith1966 and Davis (Reference Davis1987). While this additional extremum of the stability curve can become the global minimum if $ R_g$ is made sufficiently small, we do not explore the details of this transition here, for two reasons: first, our primary focus is on the interplay between elasticity and thermal effects; and second, this exchange of stabilities occurs at $ R_g$ values that are relatively small and thus lie outside the parameter range considered in this study.

As $ f_x/f_x^{\textit{sh}}$ increases, the wavenumber at the minimum point also increases, while the corresponding $ R_a$ decreases. This minimum becomes smaller than the $ R_a^{\textit{RB}}$ solution when

(4.12) \begin{equation} \frac {f_x^{{cr}}}{f_x^{\textit{sh}}}=\sqrt {1-\frac {7R_a^{\textit{RB}}}{20R_g}}. \end{equation}

If $ R_g / R_a^{\textit{RB}} \gg 1$ , then instability emerges only when the lateral compression approaches close to the wrinkling instability value. However, when $ R_g$ is of the same order as $ R_a^{\textit{RB}}$ , the instability occurs at a lower compression.

The stability diagram in the $ (f_x / f_x^{\textit{sh}}, R_a)$ plane, resulting from our example ( $ k_{\textit{sh}} \simeq 0.28$ ), is plotted in figure 8(c) and shows the transition between thermal and elastic instabilities. Below $ f_x / f_x^{\textit{sh}} \simeq 0.93$ from (4.12), instability occurs at $ R_a \simeq R_a^{\textit{RB}}$ with wavenumber $ k_{\textit{RB}}$ . However, above this value, instability occurs according to (4.11). As $ f_x / f_x^{\textit{sh}} \to 1$ , the instability emerges at vanishingly small Rayleigh numbers.

Figure 8. Analysis of the system’s stability when $ k_{\textit{sh}} \ll k_{\textit{RB}}$ . In all plots, the parameters $ R_g = 4500$ , $ \lambda = 0.05$ , $ P_r = 7$ and $ \varLambda = 10^8$ are used such that $ k_{\textit{sh}} \simeq 0.28$ . Additionally, open blue circles represent the numerical solution of (3.2) and (3.3) with $ \sigma = 0$ . (a) The marginal stability curve on the $ (k, R_a)$ plane for three different values of the lateral compression. When $ f_x / f_x^{\textit{sh}} \leqslant 0.93$ , the global minimum is close to $ k_{\textit{RB}}$ ; however, above this critical value, the global minimum occurs at $ k_{\textit{sh}}$ . (b) A zoomed-in view of the region near $ k_{\textit{sh}}$ , indicated by the grey ellipse in (a). The inset presents the perturbed temperature at the minimum wavenumber, where the solid line corresponds to (4.9), and open blue circles correspond to the numerical data. (c) Stability diagram on the $ (f_x/f_x^{\textit{sh}}, R_a)$ plane, where stable and unstable regions are denoted by light blue and light yellow regions, respectively.

5. Summary and conclusions

We investigated the effect of the interplay between thermal gradients, elastic deformations and fluid dynamics on the onset of natural convection in a system that couples a pre-compressed sheet and a fluid layer. By adding elastic effects to the classical Rayleigh–Bénard convection, we obtained new stability thresholds and emergent patterns that depend on both the system’s mechanical and thermal parameters.

Given the range of parameters investigated in this paper, we conclude that the strength of the lateral compression essentially determines the emergent pattern at the onset of instability. When the lateral compression diminishes to zero, the sheet has only a minor effect on the instability, i.e. the wavenumber of the emergent pattern is very close to $k_{\textit{RB}}$ , and the critical Rayleigh number reduces to a value only slightly below $R_a^{\textit{RB}}$ . However, when the lateral force exceeds a certain threshold, the stability diagram changes considerably as a new local minimum appears in the $(k, R_a)$ diagram, with a corresponding wavenumber that scales as $k_{\textit{sh}}$ . This minimum becomes the global minimiser when the lateral force surpasses a second threshold, typically only slightly below $f_x^{\textit{sh}}$ . The emergent pattern in this new minimum exhibits a flow field that is synchronised with the wrinkles, and its wavelength scales with that of the wrinkles. In contrast to the Rayleigh–Bénard solution, the pattern emerging from this new minimum is almost independent of the fluid depth.

Based on these observations, we derived an approximate solution for the problem using an ansatz approach for the perturbed temperature (see (4.1)). This solution allows us to investigate more quantitatively the critical threshold at which the new solution sets in. For instance, we showed that when $k_{\textit{sh}} \gg k_{\textit{RB}}$ , the instability occurs when $f_x^{{cr}} / f_x^{\textit{sh}}$ is very close to unity (see (4.6)), with deviations from this value scaling as $R_a^{\textit{RB}} / k_{\textit{sh}} R_g$ . Therefore, if the wrinkling wavelength is much smaller than the fluid depth, the new minimum pre-empts the Rayleigh–Bénard solution only very close to the static instability threshold. Conversely, when $k_{\textit{sh}} \ll k_{\textit{RB}}$ , the critical force $f_x^{{cr}} / f_x^{\textit{sh}}$ can significantly deviate from unity (see (4.12)).

We emphasise that our analysis is limited to determining the stability of the base solution, where the sheet remains flat and the fluid is at rest. However, it does not address the question of whether the wrinkle pattern emerging after the onset of instability stabilises. Recent studies on wrinkle formation in a static set-up, i.e. without considering the sheet’s inertia, have shown that wrinkles become unstable against a localised pattern known as a fold, soon after forming. More precisely, beyond the wrinkling instability, the lateral force decreases as the wrinkle amplitude increases. The wrinkles-to-fold transition occurs at $f_x^{\textit{sh}} - f_x \sim 1/L$ , where $L$ is the total length of the sheet (Diamant & Witten Reference Diamant and Witten2011; Oshri et al. Reference Oshri, Brau and Diamant2015; Oshri & Diamant Reference Oshri and Diamant2017). Consequently, the stability region for static wrinkles diminishes to zero as the system size increases. Nonetheless, it is also known that inertial effects can aid in stabilising patterns that are statically unstable, as shown, for example, in cases of a collapsed ring, sheets on viscoelastic substrates, and snapping beams (Box et al. Reference Box, Kodio, O’Kiely, Cantelli, Goriely and Vella2020; Guan et al. Reference Guan, Nguyen, Pocivavsek, Cerda and Velankar2023; Oshri et al. Reference Oshri, Goncharuk and Feldman2024). Whether wrinkles stabilise in our system due to thermal and dynamic effects remains an open question. Analysing this stabilisation requires nonlinear extensions of the linear stability theory discussed here.

The mechanical properties, such as Young’s modulus, of very thin materials at micrometric or even nanometric scales are typically not measured directly. Instead, in some cases, they are inferred indirectly through wrinkle experiments. In these experiments, a thin sheet is placed on a solid or even fluid substrate, and compressive forces are applied, causing wrinkles to form on the sheet’s surface. Theoretical models that relate the wrinkles’ wavelength to a material’s elastic properties are used to deduce these properties (Cerda & Mahadevan Reference Cerda and Mahadevan2003; Paulsen et al. Reference Paulsen, Hohlfeld, King, Huang, Qiu, Russell, Menon, Vella and Davidovitch2016). In this context, our analysis offers an alternative method for inferring the material stiffness of thin sheets by inducing a temperature difference on the underlying substrate. Nonetheless, it is important to note that this analysis is limited to the Boussinesq approximation, and in practical set-ups, non-Boussinesq effects may significantly impact the results even at relatively small temperature differences (Valori et al. Reference Valori, Elsinga, Rohde, Tummers, Westerweel and van der Hagen2017; Demou & Grigoriadis Reference Demou and Grigoriadis2019).

From a broader perspective, the analysis of our relatively simple set-up can provide valuable insights into the analytical understanding of more complex systems that require precise control over convection-driven fluid–structure interactions, such as valves, energy-harvesting devices and fluid mixers.