The effect of phase change on stability of convective flow in a layer of volatile liquid driven by a horizontal temperature gradient

TL;DR

This paper investigates how phase change influences the stability of buoyancy-thermocapillary convection in a volatile liquid layer under a horizontal temperature gradient, highlighting the role of gas phase composition and heat transport.

Contribution

It provides a linear stability analysis showing the impact of gas phase transport and phase change on convection onset, extending understanding at atmospheric conditions.

Findings

Gas phase transport significantly affects convection thresholds.

Decreasing air concentration raises the instability threshold.

Phase change effects are notable even at atmospheric conditions.

Abstract

Buoyancy-thermocapillary convection in a layer of volatile liquid driven by a horizontal temperature gradient arises in a variety of situations. Recent studies have shown that the composition of the gas phase, which is typically a mixture of vapour and air, has a noticeable effect on the critical Marangoni number describing the onset of convection as well as on the observed convection pattern. Specifically, as the total pressure or, equivalently, the average concentration of air is decreased, the threshold of the instability leading to the emergence of convective rolls is found to increase rather significantly. We present a linear stability analysis of the problem which shows that this trend can be readily understood by considering the transport of heat and vapour through the gas phase. In particular, we show that transport in the gas phase has a noticeable effect even at atmospheric conditions, when phase change is greatly suppressed.

Figures

Equations

