Exact energy stability of B\'enard-Marangoni convection at infinite Prandtl number

TL;DR

This paper uses the energy method to precisely determine the conditions under which pure conduction remains stable in Bénard-Marangoni convection at infinite Prandtl number, revealing that energy stability aligns with linear instability thresholds only for insulating boundaries.

Contribution

The study provides an exact solution to the energy stability problem for Bénard-Marangoni convection, extending previous results and clarifying boundary condition effects on stability.

Findings

Energy stability proven up to linear instability threshold for insulating boundaries.

Energy method does not exclude subcritical instabilities in non-insulating cases.

Exact stability conditions depend on boundary thermal properties.

Abstract

Using the energy method we investigate the stability of pure conduction in Pearson's model for B\'enard-Marangoni convection in a layer of fluid at infinite Prandtl number. Upon extending the space of admissible perturbations to the conductive state, we find an exact solution to the energy stability variational problem for a range of thermal boundary conditions describing perfectly conducting, imperfectly conducting, and insulating boundaries. Our analysis extends and improves previous results, and shows that with the energy method global stability can be proven up to the linear instability threshold only when the top and bottom boundaries of the fluid layer are insulating. Contrary to the well-known Rayleigh-B\'enard convection setup, therefore, energy stability theory does not exclude the possibility of subcritical instabilities against finite-amplitude perturbations.

Tables

Table 1. Table 1: Minimum critical Marangoni number for energy stability at infinite Prandtl number for selected values of the Biot number L in the extremal cases B = 0 B 0 \mbox{{B}}=0 (insulating bottom boundary) and B = ∞ B \mbox{{B}}=\infty (perfectly conducting bottom boundary). Where available, the corresponding values for linear stability  (Pearson, 1958 ) and energy stability at finite \Pran \Pran \Pran   (Davis, 1969 ) are also reported.

				\pbox4cmEnergy stability
at $\text{Pran}=\infty$		\pbox4cmLinear stability
(Pearson, 1958)		\pbox4cmEnergy stability
(Davis, 1969)
B		L		${\text{M}}_{\mathrm{cr}}$	$k_{\mathrm{cr}}$		${\text{M}}_{\mathrm{cr}}$	$k_{\mathrm{cr}}$		${\text{M}}_{\mathrm{cr}}$	$k_{\mathrm{cr}}$
0		0		48	0		48	0		—	—
0		10		210.9	2.65		383.2	2.45		—	—
0		100		769.8	2.87		3132	2.71		—	—
$\infty$		0		66.84	2.14		79.61	1.99		56.77	—
$\infty$		10		227.3	2.91		413.4	2.74		180.7	—
$\infty$		100		821.6	2.98		3304	3.10		—	—

Equations

$${\nabla }^2\mathbf{u}$$

$$\text{bnabla}\text{bcdot}\mathbf{u}$$

$$\frac{\partial \theta }{\partial t}+\mathbf{u}\text{bcdot}\text{bnabla}\theta$$

$${\left[\frac{\partial u}{\partial z}+\text{mbox}M\frac{\partial \theta }{\partial x}\right]}_{z=1}=0,$$

$${\left[\frac{\partial \theta }{\partial z}-\text{mbox}B\theta \right]}_{z=0}$$

$$\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}⟨{\theta }^2⟩=-⟨{\left|\text{bnabla}\theta \right|}^2-w\theta ⟩-\text{mbox}L\overline{{\theta }^2}(1)-\text{mbox}B\overline{{\theta }^2}(0).$$

$$⟨{\left|\text{bnabla}\theta \right|}^2-w\theta ⟩+\text{mbox}L\overline{{\theta }^2}(1)+\text{mbox}B\overline{{\theta }^2}(0)\ge 0.$$

$$f_k(z)=\frac{k\sinh k}{\sinh (2k)-2k}\left[kz\cosh (kz)-\sinh (kz)+(1-k\coth k)z\sinh (kz)\right],$$

$$⟨{\left|\text{bnabla}\theta \right|}^2-w\theta ⟩+\text{mbox}L\overline{{\theta }^2}(1)+\text{mbox}B\overline{{\theta }^2}(0)=2\sum _{k\ge 0}{\mathcal{F}}_k\{{\hat{\theta }}_k\},$$

