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 Reference Schatz and Neitzel2001), cladding processes (Kumar & Roy Reference Kumar and Roy2009), and drying of thin polymer films (Yiantsios et al. Reference Yiantsios, Serpetsi, Doumenc and Guerrier2015), and has recently received increasing attention as a paradigm for shear-driven turbulent transport processes (Boeck & Thess Reference Boeck and Thess1998, Reference Boeck and Thess2001; Hagstrom & Doering Reference Hagstrom and Doering2010).
The first mathematical model of surface-tension-driven convection was proposed by Pearson (Reference Pearson1958), 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
$M_{l}$
, independently of the fluid’s Prandtl number
$Pr$
(the ratio of the fluid’s kinematic viscosity to its thermal diffusivity). Subsequently, Davis (Reference Davis1969) used the energy method to prove that conduction is asymptotically stable against disturbances of arbitrary amplitude when
$M$
is smaller than a critical value
$M_{e}$
, also independently of
$Pr$
. In contrast to Rayleigh–Bénard convection, the linear and energy thresholds
$M_{l}$
and
$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
$Pr$
, they can be improved in the infinite Prandtl number case. This limit is an attractive model for high-Prandtl-number fluids, such as the silicone oils used in experiments (de Bruyn et al.
Reference de Bruyn, Bodenschatz, Morris, Trainoff, Hu and Cannell1996) or Earth’s mantle (Jones Reference Jones1977), because it gives accurate quantitative predictions whilst simplifying the governing equations (Boeck & Thess Reference Boeck and Thess2001). The key observation is that in the limit of infinite
$Pr$
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 Reference Hagstrom and Doering2010), allowing the formulation of an improved variational principle for energy stability. This variational principle, first considered by Hagstrom & Doering (Reference Hagstrom and Doering2010), 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
$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-
$Pr$
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
$\unicode[STIX]{x1D70C}$
, kinematic viscosity
$\unicode[STIX]{x1D708}$
, thermal diffusivity
$\unicode[STIX]{x1D705}$
, and thermal conductivity
$\unicode[STIX]{x1D706}$
, 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}/\unicode[STIX]{x1D706})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}/\unicode[STIX]{x1D705}$
, and
$Q_{0}h/\unicode[STIX]{x1D706}$
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 (Reference Hagstrom and Doering2010).
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
$\unicode[STIX]{x1D703}(\boldsymbol{x},t)$
– evolve according to
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\unicode[STIX]{x1D735}p, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}+w. & \displaystyle\end{eqnarray}$$
We assume that all variables are periodic in the horizontal (
$x$
) direction with period
$\unicode[STIX]{x1D6EC}$
, 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
$\unicode[STIX]{x1D6FE}$
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
$$\begin{eqnarray}\left[\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}+M\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}\right]_{z=1}=0,\end{eqnarray}$$
where the Marangoni number
$M=\unicode[STIX]{x1D6FE}Q_{0}h^{2}/(\unicode[STIX]{x1D706}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$
is the main governing parameter of the system. Finally, letting
$q_{bot}$
and
$q_{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
$$\begin{eqnarray}\displaystyle \left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}-B\unicode[STIX]{x1D703}\right]_{z=0}=0,\quad \left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}z}+L\unicode[STIX]{x1D703}\right]_{z=1}=0, & & \displaystyle\end{eqnarray}$$
where the Biot numbers
$B=q_{bot}h/\unicode[STIX]{x1D706}$
,
$L=q_{top}h/\unicode[STIX]{x1D706}$
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
$B,L\geqslant 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
$B=L=\infty$
. For a comprehensive discussion of the thermal boundary conditions (2.3) we refer the reader to the original work by Pearson (Reference Pearson1958).
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 \unicode[STIX]{x1D703}^{2}\rangle /2$
of a temperature perturbation, where
$\langle \cdot \rangle$
denotes the usual volume average, is found by averaging
$\unicode[STIX]{x1D703}\times$
(2.1c
) and integrating by parts using (2.1b
) and the boundary conditions (2.3) to arrive at
$$\begin{eqnarray}\frac{1}{2}\frac{\text{d}}{\text{d}t}\langle \unicode[STIX]{x1D703}^{2}\rangle =-\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}-w\unicode[STIX]{x1D703}\rangle -L\,\overline{\unicode[STIX]{x1D703}^{2}}(1)-B\overline{\unicode[STIX]{x1D703}^{2}}(0).\end{eqnarray}$$
In this equation and throughout the rest of this section, overlines denote horizontal averages. Clearly, the kinetic energy of the perturbation
$\unicode[STIX]{x1D703}$
does not increase in time if the right-hand side of (3.1) is non-positive at each instant in time, i.e.
$$\begin{eqnarray}\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}-w\unicode[STIX]{x1D703}\rangle +L\,\overline{\unicode[STIX]{x1D703}^{2}}(1)+B\overline{\unicode[STIX]{x1D703}^{2}}(0)\geqslant 0.\end{eqnarray}$$
Upon substituting the horizontal Fourier series expansion of
$\unicode[STIX]{x1D703}$
and
$w$
into (3.2) and into the boundary conditions (2.3) (we consider the case of a finite periodic domain for definiteness; 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 (Reference Hagstrom and Doering2010) that the Fourier amplitudes of the velocity perturbation
$w$
are ‘slaved’ to those of
$\unicode[STIX]{x1D703}$
according to
${\hat{w}}_{k}(z)=-M\,f_{k}(z)\hat{\unicode[STIX]{x1D703}}_{k}(1)$
, where
$$\begin{eqnarray}f_{k}(z)=\frac{k\sinh k}{\sinh (2k)-2k}[kz\cosh (kz)-\sinh (kz)+(1-k\coth k)\,z\sinh (kz)],\end{eqnarray}$$
we can rewrite
$$\begin{eqnarray}\langle \left|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\right|^{2}-w\unicode[STIX]{x1D703}\rangle +L\,\overline{\unicode[STIX]{x1D703}^{2}}(1)+B\overline{\unicode[STIX]{x1D703}^{2}}(0)=2\mathop{\sum }_{k\geqslant 0}{\mathcal{F}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\},\end{eqnarray}$$
where the sum is over all positive wavenumbers and
$$\begin{eqnarray}\displaystyle {\mathcal{F}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}:=\int _{0}^{1}\{|\hat{\unicode[STIX]{x1D703}}_{k}^{\prime }(z)|^{2}+k^{2}|\hat{\unicode[STIX]{x1D703}}_{k}(z)|^{2}+Mf_{k}(z)\text{Re}[{\hat{\unicode[STIX]{x1D703}}_{k}}^{\ast }(z)\hat{\unicode[STIX]{x1D703}}_{k}(1)]\}\,\text{d}z+L|\hat{\unicode[STIX]{x1D703}}_{k}(1)|^{2}+B|\hat{\unicode[STIX]{x1D703}}_{k}(0)|^{2}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
(Here and in the following,
$^{\ast }$
denotes complex conjugation and primes denote total differentiation with respect to
$z$
.)
Since among all possible perturbations are those defined by a single wavenumber we conclude that a necessary and sufficient condition for the global stability of Bénard–Marangoni conduction is that, for all wavenumbers
$k\geqslant 0$
,
$$\begin{eqnarray}{\mathcal{F}}_{k}\{\hat{\unicode[STIX]{x1D703}}_{k}\}\geqslant 0\end{eqnarray}$$
for all complex-valued perturbation Fourier amplitudes
$\hat{\unicode[STIX]{x1D703}}_{k}$
(hereafter simply referred to as perturbations) that satisfy
$$\begin{eqnarray}\displaystyle {\hat{\unicode[STIX]{x1D703}}_{k}}^{\prime }(0)-B\hat{\unicode[STIX]{x1D703}}_{k}(0)=0,\quad {\hat{\unicode[STIX]{x1D703}}_{k}}^{\prime }(1)+L\hat{\unicode[STIX]{x1D703}}_{k}(1)=0. & & \displaystyle\end{eqnarray}$$
In fact, we may restrict our attention to real-valued
$\hat{\unicode[STIX]{x1D703}}_{k}$
because the real and imaginary parts give identical and independent contributions to the left-hand side of (3.6). Moreover, note that (3.6) holds trivially when
$\hat{\unicode[STIX]{x1D703}}_{k}(1)=0$
, and that its left-hand side is homogeneous quadratic in
$\hat{\unicode[STIX]{x1D703}}_{k}$
. Since the boundary conditions in (3.7) are also homogeneous, we can further restrict our attention to perturbations that satisfy the normalisation condition
$\hat{\unicode[STIX]{x1D703}}_{k}(1)=1$
and, without any loss of generality, we may replace (2.3) with the boundary conditions
$$\begin{eqnarray}\displaystyle {\hat{\unicode[STIX]{x1D703}}_{k}}^{\prime }(0)-B\hat{\unicode[STIX]{x1D703}}_{k}(0)=0,\quad {\hat{\unicode[STIX]{x1D703}}_{k}}^{\prime }(1)=-L,\quad \hat{\unicode[STIX]{x1D703}}_{k}(1)=1. & & \displaystyle\end{eqnarray}$$
Putting these observations together, we define the space of admissible perturbations as
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{0}:=\left\{\!v(z):\int _{0}^{1}(|v^{\prime }|^{2}+|v|^{2})\,\text{d}z<\infty ,v^{\prime }(0)=Bv(0),v^{\prime }(1)=-L,v(1)=1\!\right\}.\end{eqnarray}$$
Finally, it is clear that (3.6) 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$
,
$$\begin{eqnarray}{\mathcal{Q}}_{k}^{\star }:=\inf _{v\in \unicode[STIX]{x1D6E4}_{0}}\int _{0}^{1}[|v^{\prime }(z)|^{2}+k^{2}|v(z)|^{2}+Mf_{k}(z)v(z)]\,\text{d}z+L+B|v(0)|^{2}\geqslant 0.\end{eqnarray}$$
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 }\geqslant 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
$$\begin{eqnarray}{\mathcal{Q}}_{k}\{v\}=\int _{0}^{1}[|v^{\prime }(z)|^{2}+k^{2}|v(z)|^{2}+Mf_{k}(z)v(z)]\,\text{d}z+L+B|v(0)|^{2}\end{eqnarray}$$
is not generally attained by any test function
$v\in \unicode[STIX]{x1D6E4}_{0}$
. In fact, a straightforward application of the calculus of variations (see e.g. Courant & Hilbert Reference Courant and Hilbert1953; Giaquinta & Hildebrandt Reference Giaquinta and Hildebrandt1996) shows that any candidate minimiser
$v_{\star }$
must satisfy the second-order Euler–Lagrange differential equation
$$\begin{eqnarray}v_{\star }^{\prime \prime }(z)-k^{2}v_{\star }(z)={\textstyle \frac{1}{2}}Mf_{k}(z),\end{eqnarray}$$
subject to the three boundary conditions
$v_{\star }^{\prime }(0)=Bv_{\star }(0)$
,
$v_{\star }^{\prime }(1)=-L$
, and
$v_{\star }(1)=1$
. This problem is overconstrained, 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
$\unicode[STIX]{x1D6E4}_{0}$
is replaced with the minimisation of
${\mathcal{F}}_{k}$
in (3.5) over all normalised test functions that satisfy the boundary conditions (3.7). As we demonstrate in appendix A for the commonly used normalisation
$\int _{0}^{1}|v(z)|^{2}\,\text{d}z=1$
, the Euler–Lagrange equations for the minimiser of
${\mathcal{F}}_{k}$
are overconstrained by so-called ‘natural conditions’ that arise when setting to zero its first variation (Courant & Hilbert Reference Courant and Hilbert1953, Chapter IV, § 5.1).
This obstacle is overcome if we can enlarge the space of test functions in such a way that (4.2) 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)=-L$
, and minimise
${\mathcal{Q}}_{k}\{v\}$
over the larger space of functions
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{1}:=\left\{\!v(z):\int _{0}^{1}(|v^{\prime }|^{2}+|v|^{2})\,\text{d}z<\infty ,v^{\prime }(0)=Bv(0),v(1)=1\!\right\}.\end{eqnarray}$$
Having removed one boundary condition, in fact, the Euler–Lagrange equation (4.2) 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
$$\begin{eqnarray}v_{\star }(z)=Mg_{k}(z)+h_{k}(z),\end{eqnarray}$$
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 }\}\leqslant {\mathcal{Q}}_{k}^{\star }$
, because
$\unicode[STIX]{x1D6E4}_{0}$
is a proper subset of
$\unicode[STIX]{x1D6E4}_{1}$
and
$v_{\star }$
minimises
${\mathcal{Q}}_{k}\{v\}$
over
$\unicode[STIX]{x1D6E4}_{1}$
. On the other hand,
${\mathcal{Q}}_{k}\{v_{\star }\}\geqslant Q_{k}^{\star }$
because we can find a sequence of functions
$(v_{n})_{n\geqslant 1}$
with
$v_{n}\in \unicode[STIX]{x1D6E4}_{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
$\unicode[STIX]{x1D709}_{n}=n/(n+1)$
and take
$$\begin{eqnarray}v_{n}(z):=\left\{\begin{array}{@{}ll@{}}v_{\star }(z)\quad & \text{if }0\leqslant z\leqslant \unicode[STIX]{x1D709}_{n},\\ v_{\star }(\unicode[STIX]{x1D709}_{n})+{\displaystyle \frac{2+L(1-\unicode[STIX]{x1D709}_{n})-2v_{\star }(\unicode[STIX]{x1D709}_{n})}{1-\unicode[STIX]{x1D709}_{n}}}(z-\unicode[STIX]{x1D709}_{n})\quad & \text{if }\unicode[STIX]{x1D709}_{n}\leqslant z\leqslant {\displaystyle \frac{1+\unicode[STIX]{x1D709}_{n}}{2}},\\ 1+L\,(1-z)\quad & \text{if }{\displaystyle \frac{1+\unicode[STIX]{x1D709}_{n}}{2}}\leqslant z\leqslant 1.\end{array}\right.\end{eqnarray}$$
Note that these test functions are simply piecewise-linear continuous functions on
$[\unicode[STIX]{x1D709}_{n},1]$
satisfying
$v_{n}(\unicode[STIX]{x1D709}_{n})=v_{\star }(\unicode[STIX]{x1D709}_{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
$\unicode[STIX]{x1D6E4}_{0}$
.
Having computed the minimiser
$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$
(4.2) by parts using the boundary conditions on
$\unicode[STIX]{x1D6E4}_{1}$
to show that
$$\begin{eqnarray}\int _{0}^{1}[|v_{\star }^{\prime }(z)|^{2}+k^{2}|v_{\star }(z)|^{2}]\,\text{d}z=v_{\star }^{\prime }(1)-B|v_{\star }(0)|^{2}-\frac{1}{2}M\int _{0}^{1}f_{k}(z)v_{\star }(z)\,\text{d}z.\end{eqnarray}$$
Upon combining this with (4.1) and (4.4) we find
$$\begin{eqnarray}{\mathcal{Q}}_{k}^{\star }=\underbrace{\left[\frac{1}{2}\int _{0}^{1}g_{k}(z)f_{k}(z)\,\text{d}z\right]}_{=:\unicode[STIX]{x1D6FC}_{k}}M^{2}+\underbrace{\left[{g_{k}}^{\prime }(1)+\frac{1}{2}\int _{0}^{1}h_{k}(z)f_{k}(z)\,\text{d}z\right]}_{=:\unicode[STIX]{x1D6FD}_{k}}M+{h_{k}}^{\prime }(1)+L.\end{eqnarray}$$
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
$\unicode[STIX]{x1D6FC}_{k}$
and
$\unicode[STIX]{x1D6FD}_{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
$\unicode[STIX]{x1D6FC}_{k}$
and
$\unicode[STIX]{x1D6FD}_{k}$
, together with
$h_{k}^{\prime }(1)$
, are plotted in figure 1 for Biot numbers
$B=0$
(corresponding to a perfectly insulating bottom boundary),
$B=1$
,
$B=10$
, and
$B=\infty$
(corresponding to a perfectly conducting bottom boundary). Note that the leading-order coefficient
$\unicode[STIX]{x1D6FC}_{k}$
is negative for all
$k$
, and it must be so because perturbations at any wavenumber
$k$
eventually become linearly unstable (Pearson Reference Pearson1958), implying that
${\mathcal{Q}}_{k}^{\star }<0$
for all sufficiently large
$M$
. Consequently, the largest Marangoni number
$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 (4.7), i.e.
$$\begin{eqnarray}M_{e}=\frac{-\unicode[STIX]{x1D6FD}_{k}-\sqrt{\unicode[STIX]{x1D6FD}_{k}^{2}-4\,\unicode[STIX]{x1D6FC}_{k}[h_{k}^{\prime }(1)+L]}}{2\unicode[STIX]{x1D6FC}_{k}}.\end{eqnarray}$$
Figure 1. Value of the coefficients
$\unicode[STIX]{x1D6FC}_{k}$
and
$\unicode[STIX]{x1D6FD}_{k}$
, defined as in (4.7), and of
$h_{k}^{\prime }(1)$
, plotted as a function of the wavenumber
$k$
for
$B=0$
(dotted line),
$B=1$
(dot-dashed line),
$B=10$
(dashed line) and
$B=\infty$
(solid line).
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$
(insulating bottom boundary) and
$B=\infty$
(perfectly conducting bottom boundary). Where available, the corresponding values for linear stability (Pearson Reference Pearson1958) and energy stability at finite
$Pr$
(Davis Reference Davis1969) are also reported.
Figure 2 illustrates the optimal stability boundary in the
$M$
–
$k$
space, given by the curve
$M_{e}(k)$
, for fixed values of the Biot numbers
$B$
and
$L$
. As in the linear stability analysis of Pearson (Reference Pearson1958), increasing the Biot number
$L$
of the upper surface raises the critical Marangoni number. This fact is obvious from (4.7), 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 Reference Davis1987). Also analogous to the linear stability problem is the fact that in the case of two insulating boundaries (
$B=L=0$
) the minimum critical Marangoni number
$M_{e}=48$
is achieved for
$k=0$
. Interestingly, this coincides with Pearson’s linear stability threshold (Pearson Reference Pearson1958), 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 (
$B=\infty$
).
Figure 2. Critical energy stability curves for Bénard–Marangoni conduction at infinite Prandtl number in the
$M$
–
$k$
space, computed using (4.7) for four different values of the Biot number of the bottom surface,
$B$
, and six values of the Biot number of the top surface,
$L$
. The extremal cases
$B=0$
and
$B=\infty$
correspond to a perfectly insulating and a perfectly conducting bottom boundary, respectively.
Table 1 presents the minimum critical Marangoni number over all wavenumbers
$k$
, denoted
$M_{cr}$
, and the critical wavenumber
$k_{cr}$
for selected Biot numbers
$L$
in the extremal cases
$B=0$
(insulating bottom boundary) and
$B=\infty$
(perfectly conducting bottom boundary). These values are compared to the corresponding linear stability results from Pearson (Reference Pearson1958) and, when available, to the energy stability results obtained by Davis (Reference Davis1969) for finite-
$Pr$
fluids: since these are actually independent of the Prandtl number, they also apply in the infinite-
$Pr$
case. As one would expect, our values are larger than those computed by Davis, because the infinite-
$Pr$
variational problem exploits the explicit coupling between the velocity and temperature fields. Moreover, for
$B=\infty$
and
$L=0$
we find
$M_{cr}=66.84$
, a 14.5 % improvement on the value 58.36 computed by Hagstrom & Doering (Reference Hagstrom and Doering2010). On the other hand, the optimal energy stability boundary is strictly smaller than the linear stability one, with the only exception of the case
$B=L=0$
(note that Pearson’s linear stability analysis is unchanged when
$Pr=\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 top and bottom boundaries. We have shown that in the infinite-
$Pr$
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
$Pr$
case, although a gap with the linear stability threshold remains in all but the insulating case
$B=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 (Reference Goulart and Chernyshenko2012) and Chernyshenko et al. (Reference Chernyshenko, Huang, Goulart, Lasagna and Tutty2013). Moreover, we remark that our results are only valid for an ‘ideal’ layer, which is infinite or periodic in the horizontal directions. How the presence of sidewalls affects stability remains an open question, the analysis being complicated by the fact that (generalised) Fourier modes cannot be decoupled to obtain independent, one-dimensional variational problems.
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. Doering & Constantin Reference Doering and Constantin1992, Reference Doering and Constantin1994, Reference Doering and Constantin1996; Constantin & Doering Reference Constantin and Doering1995a
,Reference Constantin and Doering
b
). The method, already applied to Bénard–Marangoni convection by Hagstrom & Doering (Reference Hagstrom and Doering2010), relies on the construction of a background temperature field
$\unicode[STIX]{x1D70F}(z)$
, subject to a nonlinear stability condition obtained by replacing
$Mf_{k}(z)$
with
$2M\unicode[STIX]{x1D70F}^{\prime }(z)f_{k}(z)$
in the energy stability constraint (3.10). Given a candidate background field
$\unicode[STIX]{x1D70F}$
, 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, equation (3.6) implies that the conduction solution is stable if the functional
${\mathcal{F}}_{k}$
in (3.5) satisfies
$$\begin{eqnarray}\inf _{v}{\mathcal{F}}_{k}\{v\}\geqslant 0,\end{eqnarray}$$
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 (3.7). Instead of the normalisation condition
$v(1)=1$
assumed in §§ 2–4, suppose we normalise
$v$
such that
$$\begin{eqnarray}\int _{0}^{1}|v(z)|^{2}\,\text{d}z=1.\end{eqnarray}$$
This choice of normalisation is legitimate because the problem is homogeneous, and because
${\mathcal{F}}_{k}\{0\}=0$
. Letting
$\unicode[STIX]{x1D706}$
be the Lagrange multiplier enforcing (A 2), 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
$$\begin{eqnarray}{\mathcal{L}}_{k}\{v,\unicode[STIX]{x1D706}\}:={\mathcal{F}}_{k}\{v\}+\unicode[STIX]{x1D706}\left(\int _{0}^{1}|v(z)|^{2}\,\text{d}z-1\right).\end{eqnarray}$$
Setting to zero the first variation of
${\mathcal{L}}_{k}$
with respect to
$\unicode[STIX]{x1D706}$
simply yields (A 2). Moreover, the necessary condition for
${\mathcal{L}}_{k}$
to be stationary with respect to
$v$
is that
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}{\mathcal{L}}_{k}}{\unicode[STIX]{x1D6FF}v}:=\lim _{\unicode[STIX]{x1D700}\rightarrow 0}\frac{{\mathcal{L}}_{k}\{v+\unicode[STIX]{x1D700}\,h,\unicode[STIX]{x1D706}\}-{\mathcal{L}}_{k}\{v,\unicode[STIX]{x1D706}\}}{\unicode[STIX]{x1D700}}=0\end{eqnarray}$$
for any function
$h(z)$
that satisfies
$h^{\prime }(0)-Bh(0)=0$
and
$h^{\prime }(1)+L\,h(1)=0$
. Upon integrating by parts using the boundary conditions on
$v$
, we find
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}{\mathcal{L}}_{k}}{\unicode[STIX]{x1D6FF}v} & = & \displaystyle \int _{0}^{1}[-2v^{\prime \prime }(z)+2(k^{2}+\unicode[STIX]{x1D706})v(z)+Mf_{k}(z)v(1)]h(z)\,\text{d}z\nonumber\\ \displaystyle & & \displaystyle +\,\left[Lv(1)+M\int _{0}^{1}f_{k}(z)v(z)\,\text{d}z\right]h(1)+Bv(0)h(0).\end{eqnarray}$$
After requiring the right-hand side above to vanish for all perturbations
$h$
, we conclude that the minimiser of
${\mathcal{F}}_{k}$
subject to (A 2) must satisfy the Euler–Lagrange equation
$$\begin{eqnarray}-2v^{\prime \prime }(z)+2(k^{2}+\unicode[STIX]{x1D706})v(z)+Mf_{k}(z)v(1)=0,\end{eqnarray}$$
as well as the ‘natural conditions’
$$\begin{eqnarray}\displaystyle Lv(1)+M\int _{0}^{1}f_{k}(z)v(z)\,\text{d}z=0,\quad Bv(0)=0, & & \displaystyle\end{eqnarray}$$
the normalisation condition (A 2), and the original boundary conditions
$$\begin{eqnarray}\displaystyle v^{\prime }(0)-Bv(0)=0,\quad v^{\prime }(1)+Lv(1)=0. & & \displaystyle\end{eqnarray}$$
Note that no solution exists in general: there are three equations for the two unknowns
$v$
and
$\unicode[STIX]{x1D706}$
, and furthermore there are three boundary conditions – the two original ones plus the ‘natural’ boundary condition
$Bv(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
$$\begin{eqnarray}\displaystyle & \displaystyle h_{k}(z)=\frac{B\sinh (kz)+k\cosh (kz)}{B\sinh k+k\cosh k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g_{k}(z)=\frac{P_{1}(z)\cosh (kz)+P_{2}(z)\sinh (kz)}{8k[\sinh (2k)-2k](B\sinh k+k\cosh k)}, & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle P_{1}(z) & := & \displaystyle -k\sinh k[(k\cosh k-\sinh k)z+3\sinh k]zB\nonumber\\ \displaystyle & & \displaystyle -\,k^{2}\cosh k(k\cosh k-\sinh k)z^{2}-3k^{2}z\cosh k\sinh k\nonumber\\ \displaystyle & & \displaystyle +\,k(k\cosh k\sinh k+k^{2}-2|\!\cosh k|^{2}+2),\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle P_{2}(z) & := & \displaystyle [\!k^{2}|\!\sinh k|^{2}z^{2}+\sinh k(k\cosh k-\sinh k)z\nonumber\\ \displaystyle & & \displaystyle +\,k\cosh k\sinh k+k^{2}+|\!\cosh k|^{2}-1\!]B+k^{3}z^{2}\cosh k\sinh k\nonumber\\ \displaystyle & & \displaystyle +\,k\cosh k(k\cosh k-\sinh k)z+3k\cosh k\sinh k.\end{eqnarray}$$
Appendix C. Convergence of
${\mathcal{Q}}_{k}\{v_{n}\}$
to
${\mathcal{Q}}_{k}^{\star }$
The test function
$v_{n}$
in (4.5) belongs to the functional space
$\unicode[STIX]{x1D6E4}_{0}$
because it satisfies the boundary conditions prescribed on
$\unicode[STIX]{x1D6E4}_{0}$
, it is square integrable on the interval
$[0,1]$
, and so is its (weak) derivative
$$\begin{eqnarray}v_{n}^{\prime }(z):=\left\{\begin{array}{@{}ll@{}}v_{\star }^{\prime }(z)\quad & \text{if }0\leqslant z\leqslant \unicode[STIX]{x1D709}_{n},\\ {\displaystyle \frac{2+L(1-\unicode[STIX]{x1D709}_{n})-2v_{\star }(\unicode[STIX]{x1D709}_{n})}{1-\unicode[STIX]{x1D709}_{n}}}\quad & \text{if }\unicode[STIX]{x1D709}_{n}\leqslant z\leqslant {\displaystyle \frac{1+\unicode[STIX]{x1D709}_{n}}{2}},\\ -L\quad & \text{if }{\displaystyle \frac{1+\unicode[STIX]{x1D709}_{n}}{2}}\leqslant z\leqslant 1.\end{array}\right.\end{eqnarray}$$
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}\rightarrow |v_{\star }^{\prime }(z)|^{2}$
pointwise in
$(0,1)$
as
$n\rightarrow \infty$
since
$\unicode[STIX]{x1D709}_{n}=n/(n+1)\rightarrow 1$
, and that there exists a constant
$C_{0}>0$
such that
$|v_{n}^{\prime }|^{2}\leqslant C_{0}$
for all
$z\in (0,1)$
and
$n\geqslant 1$
since (i)
$|v_{\star }^{\prime }(z)|\leqslant C_{1}$
for all
$z\in (0,\unicode[STIX]{x1D709}_{n})$
for some constant
$C_{1}>0$
because
$v_{\star }$
is a smooth function, and (ii) by virtue of Taylor’s theorem there exists
$\unicode[STIX]{x1D702}\in [\unicode[STIX]{x1D709}_{n},1]$
such that
$$\begin{eqnarray}\left|\frac{2+L(1-\unicode[STIX]{x1D709}_{n})-2v_{\star }(\unicode[STIX]{x1D709}_{n})}{1-\unicode[STIX]{x1D709}_{n}}\right|=\left|\frac{2+L(1-\unicode[STIX]{x1D709}_{n})-2[1-v_{\star }^{\prime }(\unicode[STIX]{x1D702})(1-\unicode[STIX]{x1D709}_{n})]}{1-\unicode[STIX]{x1D709}_{n}}\right|\leqslant L+2C_{1}.\end{eqnarray}$$
Lebesgue’s dominated convergence theorem then implies that
$\int _{0}^{1}|v_{n}^{\prime }|^{2}\,\text{d}z\rightarrow \int _{0}^{1}|v_{\star }^{\prime }|^{2}\,\text{d}z$
as
$n\rightarrow \infty$
. Similar arguments can be applied to
$|v_{n}|^{2}$
and
$f_{k}v_{n}$
.