$$\rho \left({\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho \left(T,c_a\right)\mathbf{g}$$

$$\nabla \cdot \mathbf{u}=0.$$

$${\rho }_l={\rho }_l^0[1-{\beta }_l\left(T-T_0\right)],$$

$${\rho }_g={\rho }_a+{\rho }_v,$$

$${\rho }_{a,v}=\frac{p_{a,v}}{{\bar{R}}_{a,v}T},$$

$$p_g=p_a+p_v.$$

$${\partial }_t{n}_{a,v}+\nabla \cdot {\mathbf{j}}_{a,v}=0,$$

$$n_{a,v}=\frac{{\rho }_{a,v}}{m_{a,v}},$$

$${\mathbf{j}}_{a,v}=n_{a,v}\mathbf{u}-n_{g}D\nabla c_{a,v}=n_g(\mathbf{u}{c}_{a,v}-D\nabla c_{a,v}),$$

$$c_{a,v}=\frac{n_{a,v}}{n_g}=\frac{p_{a,v}}{p_g}$$

$${\partial }_t(n_g{c}_{a,v})+n_g\mathbf{u}\cdot \nabla c_{a,v}=\nabla \cdot (n_{g}D\nabla c_{a,v}).$$

$$n_g=\frac{p_g}{k_{B}T},$$

$${\partial }_t{c}_{a,v}+\mathbf{u}\cdot \nabla c_{a,v}=\nabla \cdot (D\nabla c_{a,v}).$$

$${\partial }_{t}T+\mathbf{u}\cdot \nabla T=\alpha {\nabla }^{2}T,$$

$$J=\frac{2\chi }{2-\chi }{\rho }_v\sqrt{\frac{{\bar{R}}_v{T}_i}{2\pi }}\left[\frac{p_l-p_g}{{\rho }_l{\bar{R}}_v{T}_i}+\frac{\mathcal{L}}{{\bar{R}}_v{T}_i}\frac{T_i-T_s}{T_s}\right],$$

$$\ln \frac{p_v}{p_v^0}=\frac{\mathcal{L}}{{\bar{R}}_v}\left[\frac{1}{T_0}-\frac{1}{T_s}\right],$$

$$\frac{J}{m_v}=\hat{\mathbf{n}}\cdot {\mathbf{j}}_v=n_g\hat{\mathbf{n}}\cdot ([{\mathbf{u}}_g-{\mathbf{u}}_i]c_v-D\nabla c_v),$$

$$0=\hat{\mathbf{n}}\cdot {\mathbf{j}}_a=n_g\hat{\mathbf{n}}\cdot ([{\mathbf{u}}_g-{\mathbf{u}}_i]c_a-D\nabla c_a).$$

$$\hat{\mathbf{n}}\cdot \nabla c_v=-\frac{c_{a}J}{m_v{n}_{g}D}$$

$$\hat{\mathbf{n}}\cdot ({\mathbf{u}}_g-{\mathbf{u}}_i)=\frac{J}{m_v{n}_g}.$$

$$\mathcal{L}J=\hat{\mathbf{n}}\cdot k_g\nabla T_g-\hat{\mathbf{n}}\cdot k_l\nabla T_l,$$

$$\hat{\mathbf{n}}\cdot \left({\mathbf{u}}_l-{\mathbf{u}}_i\right)=\frac{J}{{\rho }_l}.$$

$$T_l=T_i=T_v$$

$$(\mathbb{I}-\hat{\mathbf{n}}\hat{\mathbf{n}})\cdot ({\mathbf{u}}_{\mathbf{l}}-{\mathbf{u}}_{\mathbf{g}})=0.$$

$$({\Sigma }_l-{\Sigma }_g)\cdot \hat{\mathbf{n}}=\hat{\mathbf{n}}\kappa \sigma -\gamma {\nabla }_s{T}_i$$

$$Ma=\frac{\gamma d_l^2}{{\mu }_l{\alpha }_l}\tau$$

$$Ra=\frac{g{\beta }_l{d}_l^4\tau }{{\nu }_l{\alpha }_l},$$

$$B{o}_D=\frac{Ra}{Ma}=\frac{{\rho }_{l}g{\beta }_l{d}_l^2}{\gamma }$$

$$\tilde{\nabla }\cdot {\tilde{\mathbf{u}}}_l$$

$${\partial }_{\tilde{t}}{\tilde{\mathbf{u}}}_l+{\tilde{\mathbf{u}}}_l\cdot \tilde{\nabla }{\tilde{\mathbf{u}}}_l$$

$${\partial }_{\tilde{t}}{\tilde{T}}_l+{\tilde{\mathbf{u}}}_l\cdot \tilde{\nabla }{\tilde{T}}_l$$

Full text

1 Introduction

Convection in fluids with a free surface driven by a horizontal temperature gradient has been studied extensively due to applications in crystal growth and thermal management. The first systematic study of convection in nonvolatile fluids appears to be due to Birikh (1966) who derived an analytic solution for a planar return flow in a uniform laterally unbounded layer due to a combination of buoyancy and thermocapillary stresses. This solution also describes the flow away from the end walls in a laterally bounded geometry: thermocapillary stresses drive the flow from the hot end towards the cold end near the free surface, with a return flow near the bottom of the layer. Kirdyashkin (1984) repeated Birikh’s theoretical analysis and validated the analytical solutions experimentally.

Smith & Davis (1983a, b) performed a linear stability analysis of such flows in the limit of vanishing dynamic Bond number, $Bo_{D}$ (i.e., ignoring buoyancy effects). They predicted that, depending on the Prandtl number of the liquid, the uniform base flow would undergo an instability towards either surface waves (for $Pr<0.15$, which corresponds to liquid metals) or hydrothermal waves (for $Pr>0.15$, which corresponds to gases and nonmetallic liquids) above a critical Marangoni number $Ma$, which characterizes the magnitude of thermocapillary stresses. In particular, hydrothermal waves were predicted to form at an angle to the direction of the thermal gradient and travel in the direction of thermal gradient. As $Pr$ increases, the angle changes smoothly from nearly transverse to nearly parallel to the thermal gradient. The theoretical predictions have since been thoroughly tested and verified both in microgravity and for thin films in terrestrial conditions. A thorough overview of these experiments is presented in a review paper by Schatz & Neitzel (2001).

A different type of instability is found at $Bo_{D}=O(1)$, when buoyancy is nonnegligible. Villers & Platten (1992) studied buoyancy-thermocapillary convection in a rectangular cavity for acetone ($Pr=4.24$) experimentally and numerically. Although acetone is volatile, reasonable agreement was found between the experimental observations at atmospheric conditions and the numerical simulations based on a one-sided model that ignored heat and mass transfer in the gas phase. For low $Ma$ a featureless planar return flow was found, which is well-described by Birikh’s solution. At higher $Ma$ a steady cellular pattern featuring multiple convection rolls was found instead of hydrothermal waves. The (transverse) convection rolls were found to rotate in the same direction, unlike the case of pure buoyancy (or Rayleigh-Bénard) convection driven by a vertical temperature gradient. At even higher $Ma$ the steady state was found to be replaced by an oscillatory pattern that was also unlike a hydrothermal wave: the convection rolls were observed to travel in the direction opposite to that of the thermal gradient. Similar results were obtained later by De Saedeleer et al. (1996) for decane ($Pr=15$) and Garcimartín et al. (1997) for 0.65 cSt and 2.0 cSt silicone oil ($Pr=10$ and 30, respectively) in rectangular cavities with strong confinement in the spanwise direction.

Riley & Neitzel (1998) performed one of the most extensive and detailed experimental studies of convection in 1 cSt silicone oil ($Pr=13.9$) in a rectangular cavity with a spanwise dimension $L_{y}$ comparable to the streamwise dimension $L_{x}$. They discovered that a direct transitions from steady, unicellular flow to hydrothermal waves takes place for small values of the dynamic Bond number ($Bo_{D}\lesssim 0.2$), while for $Bo_{D}\gtrsim 0.2$ the results are similar to those of earlier studies with spanwise confinement: the featureless return flow first transitions to steady co-rotating convection cells and, upon further increase in $Ma$, to an oscillatory multicellular pattern. Riley and Neitzel also determined the critical values of $Ma$ and the wavelength $\lambda$ of the convective pattern as a function of $Bo_{D}$. Burguete et al. (2001) performed experiments on convection in 0.65 cSt silicone oil ($Pr=10.3$) in a rectangular cavity with different aspect ratios where spanwise dimension was greater than the streamwise dimension. Similarly, they found that the base return flow destabilizes into either oblique travelling waves or longitudinal stationary rolls, respectively, for low and high $Bo_{D}$.

Convective patterns have also been studied extensively using numerical simulations. Most of the numerical studies (Ban Hadid & Roux, 1990; Villers & Platten, 1992; Ben Hadid & Roux, 1992; Mundrane & Zebib, 1994; Lu & Zhuang, 1998; Shevtsova et al., 2003) were based on one-sided models which ignore the transport in the gas phase, assume that the free surface is flat and nondeformable, the bottom wall and the interface are adiabatic, and phase change is negligible. These numerical simulations were able to reproduce some features of the experimental studies. For example, Shevtsova et al. (2003) and Shevtsova & Legros (2003) performed numerical simulations for decane ($Pr=14.8$) in a rectangular layer at different $Bo_{D}$. They found that as $Ma$ increases, the primary instability leads to hydrothermal waves for $Bo_{D}\leq 0.25$, while for $Bo_{D}\geq 0.32$ the primary (secondary) instability produces a steady (oscillatory) multicellular flow.

Since it does not account for buoyancy, the linear stability analysis of Smith & Davis (1983a, b) fails to predict the stationary patterns that emerge for $Bo_{D}=O(1)$. However, most of the linear stability analyses accounting for buoyancy also failed to predict the correct pattern, i.e., stationary (transverse) rolls that were observed in the experiments. Using adiabatic boundary conditions at the top and bottom of the liquid layer, Parmentier et al. (1993) predicted transition to travelling waves rather than steady multicellular pattern for a range of $Pr$ from 0.01 to 10, regardless of the value of $Bo_{D}$. Chan & Chen (2010), who used similar assumptions, also predicted transition to travelling waves for a $Pr=13.9$ fluid. Moreover, their predicted critical $Ma$ and wavelength $\lambda$ do not match the experiment (Riley & Neitzel, 1998). In both cases the predicted travelling waves are oblique for smaller $Bo_{D}$ and become transverse for $Bo_{D}$ greater than some critical $O(1)$ value.

Mercier & Normand (1996) showed that transition to a stationary convective pattern can take place if the adiabatic boundary conditions are replaced with Newton’s cooling law, although that requires an unrealistically large surface Biot number ($Bi\gtrsim 185/Bo_{D}$). Moreover, the predicted pattern corresponds to longitudinal convection rolls, while in most experiments transverse rolls were observed. In a subsequent paper Mercier & Normand (2002) considered the effects of the end walls, which they described as spatial disturbances superimposed on the uniform base flow. Their analysis predicted that, depending upon the Prandtl number, convection rolls would develop near the hot end (for $Pr>4$), near the cold end (for $Pr<0.01$), or at both end walls (for $0.01<Pr<4$).

To our knowledge, the study by Priede & Gerbeth (1997) is the only one to date which correctly predicts the formation of a stationary pattern at $Bo_{D}=O(1)$. They argued that travelling waves, being convectively unstable, cannot get sufficiently amplified by linear instability in a laterally bounded system. At the same time, the effect of lateral walls extends far into the bulk of the liquid layer, resulting in a stationary pattern of transverse convection rolls. Their predicted critical values of $Ma$ are in reasonable agreement with the threshold values found by Riley & Neitzel (1998), although their is a systematic deviation, suggesting that some important effects have not been considered.

The volatility of the fluids used in the experiments is one source of the discrepancy with analytical (and most numerical) predictions. Although at atmospheric conditions phase change is usually strongly suppressed, it can still play a role. The latent heat associated with phase change can significantly modify the interfacial temperature, and hence the thermocapillary stresses. However, there are very few studies that investigated this effect. Li et al. (2012) have studied nonadiabatic effects by using Newton’s cooling law with a small Biot number. Their numerical simulations ignored phase change, but were able to reproduce many features of the experimental observations at atmospheric conditions. Ji et al. (2008) considered phase change, but ignored buoyancy, so their analysis is only applicable for thin films or under microgravity, i.e., when $Bo_{D}\approx 0$.

A few recent studies investigated the role played by the gas phase, which is generally a binary mixture of a noncondensable gas (typically air) and the vapour, in more detail. In particular, the experimental study of Li et al. (2014), which used a volatile 0.65 cSt silicone oil ($Pr=9.2$), showed that for $Bo_{D}=O(1)$ the threshold values of $Ma$ increase rather dramatically as the air is removed from the experimental apparatus. The numerical simulations of this experimental setup (Qin et al., 2014, 2015; Qin & Grigoriev, 2015) based on a two-sided model, which took phase change into account and explicitly described the transport of heat and vapour through the gas phase, reproduced the experimental results. The main effect of air is to suppress phase change by impeding the transport of vapour towards, or away from, the interface, but its presence also affects thermal conductivity of the gas layer. This suggests that phase change and transport in the gas phase play an important role in this problem, especially at reduced pressures (reduced concentrations of air).

In order to better understand the mechanism of the instability and the effect of the gas phase we formulated a two-sided model of the flow valid for arbitrary composition of the gas phase; it is described in Section 2. The analytical solution describing the uniform return flow (in both phases) is derived in Section 3. The stability of that solution is investigated in Section 4. The results of linear stability analysis are compared with previous analytical, experimental, and numerical studies in Section 5 and are discussed in Section 6. Conclusions are presented in Section 7.

2 Mathematical Model

2.1 Governing Equations

We will assume that the liquid layer has depth $d_{l}$ and lateral dimensions much larger than $d_{l}$. Since in most experiments the system is covered by a horizontal plate to limit evaporation, we will assume that the gas layer also has a finite depth $d_{g}$. To describe convection in this two-layer system, we will use a variation of the two-sided model originally introduced by Qin et al. (2014) for near-atmospheric conditions and later extended by Qin & Grigoriev (2015) to the limit when the gas phase is dominated by the vapour, rather than air. A version of the model interpolating between the two limits is summarized below. The momentum transport in the bulk is described, for both the liquid and the gas phase, by the Navier-Stokes equation

[TABLE]

where $p$ is the fluid pressure, $\rho$ and $\mu$ are the fluid’s density and viscosity, respectively, $c_{a}$ is the concentration of air, and ${\bf g}$ is the gravitational acceleration. (The air is noncondensable, so $c_{a}=0$ in the liquid phase.) Following standard practice, we use the Boussinesq approximation, where the density is considered constant everywhere except in the last term on the right-hand-side representing the buoyancy force. In particular, the mass conservation equation in both phases reduces to

[TABLE]

To account for buoyancy, the density of the liquid phase is assumed to depend linearly on the temperature

[TABLE]

where $\rho^{0}_{l}$ is the density at the temperature $T_{0}$ describing global thermodynamic equilibrium, and $\beta_{l}=-\rho_{l}^{-1}\partial\rho_{l}/\partial T$ is the coefficient of thermal expansion. Here and below, subscripts $l$, $g$, $v$, $a$, and $i$ denote properties of the liquid and gas phase, vapour and air component, and the liquid-vapour interface, respectively. The index [math] will be used throughout the paper to distinguish the equilibrium values from the nonequilibrium ones, where the choice is not obvious. In the gas phase

[TABLE]

where both vapour and air are considered to be ideal gases

[TABLE]

${\bar{R}}_{a,v}=R/M_{a,v}$, $R$ is the universal gas constant, and $M_{a,v}$ is the molar mass of air/vapour. Correspondingly, the total pressure in the gas is

[TABLE]

Mass transport in the gas phase is described by the standard conservation equation(s)

[TABLE]

for the number density of the two components

[TABLE]

where $m_{a,v}=M_{a,v}/N_{A}$ is the mass of one air/vapour molecule. The number density flux is given by

[TABLE]

where the first and the second term on the right-hand-side represent the contributions due to advection and diffusion, respectively, $D$ is the binary diffusion coefficient, and

[TABLE]

are the concentrations (or, more precisely, the molar fractions) of the two components. Therefore, the conservation equation(s) (7) can be rewritten as

[TABLE]

From the ideal gas law, the total number density in the gas phase

[TABLE]

where $k_{B}=R/N_{A}$ is the Boltzman constant. As we showed previously (Qin et al., 2015), spatial variation of $p_{g}$ can in practice be neglected. Furthermore, the largest temperature variation in relevant experimental and numerical studies is around $5\%$ (and typically much smaller than that), such that $n_{g}$ can be assumed constant. Consequently (11) reduces to an advection-diffusion equation for the two concentrations

[TABLE]

These two equations are equivalent, so either one can be used, since $c_{a}+c_{v}=1$.

Finally, the transport of heat is also described using an advection-diffusion equation

[TABLE]

where $\alpha=k/(\rho C_{p})$ is the thermal diffusivity, $k$ is the thermal conductivity, and $C_{p}$ is the heat capacity, of the fluid.

2.2 Boundary Conditions at the Interface

The system of coupled evolution equations for the velocity, pressure, temperature, and concentration fields should be solved in a self-consistent manner, subject to the boundary conditions describing the balance of momentum, heat, and mass fluxes between the two phases. The phase change at the interface can be described using Kinetic Theory (Schrage, 1953). As we have shown previously (Qin et al., 2015), the choice of the phase change model has negligible effect on the results. The mass flux across the interface is given by (Klentzman & Ajaev, 2009)

[TABLE]

where $\chi$ is the accommodation coefficient, $\mathcal{L}$ is the latent heat of vaporization, and subscript $s$ denotes saturation values for the vapour. The dependence of the local saturation temperature $T_{s}$ on the partial pressure of vapour $p_{v}$ is described using the Clausius-Clapeyron equation for phase equilibrium

[TABLE]

where $p_{v}^{0}$ is the saturation pressure of the vapour at the equilibrium temperature $T_{0}$. The first term in (15) is proportional to the Young-Laplace pressure and can be ignored in this problem, since the interface is considered flat.

The mass flux balance at the interface can be expressed with the help of (7). In the reference frame of the interface, the mass flux of the vapour is given by

[TABLE]

where ${\bf u}_{i}$ is the velocity of the liquid at the interface. Since air is noncondensable, its mass flux across the interface is zero:

[TABLE]

Since $c_{a}+c_{v}=1$, these two relations can be solved yielding two of the boundary conditions for (1) and (7) in the gas phase

[TABLE]

and

[TABLE]

The heat flux balance is given by

[TABLE]

where the advective contribution to the heat flux is negligible on both sides of the interface. Indeed, in the gas phase, conduction is the dominant contribution (Qin & Grigoriev, 2015), while on the liquid side

[TABLE]

Since the liquid density is much greater than that of the gas, the left-hand-side of (22) is very small compared with $\hat{\bf n}\cdot({\bf u}_{g}-{\bf u}_{i})$ and can be ignored.

The remaining boundary conditions for $\bf u$ and $T$ at the liquid-vapour interface are standard: the temperature is continuous

[TABLE]

as are the tangential velocity components

[TABLE]

The stress balance

[TABLE]

incorporates both the viscous drag between the two phases and the thermocapillary effect. Here $\Sigma=\mu\left[\nabla{\bf u}+{(\nabla{\bf u})}^{T}\right]-p$ is the stress tensor, \kappa$$(=0) is the interfacial curvature, ${\nabla}_{s}=\left(\mathbb{I}-{\hat{\bf n}}\hat{\bf n}\right)\cdot\nabla$

is the surface gradient, and $\gamma=-\partial\sigma/\partial T$ is the temperature coefficient of surface tension.

We will further assume that side/top/bottom walls are adiabatic, $\hat{\bf n}\cdot\nabla T=0$, which is a good approximation for most experimental setups. The heat transport through the end walls at $x=0$ and $x=L_{x}$ is not treated explicitly in this study, but in principle a variety of boundary conditions, from constant temperature to constant flux to mixed boundary conditions can be accommodated. Standard no-slip boundary conditions ${\bf u}=0$ for the velocity and no-flux boundary conditions $\hat{\bf n}\cdot\nabla c_{a,v}=0$ for concentration apply on all the walls.

3 The Base Flow

In liquid layers that are not too thin, under normal gravity, convection is driven by both buoyancy and thermocapillarity. The strength of these two effects can be quantified, respectively, in terms of the Marangoni number

[TABLE]

and the Rayleigh number

[TABLE]

where $\nu_{l}=\mu_{l}/\rho_{l}$ is the kinematic viscosity of the liquid and $\tau$ is the horizontal component of the temperature gradient, assumed to be in the positive $x$ direction. In the numerics and experiment we set $\tau=\hat{\bf x}\cdot\nabla_{s}T_{i}$, which is found to be nearly independent of the location in the central region of the flow (Qin & Grigoriev, 2015). The dynamic Bond number

[TABLE]

is independent of $\tau$ and quantifies the relative strength of buoyancy and thermocapillarity. In defining nondimensional combinations, as well as various scales, we will use the values of material parameters at the equilibrium temperature $T_{0}$.

As discussed in the introduction, at sufficiently low $Ma$, a steady return flow is found in the liquid layer. The analytical solution of Birikh (1966) describes such a flow in laterally unbounded layers of nonvolatile liquids. This solution also describes the flow observed in laterally bounded layers away from the end walls (cf. Fig. 1) even for volatile fluids under atmospheric conditions when the phase change is strongly suppressed (Qin et al., 2014). For volatile fluids at reduced pressures, Birikh’s solution becomes invalid due to the increasing role of phase change (Li et al., 2014; Qin et al., 2014). Instead, we should look for a solution to the two-sided model described in the previous section, which describes the flow in both layers. In order to reduce the number of parameters, the governing equations (1), (2), (13), and (14) are nondimensionalized by introducing the length scale $d_{l}$, time scale $d_{l}^{2}/\nu_{l}$, velocity scale $\nu_{l}/d_{l}$, density scale $\rho_{l}$, pressure scale $\rho_{l}(\nu_{l}/d_{l})^{2}$, and temperature scale $\tau d_{l}=\mu_{l}\alpha_{l}Ma/(\gamma d_{l})$.

The dimensionless governing equations for the liquid layer become

[TABLE]

where $\tilde{\nabla}=d_{l}\nabla$, $\tilde{\bf u}_{l}={\bf u}_{l}d_{l}/\nu_{l}$, $\tilde{T}_{l}=(T_{l}-T_{0})/(\tau d_{l})$, and

[TABLE]

is the Grashof number. We will use a coordinate system defined such that the liquid-vapour interface corresponds to the plane $\tilde{z}=0$ (so that the liquid layer corresponds to $-1<\tilde{z}<0$). Recall that the $x$ axis points in the direction of the applied temperature gradient, with the origin chosen such that $T_{i}=T_{0}$ at $\tilde{x}=0$.

The dimensionless governing equations for the gas layer ($0<\tilde{z}<A$) are

[TABLE]

where $\tilde{\bf u}_{g}={\bf u}_{g}d_{l}/\nu_{l}$, $\tilde{T}_{g}=(T_{g}-T_{0})/(\tau d_{l})$, $\tilde{c}_{v}=c_{v}-c_{v}^{0}$, $c_{v}^{0}=1-c_{a}^{0}$, $A=d_{g}/d_{l}$, $K_{\nu}={\nu_{g}}/{\nu_{l}}$, $K_{\alpha}={\alpha_{g}}/{\nu_{l}}$, and $K_{D}=D/\nu_{l}$. The nondimensional combinations

[TABLE]

describe the contributions to the buoyancy force in the gas layer due to perturbations in the temperature and composition of the gas, respectively. Note that $\Xi_{T}$ depends on the imposed temperature gradient $\tau$, but not $c_{a}^{0}$, while $\Xi_{c}$ depends on $c_{a}^{0}$, but not $\tau$.

Both the imposed temperature gradient $\tau$ and the composition of the gas phase, parametrized by the equilibrium concentration of air $c_{a}^{0}$, play a key role in this problem. In experiment (Li et al., 2014), $c_{a}^{0}$ was controlled indirectly by varying the net gas pressure

[TABLE]

In the following analysis, it will be more convenient to describe the composition directly in terms of $c_{a}^{0}$. In terms of $p_{g}^{0}$, we can write $n_{g}\approx p_{g}^{0}/(k_{B}T_{0})$, while $\tilde{p}$ in (29) and (31) represents the nondimensional form of the difference $p_{g}-p_{g}^{0}$.

In order to satisfy the incompressibility condition, we will assume that the flow is strictly two-dimensional and introduce a stream function for each layer, such that

[TABLE]

Eliminating the pressure, the governing equations (29) for the liquid layer can be rewritten as

[TABLE]

For the gas layer we have

[TABLE]

3.1 Boundary Conditions

At the bottom of the liquid layer ($\tilde{z}=-1$) and the top of the gas layer ($\tilde{z}=A$), no-slip and adiabatic boundary conditions apply

[TABLE]

At the interface ($\tilde{z}=0$), the temperature and velocity fields are continuous

[TABLE]

where for our choice of the origin of coordinate system $\tilde{T}_{i}=\tilde{x}$. Since $\mu_{l}\gg\mu_{g}$ and typically $d_{g}>d_{l}$, the viscous stress in the gas layer can be ignored, yielding a simplified expression for the shear stress balance at the interface

[TABLE]

where

[TABLE]

is the Reynolds number.

The heat flux balance (21) at the interface reduces to

[TABLE]

where $k_{g}=c_{v}^{0}k_{v}+c_{a}^{0}k_{a}$, and $\tilde{J}=Jd_{l}/(Dm_{v}n_{g})$ is the dimensionless mass flux. The dimensionless combination

[TABLE]

or more precisely the ratio $V/Ma$, describes the relative magnitude of the latent heat released (absorbed) at the interface due to condensation (evaporation) compared with the vertical heat flux in the liquid layer due to conduction. It should be noted that, although the product $Dp_{g}^{0}$, and consequently $V$, is a function of $T_{0}$, it is independent of the gas pressure, and therefore $c_{a}^{0}$.

Nondimensionalizing (19) and (20), we obtain

[TABLE]

and

[TABLE]

The mass flux vanishes at the top of the gas layer

[TABLE]

The base uniform flow corresponds to a vanishing mass flux at the interface, $\tilde{J}=0$, which is also an assumption made by all one-sided models. The base flow is spatially uniform in the lateral direction(s) which, coupled with incompressibility, requires ${\hat{\bf n}}\cdot{\bf u}_{i}=0$. This leads to a number of simplifications. In particular, (43) and (44) require that $\tilde{u}_{l,z}=\tilde{u}_{g,z}=0$ and $\partial_{\tilde{z}}\tilde{c}_{v}=0$ at the interface. Furthermore, the heat flux at the interface should also vanish, so that (21) gives $\partial_{\tilde{z}}\tilde{T}_{l}=\partial_{\tilde{z}}\tilde{T}_{g}=0$.

The boundary conditions for the stream function can be easily obtained from those for the velocities. At the bottom of the liquid layer and the top of the gas layer

[TABLE]

For a uniform flow, the net flux through any vertical plane vanishes. Combined with conditions (24) and (25) this requires that at the interface

[TABLE]

3.2 Fluid Flow and Temperature in the Liquid Layer

For a uniform horizontal flow, where both $\tilde{\psi}_{l}$ and $\tilde{\psi}_{g}$ are functions of $\tilde{z}$ alone, we can look for solutions to (35) in the form

[TABLE]

where $\tilde{\theta}_{l}(0)=0$. With this choice, the system (35) reduces to coupled ordinary differential equations (ODEs)

[TABLE]

where prime stands for the derivatives with respect to the $\tilde{z}$ coordinate.

Solving the system (49) subject to the boundary conditions at the bottom and the free surface of the liquid layer, we find the steady state solutions for the stream function

[TABLE]

velocity

[TABLE]

and temperature field

[TABLE]

describing the base flow. They agree with the analytical solutions originally obtained by Birikh (Birikh, 1966) and later rederived by Kirdyashkin (Kirdyashkin, 1984) and Villers and Platten (Villers & Platten, 1987) using a one-sided model that ignores the effects of the gas phase.

The assumption that the interfacial temperature varies linearly has been widely used in previous studies, without much justification, for deriving the solutions (51) and (52) for the return flow underlying the stability analyses (Parmentier et al., 1993; Mercier & Normand, 1996; Priede & Gerbeth, 1997) as well as in models of the adiabatic section of heat pipes (Ha & Peterson, 1994; Suman & Kumar, 2005; Markos et al., 2006), which assume unidirectional flow in the liquid phase. However, the validity of this assumption cannot be established by a one-sided model which ignores heat and mass transport in the gas phase. In fact, when $c_{a}^{0}$ becomes sufficiently low, the interfacial temperature profile becomes nonlinear (Li et al., 2014; Qin & Grigoriev, 2015). Proper justification of the linearity assumption requires showing that it is consistent with a steady-state solution of the transport equations in the gas phase, which satisfies all of the boundary conditions at the free surface. We turn to this next.

3.3 Fluid Flow, Temperature, and Composition in the Gas Layer

The solutions for the velocity, temperature, and composition of the gas phase can be found in the same way the solutions (51) and (52) were obtained for the liquid phase. We will start by finding the solution for the vapour concentration at the interface before deriving the solution in the bulk. Since there is no phase change, $T_{s}=T_{i}$, according to (15), and $\partial_{x}T_{s}=\tau$. Using Clausius-Clapeyron relation (16) and neglecting the deviation of $p_{g}$ from the equilibrium value (33) we therefore find

[TABLE]

where for $\tau x\ll T_{0}$ we can replace $T_{i}$ with the equilibrium temperature $T_{0}$. In general, the solution to this equation yields an exponential concentration profile for $c_{v}$ and hence a nonlinear interfacial temperature profile (Qin & Grigoriev, 2015). However, when the deviation of $c_{v}$ from $c_{v}^{0}$ is small, the concentration profile at the interface becomes approximately linear

[TABLE]

where

[TABLE]

and

[TABLE]

is another nondimensional parameter, the meaning of which will become clear later. Again, since the product $Dp_{g}^{0}$ is independent of $c_{a}^{0}$, so is $H$.

Given the boundary conditions for $\tilde{T}_{g}$ and $\tilde{c}_{v}$, we can look for solutions to these two fields in the gas layer of the form

[TABLE]

where

[TABLE]

For a uniform flow, $\tilde{\psi}_{g}$ should only depends on $\tilde{z}$, so the system (36) reduces to

[TABLE]

Solving these equations subject to the boundary conditions (45), (46), (47), and (3.3) at the interface and the top of the gas layer yields the steady state solutions for the stream function

[TABLE]

velocity

[TABLE]

temperature

[TABLE]

and vapour concentration in the gas phase

[TABLE]

The parameters $\mathcal{R}$ and $\mathcal{B}$ are analogous to the Reynolds number and the dynamic Bond number, but incorporate the properties of both fluid layers:

[TABLE]

where we defined

[TABLE]

Note that the form of the analytical solutions (51)-(52) and (61)-(63) is different from that of the solutions derived by Qin et al. (2014) because (i) the buoyancy force caused by the variation in the composition of the gas is explicitly taken into account in the present analysis, (ii) the mass transport in the gas phase is described using the concentration rather than the mass density of the vapour, and (iii) a different choice of the origin of the coordinate system is made. Figure 2 compares the analytical solutions with the numerical ones (cf. Fig. 1) computed by Qin & Grigoriev (2015). We find that the analytical solutions accurately describe the core region of the flow away from the end walls for sufficiently low $\tau$ at which the base flow is stable.

4 Linear Stability Analysis

The main focus of our study is on relatively thick layers where the effect of buoyancy is nonnegligible, such that $Bo_{D}=O(1)$. In this parameter regime gravity is strong enough to keep the liquid interface nearly flat, so in our analysis we will assume that the thickness of the liquid layer remains uniform even after the onset of convection pattern. To date, only one study (Priede & Gerbeth, 1997) of pattern formation in buoyancy-thermocapillary convection correctly predicted the formation of a stationary pattern of convection rolls observed in experiment at $Bo_{D}=O(1)$. This study was based on a one-sided model which completely ignored phase change and assumed adiabatic conditions at the free surface. While such description may be acceptable for nonvolatile liquids or at high concentrations of air, it fails to describe volatile liquids at lower concentrations of air. Below a generalized analysis is presented, which accounts for both heat and mass flux across the interface associated with phase change and which is applicable regardless of the amount of air present in the gas layer (quantified by $c_{a}^{0}$).

4.1 Diffusion-Dominated Case

While the nondimensional equations (29) for the liquid layer can be used without modification, in order to obtain a simplified problem, which can generate useful physical insight, we will start by making several approximations in the treatment of heat and mass transport on the gas side. The relative contribution of advection and diffusion to the mass and heat transport in the gas layer are described by the Péclet numbers $Pe_{m}$ and $Pe_{t}$, respectively. Since the relevant length scale here is the wavelength of the convective pattern, which is comparable to the depth of the liquid layer, and the velocity scale is determined by the interfacial velocity $u_{i}=\nu_{l}\tilde{u}_{l}(0)/d_{l}$ (cf. Eq. (51)), we have

[TABLE]

and

[TABLE]

In typical experiments (Villers & Platten, 1992; Riley & Neitzel, 1998; Li et al., 2014) both Péclet numbers are at most $O(1)$, so we can drop the advective terms in the transport equations (13) and (14) in the gas phase. Furthermore, if we are only interested in the transition thresholds for stationary instability, we can also ignore the time derivatives, such that the system (36) simplifies to

[TABLE]

Note that neglecting the transient dynamics of the gas phase changes neither the critical Marangoni number nor the critical wavelength of a (stationary) pattern.

The base solution $\tilde{\psi}_{l0}$, $\tilde{T}_{l0}$ describing the uniform return flow in the liquid layer is given by (50)-(52). The base solution to (68) describing the gas layer in the diffusion-dominate case can be obtained by setting $\mathcal{R}=0$ in (62)-(63) which is equivalent to setting $\tilde{\bf u}_{g}=0$ ($\tilde{\psi}_{g}=0$) and gives:

[TABLE]

Although the base state satisfies the adiabatic boundary condition at the interface, perturbations in the temperature will give rise to heat and mass flux across the interface and, consequently, through the gas phase. Perturbed solutions can be written in the form of Fourier integrals

[TABLE]

where $m=\{l,g\}$ denotes the phase, $q$ is the wavenumber describing the variation in the horizontal direction, and $\sigma_{q}$ is the temporal growth rate. The functions $\tilde{\psi}_{mq}(\tilde{z})$, $\tilde{\theta}_{mq}(\tilde{z})$, and $\tilde{\varsigma}_{vq}(\tilde{z})$ define the vertical profile for the perturbations in, respectively, the stream function $\tilde{\psi}$, temperature $\tilde{T}$, and vapour concentration $\tilde{c}_{v}$ with wavenumber $q$.

Rewriting (68) in terms of the perturbations, we obtain

[TABLE]

Temperature continuity at the interface requires

[TABLE]

At the top of the gas layer we have

[TABLE]

The solution of (71) satisfying these boundary conditions is

[TABLE]

The solution for the perturbation in the vapour concentration is analogous,

[TABLE]

where $\tilde{\varsigma}_{vq}(0)$ can be related to the perturbation in the interfacial temperature via the Clausius-Clapeyron relation (16)

[TABLE]

Fourier transforming (43) yields

[TABLE]

where $\tilde{J}_{q}$ is the Fourier coefficient of $\tilde{J}(\tilde{x})$. Subtracting the base solutions from the heat balance (41), Fourier transforming the result, and using (77) gives the following relation

[TABLE]

With the help of (74), (75), and (76) this can rewritten in the form of Newton’s cooling law

[TABLE]

where we introduced a wavenumber-dependent analogue of the Biot number

[TABLE]

which will be referred to as the Biot coefficient below to highlight the fact that it is a function of $q$ and $c_{a}^{0}$, not a constant. The prefactor $\tanh(qA)$ describes the effect of finite thickness of the gas layer. For $d_{g}\gg d_{l}$ ($A\gg 1$), it approaches unity and (80) reduces to the expression derived by Chauvet et al. (Chauvet et al., 2012) in the context of stability of a volatile liquid layer in the presence of a vertical temperature gradient. The first term describes the effect of thermal conduction through the gas layer, while the second term describes the effect of latent heat released/absorbed at the interface as a result of condensation/evaporation.

By linearizing the governing equations (35) around the base state (50)-(52), we obtain the evolution equations for the perturbations $\tilde{\psi}_{lq}$ and $\tilde{\theta}_{lq}$ in the liquid layer

[TABLE]

where we defined $\tilde{\nabla}_{q}^{2}=\partial_{\tilde{z}}^{2}-q^{2}$ and

[TABLE]

This is a system of ODEs which is fourth-order with respect to $\tilde{\psi}_{lq}$ and second-order with respect to $\tilde{\theta}_{lq}$ and hence needs a total of six boundary conditions. These boundary conditions are:

[TABLE]

at the bottom of the liquid layer and

[TABLE]

and

[TABLE]

at the free surface. In the subsequent discussion we will refer to the system (81) with the boundary conditions (83)–(85) as the enhanced one-sided model, since it incorporates the effect of the gas phase entirely through boundary conditions at the interface.

4.2 Transient Dynamics in the Gas Layer

A more accurate description of the instability, stationary or oscillatory, can be obtained by restoring the time dependence of the temperature and composition of the gas phase. This corresponds to replacing the Laplace equations (68) with

[TABLE]

The corresponding equations for the perturbations are

[TABLE]

These equations should be solved subject to the boundary conditions (72), (73), (76), (83), (84) and

[TABLE]

which follows from the heat flux balance (41) and the mass flux balance (77) and replaces (79).

Since it involves an incomplete description of the gas layer (the effects of advection are not included), this model based on the governing equations (81) and (87) can be called “one-and-a-half-sided.” Note that, since both $K_{\alpha}$ and $K_{D}$ are typically large compared with unity, the terms on the right-hand-side of (87) are small. Dropping them would lead to the enhanced one-sided model derived in Section 4.1.

4.3 The Effect of Advection in the Gas Layer

The effect of advection can also be incorporated for the range of $c_{a}^{0}$ where the base solution (60)-(63) is valid. By setting $\tilde{\psi}_{gq}\neq 0$ we ensure that all of the transport mechanisms in the gas phase are accounted for and both layers are treated using an equally comprehensive model. The linearized evolution equations for the perturbations $\tilde{\psi}_{gq}$, $\tilde{\varsigma}_{vq}$, and $\tilde{\theta}_{gq}$ in the gas phase are

[TABLE]

where

[TABLE]

The boundary conditions are given by (72), (73), (76), (83), (84), and (88). In addition, the boundary conditions for $\bar{\psi}_{gq}(\bar{z})$ follow from (46) and (47):

[TABLE]

at the top of the gas layer and

[TABLE]

at the free surface. In the subsequent discussion we will refer to the system (81) and (89) with the appropriate boundary conditions as the two-sided model, since it treats both the liquid and the gas phase with the same level of detail.

4.4 Different Modes of Instability

The boundary value problem, whether it is described by the system (81) by itself or combined with (87) or (89), has a spectrum of eigenvalues $\sigma_{q}^{n}$ and the corresponding eigenfunctions $\{\psi_{lq}^{n},\theta_{lq}^{n},\cdots\}$, $n=0,1,\cdots$, for a given complex wavenumber $q=k+is$ and Marangoni number $Ma$. The stability of the base flow is determined by the eigenvalue $\sigma_{q}^{0}=\kappa_{q}+i\omega_{q}$ with the largest real part $\kappa_{q}$ and the character of the instability (oscillatory or stationary) is determined by the imaginary part $\omega_{q}$.

So far we have only focused on the boundary conditions in the vertical (confined) direction. In both experiment and simulations the fluid layers are also bounded in the horizontal (extended) directions $0<x<L_{x}$ and $0<y<L_{y}$. Comparison of various experimental results suggests that the spanwise aspect ratio $\Gamma_{y}=L_{y}/d_{l}$ is not particularly important: aside from the orientation of hydrothermal waves at low $Bo_{D}$, the pattern that emerges above the threshold of primary instability is basically the same for both $\Gamma_{y}=O(1)$ (Villers & Platten, 1992; De Saedeleer et al., 1996; Garcimartín et al., 1997; Li et al., 2014) and $\Gamma_{y}\gg 1$ (Riley & Neitzel, 1998; Burguete et al., 2001). This is not very surprising since the boundaries in the spanwise direction play a passive role – they tend to suppress, rather than enhance convection.

The boundaries in the streamwise direction are less trivial and play a crucial role in selecting the convection pattern for any $\Gamma_{x}=L_{x}/d_{l}$. The corresponding boundary conditions are of the enhancing type (Cross & Greenside, 2009): buoyancy, which plays a significant role for $Bo_{D}=O(1)$, drives an upward flow near the hot end wall, generating a single localized convection roll. Similarly, buoyancy drives a downward flow near the cold end wall, generating another, much weaker, convection roll. These rolls appear for $Ma$ well below the critical value $Ma_{c}$ at which convection pattern emerges away from the end walls. The rolls are steady over a wide range of $Ma$ and have a horizontal extent of $O(d_{l})$. For $\Gamma_{x}\gg 1$ and sufficiently low $Ma$, the flow in the core region between these two convection rolls, i.e., for $O(d_{l})<x<L_{x}-O(d_{l})$, is described well by the base flow solution (51) and (61), as mentioned previously (cf. Fig. 1). Riley & Neitzel (1998) refer to this state as a steady unicellular flow (SUF). Since the perturbation strength is always finite at the end walls, defining the threshold of the primary instability at finite $\Gamma_{x}$ requires some care.

4.4.1 Convective Instability

The dominant mode of primary instability in laterally unbounded systems (infinite $\Gamma_{x}$) tends to be convective. The threshold of convective instability corresponds to $\kappa_{q}$ being negative over all $q$-real, except for $q=k_{c}$, where $k_{c}$ is the (real) critical wavenumber which defines the wavelength $\lambda/d_{l}=2\pi/k_{c}$ of the pattern that develops above threshold. Therefore, the critical wavenumber $k_{c}$ and the critical Marangoni number $Ma_{c}$ are given by the solution of the two equations

[TABLE]

over $q$-real. Generally, the critical frequency $\omega_{c}=\omega_{q}|_{q_{c}}\neq 0$, so the resulting pattern will be oscillatory, with disturbances amplified in the frame of reference moving with the group velocity (Huerre, 2000)

[TABLE]

In laterally bounded systems, such as the fluid volumes with finite aspect ratios $\Gamma_{x}$ discussed in this paper, convective instabilities should not be observed in the absence of a localized oscillatory source of the right temporal frequency ($\omega\approx\omega_{c}$). An example of such a source could be an oscillating convection roll near the left (right) end wall for $v_{g}>0$ ($v_{g}<0)$ driven by a completely different mechanism (e.g., buoyancy due to cooling or heating at the end wall). While such oscillating convection rolls have indeed been observed in experiment (De Saedeleer et al., 1996; Garcimartín et al., 1997; Li et al., 2014) and simulations (Qin et al., 2014), they arise after convection pattern in the core region of the flow has already been established. Hence, this is a secondary, rather than a primary, instability. We should only expect to see a primary convective instability in systems with very large streamwise aspect ratio $\Gamma_{x}$ where infinitesimal disturbances have sufficient time $\Gamma_{x}/v_{g}$ to develop.

4.4.2 Absolute Instability

Absolute instabilities arise when infinitesimal disturbances grow in the laboratory frame. Since they do not need a finite-amplitude source, they can arise both in laterally unbounded systems as well as in bounded systems with moderate $\Gamma_{x}$. Absolute instability requires that the group velocity $v_{g}=0$ at onset, yielding a system of three real equations for three real unknown $k_{c}=\mathrm{Re}(q_{c})$, $s_{c}=\mathrm{Im}(q_{c})$, and $Ma_{c}$, which can be rewritten as a system of two equations (Huerre, 2000)

[TABLE]

the second of which is complex, for $q_{c}$-complex and $Ma_{c}$.

The resulting pattern will generally be oscillatory, $\omega_{c}\neq 0$, just like in the case of convective instability. Although its group velocity vanishes, the phase velocity $v_{p}=-\omega_{c}/k_{c}$ of the pattern will be nonzero. A good example of such an absolute instability is hydrothermal waves (Smith & Davis, 1983a, b) that are observed in the limit of low $Bo_{D}$ (Daviaud & Vince, 1993; Riley & Neitzel, 1998; Burguete et al., 2001). In the following, where necessary to avoid confusion, we will use the superscript $a$ to denote the critical values corresponding to absolute instability.

4.4.3 Stationary Instability

The threshold of a stationary instability corresponds to both the real and imaginary part of $\sigma_{q}^{0}$ vanishing

[TABLE]

These two equations should be solved for the critical wavenumber $k_{c}$ and the critical Marangoni number $Ma_{c}$, which both depend on the spatial attenuation rate $s$. In the following, where necessary to avoid confusion, we will use the superscript $s$ to denote the critical values corresponding to stationary instability. In particular, we will write $Ma_{c}^{s}(s)$ to explicitly denote the dependence on $s$.

In laterally infinite systems the transition threshold is well-defined and corresponds to a spatially uniform stationary pattern ($s=0$). Such a pattern will be referred to as the steady multicellular (SMC) state below. For $Ma<Ma_{c}^{s}(0)$ the convection pattern arises due to forcing at one of the boundaries (convection rolls driven by buoyancy, cf. Fig. 1(c) and (d)), rather than spontaneously, and generally does not extend across the entire length of the system. In this case the absolute value of $s=\mathrm{Im}(q_{c})$ defines the spatial attenuation of the perturbation (or the number of convection rolls that can be detected in the liquid layer), while $k_{c}=\mathrm{Re}(q_{c})$ determines the wavelength $\lambda$ of the pattern (or the distance between convection rolls). For $s\neq 0$ the pattern appears below $Ma_{c}^{s}(0)$ and is localized near the right (hot) end wall for $s<0$ or near the left (cold) end wall for $s>0$. As we will see below, $dMa_{c}^{s}/ds>0$, which means that the pattern always appears near the hot end wall (cf. Fig. 1(d)) and we can only consider $s\leq 0$.

When $e^{-|s|\Gamma_{x}}\ll 1$, the pattern extends over a region of length $O(2d_{l}/|s|)<L_{x}$. We will refer to this as partial multicellular (PMC) state following Qin & Grigoriev (2015) and Li et al. (2014). In this analysis, we set $s_{\mathrm{PMC}}=-1$ to define the transition from SUF to PMC at $Ma_{c}^{s}(-1)<Ma_{c}^{s}(0)$, which roughly corresponds to one extra roll nucleating near the (hot) end wall. As $Ma$ is increased beyond this value, the convection pattern extends further away from the end wall, until it eventually covers the entire fluid layer. In practice this roughly corresponds to $e^{-|s|\Gamma_{x}}\approx 0.1$ or $s_{\mathrm{SMC}}=\pm 2/\Gamma_{x}$. This is the definition that has to be used for $Pr\lesssim 1$ (i.e., for gases or metallic liquids featuring sensitive dependence of $Ma_{c}^{s}$ on $s$ for small $|s|$, as shown below). For typical (nonmetallic) liquids with $Pr=O(10)$, we can set $s_{\mathrm{SMC}}=0$, as in the case of laterally unbounded layers, as long as $\Gamma_{x}\gtrsim 20$. This is the case in both the experimental (Riley & Neitzel, 1998; Li et al., 2014) and numerical studies (Qin et al., 2014; Qin & Grigoriev, 2015) used for comparison below.

Note that, for $\sigma_{q}^{0}=0$, equations (87) reduce to (71) even if $K_{\alpha}$ and $K_{D}$ are not large. Hence the two formulations always give identical predictions for the threshold of stationary instability (even though they predict different growth rates for modes with wavenumbers $q$ different from the critical one). Therefore, in what follows we will only focus on the limits when either both $Pe_{m}$ and $Pe_{t}$ are vanishingly small and stability is described by (81) alone, or when they are both $O(1)$, so (81) and (89) have to be solved simultaneously.

5 Results

Previous theoretical studies of stability focused mainly on the dependence of the critical Marangoni number $Ma$ on the dynamic Bond number $Bo_{D}$ and Prandtl number $Pr$ characterizing the liquid layer. We have also investigated how both the critical $Ma$ and the critical wavelength $\lambda$ depend on the equilibrium concentration $c_{a}^{0}$ of air in the gas layer. Where possible, the predictions of linear stability analysis have been compared with the available experimental and numerical data.

The boundary value problems, which determine $\sigma_{q}^{0}$ as a function of the wavenumber $q$ and various nondimensional parameters, were solved using the function bvp5c in Matlab 2016a (cf. Appendix A for details) and the equations (93), (95), and (96), defining the conditions for various types of instabilities where solved using the Newton method.

5.1 Convection at atmospheric conditions

The most comprehensive experimental study to date that investigated convection patterns in the problem considered here is due to Riley & Neitzel (1998), who determined the dependence of both the critical $Ma$ and the critical wavelength $\lambda$ corresponding to the onset of various instabilities on the dynamic Bond number $Bo_{D}$ for the 1 cSt silicone oil with $Pr=13.9$. Since the authors have not reported the values of all material parameters of the liquid, we took the missing ones from Kavehpour et al. (2002). No data has been found for the material properties of the vapour; these were assumed to be the same as those of the 0.65 cSt silicone oil (Yaws, 2003, 2009). Any potential uncertainty in the values of these material parameters, however, will have a negligible effect on the results. The 1.0 cSt silicone oil is not particularly volatile (the saturation vapour pressure is $p_{v}=0.5$ kPa at $T=298$ K, which corresponds to $c_{a}^{0}=0.995$ at atmospheric conditions), therefore the properties of the gas phase are essentially those of air.

The predictions of the linear stability analysis for this particular fluid are shown in Figs. 3 and 4. The enhanced one-sided model and the two-sided model produce essentially identical results at atmospheric conditions; we only show the results of the latter. In particular, the thresholds of different types of instabilities for a laterally unbounded system ($\Gamma_{x}=\infty$) are shown as solid lines in Fig. 3(a). These are in decent agreement with the thresholds of convective, absolute, and stationary instabilities predicted by Priede & Gerbeth (1997) using an adiabatic one-sided model (which is equivalent to our enhanced one-sided model with $Bi_{q}=0$). As expected, convective instability is predicted to occur at the lowest value of $Ma$ regardless of $Bo_{D}$. Fig. 3(b) shows the group velocity $v_{g}$ associated with convective instability (along with the phase velocity $v_{p}$ associated with absolute instability). In the absence of end walls, one should expect a pattern arising from convective instability to develop sooner on one side than the other, with the interface between patterned and unpatterned flow propagating in the direction controlled by the sign of $v_{g}$, i.e., in the direction of the thermal gradient for $Bo_{D}<0.18$ and opposite the thermal gradient for $Bo_{D}>0.18$.

In practice, the aspect ratio $\Gamma_{x}$ is too low for convective instability to amplify tiny disturbances that are naturally present in the flow sufficiently so they could be observed before they collide with one of the walls and disappear. For instance, in the study of Riley & Neitzel (1998) $12<\Gamma_{x}<40$, while in most other experimental and numerical studies $\Gamma_{x}\leq 20$. As a result, convective instability would not be observed (and will not be considered in detail here). Instead an absolute instability, a stationary instability, or a combination of both instabilities would occur, and the order in which these instabilities occur would determine the type of convection pattern. In particular, hydrothermal waves (HTW) are usually found for $Bo_{D}\to 0$ and a stationary pattern is found for $Bo_{D}=O(1)$. However, we are not aware of any studies that predict what the pattern might look like at the intermediate values of $Bo_{D}$. The analysis for large, but finite $\Gamma_{x}$ (so there are end walls) is presented below.

Comparison of the theoretical predictions with the results of Riley & Neitzel (1998) is complicated by the fact that in the experimental study the transition threshold for SMC is identified as the instant when multiple convection rolls appear (near the hot end wall), but do not extend all the way to the cold end wall, which is not a very precise definition. However, the experimental data, which corresponds to an unknown value $-1<s<0$, should lie between the analytical prediction for the onset of PMC ($s=-1$) and SMC ($s=0$). As shown in Fig. 4, both the critical $Ma$ and the critical wavelength $\lambda$ observed in the experiments are indeed bracketed by the theoretical values corresponding to the onset of PMC and SMC, and lie closer to the SMC boundary in the entire range of $Bo_{D}$ investigated in the experiment, where stationary convective patterns are observed. Riley and Neitzel’s supporting figures show the presence of numerous convection rolls at what they define to be the critical $Ma$, so their data points should indeed be closer to the SMC threshold then the PMC threshold.

For $Bo_{D}<0.04$ the absolute instability boundary lies below that of PMC, $Ma_{c}^{a}<Ma_{c}^{s}(-1)$, so a pattern should appear at $Ma_{c}^{a}$ over the entire extent of the fluid layer and feature traveling convection rolls, since $\omega_{c}\neq 0$. The predicted phase velocity is positive, which is consistent with hydrothermal waves that are known to propagate from the cold end wall to the hot end wall. All of this is in perfect agreement with the analysis of Smith & Davis (1983a) applicable in the limit $Bo_{D}\to 0$. Riley and Neitzel’s experimental observations for $0.04<Bo_{D}<0.29$ deserve a separate discussion, since this is the range where the threshold of absolute instability lies between the PMC and SMC thresholds, $Ma_{c}^{s}(-1)>Ma_{c}^{a}>Ma_{c}^{s}(0)$. In this range the base flow can undergo an instability with respect to two different patterns in different regions. Consider the experimental protocol in which $Ma$ is gradually increased at a fixed $Bo_{D}$. Our analysis predicts that a stationary pattern first develops near the hot end wall at the SMC threshold $Ma_{c}^{s}(-1)$. As $Ma$ is increased beyond this value, the pattern spreads towards the cold end wall until $Ma$ reaches the threshold $Ma_{c}^{a}$ of absolute instability, which should cause hydrothermal waves to appear near the cold end wall. Hence the fraction of the fluid layer occupied by the two patterns should then depend on $Bo_{D}$ and $\Gamma_{x}$, with HTW dominating for smaller $Bo_{D}$ and PMC at higher $Bo_{D}$, which is qualitatively consistent with experimental observations. As Fig. 4(a) illustrates, our quantitative predictions for $Ma_{c}$ and $\lambda$ are also consistent with the experimental results in this range of $Bo_{D}$.

If $Ma$ is increased further beyond $Ma_{c}^{s}(0)$, the SMC pattern itself eventually becomes unstable and is replaced by the oscillatory multicellular convection (OMC) state (Villers & Platten, 1992; De Saedeleer et al., 1996; Garcimartín et al., 1997; Riley & Neitzel, 1998; Qin et al., 2014). OMC features convection rolls traveling in the direction opposite to that of HTW, i.e., against thermal gradient, with the right-most convection roll acting as the oscillatory source of disturbance (Garcimartín et al., 1997; Qin et al., 2014; Li et al., 2014). While it may be tempting to associate the extension of the curve $Ma=Ma_{c}^{a}$ for $Bo_{D}>0.29$ in Fig. 4(a) with the boundary between SMC and OMC, it is incorrect to do so, since the absolute instability describes destabilization of the uniform base flow rather than the SMC pattern. On the other hand, the group velocity $v_{g}$ associated with convective instability appears to predict the location of the boundary between HTW and OMC (the dashed line in Fig. 4(a)). As mentioned previously, $v_{g}>0$ for $Bo_{D}<0.18$, such that disturbances originating near the cold end wall should amplify as they propagate towards the hot end wall, while $v_{g}<0$ for $Bo_{D}>0.18$, so the opposite is true, which is consistent with the location of the wave sources for HTW and OMC. For comparison, in the experiment the boundary is at $Bo_{D}=0.22$.

Next we compare the analytical predictions with the results of numerical and experimental studies (Qin et al., 2014; Qin & Grigoriev, 2015; Li et al., 2014), which used a more volatile liquid (0.65 cSt silicone oil with $Pr=9.2$) and characterized both the SMC$\to$PMC and the PMC$\to$SMC transition. The predicted dependence of the critical Marangoni number and wavelength on the dynamic Bond number at atmospheric conditions ($c_{a}^{0}=0.96$) is shown in Fig. 5. Again, the predictions of the enhanced one-sided and two-sided model are essentially indistinguishable and are in reasonable agreement with both experimental and numerical data. Just as in the case of 1 cSt silicone oil, the critical values of $Ma$ and $\lambda$ increase monotonically with $Bo_{D}$ over the range where a stationary pattern is found ($Bo_{D}>0.31$). This increase in $Ma_{c}$ can be easily understood by noting that buoyancy has a stabilizing effect on the flow. The base flow solution (52) shows that the temperature increases, and hence the density of the liquid decreases, with height. An increase in $Bo_{D}$ therefore corresponds to an increasingly strong effect of the stable density stratification.

At low $Bo_{D}$ we again find that the absolute instability precedes the stationary instability, and we should expect to find HTW, at least near the cold end wall. According to Fig. 5(a), for $Bo_{D}<0.05$, the threshold of absolute instability lies below the PMC threshold, $Ma_{c}^{a}<Ma_{c}^{s}(-1)$, so HTW should be observed over nearly the entire extent of the fluid layer. For $0.05<Bo_{D}<0.31$, the threshold of absolute instability lies between the PMC and SMC thresholds, $Ma_{c}^{s}(-1)<Ma_{c}^{a}<Ma_{c}^{s}(0)$. As we discussed previously, in this range we expect the two patterns to coexist for sufficiently large $\Gamma_{x}$, with the boundary between HTW and PMC at onset of the absolute instability gradually shifting towards the cold end wall as $Bo_{D}$ increases. The competition between these patterns above the threshold should be described by nonlinear theory and hence is not addressed here.

The dependence of the threshold values of $Ma$ and $\lambda$ on the Prandtl number of the liquid has previously been investigated only in the context of linear stability analysis. In particular, Priede & Gerbeth (1997), who only considered uniform disturbances ($s=0$) and ignored heat and mass flux across the interface, found that the critical values of $Ma$ for the stationary instability diverges at $Pr_{0}=O(1)$ and reaches a minimum at $Pr_{1}$ just slightly higher than $Pr_{0}$. (Oddly enough, they only show the data for $Bo_{D}=0$, where HTW are expected instead of SMC.) Our enhanced one-sided model predicts that this result also holds for $Bo_{D}=O(1)$, once transport in the gas phase has been taken into account (cf. Fig. 6). The corresponding data for the absolute instability at lower $Bo_{D}$ is shown in Fig. 7. Predictions of the two-sided model are not shown, since it involves too many parameters which will vary with $Pr$, as different $Pr$ correspond to different liquids.

Qualitatively similar trends are found for the threshold of SMC for all $Bo_{D}=O(1)$, with both $Pr_{0}$ and $Pr_{1}$ slowly increasing with $Bo_{D}$. Over the entire range of $Pr$ corresponding to nonmetallic liquids, $Ma_{c}$ is a monotonically increasing function of $Pr$. The dependence of the critical wavelength on the Prandtl number is nonmonotonic even for nonmetallic liquids; $\lambda$ has a minimum at $Pr_{2}=O(10)$ and also diverges at $Pr_{0}$. The divergence of the critical values of $Ma$ and $\lambda$ at $Pr_{0}$, however, does not imply that no stationary pattern emerges at $Pr<Pr_{0}$. Indeed, the critical $Ma$ for PMC ($s=-1$) decreases monotonically as $Pr$ decreases, with $Ma_{c}\to 0$ as $Pr\to 0$. For instance, for a moderate value of $Ma\approx 600$, a stationary (spatially nonuniform) pattern is predicted to appear near the hot end wall for any $Pr<Pr_{0}$. As $Ma$ increases, the spatial extent of the pattern should increase as well.

As noted in the introduction, the analysis of Mercier & Normand (2002) predicted that, steady convection rolls should develop near the hot end for $Pr>4$, near the cold end for $Pr<0.01$, or at both end walls for $0.01<Pr<4$. Our results unequivocally contradict that prediction for $Bo_{D}=O(1)$. According to Fig. 6(b), for $Pr\lesssim 0.1$ the wavelength of the pattern exceeds typical values of $\Gamma_{x}$ for any value of the the spatial attenuation rate $s$, so it makes no sense to discuss spatial localization. The dependence of the critical $Ma$ and $\lambda$ on $s$ for the stationary instability at $Pr=0.1$, 1, and 10 are shown in Fig. 8). We find that $Ma_{c}^{s}$ is a monotonically increasing function of $s$, regardless of the values of $Pr$ for $Bo_{D}=O(1)$, indicating that a stationary pattern always emerges first near the hot end wall.

The analysis of Mercier & Normand (2002) was performed mainly to understand the results of numerical simulations by Ban Hadid & Roux (1990) in the limit of low $Pr$ which indicated the formation of a convective pattern near the cold end wall. Both of these studies made two contradicting assumptions: (1) that the dynamic Bond number $Bo_{D}\ll 1$ and (2) that the liquid-gas interface is flat. As Smith & Davis (1983b) showed, in the limit $Bo_{D}\to 0$ the uniform base flow undergoes an instability towards surface waves for $Pr<0.15$, violating the assumption of flat interface before a stationary instability can occur. Indeed, as Fig. 7(b) shows, even for finite $Bo_{D}$ the critical wavelength associated with the absolute instability diverges as $Pr$ approaches $Pr_{0}=O(1)$ from above, indicating that for $Pr<Pr_{0}$ a long-wavelength instability will arise that necessarily causes substantial deformation of the interface. Even for $Pr>Pr_{0}$, no stationary pattern will appear near the cold end wall for small $Bo_{D}$, since the absolute instability precedes the stationary one, as Figs. 4 and 5 show.

In conclusion of this section, let us point out that even though phase change is greatly suppressed at atmospheric conditions, this does not mean that the effect of the gas phase is negligible. Fig. 3 clearly shows the difference in the thresholds of all three instability types predicted by the adiabatic one-sided model (which completely ignores the effect of the gas layer) and the two-sided model (which accounts for heat, mass, and momentum transport in the gas) for a fluid of relatively low volatility (1 cSt silicone oil). The difference is even more significant for the more volatile 0.65 cSt silicone oil, as Fig. 5 illustrates. We find that the one-sided model underestimates the critical values of both $Ma$ and $\lambda$ which correspond to the PMC/SMC boundary by between 4% at $Bo_{D}\approx 0.3$ and 11% at $Bo_{D}\approx 1$. The corresponding discrepancy for the absolute stability boundary ranges from 3% at $Bo_{D}=0$ to 10% for $Bo_{D}\approx 0.3$. Not surprisingly, we find that transport of heat and mass through the gas phase delays all of the transitions. The effect of the gas phase increases rather dramatically, however, once the pressure is reduced below the atmospheric value, as we show below.

5.2 Convection at reduced pressures

In this section we will discuss the dependence of the critical $Ma$ on the average concentration $c_{a}^{0}$ of air. We will focus exclusively on the 0.65 cSt silicone oil and $Bo_{D}=0.7$, since relevant numerical and experimental data tend to cluster around that value. As we discussed previously, at this $Bo_{D}$ the stationary instability dominates, so a time-independent pattern emerges. The one-sided model of Priede & Gerbeth (1997) does not account for phase change (recall that it corresponds to our enhanced one-sided model with $Bi_{q}=0$), so its predictions are independent of $c_{a}^{0}$. The predictions of the enhanced one-sided and the two-sided model are compared with available data in Fig. 9. The results of the two models are essentially identical, except for low $c_{a}^{0}$, where the validity of the assumptions underlying both models is questionable, as discussed below.

The predicted $Ma_{c}$ increases rather significantly as $c_{a}^{0}$ is decreased from the atmospheric values down to the regime when the gas is dominated by the vapour, rather than air. The results of linear stability analysis are generally consistent with the experimental data, although at atmospheric conditions the predicted thresholds for the onset of PMC and SMC states are higher by about 20% compared with those found in experiment (cf. Fig. 9(a)). Numerical simulations have only been performed in the limits when either air or vapour dominate. When the gas phase is dominated by air ($c_{a}^{0}\geq 0.85$), the predictions of linear stability analysis are consistent with the results of numerical simulations. When the vapour is dominant ($c_{a}^{0}\leq 0.16$), the numerical simulations only find the SUF flow over the limited range of the imposed temperature gradients considered, so no prediction of $Ma_{c}$ can be made. However, the corresponding lower bound supports the theoretical prediction that the threshold increases considerably.

The dependence of the critical wavelength on $c_{a}^{0}$ shows a trend similar to $Ma_{c}$: the predicted $\lambda$ increases as $c_{a}^{0}$ decreases (cf. Fig. 9(b)). These predictions are in reasonable agreement with the available numerical and experimental results (Qin & Grigoriev, 2015; Li et al., 2014). When the gas phase is dominated by air ($c_{a}^{0}\geq 0.85$), the predicted critical wavelength for both the PMC and SMC state is in good agreement with the values found in the numerics and experiment. Only experimental data is available for lower $c_{a}^{0}$, and again it is consistent with theoretical predictions for $c_{a}^{0}>0.3$: the wavelength of the pattern is found to lie between the PMC and SMC curves, as it should. For $c_{a}^{0}<0.16$ the critical wavelength becomes comparable to the length $L$ of the cavity studied in Li et al. (2014), so some discrepancy between the theoretical predictions and the experimental results is expected.

In general, our linear stability analysis assumes a flat interface and will break down for $\lambda\gg d_{l}$ (i.e., at low $c_{a}^{0}$ and/or low $Pr$), when the thickness of the liquid layer is expected to vary noticeably. Furthermore, for sufficiently low $c_{a}^{0}$ and/or large $\Gamma_{x}$ the interfacial temperature profile in the SUF state deviates significantly from linear, resulting in the breakdown of the analytical solution for the base flow (Qin & Grigoriev, 2015). Therefore, the theoretical predictions in these limits are not expected to be particularly accurate.

6 Discussion

As our results illustrate, the predictions of the linear stability analysis based on the two-sided model are in general agreement with the reported numerical and experimental data. In particular, this analysis correctly predicts that the transitions between different flow regimes (SUF, PMC, SMC) are delayed (i.e., shifted towards higher $Ma$) when the concentration of air decreases, although quantitative agreement is not expected at low $c_{a}^{0}$ when the assumptions underlying our analysis become invalid. The general agreement found at higher $c_{a}^{0}$ confirms the validity of the assumptions made in the linear stability analysis, giving us confidence that it captures all the essential physical processes governing the flow stability for volatile liquids driven by a horizontal temperature gradient.

The predictions of linear stability analysis are consistent with the available numerical results. The discrepancy with the experiment can be due to a number of different reasons. One is the effect of transverse confinement: the liquid is wetting and rises at the side walls, increasing the thickness of the liquid layer there by up to 60%. The experimental study by Li et al. (2014) used the (smaller) thickness $d_{l}$ of the liquid layer at the symmetry mid-plane $y=L_{y}/2$ of the cavity to compute $Ma$. Furthermore, the experimental study indirectly deduced the interfacial temperature gradient by comparing the experimental velocity profile with the analytical solution as in Fig. 2(a) rather than measure it directly, which introduces additional uncertainty into the reported values of $Ma_{c}$.

In this problem, buoyancy plays a stabilizing role, hence the instability leading to the formation of convection rolls is driven primarily by thermocapillary stresses, which depend on the interfacial temperature gradient. The fluctuations in the interfacial temperature are significantly affected by the heat and mass transport through the gas phase described – in the enhanced one-sided model – by the Biot coefficient $Bi_{q}$, which is a function of both the wavenumber $q$ and the concentration of air $c_{a}^{0}$. In particular, the first term in the square brackets in (80)

[TABLE]

describes the effect of conductive heat transport through the gas layer and reflects the dependence of thermal conductivity of the gas phase on its composition, while the second term

[TABLE]

where $H$ is the nondimensional parameter defined by (56), describes the effect of the latent heat released or absorbed at the interface as a result of phase change and reflects the variation in the diffusive transport of vapour through the gas phase, which also depends on the gas composition. Both heat and mass transport through the gas phase suppress fluctuations in the interfacial temperature and hence play a stabilizing role.

To see the relative impact of the heat and mass transport on stability, it is instructive to compare the trends and characteristic magnitudes of the two terms. For the 0.65 cSt silicone oil, $B_{1}$ decreases slightly as $c_{a}^{0}$ decreases (thermal conductivity of the vapour is somewhat smaller than that of the air), but $B_{2}$ increases significantly, reflecting the enhancement of phase change (cf. Fig. 10). At atmospheric conditions, the heat flux in the gas phase has a slightly larger effect ($B_{1}=0.23$ vs. $B_{2}=0.17$). The effect of phase change becomes dominant for $c_{a}^{0}<0.95$ and increases as $c_{a}^{0}$ decreases. In fact, $B_{2}$ diverges as $c_{a}^{0}\rightarrow 0$ which, according to (79), implies that the temperature fluctuations at the interface vanish and the critical $Ma$ becomes infinite. The trends are similar for the less volatile 1 cSt silicone oil, although the effect of phase change is clearly weaker.

The increase in the critical $Ma$ with decreasing $c_{a}^{0}$ is due primarily to the enhancement of phase change which increases the amount of latent heat released/absorbed by the warmer/cooler regions of the interface. Interestingly, the composition of the gas phase has a more significant effect on flow stability than the base flow itself. As discussed by Qin & Grigoriev (2015), the concentration of air has a relatively weak effect on the base flow, since the average interfacial temperature gradient, which determines the thermocapillary stresses and hence the speed of the base flow, is insensitive to $c_{a}^{0}$ above 10% or so.

The one-sided model with the adiabatic boundary condition is ill-suited for describing the flow of volatile fluids. As discussed previously, even at atmospheric conditions, when phase change is strongly suppressed, it underestimates the critical $Ma$ corresponding to all three instability types. At reduced pressures, phase change at the interface and heat and mass transport through the gas phase are all significantly enhanced, and predictions of the one-sided model become very inaccurate, since it completely fails to describe these effects. To illustrate the role of the gas phase, it is useful to focus on the dimensionless combinations that have been introduced over the years to describe the deviation from thermodynamic equilibrium, which gives rise to convection in the presence of evaporation. In particular, in their analysis of pure Marangoni (thermocapillary) instability in volatile fluids subjected to a vertical temperature gradient, Burelbach et al. (1988) introduced two nondimensional parameters: the evaporation number

[TABLE]

which defines the ratio of the evaporative time scale (how long it would take for a liquid layer to evaporate completely) to the viscous time scale, and the “nonequilibrium parameter”

[TABLE]

which defines the ratio of the latent heat flux at the interface to the conductive heat flux in the liquid. The stability analysis of that problem was performed correctly only recently by Chauvet et al. (2012) who found that a version of the Biot coefficient (80) appears there as well. However, neither $E$ nor $K$ appears in the stability analysis, regardless of the direction of the temperature gradient, because both parameters fail to account for the transport (of either heat or mass) in the gas phase. Instead, the Biot coefficient involves a nondimensional combination $H$ which describes the effect of the latent heat associated with phase change and explicitly takes into account the mass transport in the gas phase.

The study by Normand et al. (1977) was probably the first to note that the Biot “number” should depend on the wavenumber of the pattern for convective flows. In their analysis of buoyancy-thermocapillary convection driven by a horizontal temperature gradient, Mercier & Normand (1996) repeated that statement, but proposed to use the definition

[TABLE]

based on conduction in the two layers instead. Although wavenumber-independent, this definition is similar in form to the correct one in the combined limits of an infinitely thick gas layer and vanishing volatility. Indeed, if we ignore phase change, for $d_{g}\gg d_{l}$, the expression (80) would reduce to

[TABLE]

Here the wavelength of the pattern effectively plays the role of the gas layer thickness, the result that follows immediately from the second equation in (68) which controls the heat transport in the gas phase.

The wavenumber dependence of the Biot coefficient (80) implies that, once transformed to the real space, Newton’s cooling law

[TABLE]

becomes invalid and has to be replaced by a spatially nonlocal relation between the interfacial temperature and the normal temperature gradient (or heat flux) in the liquid phase. Given that the Biot number (101) describes one-dimensional conduction, it is not surprising that (103) fails when spatial variation of the solution in the extended direction(s) is taken into account.

Finally, by comparing the results of linear stability analysis based on three models of the gas layer that incorporate varying levels of detail, we found that advection (of heat, mass, or momentum) plays a relatively minor role and can be ignored except for low $c_{a}^{0}$ (e.g., when vapor dominates in the gas phase). In this limit, the two-sided model, which accounts for advection in the gas layer predicts a slightly higher critical $Ma$ for the threshold of both PMC and SMC. The increased role of advection can be readily understood by considering the associated changes in the Péclet numbers (66) and (67). Both $\alpha_{g}$ and $D$ are proportional to $p_{g}^{-1}\propto 1-c_{a}^{0}$, so both $Pe_{m}$ and $Pe_{t}$ should increase with decreasing $c_{a}^{0}$ (even if we ignore the increase in the characteristic velocity of the gas associated with enhanced phase change). We should also expect advection to play a more important role for low $Pr$, since both $Pe_{m}$ and $Pe_{t}$ are proportional to $Re=Ma/Pr$. The enhanced advective transport in the gas phase at low $c_{a}^{0}$ or low $Pr$ tends to suppress the variation in the temperature at the interface, the thermocapillary stresses, and hence the onset of convection.

7 Conclusions

As several recent experimental (Li et al., 2014) and numerical (Qin & Grigoriev, 2015) studies have demonstrated, phase change has a rather dramatic effect on the stability of the flow in a layer of a volatile liquid (0.65 cSt silicone oil) driven by a horizontal temperature gradient. These studies showed that the onset of instability is delayed rather significantly when the phase change is enhanced by removing noncondensable gases, such as air, from the system. Linear stability analysis based on a model that takes heat and mass transport in the gas phase into account produces results that are overall in good quantitative agreement with these studies. In particular, it confirms that the instability that is responsible for generating a pattern of convection rolls is caused by the thermocapillary stresses, while buoyancy in the liquid layer has a stabilizing effect.

The analysis also shows that phase change plays a very important role, with the magnitude of the effect described by the nondimensional combination $H$ defined by (56), the same combination that appears in the convection problem with a vertical, rather than horizontal, temperature gradient (Chauvet et al., 2012). The similarity of the two problems is not coincidental: due to the weakness of advective fluxes, the transport of heat and mass in the gas layer is primarily due to diffusion of the perturbations about the base state. This diffusive transport is independent of the direction of thermal gradient, which mainly affects the base state.

Even when phase change is suppressed (e.g., at atmospheric conditions or for a fluid with relatively low volatility), we find that the gas phase can have a noticeable effect on the stability of the flow. Our analysis shows that the transport of heat in the gas layer can account for an increase of more than ten percent in the critical Marangoni number. The effect of heat transport through the gas layer has been traditionally described using Newton cooling law, with the magnitude of the heat flux characterized by the Biot number. While this law may be a reasonable approximation for the base flow, it breaks down completely the moment convection rolls appear. In fact, heat transport in the gas phase has to be treated explicitly in order to obtain quantitatively accurate results in the entire range of parameters.

Acknowledgements

This work was supported in part by the National Science Foundation under Grant No. CMMI-1511470.

Appendix A Numerical Solution of the Boundary Value Problem

The boundary value solver bvp5c is designed to deal with systems of first-order ODEs, so the higher-order ODEs arising in the linear stability analysis were converted to this form.

A.1 Diffusion-Dominated Case

For the diffusion-dominated case, we converted the boundary value problem (81) to a system of six first-order ODEs

[TABLE]

where $y_{1}=\tilde{\psi}_{lq}$, $y_{5}=\tilde{\theta}_{lq}$, and

[TABLE]

Respectively, the boundary conditions become

[TABLE]

and

[TABLE]

A.2 Transient Dynamics in the Gas Layer

To avoid solving coupled boundary value problems (81) and (89) defined on adjacent domains, we remapped the domain of $\tilde{\theta}_{gq}$ and $\tilde{\varsigma}_{vq}$ from $[0,A]$ to $[-1,0]$ by defining $\bar{z}=-\tilde{z}/A$ and introducing two new functions $\bar{\theta}_{gq}(\bar{z})=\tilde{\theta}_{gq}(-\tilde{z})$ and $\bar{\varsigma}_{vq}(\bar{z})=\tilde{\varsigma}_{vq}(-\tilde{z})$. In terms of these new functions, the system (87) becomes

[TABLE]

where the prime now denotes the derivative with respect to $\bar{z}$, and the boundary condition (88) is replaced with

[TABLE]

Furthermore, the no-flux boundary conditions at the top of the gas layer become

[TABLE]

Temperature continuity and the local phase equilibrium (76) at the free surface yield

[TABLE]

Converting (108) into a system of first-order ODEs yields four additional equations

[TABLE]

where $y_{7}=\bar{\theta}_{gq}$, and $y_{9}=\bar{\varsigma}_{vq}$, which should be added to the system (104). The boundary conditions (106) remain and (107) is replaced with

[TABLE]

A.3 The Effect of Advection in the Gas Layer

Just like in the previous case, we introduce new functions $\bar{\theta}_{gq}(\bar{z})=\tilde{\theta}_{gq}(-\tilde{z})$, $\bar{\varsigma}_{vq}(\bar{z})=\tilde{\varsigma}_{vq}(-\tilde{z})$, and $\bar{\psi}_{gq}(\bar{z})=\tilde{\psi}_{gq}(-\tilde{z})$, so the system (89) becomes

[TABLE]

where we have defined $\bar{\nabla}_{q}^{2}=A^{-2}\partial_{\bar{z}}^{2}-q^{2}$ and

[TABLE]

The functions $\bar{\varsigma}_{vq}(\bar{z})$ and $\bar{\theta}_{gq}(\bar{z})$ again satisfy the boundary conditions (110) and (111). The boundary conditions for the function $\bar{\psi}_{gq}(\bar{z})$ follow from (46) and (47):

[TABLE]

Converting (114) into a system of first-order ODEs yields eight additional equations

[TABLE]

which augment the system (104) describing the liquid layer. Here $y_{7}=\bar{\theta}_{gq}$, $y_{9}=\bar{\varsigma}_{vq}$, $y_{11}=\bar{\psi}_{gq}$, and

[TABLE]

The boundary conditions are given by (106), (113), and

[TABLE]

In all of the expressions given in the Appendix, the nondimensional parameters $Re$, $Gr$, $\mathcal{R}$, $\mathcal{B}$, $\Xi_{T}$, and $\Xi_{\varsigma}$ were expressed in terms of $Ma$, $Pr$, $Bo_{D}$, and $c_{a}^{0}$ using the definitions (28), (30), (3), (40), (55), (3.3), and (65). Similarly, parameters $K_{D}$, $K_{\alpha}$ and $K_{\nu}$ were evaluated as functions of $c_{a}^{0}$.