$$\begin{gathered} {\mathcal{F}}_k\{{\hat{\theta }}_k\}:={\int }_0^1\left\{|{\hat{\theta }}_k^{\prime }(z){|}^2+k^2|{\hat{\theta }}_k(z){|}^2+\text{mbox}M{f}_k(z)\text{Real}\left[{{\hat{\theta }}_k}^{\ast }(z){\hat{\theta }}_k(1)\right]\right\}\mathrm{d}z \\ +\text{mbox}L|{\hat{\theta }}_k(1){|}^2+\text{mbox}B|{\hat{\theta }}_k(0){|}^2. \end{gathered}$$

$${\mathcal{F}}_k\{{\hat{\theta }}_k\}\ge 0$$

$${{\hat{\theta }}_k}^{\prime }(0)-\text{mbox}B{\hat{\theta }}_k(0)$$

$${\Gamma }_0:=\left\{v(z):{\int }_0^1\left(|v^{\prime }{|}^2+|v{|}^2\right)\mathrm{d}z<\infty ,v^{\prime }(0)=\text{mbox}Bv(0),v^{\prime }(1)=-\text{mbox}L,v(1)=1\right\}.$$

$${\mathcal{Q}}_k^{\star }:=\underset{v\in {\Gamma }_0}{\inf }{\int }_0^1\left[|v^{\prime }(z){|}^2+k^2|v(z){|}^2+\text{mbox}M{f}_k(z)v(z)\right]\mathrm{d}z+\text{mbox}L+\text{mbox}B|v(0){|}^2\ge 0.$$

$${\mathcal{Q}}_k\{v\}={\int }_0^1\left[|v^{\prime }(z){|}^2+k^2|v(z){|}^2+\text{mbox}M{f}_k(z)v(z)\right]\mathrm{d}z+\text{mbox}L+\text{mbox}B|v(0){|}^2$$

$$v_{\star }^{\prime \prime }(z)-k^2{v}_{\star }(z)=\frac{1}{2}\text{mbox}M{f}_k(z),$$

$${\Gamma }_1:=\left\{v(z):{\int }_0^1\left(|v^{\prime }{|}^2+|v{|}^2\right)\mathrm{d}z<\infty ,v^{\prime }(0)=\text{mbox}Bv(0),v(1)=1\right\}.$$

$$v_{\star }(z)=\text{mbox}M{g}_k(z)+h_k(z),$$

$$v_n(z):=\{\begin{array}{ll}{v}_{\star }(z)& \text{if}0\le z\le {\xi }_n,\\ v_{\star }({\xi }_n)+\frac{2+\text{mbox}L(1-{\xi }_n)-2v_{\star }({\xi }_n)}{1-{\xi }_n}\left(z-{\xi }_n\right)& \text{if}{\xi }_n\le z\le \frac{1+{\xi }_n}{2},\\ 1+\text{mbox}L(1-z),& \text{if}\frac{1+{\xi }_n}{2}\le z\le 1.\end{array}$$

$${\int }_0^1\left[|v_{\star }^{\prime }(z){|}^2+k^2|v_{\star }(z){|}^2\right]\mathrm{d}z=v_{\star }^{\prime }(1)-\text{mbox}B|v_{\star }(0){|}^2-\frac{1}{2}\text{mbox}M{\int }_0^1{f}_k(z)v_{\star }(z)\mathrm{d}z.$$

$${\mathcal{Q}}_k^{\star }=\underset{=:{\alpha }_k}{\underbrace{\left[\frac{1}{2}{\int }_0^1{g}_k(z)f_k(z)\mathrm{d}z\right]}}\text{mbox}{M}^2+\underset{=:{\beta }_k}{\underbrace{\left[{g_k}^{\prime }(1)+\frac{1}{2}{\int }_0^1{h}_k(z)f_k(z)\mathrm{d}z\right]}}\text{mbox}M+{h_k}^{\prime }(1)+\text{mbox}L.$$

$$\text{mbox}{M}_e=\frac{-{\beta }_k-\sqrt{{\beta }_k^2-4{\alpha }_k\left[h_k^{\prime }(1)+\text{mbox}L\right]}}{2{\alpha }_k}.$$

$$\underset{v}{\inf }{\mathcal{F}}_k\{v\}\ge 0,$$

$${\int }_0^1|v(z){|}^2\mathrm{d}z=1.$$

$${\mathcal{L}}_k\{v,\lambda \}:={\mathcal{F}}_k\{v\}+\lambda \left({\int }_0^1|v(z){|}^2\mathrm{d}z-1\right).$$

$$\frac{\delta {\mathcal{L}}_k}{\delta v}:=\underset{\epsilon \to 0}{\lim }\frac{{\mathcal{L}}_k\{v+\epsilon h,\lambda \}-{\mathcal{L}}_k\{v,\lambda \}}{\epsilon }=0$$

$$\begin{gathered} \frac{\delta {\mathcal{L}}_k}{\delta v}={\int }_0^1\left[-2v^{\prime \prime }(z)+2\left(k^2+\lambda \right)v(z)+\text{mbox}M{f}_k(z)v(1)\right]h(z)\mathrm{d}z \\ +\left[\text{mbox}Lv(1)+\text{mbox}M{\int }_0^1{f}_k(z)v(z)dz\right]h(1)+\text{mbox}Bv(0)h(0). \end{gathered}$$

$$-2v^{\prime \prime }(z)+2\left(k^2+\lambda \right)v(z)+\text{mbox}M{f}_k(z)v(1)=0,$$

$$\text{mbox}Lv(1)+\text{mbox}M{\int }_0^1{f}_k(z)v(z)dz$$

$$v^{\prime }(0)-\text{mbox}Bv(0)$$

Full text

1 Introduction

Bénard–Marangoni convection describes the motion of a layer of fluid driven by shear stresses due to gradients in surface tension at the interface between the fluid and its surroundings. This type of convection arises in numerous engineering applications, including the growth of crystals in semiconductors (Schatz & Neitzel, 2001), cladding processes (Kumar & Roy, 2009), and drying of thin polymer films (Yiantsios et al., 2015), and has recently received increasing attention as a paradigm for shear-driven turbulent transport processes (Boeck & Thess, 1998, 2001; Hagstrom & Doering, 2010).

The first mathematical model of surface-tension-driven convection was proposed by Pearson (1958), who showed that pure conduction is linearly unstable when the Marangoni number M, a non-dimensional measure of the surface tension effects, exceeds a critical threshold $\mbox{{M}}_{\!l}$, independently of the fluid’s Prandtl number \Pran (the ratio of the fluid’s kinematic viscosity to its thermal diffusivity). Subsequently, Davis (1969) used the energy method to prove that conduction is asymptotically stable against disturbances of arbitrary amplitude when M is smaller than a critical value $\mbox{{M}}_{\!e}$, also independently of \Pran. In contrast to Rayleigh–Bénard convection, the linear and energy thresholds $\mbox{{M}}_{\!l}$ and $\mbox{{M}}_{\!e}$ do not coincide, allowing the possibility of subcritical instabilities.

Although Davis’s analysis and computations yield the best global stability boundary that can be attained with the energy method for a fluid of finite \Pran, they can be improved in the infinite Prandtl number case. This limit is an attractive model for high-Prandtl-number fluids, such as the silicon oils used in experiments (de Bruyn et al., 1996) or Earth’s mantle (Jones, 1977), because it gives accurate quantitative predictions whilst simplifying the governing equations (Boeck & Thess, 2001). The key observation is that in the limit of infinite $\Pran$ the inertial term in the momentum equation can be dropped, and as a result the velocity field can be “slaved” to the temperature field (see e.g. Hagstrom & Doering, 2010), allowing the formulation of an improved variational principle for energy stability. This variational principle, first considered by Hagstrom & Doering (2010), is interesting from a mathematical perspective because it requires the minimisation of a quadratic functional that depends explicitly on the boundary values of the argument function and, as we will show, whose Euler–Lagrange differential equation is overconstrained. Hagstrom & Doering bypassed this difficulty by applying elementary functional estimates to the quadratic functional directly, and raised the energy stability boundary $\mbox{{M}}_{\!e}$ from $56.77$ to $58.36$ in the case of a perfectly conducting bottom boundary.

The main contribution of this work is to show that further improvements are possible. We consider an extended version of Hagstrom & Doering’s infinite-\Pran energy stability variational principle—which applies to Pearson’s model of Bénard–Marangoni convection with perfectly conducting, imperfectly conducting, and insulating bottom boundaries—and compute the optimal energy stability boundary by extending the domain of the variational problem in such a way that the Euler–Lagrange equation admits a unique solution.

The rest of this work is organised as follows. Section 2 reviews Pearson’s model for Bénard–Marangoni convection at infinite Prandtl number. We formulate the variational principle for energy stability in §3, and derive an exact solution in §4. Further remarks and suggestions for future investigations are offered in §5.

2 Pearson’s model

Consider a layer of fluid of density $\rho$, kinematic viscosity $\nu$, thermal diffusivity $\kappa$, and thermal conductivity $\lambda$, bounded by two non-deformable surfaces at $z=0$ and $z=h$. When the fluid is at rest and heat is transported by conduction alone, the temperature of the fluid is given by $T(z)=T_{0}-(Q_{0}/\lambda)z$, where $T_{0}$ is the temperature of the bottom boundary and $Q_{0}$ is the imposed heat flux through the layer per unit area. We choose the temperature scale such that $T_{0}=0$, and make the system non-dimensional by using $h$, $h^{2}/\kappa$, and $Q_{0}h/\lambda$ as the characteristic length, time, and temperature units. For simplicity, we work in two dimensions and denote the non-dimensional position vector by $\boldsymbol{x}=x\boldsymbol{i}+z\boldsymbol{k}$; the model and all our results extend with no modifications to three dimensions as described by Hagstrom & Doering (2010).

At infinite Prandtl number, non-dimensional velocity, pressure, and temperature disturbances to the conductive state—denoted by $\boldsymbol{u}(\boldsymbol{x},t)=u(\boldsymbol{x},t)\boldsymbol{i}+w(\boldsymbol{x},t)\boldsymbol{k}$, $p(\boldsymbol{x},t)$, and $\theta(\boldsymbol{x},t)$—evolve according to

[TABLE]

We assume that all variables are periodic in the horizontal ($x$) direction with period $\Lambda$, or that their Fourier transform exists. We impose the no-slip condition $\boldsymbol{u}|_{z=0}=0$ at the bottom boundary, and the impenetrability condition $w|_{z=1}=0$ at the top boundary. Moreover, if $\gamma$ denotes the negative of the derivative of surface tension with respect to surface temperature, the balance of surface stresses and tension forces is expressed by the boundary condition

[TABLE]

where the Marangoni number $\mbox{{M}}=\gamma Q_{0}h^{2}/(\lambda\rho\nu\kappa)$ is the main governing parameter of the system. Finally, letting $q_{\rm bot}$ and $q_{\rm top}$ denote the derivative of the outward heat fluxes through the top and bottom surfaces with respect to the surface temperature, balancing the heat fluxes through boundaries requires that

[TABLE]

where the Biot numbers $\mbox{{B}}=q_{\rm bot}h/\lambda$, $\mbox{{L}}=q_{\rm top}h/\lambda$ describe the conductivity of the boundaries. (The sign difference between the two boundaries is due to the convention that outward heat flux is positive). We consider $\mbox{{B}},\mbox{{L}}\geq 0$, a reasonable assumption because an increase in the fluid’s surface temperature should raise the heat flux to the surroundings; the case $B=L=0$ corresponds to perfectly insulating boundaries, while the perfectly conducting case corresponds to the (formal) choice $\mbox{{B}}=\mbox{{L}}=\infty$. For a comprehensive discussion of the thermal boundary conditions (3) we refer the reader to the original work by Pearson (1958).

3 Energy stability analysis

Stability analysis via the energy method relies on the simple observation that stationary conduction (i.e. when the fluid is at rest) is stable if the kinetic energy of a temperature perturbation does not increase in time, irrespective of the perturbation’s initial amplitude. The evolution equation for the average kinetic energy $\langle\theta^{2}\rangle/2$ of a temperature perturbation, where $\langle\cdot\rangle$ denotes the usual volume average, is found by averaging $\theta\times$(1c) and integrating by parts using (1b) and the boundary conditions (3) to arrive at

[TABLE]

In this equation and throughout the rest of this section, overlines denote horizontal averages. Clearly, the kinetic energy of the perturbation $\theta$ does not increase in time if the right-hand side of (4) is non-positive at each instant in time, i.e.,

[TABLE]

Upon substituting the horizontal Fourier series expansion of $\theta$ and $w$ into (5) and into the boundary conditions (3) (we consider the case of a finite periodic domain for definitess; similar arguments hold for the infinite domain if the Fourier series is replaced by the Fourier transform), dropping the time dependence, and recalling from Hagstrom & Doering (2010) that the Fourier amplitudes of the velocity perturbation $w$ are “slaved” to those of $\theta$ according to $\hat{w}_{k}(z)=-\mbox{{M}}\,f_{k}(z)\,\hat{\theta}_{k}(1)$, where

[TABLE]

we can rewrite

[TABLE]

where the sum is over all positive wavenumbers and

[TABLE]

(Here and in the following, ∗ denotes complex conjugation and primes denote total differentiation with respect to $z$.)

Since among all possible perturbations are those defined by a single wavenumbers, we conclude that a necessary and sufficient condition for the global stability of Bénard–Marangoni conduction is that, for all wavenumbers $k\geq 0$,

[TABLE]

for all complex-valued perturbation Fourier amplitudes $\hat{\theta}_{k}$ (hereafter simply referred to as perturbations) that satisfy

[TABLE]

In fact, we may restrict our attention to real-valued $\hat{\theta}_{k}$ because the real and imaginary parts give identical and independent contributions to the left-hand side of (9). Moreover, note that (9) holds trivially when $\hat{\theta}_{k}(1)=0$, and that its left-hand side is homogeneous quadratic in $\hat{\theta}_{k}$. Since the boundary conditions in (10) are also homogeneous, we can further restrict our attention to the perturbations that satisfy the normalisation condition $\hat{\theta}_{k}(1)=1$ and, without any loss of generality, we may replace (3) with the boundary conditions

[TABLE]

Putting these observations together, we define the space of admissible perturbations as

[TABLE]

Finally, it is clear that (9) holds if and only if the infimum of its left-hand side over all admissible perturbation fields is non-negative. The stability of Bénard–Marangoni conduction at given Marangoni and Biot numbers M, B, and L is then established if we can prove that, for all wavenumbers $k$,

[TABLE]

In particular, for fixed values of the Biot numbers B and L we can compute the energy stability boundary in the M–$k$ space—i.e., the largest Marangoni number for which a perturbation of wavenumber $k$ is stable—by solving the variational problem for the infimum $\mathcal{Q}_{k}^{\star}$ as a function of the Marangoni number M for each $k$ and choosing the largest M for which $\mathcal{Q}_{k}^{\star}\geq 0$ .

4 Solution of the variational problem

As we have anticipated in §1, the variational problem for $\mathcal{Q}_{k}^{\star}$ is interesting from the mathematical point of view because the infimum of the quadratic form

[TABLE]

is not generally attained by any test function $v\in\Gamma_{0}$. In fact, a straightforward application of the calculus of variations (see e.g. Courant & Hilbert, 1953; Giaquinta & Hildebrandt, 1996) shows that any candidate minimiser $v_{\star}$ must satisfy the second-order Euler–Lagrange differential equation

[TABLE]

subject to the three boundary conditions $v_{\star}^{\prime}(0)=\mbox{{B}}v_{\star}(0)$, $v_{\star}^{\prime}(1)=-\mbox{{L}}$, and $v_{\star}(1)=1$. This problem is over-constrained, and admits no solution (with the possible exception of selected values of B, L, M and $k$).

It should be noted that the lack of a minimiser for $\mathcal{Q}_{k}$ is not due to our normalisation convention for the test functions, which is the source of the extra boundary condition $v_{\star}(1)=1$. When a different normalisation is used, in fact, the minimisation of $\mathcal{Q}_{k}$ over $\Gamma_{0}$ is replaced with the minimisation of $\mathcal{F}_{k}$ in (8) over all normalised test functions that satisfy the boundary conditions (10). As we demonstrate in appendix A for the commonly used normalisation $\int_{0}^{1}|v(z)|^{2}\,\mathrm{d}z=1$, the Euler–Lagrange equations for the minimiser of $\mathcal{F}_{k}$ are over-constrained by so-called “natural conditions” that arise when setting to zero its first variation (Courant & Hilbert, 1953, Chapter IV, Section 5.1).

This obstacle is overcome if we can enlarge the space of test functions in such a way that (15) has a unique solution $v_{\star}$, and moreover $\mathcal{Q}_{k}\{v_{\star}\}=\mathcal{Q}_{k}^{\star}$. This is indeed the case if we drop the boundary condition $v^{\prime}(1)=-\mbox{{L}}$, and minimise $\mathcal{Q}_{k}\{v\}$ over the larger space of functions

[TABLE]

Having removed one boundary condition, in fact, the Euler–Lagrange equation (15) becomes a standard second-order inhomogeneous ordinary differential equation and it can be solved analytically. For each wavenumber $k$ the solution can be written in the form

[TABLE]

where $g_{k}(z)$ and $h_{k}(z)$ are two known smooth functions whose expressions, given in Appendix B, depend on the Biot number B. Furthermore, to see that $\mathcal{Q}_{k}^{\star}=\mathcal{Q}_{k}\{v_{\star}\}$ we note that on one hand we must have $\mathcal{Q}_{k}\{v_{\star}\}\leq\mathcal{Q}_{k}^{\star}$, because $\Gamma_{0}$ is a proper subset of $\Gamma_{1}$ and $v_{\star}$ minimises $\mathcal{Q}_{k}\{v\}$ over $\Gamma_{1}$. On the other hand, $\mathcal{Q}_{k}\{v_{\star}\}\geq Q_{k}^{\star}$ because we can find a sequence of functions $(v_{n})_{n\geq 1}$ with $v_{n}\in\Gamma_{0}$ such that $\mathcal{Q}_{k}\{v_{n}\}$ converges to $\mathcal{Q}_{k}\{v_{\star}\}$; for example, in Appendix C we show that this is the case if we let $\xi_{n}=n/(n+1)$ and take

[TABLE]

Note that these test functions are simply piecewise-linear continuous functions on $[\xi_{n},1]$ satisfying $v_{n}(\xi_{n})=v_{\star}(\xi_{n})$, $v_{n}(1)=1$, and $v_{n}^{\prime}(1)=-L$, and that they could be smoothed around the corner points without changing their boundary values and derivatives to meet any regularity requirements prescribed on the space $\Gamma_{0}$.

Having computed the minimizer $v_{\star}$, we now turn to the computation of the minimum $\mathcal{Q}_{k}^{\star}=\mathcal{Q}_{k}\{v_{\star}\}$. To simplify the analysis, we integrate $v_{\star}\times$(15) by parts using the boundary conditions on $\Gamma_{1}$ to show that

[TABLE]

Upon combining this with (14) and (17) we find

[TABLE]

For each wavenumber $k$ and given Biot numbers B and L, the infimum $\mathcal{Q}_{k}^{\star}$ is a quadratic form of the Marangoni number M, and the coefficients $\alpha_{k}$ and $\beta_{k}$ have explicit expressions ($f_{k}$, $g_{k}$ and $h_{k}$ are known, and their products can be integrated analytically). These are too long to be reported, but $\alpha_{k}$ and $\beta_{k}$, together with $h_{k}^{\prime}(1)$, are plotted in figure 1 for Biot numbers $\mbox{{B}}=0$ (corresponding to a perfectly insulating bottom boundary), $\mbox{{B}}=1$, $\mbox{{B}}=10$, and $\mbox{{B}}=\infty$ (corresponding to a perfectly conducting bottom boundary). Note that the leading order coefficient $\alpha_{k}$ is negative for all $k$, and it must be so because perturbations at any wavenumber $k$ eventually become linearly unstable (Pearson, 1958), implying that $\mathcal{Q}_{k}^{\star}<0$ for all sufficiently large M. Consequently, the largest Marangoni number $\mbox{{M}}_{\!e}$ at which Bénard–Marangoni conduction at infinite Prandtl number is stable against perturbations of wavenumber $k$ and arbitrary amplitude is given by the largest root of the quadratic form in (20), i.e.

[TABLE]

Figure 2 illustrates the optimal stability boundary in the M–$k$ space, given by the curve $\mbox{{M}}_{\!e}(k)$, for fixed values of the Biot numbers B and L. As in the linear stability analysis of Pearson (1958), increasing the Biot number L of the upper surface raises the critical Marangoni number. This fact is obvious from (20), and it corresponds to the physical observation that improving the conductivity of the upper boundary reduces the surface temperature gradients and, consequently, the surface tension driving the flow (Davis, 1987). Also analogous to the linear stability problem is the fact that in the case of two insulating boundaries ($\mbox{{B}}=\mbox{{L}}=0$) the minimum critical Marangoni number $\mbox{{M}}_{\!e}=48$ is achieved for $k=0$. Interestingly, this coincides with Pearson’s linear stability threshold (Pearson, 1958), i.e., Bénard–Marangoni conduction between insulating boundaries is globally stable until an “infinite wavelength” linear instability occurs. This instability is suppressed by any increase in B, and the qualitative distribution of the energy stability boundaries for an imperfectly conducting bottom boundary (B finite) is the same as in the perfectly conducting case ($\mbox{{B}}=\infty$).

Table 1 presents the minimum critical Marangoni number over all wavenumbers $k$, denoted $\mbox{{M}}_{\rm cr}$, and the critical wavenumber $k_{\rm cr}$ for selected Biot numbers L in the extremal cases $\mbox{{B}}=0$ (insulating bottom boundary) and $\mbox{{B}}=\infty$ (perfectly conducting bottom boundary). These values are compared to the corresponding linear stability results from Pearson (1958) and, when available, to the energy stability results obtained by Davis (1969) for finite-\Pran fluids: since these are actually independent of the Prandtl number, they also apply in the infinite-\Pran case. As one would expect, our values are larger than those computed by Davis, because the infinite-$\Pran$ variational problem exploits the explicit coupling between the velocity and temperature fields. Moreover, for $\mbox{{B}}=\infty$ and $\mbox{{L}}=0$ we find $\mbox{{M}}_{\rm cr}=66.84$, a 14.5% improvement on the value $58.36$ computed by Hagstrom & Doering (2010). On the other hand, the optimal energy stability boundary is strictly smaller than the linear stability one, with the only exception of the case $\mbox{{B}}=\mbox{{L}}=0$ (note that Pearson’s linear stability analysis is unchanged when $\Pran=\infty$). This means that, unlike in Rayleigh–Bénard convection, there generally exists a finite range of Marangoni numbers for which the flow is linearly stable, but subcritical instabilities due to perturbations of finite amplitude may occur.

5 Conclusion

To summarise, we have studied the global stability of the purely conductive state of infinite-Prandtl-number Bénard–Marangoni convection using the method of energy, and we have computed the exact critical Marangoni number in wavenumber space for thermal boundary conditions corresponding to perfectly conducting, imperfectly conducting, and perfectly insulating boundaries. We have shown that in the infinite-$\Pran$ limit, the explicit slaving of the velocity field to the temperature field can be exploited to raise the energy stability boundary compared to the finite \Pran case, although a gap with the linear stability threshold remains in all but the insulating case $\mbox{{B}}=\mbox{{L}}=0$. Whether global stability attains up to the linear stability boundary or subcritical instabilities exist, should be determined by bifurcation analysis, numerical simulations, or alternative techniques for global stability analysis, such as those of Goulart & Chernyshenko (2012) and Chernyshenko et al. (2013).

Finally, we note that the analysis presented in this work may be of use in the computation of upper bounds on the convective heat transport using the background field method (see e.g. Constantin & Doering, 1995a, b; Doering & Constantin, 1992, 1994, 1996). The method, already applied to Bénard–Marangoni convection by Hagstrom & Doering (2010), relies on the construction of a background temperature field $\tau(z)$, subject to a nonlinear stability condition obtained by replacing $\mbox{{M}}\,f_{k}(z)$ with $2\,\mbox{{M}}\,\tau^{\prime}(z)\,f_{k}(z)$ in the energy stability constraint (13). Given a candidate background field $\tau$, this nonlinear condition can be analysed using the same ideas presented in §4, and the corresponding the Euler–Lagrange equation has an analytic solution. Whether this allows one to lower Hagstrom & Doering’s original bound, in the same way that their energy stability result was improved in this work, remains an intriguing open question for future work.

Appendix A On the issue of “natural conditions”

Upon restricting attention to real-valued perturbations, (9) implies that the conduction solution is stable if the functional $\mathcal{F}_{k}$ in (8) satisfies

[TABLE]

the infimum being taken over all square-integrable test functions $v$ with a square-integrable (weak) first derivative, and that satisfy the boundary conditions in (10). Instead of the normalisation condition $v(1)=1$ assumed in §§2-4, suppose we normalise $v$ such that

[TABLE]

This choice of normalisation is legitimate because the problem is homogeneous, and because $\mathcal{F}_{k}\{0\}=0$. Letting $\lambda$ be the Lagrange multiplier enforcing (23), the Euler–Lagrange equations for a candidate minimiser of $\mathcal{F}_{k}$ are found by setting to zero the first variations of the augmented functional

[TABLE]

Setting to zero the first variation of $\mathcal{L}_{k}$ with respect to $\lambda$ simply yields (23). Moreover, the necessary condition for $\mathcal{L}_{k}$ to be stationary with respect to $v$ is that

[TABLE]

for any function $h(z)$ that satisfies $h^{\prime}(0)-\mbox{{B}}\,h(0)=0$ and $h^{\prime}(1)+\mbox{{L}}\,h(1)=0$. Upon integrating by parts using the boundary conditions on $v$, we find

[TABLE]

After requiring the right-hand side above to vanish for all perturbations $h$, we conclude that the minimiser of $\mathcal{F}_{k}$ subject to (23) must satisfy the Euler–Lagrange equation

[TABLE]

as well as the “natural conditions”

[TABLE]

the normalisation condition (23), and the original boundary conditions

[TABLE]

Note that no solution exists in general: there are three equations for the two unknowns $v$ and $\lambda$, and furthermore there are three boundary conditions—the two original ones plus the “natural” boundary condition $\mbox{{B}}\,v(0)=0$—but only two integration constants for $v$.

Appendix B Expressions for $g_{k}(z)$ and $h_{k}(z)$

The expressions for $h_{k}(z)$ and $g_{k}(z)$ are

[TABLE]

where

[TABLE]

and

[TABLE]

Appendix C Convergence of $\mathcal{Q}_{k}\{v_{n}\}$ to $\mathcal{Q}_{k}^{\star}$

The test function $v_{n}$ in (18) belongs to the functional space $\Gamma_{0}$ because it satisfies the boundary conditions prescribed on $\Gamma_{0}$, it is square integrable on the interval $[0,1]$, and so is its (weak) derivative

[TABLE]

The convergence of $\mathcal{Q}_{k}\{v_{n}\}$ to $\mathcal{Q}_{k}^{\star}$ follows from a relatively straightforward application of Lebesgue’s dominated convergence theorem. For example, note that $|v_{n}^{\prime}(z)|^{2}\to|v_{\star}^{\prime}(z)|^{2}$ pointwise in $(0,1)$ as $n\to\infty$ since $\xi_{n}=n/(n+1)\to 1$, and that there exists a constant $C_{0}>0$ such that $|v_{n}^{\prime}|^{2}\leq C_{0}$ for all $z\in(0,1)$ and $n\geq 1$ since (i) $|v_{\star}^{\prime}(z)|\leq C_{1}$ for all $z\in(0,\xi_{n})$ for some constant $C_{1}>0$ because $v_{\star}$ is a smooth function, and (ii) by virtue of Taylor’s theorem there exists $\eta\in[\xi_{n},1]$ such that

[TABLE]

Lebesgue’s dominated convergence theorem then implies that $\int_{0}^{1}|v_{n}^{\prime}|^{2}\,\mathrm{d}z\to\int_{0}^{1}|v_{\star}^{\prime}|^{2}\,\mathrm{d}z$ as $n\to\infty$. Similar arguments can be applied to $|v_{n}|^{2}$ and $f_{k}v_{n}$.