Instability of axisymmetric flow in thermocapillary liquid bridges: Kinetic and thermal energy budgets for two-phase flow with temperature-dependent material properties

Abstract

In numerical linear stability investigations, the rates of change of the kinetic and thermal energy of the perturbation flow are often used to identify the dominant mechanisms by which kinetic or thermal energy is exchanged between the basic and the perturbation flow. Extending the conventional energy analysis for a single-phase Boussinesq fluid, the energy budgets of arbitrary infinitesimal perturbations to the basic two-phase liquid–gas flow are derived for an axisymmetric thermocapillary bridge when the material parameters in both phases depend on the temperature. This allows identifying individual transport terms and assessing their contributions to the instability if the basic flow and the critical mode are evaluated at criticality. The full closed-form energy budgets of linear modes have been derived for thermocapillary two-phase flow taking into account the temperature dependence of all thermophysical parameters. The influence of different approximations to the temperature dependence on the linear stability boundary of the axisymmetric flow in thermocapillary liquid bridges is tested regarding their accuracy. The general mechanism of symmetry breaking turns out to be very robust.

1. Introduction

Thermocapillary flow in axisymmetric liquid bridges which are heated differentially represents one of the most popular paradigms of thermocapillary flow [8]. It originated from the desire to better understand the formation of striations in crystals grown by the floating-zone method [18]. One intriguing aspect is the spontaneous breaking of the steady axisymmetric flow and the relevant physical mechanisms at work. The onset of the three-dimensional flow in high-Prandtl-number liquids is characterised by a thermocapillary Reynolds number ${\textit{Re}} \sim \Delta T d/\mu ^2$ which scales linearly with the length of the liquid bridge $d$ and the temperature difference $\Delta T$ applied between the supporting rods or, equivalently, with the total variation of the surface tension. Since the length $d$ under normal gravity is limited by the Rayleigh–Plateau instability causing a mechanical breaking of the bridge, driving the system into a three-dimensional flow state can require a relatively large temperature difference $\Delta T$ , in particular if the dynamic viscosity $\mu$ is large. Hence, the temperature dependence of the material properties, like the dynamic viscosity, may not be negligible. For that reason, a temperature-dependent viscosity has been taken into account in stability analyses [7] and numerical simulations [22, 23]. However, most numerical results have been obtained for constant material properties [4, 9, 12, 14, 15, 16, 29].

Since the work of Reynolds [19], Orr [17] and, for thermocapillary flows, Smith [25], the balance of kinetic and thermal perturbation energies of the (infinitesimal) perturbation flow has proven a valuable tool to determine the total amount and the spatial distribution of the energy production and dissipation. Knowledge of these properties can be key to understanding the physical mechanisms of linear instability of the flow. For instance, the thermocapillary instability for low-Prandtl-number liquids is caused by the lift-up mechanism in the free surface shear layer [10, 33], similar to the vortex-ring instability [34], while the temperature field is passive with respect to the instability. The temperature field is only required to drive the basic flow. For large Prandtl numbers, on the other hand, the temperature field is important and the flow instability arises in form of a pair of azimuthally propagating hydrothermal waves [33], similar as for plane layers [26]. The hydrothermal waves draw their energy from temperature gradients of the basic flow in the bulk. The strong internal temperature gradients arise due to the basic recirculation, driven by thermocapillary forces, which transports hot fluid from the free surface over the cold wall deep into the bulk.

Figure 1. Sketch of the thermocapillary liquid bridge held in place between the hot rod at temperature $\bar T + \Delta T/2$ and the cold rod at temperature $\bar T -\Delta T/2$ . The flow is driven by (i) the thermocapillary effect, (ii) buoyancy forces in the gravity field $\boldsymbol{{g}}$ and (iii) a gas flow with a given inlet velocity $w_{\text{G,in}}$ . $A_{\text{fs}}$ and $A_{\text{out}}$ denote the liquid–gas interface and the outlet section, respectively. Polar coordinates are indicated.

Uncertainties in the computation of the critical Reynolds number are mainly caused by (a) discretisation errors, (b) neglect of the temperature dependence of the material properties and (c) simplifying assumptions about the ambient conditions like the assumption of Newtons’s law of heat transfer. As numerical capabilities have improved, the discretisation errors can be well controlled. Moreover, it has become feasible to take into account temperature-dependent material properties as well as the flow in the ambient atmosphere. Due to its usefulness for the understanding of the physical instability mechanics as well as for a check of the energy preservation of the numerically computed critical mode, we establish the Reynold–Orr equations governing the temporal evolution of the kinetic and thermal energy budgets of the critical perturbation mode of the linear theory for an axisymmetric liquid bridge surrounded by a gas. The energy balances obtained take into account the temperature dependence of all thermophysical material properties and are valid both in the liquid and in the gas phase, of which the latter is confined to a concentric tube surrounding the liquid bridge.

2. Geometrical configuration

The flow in a liquid droplet suspended between two coaxial cylindrical rods is considered, assuming that the physical properties of the liquid and gas phases are temperature-dependent (Figure 1). The two rods are kept at constant temperatures with $T_{\text{hot}} = \bar T + \Delta T/2$ (hot rod) and $T_{\text{cold}} = \bar T-\Delta T/2$ (cold rod), where $\bar T=(T_{\text{hot}}+T_{\text{cold}})/2$ denotes the arithmetic mean temperature which is assumed as reference temperature. The liquid bridge is surrounded by a gas confined to a coaxial cylindrical tube, intended to prevent uncontrollable circulations from the ambience in experiments. The flow is driven (i) by thermocapillary forces acting on the interface due to a variable surface tension $\sigma (T)$ , (ii) buoyancy forces in the presence of the acceleration of gravity which is assumed to be directed in the negative axial direction (Figure 1) and (iii) by an externally imposed gas flow ( $w_{\text{G,in}}$ in Figure 1), assumed to be axisymmetric and non-swirling. The external gas flow affects the flow in the liquid phase via viscous shear stresses acting on the liquid–gas interface and by the thermal coupling between liquid and gas.

On all solid walls of the support rods and the tube, we assume no-slip boundary conditions. On the surfaces of the rods, the temperatures $T_{\text{hot}}$ and $T_{\text{cold}}$ are imposed as indicated by colour in Figure 1, while the cylindrical tube is assumed adiabatic. Gas may enter the system with a given axial velocity profile $w_{\text{G,in}}(r)$ and a given temperature, also satisfying outflow conditions at the outlet (denoted $A_{\text{out}}$ in Figure 1). Alternatively, the gas tube may be closed ( $w_{\text{G,in}}\equiv w_{\text{G,out}}\equiv 0$ ) and confined by either adiabatic or conductive walls. The thermocapillary flow is driven along the interface $A_{\text{fs}}$ by the thermocapillary effect which creates the tangential stress [11]

(2.1) \begin{align} \nabla _\| \sigma (T) &= \frac{\partial \sigma }{\partial T} \nabla _\| T = \frac{\partial }{\partial T}\left [\sigma (\bar T) -\gamma (T - \bar T) + \delta (T - \bar T)^2 + \ldots \right ]\nabla _\| T \nonumber \\ &= \left [-\gamma + 2 \delta (T - \bar T) + \ldots \right ] \nabla _\| T, \end{align}

where $\nabla _\|$ denotes the tangential Nabla operator. For the overwhelming majority of liquids, the surface tension decreases with temperature $(\gamma \gt 0)$ . Therefore, the thermocapillary effect will typically generate a flow which is directed along the interface away from the hot rod (low surface tension) and towards the cold rod (high surface tension). The usual approximation is to neglect quadratic terms in the Taylor expansion of the surface tension leading to $\nabla _\| \sigma (T) \approx -\gamma \nabla _\| T$ . We note that the present analysis is independent of the exact functional dependence of $\sigma (T)$ . While the Taylor coefficients $\gamma$ , $\delta$ , etc. in (2.1) crucially affect the flow fields by coupling the velocity and temperature fields on the interface, they do not explicitly appear in the energy budgets of a given linear instability mode.

3. Governing equations

The general flow problem is governed by the Navier–Stokes and energy equations in both the liquid and the gas phase. We consider the strong conservative forms

(3.1a) \begin{align} \partial _t\rho +\nabla \cdot \big (\rho \hat{\boldsymbol{{U}}}\big )&=0, \end{align}

(3.1b) \begin{align} \partial _t\big (\rho \hat{\boldsymbol{{U}}}\big ) + \nabla \cdot \big (\rho \hat{\boldsymbol{{U}}} \hat{\boldsymbol{{U}}}\big ) &= -\nabla \hat P +\rho \boldsymbol{{g}}+ \nabla \cdot \boldsymbol{\mathcal{T}}, \end{align}

(3.1c) \begin{align} \partial _t \big (\rho c_{p} \hat T\big )+\nabla \cdot \big (\rho c_{p} \hat T \hat{\boldsymbol{{U}}}\big )&=\nabla \cdot \big (\lambda \nabla \hat T\big ), \end{align}

where $\hat{\boldsymbol{{U}}}$ , $\hat P$ and $\hat T$ are the velocity, pressure and temperature fields. Henceforth, the hat ( $\hat{}$ ) indicates the total flow fields, including a basic flow and a perturbation. Equations (3.1a)–(3.1c) describe the transport in both the liquid and the gas phase. As long as the formulation for both phases is the same, we do not distinguish between them. We assume the fluids are Newtonian with stress tensor

(3.2) \begin{equation} \boldsymbol{\mathcal{T}} = \mu \left [\nabla \hat{\boldsymbol{{U}}} + \big(\nabla \hat{\boldsymbol{{U}}}\big)^{\text{T}} \right ]-\frac{2}{3}\mu \big(\nabla \cdot \hat{\boldsymbol{{U}}} \big)\boldsymbol{{I}}, \end{equation}

where $\boldsymbol{{I}}$ is the identity matrix. The temperature-dependent density, dynamic viscosity, thermal conductivity and specific heat at constant pressure are denoted $\rho (\hat T)$ , $\mu (\hat T)$ , $\lambda (\hat T)$ and $c_p(\hat T)$ , respectively. Since the density is treated as temperature-dependent, the velocity field is not solenoidal. The formulation used for the temperature equation (3.1c) neglects the pressure work and the viscous dissipation. These assumptions are justified, respectively, if the conditions

(3.3) \begin{equation} \chi \frac{\bar{T}}{\Delta T}\leq 0.1\quad \mathrm{and}\quad \chi{\textit{Pr}}\leq 0.1\quad \mathrm{with}\quad \chi =\frac{\bar{\beta } g d}{\bar{c}_p} \end{equation}

are satisfied, where ${\textit{Pr}}=\bar{\mu }\bar{c}_p/\bar{\lambda }$ is the Prandtl number and $\beta =-\rho ^{-1}({\partial } \rho/{\partial } T)_p$ is the thermal expansion coefficient [3]. The overbar indicates reference values at the reference temperature $\bar{T}$ . Similarly, we disregard the pressure contribution to the enthalpy in (3.1c) assuming $p/\rho \ll \vert c_p T\vert$ . In Section 6, we verify the conditions (3.3) for two different cases, confirming the validity of (3.1c).

For equations (3.1a)–(3.1c) and for the assumed steady axisymmetric boundary conditions a steady axisymmetric basic flow $(\boldsymbol{{u}}_0,T_0,p_0,h_0)$ exists. The axisymmetric shape function $h_0 (z)$ marks the radial coordinate of the location of the liquid–gas interphase in the axisymmetric steady flow which is assumed to be pinned to the sharp circular edges of the supporting rods. The shape $h_0(z)$ is determined by the flow-induced normal stresses and the Laplace pressure, which also depends on the full surface tension $\sigma (\hat T)$ and the hydrostatic pressure difference.

Here we are not concerned with computing the basic flow. We assume it has been obtained numerically, taking into account the full temperature dependence of the material properties. Therefore, the exact form of the boundary conditions and forcing terms for the basic flow does not enter the present problem. Furthermore, we assume a linear stability analysis has been carried out by solving the associated eigenvalue problem (see e.g. Stojanović et al. [27]) such that the neutrally stable linear mode $(\boldsymbol{{u}},T,p,h)$ is available as well. We consider the formal decomposition

(3.4a) \begin{align} \hat{\boldsymbol{{U}}}&=\boldsymbol{{u}}_0(r,z)+\boldsymbol{{u}}(r,\varphi,z,t), \end{align}

(3.4b) \begin{align} \hat T&=T_0(r,z)+T(r,\varphi,z,t), \end{align}

(3.4c) \begin{align} \hat P&=p_0(r,z)+p(r,\varphi,z,t), \end{align}

(3.4d) \begin{align} \hat H&=h_0(z)+h(\varphi,z,t), \end{align}

of the total flow ( $\,\hat{ }\,$ ) into the basic state (index 0) and a perturbation $(\boldsymbol{{u}},T,p,h)$ . All flow fields are described using cylindrical coordinates $(r,\varphi,z)$ and associated velocity components $(u,v,w)$ such that $\boldsymbol{{u}} = u {\boldsymbol{{e}}_\boldsymbol{{r}}} + v\boldsymbol{{e}}_\varphi + w{\boldsymbol{{e}}_\boldsymbol{{z}}}$ , where $\boldsymbol{{e}}$ denotes a unit vector.

The neutral mode is typically obtained by a linearisation of the governing equations which requires the perturbation quantities to be asymptotically small. We do not explicitly introduce a smallness parameter $\epsilon$ , but keep in mind that the perturbation quantities $(\boldsymbol{{u}},T,p,h)$ are all of the order of $O(\epsilon )$ in the sense of the linearisation required for the linear stability analysis. For convenience, we shall not express the perturbation quantities by normal modes. This is easily accomplished a posteriori.

In order to keep the effort required in deriving the energy budgets for the neutral mode at a meaningful level, we make the ad hoc assumption that the perturbation flow does not affect the interfacial shape. This is motivated by the experimental observations that the interfacial deformations due to the supercritical three-dimensional flow are very small under typical laboratory conditions [36]. We note, however, that this simplification precludes surface waves from the range of critical modes. With the assumption $h=0$ , the outward pointing unit vector normal to the free surface is

(3.5) \begin{equation} \boldsymbol{{n}} = \frac{1}{N}\big(\boldsymbol{{e}}_r-h_{0z}\boldsymbol{{e}}_z\big), \end{equation}

with the normalising denominator $N=\sqrt{1+h_{0z}^2}$ , where we use the notation $h_{0z} = \text{d} h_0/\text{d}z$ . Likewise, we define $h_{0zz} = \text{d}^2 h_0/\text{d}z^2$ .

Let us assume the basic flow has been computed from the axisymmetric steady version of (3.1a)–(3.1c). Furthermore, we assume the perturbation flow has been computed by solving (3.1a)–(3.1c) linearised about the basic state. Within a postprocessing step, we are then interested in the temporal evolution of the kinetic and thermal perturbation energy densities

(3.6a) \begin{align} \varepsilon _{\text{kin}}({\boldsymbol{{x}},\boldsymbol{{t}}}) &= \frac{1}{2} \rho (T_0) \boldsymbol{{u}}^2, \end{align}

(3.6b) \begin{align} \varepsilon _{\text{therm}}({\boldsymbol{{x}},\boldsymbol{{t}}}) &= \frac{1}{2} \rho (T_0) c_{p}(T_0) T^2, \end{align}

in the liquid and the gas phase. These energy densities must be considered measures of the perturbation flow in the sense of Joseph [5]. With the perturbations being of $O(\epsilon )$ , the above energy densities are of $O(\epsilon ^2)$ . The aim is to express $\partial _t \varepsilon _{\text{kin}}({\boldsymbol{{x}},\boldsymbol{{t}}})$ and $\partial _t \varepsilon _{\mathrm{therm}}({\boldsymbol{{x}},\boldsymbol{{t}}})$ by a sum of contributions which describe individual transport processes and can be interpreted in physical terms. While the local rates of change of the energy densities (3.6a) and (3.6b) are generally non-zero and depend on $\boldsymbol{{x}}$ , the total change rates obtained by integration over the volume occupied by the respective fluid must vanish, if the perturbation flow field represents a critical or a neutral mode for which the growth rate vanishes.

In the following, we assume that the fluid properties depend solely on the temperature and not on the pressure. This simplifying assumption is commonly made by the manufacturers of batch liquids employed for silicone oil liquid bridges [24] and for liquids in general far from their phase-change critical points. To take into account the temperature dependence of the material parameters, we assume the parameters, as well as their first and second derivatives ${\rho^{\prime}} (\hat T)$ , ${\mu^{\prime}} (\hat T)$ , ${\lambda^{\prime}} (\hat T)$ , ${c^{\prime}_p} (\hat T)$ and ${\rho^{\prime\prime}} (\hat T)$ , ${\mu^{\prime\prime}} (\hat T)$ , ${\lambda^{\prime\prime}} (\hat T)$ , ${c^{\prime\prime}_p} (\hat T)$ , respectively, are available in closed-form expressions as functions of the temperature $\hat T$ . The functional dependence on $\hat T$ could be established, for instance, by fitting discrete data by suitable ansatz functions (polynomials, exponentials, etc.) or spline functions. With this information available, all material parameters can be expanded about the local temperature $T_0(r,z)$ of the basic flow and up to second order

(3.7a) \begin{align} \rho \big(\hat T\big) &= \rho (T_0+T)\approx \rho (T_0)+\left .{\partial }_T\rho \right |_{T_0}T +\frac{1}{2}\left .{\partial }^2_T\rho \right |_{T_0}T^2 \,:\!=\rho _{0}+\rho^{\prime}_{0}T+\frac{1}{2}\rho^{\prime\prime}_{0}T^2, \end{align}

(3.7b) \begin{align} \mu \big(\hat T\big) &= \mu (T_0+T)\approx \mu (T_0)+\left .{\partial }_T\mu \right |_{T_0}T+\frac{1}{2}\left .{\partial }^2_T\mu \right |_{T_0}T^2\,:\!=\mu _{0}+\mu^{\prime}_{0}T+\frac{1}{2}\mu^{\prime\prime}_{0}T^2, \end{align}

(3.7c) \begin{align} \lambda \big(\hat T\big) &= \lambda (T_0+T)\approx \lambda (T_0)+\left .{\partial }_T\lambda \right |_{T_0}T+\frac{1}{2}\left .{\partial }^2_T\lambda \right |_{T_0}T^2\,:\!=\lambda _{0}+\lambda^{\prime}_{0}T+\frac{1}{2}\lambda^{\prime\prime}_{0}T^2, \end{align}

(3.7d) \begin{align} c_p\big(\hat T\big) &= c_p(T_0+T)\approx c_p(T_0)+\left .{\partial }_T c_p\right |_{T_0}T +\frac{1}{2}\left .{\partial }^2_T c_p\right |_{T_0}T^2\,:\!=c_{p0} + c^{\prime}_{p0}T + \frac{1}{2}c^{\prime\prime}_{p0}T^2. \end{align}

The Taylor coefficients $(\rho _{0},\mu _{0},\lambda _{0},c_{p0})$ , $(\rho^{\prime}_{0},\mu^{\prime}_{0},\lambda^{\prime}_{0},c^{\prime}_{p0})$ and $(\rho^{\prime\prime}_{0},\mu^{\prime\prime}_{0},\lambda^{\prime\prime}_{0},c^{\prime\prime}_{p0})$ are scalar fields which depend continuously on the basic temperature $T_0$ . Note all the above thermophysical properties and coefficients depend on the phase.

Before deriving the rates of change of the kinetic energies, it is useful to consider the continuity equation. Inserting expansions (3.7a)–(3.7d) into the continuity equation (3.1a) and neglecting quadratic terms yields

(3.8) \begin{equation} \partial _t(\rho^{\prime}_0 T)+\nabla \cdot (\rho _0 \boldsymbol{{u}}_0)+\nabla \cdot (\rho _0 \boldsymbol{{u}})+\nabla \cdot (\rho^{\prime}_0 T \boldsymbol{{u}}_0) =0. \end{equation}

As this equation involves different orders of magnitude, the terms of each order of magnitude in this equation must vanish separately,

(3.9a) \begin{align} O\big(\epsilon ^0\big):&& \nabla \cdot (\rho _0 \boldsymbol{{u}}_0)&=0, \end{align}

(3.9b) \begin{align} O\big(\epsilon ^1\big):&&\rho^{\prime}_0 \partial _t T+\nabla \cdot (\rho _0 \boldsymbol{{u}})+\nabla \cdot (\rho^{\prime}_0 T \boldsymbol{{u}}_0)&=0. \end{align}

The terms of order $O(\epsilon ^0)$ arise in the equations for the basic flow and thus do not enter the equations of order $O(\epsilon ^1)$ for the perturbation flow. Equation (3.9b), on the other hand, balances the terms of $O(\epsilon ^1)$ and thus represents the continuity equation in linear order entering the linear stability analysis.

4. Thermal energy budget

The basic flow $(\boldsymbol{{u}}_0,T_0,p_0,h_0)$ is assumed stationary and of order $O(\epsilon ^0)$ . Since the perturbation flow is of order $O(\epsilon ^1)$ , it has no effect on the energy budget of the basic flow which is also $O(\epsilon ^0)$ . Nevertheless, the rates of change of the perturbation energies (3.6a) and (3.6b), which are of order $O(\epsilon ^2)$ , contain terms which describe an energy exchange with the basic flow.

To derive the local thermal energy budget, we first derive the linear equation governing the evolution of the perturbation temperature. To that end, the flow decomposition (3.4a)–(3.4d) is inserted into the temperature equation (3.1c) to obtain

(4.1) \begin{equation} \partial _t(\rho c_p T_0)+\partial _t(\rho c_p T)+\nabla \cdot \left [\rho c_p (T_0+T)(\boldsymbol{{u}}_0+\boldsymbol{{u}})\right ] = \nabla \cdot (\lambda \nabla T_0)+\nabla \cdot (\lambda \nabla T), \end{equation}

which contains both the basic state and the perturbation flow. Neglecting terms of order $O(\epsilon ^2)$ yields

(4.2) \begin{equation} \partial _t(\rho c_p T_0)+\partial _t(\rho c_p T)+\nabla \cdot (\rho c_p T_0\boldsymbol{{u}}_0)+\nabla \cdot (\rho c_p T_0\boldsymbol{{u}})+\nabla \cdot (\rho c_p T\boldsymbol{{u}}_0)=\nabla \cdot (\lambda \nabla T_0)+\nabla \cdot (\lambda \nabla T). \end{equation}

Inserting the Taylor expansions (3.7a)–(3.7d) of $\rho$ , $c_p$ and $\lambda$ we obtain, after linearisation,

(4.3) \begin{align} &T_0\left .{\partial }_T(\rho c_p)\right |_{T_0}\partial _t T+\rho _0 c_{p0}\partial _t T+\nabla \cdot (\rho _0 c_{p0} T_0\boldsymbol{{u}}_0)+\nabla \cdot (\rho^{\prime}_0 c_{p0} T_0\boldsymbol{{u}}_0 T)+\nabla \cdot (\rho _0 c^{\prime}_{p0} T_0\boldsymbol{{u}}_0 T)\\ &\quad +\nabla \cdot (\rho _0 c_{p0} T\boldsymbol{{u}}_0)+\nabla \cdot (\rho _0 c_{p0} T_0\boldsymbol{{u}})=\nabla \cdot (\lambda _0 \nabla T_0)+\nabla \cdot (\lambda^{\prime}_0 T \nabla T_0)+\nabla \cdot (\lambda _0 \nabla T), \nonumber\end{align}

where the coefficients, like $\lambda _0=\lambda (T_0)$ , are functions of the basic state temperature field $T_0$ . Separating the orders of magnitude, the terms of order $O(\epsilon ^0)$ enter the basic state equation for $T_0$

(4.4) \begin{equation} \nabla \cdot (\rho _0 c_{p0} T_0\boldsymbol{{u}}_0)=\nabla \cdot (\lambda _0 \nabla T_0), \end{equation}

while the terms of order $O(\epsilon ^1)$

(4.5) \begin{align} & \left [T_0 (\rho^{\prime}_0 c_{p0}+\rho _0 c^{\prime}_{p0})+\rho _0 c_{p0}\right ]\partial _t T+\nabla \cdot (\rho^{\prime}_0 c_{p0} T_0\boldsymbol{{u}}_0 T)+\nabla \cdot (\rho _0 c^{\prime}_{p0} T_0\boldsymbol{{u}}_0 T)\nonumber \\ &\quad +\nabla \cdot (\rho _0 c_{p0} T\boldsymbol{{u}}_0)+\nabla \cdot (\rho _0 c_{p0} T_0\boldsymbol{{u}})=\nabla \cdot (\lambda^{\prime}_0 T \nabla T_0)+\nabla \cdot (\lambda _0 \nabla T) \end{align}

constitute the linear perturbation equation for $T$ .

To obtain the rate of change of the thermal energy density (3.6b), equation (4.5) is multiplied by $T$ to yield

(4.6) \begin{align}{\underbrace{\vphantom{\bigl [} \rho _0 c_{p0} T\partial _t T}_{\texttt{T1}}}&= -{\underbrace{\vphantom{\bigl [} (\rho^{\prime}_0 c_{p0}+\rho _0 c^{\prime}_{p0}) T_0 T \partial _t T}_{\texttt{T2}}} -{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\rho^{\prime}_0 c_{p0} T_0\boldsymbol{{u}}_0 T)}_{\texttt{T3}}} -{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\rho _0 c^{\prime}_{p0} T_0\boldsymbol{{u}}_0 T)}_{\texttt{ T4}}}\nonumber \\ &\quad -{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\rho _0 c_{p0} T\boldsymbol{{u}}_0)}_{\texttt{T5}}} -{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\rho _0 c_{p0} T_0\boldsymbol{{u}})}_{\texttt{T6}}} +{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\lambda^{\prime}_0 T \nabla T_0)}_{\texttt{T7}}} +{\underbrace{\vphantom{\bigl [} T\nabla \cdot (\lambda _0 \nabla T)}_{\texttt{T8}}}. \end{align}

As far as the transport mechanisms are concerned, we recognise that the term T1 represents the rate of change of thermal perturbation energy density $\partial _t\varepsilon _{\text{therm}}$ . Moreover, the term $\texttt{T2}$ describes an additional rate of change of thermal perturbation energy density due to the dependence of $\rho$ and $c_p$ on the temperature $T_0$ .

The remaining divergence terms describe the rates of change of thermal perturbation energy density due to the divergence of thermal perturbation energy flux densities. These thermal perturbation energy flux densities are caused by the basic state velocity $\boldsymbol{{u}}_0$ and the thermal energy densities given by $(\rho^{\prime}_0 T_0)c_{p0} T$ in $\texttt{T3}$ , $(c^{\prime}_{p0} T_0)\rho _0 T$ in $\texttt{T4}$ and $\rho _0 c_{p0} T$ in $\texttt{T5}$ , where the first and second terms are due to the variation with $T_0$ of $\rho$ and $c_p$ , respectively. The term $\texttt{T6}$ is due to the thermal perturbation energy flux density which is caused by the transport of basic state thermal energy $\rho _0 c_{p0} T_0$ by the perturbation velocity $\boldsymbol{{u}}$ .

Finally, the term $\texttt{T8}$ describes the dissipation of thermal energy due to Fourier’s diffusive perturbation heat flux $-\lambda _0 \nabla T$ , and $\texttt{T7}$ is due to the diffusive heat flux caused by gradients of the basic temperature in combination with the temperature dependence of $\lambda$ which is not taken care of in the $O(\epsilon ^0)$ equations for the basic state.

To arrive at the total, that is integral, thermal energy budget for each fluid phase, the rate of change of the thermal energy density (4.6) must be integrated over the volume $V_i$ occupied by the respective phase (liquid or gas), where the subscript $i\in [\text{L, G}]$ indicates the phase. Moreover, we define the coefficient

(4.7) \begin{equation} \alpha _i = \begin{cases} 1 & i = \text{L},\\ -1 & i = \text{G}. \end{cases} \end{equation}

Since the integration is rather technical, the derivation of the integral thermal energy budget is made in Appendix A. As a result, we obtain the total rate of change of thermal energy $\partial _t E_T = \partial _t \int _{V_i} \varepsilon _{\text{therm}}\,\text{d} V$ in the phase $i$

(4.8) \begin{equation} {\partial }_t E_T= -D_{\text{th}}+J+H_{\text{fs}}+K_{\text{G,th}}+\Pi _\rho +\Pi _{c_p}+\Pi _\lambda -{\partial }_t E^{\prime}_T, \end{equation}

with the abbreviations

(4.9a) \begin{align} D_{\text{th}} &=\int _{V_i} \lambda _0 (\nabla T)^2\,\text{d} V, \end{align}

(4.9b) \begin{align} J=\sum _{j=1}^2 J_j &=-\int _{V_i} \rho _0 T \left [ u{\partial }_r (c_{p0} T_0)+w{\partial }_z (c_{p0} T_0) \right ]\,\text{d} V, \end{align}

(4.9c) \begin{align} H_{\text{fs}} &=\alpha _i\int _{A_{\text{fs}}} \lambda _0 T \nabla T \cdot \boldsymbol{{n}}\,\text{d} S, \end{align}

(4.9d) \begin{align} K_{\text{G,th}} &=-\frac{1-\alpha _i}{4}\int _{A_{\text{out}}}\rho _0 c_{p0} T^2w_0 \,\text{d} S, \end{align}

(4.9e) \begin{align} \Pi _\rho &=\int _{V_{i}} \rho^{\prime}_{0} c_{p0} T_{0} \boldsymbol{{u}}_{0}\cdot \nabla{T}^{2} \,\text{d} V-{\frac{1}{2}}\int _{V_{i}} \left (\frac{\rho^{\prime 2}_{0}}{\rho _{0}}-\rho^{\prime\prime}_{0}\right ){c}_{p0} T^{2} \boldsymbol{{u}}_{0}\cdot \nabla T_0^{2}\,\text{d} V\nonumber \\ &\quad -\frac{1-\alpha _i}{2}\int _{{A}_{\text{out}}}\rho^{\prime}_{0} c_{p0} T_{0} T^{2} w_{0} \,\text{d} S, \end{align}

(4.9f) \begin{align} \Pi _{c_p}&=\frac{1}{2}\int _{V_i} \rho _0 c^{\prime}_{p0} T_0\boldsymbol{{u}}_0\cdot \nabla T^2 \,\text{d} V-\frac{1}{2} \int _{V_i} \rho _0 c^{\prime}_{p0}T^2 \boldsymbol{{u}}_0\cdot \nabla T_0 \,\text{d} V\nonumber \\ &\quad -\frac{1-\alpha _i}{2}\int _{A_{\text{out}}}\rho _0 c^{\prime}_{p0} T_0 T^{2} w_0 \,\text{d} S, \end{align}

(4.9g) \begin{align} \Pi_{\lambda} & ={\alpha_i} \int_{A_{\text{fs}}} \lambda^{\prime}_0 {T^2} \nabla {T_0} \cdot \boldsymbol{{n}} \,\text{d} S - \frac{1}{2}\int_{{V_i}} {\lambda^{\prime}_0} \nabla {T_0} \cdot \nabla {T^2} \,\text{d} V, \end{align}

(4.9h) \begin{align}{\partial }_t E^{\prime}_T &=\frac{1}{2} \int _{V_i} T_0 \rho _0 c^{\prime}_{p0}{\partial }_t T^2 \,\text{d} V, \end{align}

where the index $i$ has been suppressed for the thermophysical properties of the two phases.

The terms (4.9a)–(4.9c) are well known from the energy budget for constant material properties (see e.g. Nienhuser and Kuhlmann [16]). The surface integrals in (4.9d)–(4.9f) represent rates of change of thermal energy of the gas phase $(\alpha _i=-1)$ due to convective heat fluxes through the outlet boundary $A_{\text{out}}$ of the gas container. These fluxes vanish if the gas container is closed $(\left. \!\!w_0\right |_{A_{\text{out}}}=0)$ or if the temperature is prescribed at the outlet $(\left. \!\!T\right |_{A_{\text{out}}}=0)$ . Since we assume the gas enters the container with a prescribed (basic state) temperature, no perturbation energy is introduced through the inlet. The terms $\Pi _\rho$ , $\Pi _{c_p}$ and $\Pi _\lambda$ arise due to the temperature dependence of the material parameters. They vanish, respectively, if $\rho =\text{const.}$ , $c_p=\text{const.}$ or $\lambda =\text{const}$ . Similar to the thermal perturbation energy density (3.6b), the term (4.9h) also depends on the temporal evolution of the perturbation temperature.

5. Kinetic energy budget

The rate of change of the kinetic energy density $\partial _t\varepsilon _{\text{kin}}$ is derived by linearising the momentum equation with respect to the perturbation quantities followed by a scalar multiplication of the linearised momentum equation with the perturbation velocity $\boldsymbol{{u}}$ . Inserting (3.4a) and (3.4c) in (3.1b), we obtain

(5.1) \begin{align} &\partial _t(\rho \boldsymbol{{u}}_0) + \partial _t(\rho \boldsymbol{{u}}) + \nabla \cdot \left [\rho (\boldsymbol{{u}}_0+\boldsymbol{{u}})(\boldsymbol{{u}}_0+\boldsymbol{{u}})\right ] = -\nabla p_0-\nabla p+\rho _0\boldsymbol{{g}}+\rho^{\prime}_0 T\boldsymbol{{g}} \\ & \quad+\nabla \cdot \big\{\mu \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ]+\nabla \cdot \big\{\mu \left [\nabla \boldsymbol{{u}} + (\nabla \boldsymbol{{u}})^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu (\nabla \cdot \boldsymbol{{u}} )\boldsymbol{{I}}\right ]. \nonumber \end{align}

Linearising this equation with respect to the perturbation quantities by neglecting terms of $O(\epsilon ^2)$ yields

(5.2) \begin{align} &\partial _t(\rho \boldsymbol{{u}}_0) + \partial _t(\rho \boldsymbol{{u}}) + \nabla \cdot (\rho \boldsymbol{{u}}_0\boldsymbol{{u}}_0)+\nabla \cdot \left [\rho (\boldsymbol{{u}}_0\boldsymbol{{u}}+\boldsymbol{{u}}\boldsymbol{{u}}_0)\right ] = -\nabla p_0-\nabla p+\rho _0\boldsymbol{{g}}+\rho^{\prime}_0 T\boldsymbol{{g}} \\ & \quad+\nabla \cdot \big\{\mu \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ]+\nabla \cdot \big\{\mu \left [\nabla \boldsymbol{{u}} + (\nabla \boldsymbol{{u}})^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu (\nabla \cdot \boldsymbol{{u}} )\boldsymbol{{I}}\right ]. \nonumber \end{align}

Now the Taylor expansion of the material parameters (3.7a)–(3.7d) is inserted in (5.2) to obtain, after linearisation,

(5.3) \begin{align} &\boldsymbol{{u}}_0\left .{\partial }_T\rho \right |_{T_0}\partial _t T +\rho _0\partial _t\boldsymbol{{u}} + \nabla \cdot (\rho _0\boldsymbol{{u}}_0\boldsymbol{{u}}_0)+ \nabla \cdot (\rho^{\prime}_0 T\boldsymbol{{u}}_0\boldsymbol{{u}}_0)+\nabla \cdot \left [\rho _0(\boldsymbol{{u}}_0\boldsymbol{{u}}+\boldsymbol{{u}}\boldsymbol{{u}}_0)\right ] \\ & \quad= -\nabla p_0-\nabla p +\nabla \cdot \big\{\mu _0 \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \} -\frac{2}{3}\nabla \cdot \left [\mu _0 (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ]+\nabla \cdot \big\{\mu^{\prime}_0 T \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \}\nonumber \\ & \qquad+\rho _0\boldsymbol{{g}}+\rho^{\prime}_0 T\boldsymbol{{g}} -\frac{2}{3}\nabla \cdot \left [\mu^{\prime}_0 T (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ]+\nabla \cdot \big\{\mu _0 \left [\nabla \boldsymbol{{u}} + (\nabla \boldsymbol{{u}})^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu _0 (\nabla \cdot \boldsymbol{{u}} )\boldsymbol{{I}}\right ]. \nonumber \end{align}

Separating again the orders of magnitude yields the basic state momentum equation at $O(\epsilon ^0)$

(5.4) \begin{equation} \nabla \cdot (\rho _0\boldsymbol{{u}}_0\boldsymbol{{u}}_0)= -\nabla p_0+\rho _0\boldsymbol{{g}}+\nabla \cdot \big\{\mu _0 \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \} -\frac{2}{3}\nabla \cdot \left [\mu _0 (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ], \end{equation}

and the momentum perturbation equation at $O(\epsilon ^1)$

(5.5) \begin{align} &\boldsymbol{{u}}_0\rho^{\prime}_0\partial _t T +\rho _0{\partial }_t \boldsymbol{{u}}+ \nabla \cdot (\rho^{\prime}_0 T\boldsymbol{{u}}_0\boldsymbol{{u}}_0)+\nabla \cdot \left [\rho _0(\boldsymbol{{u}}_0\boldsymbol{{u}}+\boldsymbol{{u}}\boldsymbol{{u}}_0)\right ] \\ & \quad= -\nabla p+\rho^{\prime}_0 T\boldsymbol{{g}}+\nabla \cdot \big\{\mu^{\prime}_0 T \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu^{\prime}_0 T (\nabla \cdot \boldsymbol{{u}}_0 )\boldsymbol{{I}}\right ]\nonumber \\ &\quad+\nabla \cdot \big\{\mu _0 \left [\nabla \boldsymbol{{u}} + (\nabla \boldsymbol{{u}})^{\mathrm{T}} \right ]\big \}-\frac{2}{3}\nabla \cdot \left [\mu _0 (\nabla \cdot \boldsymbol{{u}} )\boldsymbol{{I}}\right ]. \nonumber \end{align}

Finally, the scalar product between the momentum perturbation equation (5.5) is taken with the perturbation velocity field $\boldsymbol{{u}}$ , yielding

(5.6) \begin{align} {\underbrace{\vphantom{\bigl [} \rho _0\boldsymbol{{u}}\cdot \partial _t\boldsymbol{{u}}}_{\texttt{K1}}} & = -{\underbrace{\vphantom{\bigl [} \rho^{\prime}_0 \boldsymbol{{u}}_0\cdot \boldsymbol{{u}}\partial _t T}_{\texttt{K2}}} -{\underbrace{\vphantom{\bigl [} \boldsymbol{{u}}\cdot \{\nabla \cdot \left [\rho _0(\boldsymbol{{u}}_0\boldsymbol{{u}}+\boldsymbol{{u}}\boldsymbol{{u}}_0)\right ]\}}_{\texttt{ K3}}} -{\underbrace{\vphantom{\bigl [} \boldsymbol{{u}}\cdot \left [\nabla \cdot (\rho^{\prime}_0 T\boldsymbol{{u}}_0\boldsymbol{{u}}_0)\right ]}_{\texttt{K4}}} \\ & -{\underbrace{\vphantom{\frac{2}{3}} \boldsymbol{{u}}\cdot \nabla p}_{\texttt{K5}}} +{\underbrace{\vphantom{\frac{2}{3}} \rho^{\prime}_0 T\boldsymbol{{u}}\cdot \boldsymbol{{g}}}_{\texttt{K6}}} +{\underbrace{\vphantom{\frac{2}{3}} \boldsymbol{{u}}\cdot \big\{\nabla \cdot \big\{\mu _0 \left [\nabla \boldsymbol{{u}} + (\nabla \boldsymbol{{u}})^{\mathrm{T}} \right ]\big \}\big\}}_{\texttt{K7}}} -{\underbrace{\frac{2}{3}\boldsymbol{{u}}\cdot \big\{\nabla \cdot \left [\mu _0 (\nabla \cdot \boldsymbol{{u}})\boldsymbol{{I}}\right ]\big\}}_{\texttt{K8}}}\\ & \quad \quad \quad \quad \quad \quad +{\underbrace{\vphantom{\frac{2}{3}} \boldsymbol{{u}}\cdot \big\{\nabla \cdot \big\{\mu^{\prime}_0 T \left [\nabla \boldsymbol{{u}}_0 + (\nabla \boldsymbol{{u}}_0)^{\mathrm{T}} \right ]\big\}\big\}}_{\texttt{K9}}} -{\underbrace{\frac{2}{3}\boldsymbol{{u}}\cdot \big\{\nabla \cdot \left [\mu^{\prime}_0 T (\nabla \cdot \boldsymbol{{u}}_0)\boldsymbol{{I}}\right ]\big\}}_{\texttt{K10}}}. \end{align}

Equation (5.6) represents the balance of kinetic energy density at order $O(\epsilon ^2)$ . The first term K1 is recognised as the rate of change of the kinetic energy density of the perturbation flow $\partial _t \varepsilon _{\text{kin}}$ . The physical processes leading to the change of energy density appear on the right-hand side of (5.6). Similar to the thermal budget, the term K2 describes a rate of change of the kinetic perturbation energy density due to the temperature dependence of the density. This term is conservative in the sense that it vanishes when integrated over the volume, as explained in Appendix B.

The terms K3 and K4 describe the rate of change of kinetic perturbation energy density due to the divergence of kinetic perturbation energy flux densities. These fluxes arise due to the transfer of momentum between the basic and the perturbation flow (K3) and due to the second-order density dependence on the temperature (K4), after evaluation of the divergence ( $\nabla \rho^{\prime}_0=\rho^{\prime\prime}_0 \nabla T$ ). The term K5 describes the work per volume and time done by pressure forces which is enabled by the weak compressibility of the perturbation flow due to spatial variation of $\rho$ . The term K6 represents the work done by buoyancy forces. It also arises in the framework of the Oberbeck–Boussinesq approximation.

The remaining terms K7 to K10 describe the rate of change of kinetic perturbation energy density due to viscous dissipation of the perturbation flow (K7), corrected by the effects due to the spatial variation of the density (K8), the spatial variation of the dynamic viscosity (K9) and the spatial variation of both, density and dynamic viscosity (K10).

As for the thermal energy budget, the integral kinetic energy budget is detailed in Appendix B. Integration over the volume occupied by the liquid and the gas separately yields the total rate of change of kinetic energy $\partial _t E_{\text{kin}} = \partial _t \int _{V_i} \varepsilon _{\text{kin}}\,\text{d} V$ in the phase $i$

(5.7) \begin{equation} {\partial }_t E_{\mathrm{kin}}=-D_{\text{kin}}+M_r+M_\varphi +M_z+I+B+K_{\text{G}}+\Lambda _\rho +\Lambda _\mu +\Lambda _{\rho \mu }, \end{equation}

where we introduced the abbreviations

(5.8a) \begin{align} D_{\text{kin}}& =\int _{V_i} \mu _0(\nabla \boldsymbol{{u}}):(\nabla \boldsymbol{{u}}) \,\text{d} V+\alpha _i\int _{A_{\text{fs}}} \mu _0(h_0 h_{0zz} w^2-v^2)\,\text{d} \varphi \,\text{d} z, \end{align}

(5.8b) \begin{align} M_r & = \alpha _i\int _{A_{\text{fs}}}\mu _0 h_0 h_{0z} u\left (\partial _r w - \partial _z u\right ) \,\text{d} \varphi \, \,\text{d} z, \end{align}

(5.8c) \begin{align} M_\varphi & = \alpha _i\int _{A_{\text{fs}}}\mu _0 h_0 v \left (\partial _r v-\frac{v}{h_0}- h_{0z} \partial _z v\right ) \,\text{d} \varphi \, \,\text{d} z, \end{align}

(5.8d) \begin{align} M_z & = \alpha _i\int _{A_{\text{fs}}}\mu _0 h_0 w \left (\partial _r w+h_{0zz}w-h_{0z} \partial _z w\right ) \,\text{d} \varphi \, \,\text{d} z, \end{align}

(5.8e) \begin{align} I & = \sum _{j=1}^5 I_j = -\int _{V_i} \rho _0\left ( u_0\frac{v^2}{r}+u^2\partial _r u_0+u w \partial _z u_0+u w \partial _r w_0+w^2\partial _z w_0 \right )\,\text{d} V, \end{align}

(5.8f) \begin{align} B &= -\int _{V_i}\rho^{\prime}_0 T g w\,\text{d} V, \end{align}

(5.8g) \begin{align} K_{\text{G}}& =-\frac{1-\alpha _i}{4}\int _{A_{\mathrm{out}}} \rho _0w^2w_0\,\text{d} S, \end{align}

(5.8h) \begin{align} \Lambda _\rho &= -\int _{V_i}\rho^{\prime}_0 T\boldsymbol{{u}}\cdot (\boldsymbol{{u}}_0\cdot \nabla \boldsymbol{{u}}_0)\,\text{d} V+\int _{V_i}\zeta \left (p- \frac{1}{3}\mu _0\zeta \right ) \,\text{d} V, \end{align}

(5.8i) \begin{align} \Lambda _\mu &=\int _{V_i} \mu^{\prime}_0 \boldsymbol{{u}}\cdot [\mathcal{S}+(\nabla \boldsymbol{{u}})^{\text{T}}]\cdot \nabla T_0\,\text{d} V+\int _{V_i} (\mu^{\prime}_0+\mu^{\prime\prime}_0 T_0) \boldsymbol{{u}}\cdot [\mathcal{S}_0+(\nabla \boldsymbol{{u}}_0)^{\text{T}}]\cdot \nabla T\,\text{d} V\nonumber \\ &- \int _{V_i} \mu^{\prime}_0 T(\nabla \boldsymbol{{u}}_0):(\nabla \boldsymbol{{u}})\,\text{d} V+\alpha _i\int _{A_{\text{fs}}} \mu^{\prime}_0 w T\left (N^2 \partial _r w_0 - N^2 h_{0z} \partial _z w_0 - h_{0z}^2 h_{0zz} w_0\right )\,\text{d} \varphi \,\text{d} z, \end{align}

(5.8j) \begin{align} \Lambda _{\rho \mu }&=-\int _{V_i}\mu^{\prime}_0 \zeta _0\left (\frac{1}{3} T\zeta +\boldsymbol{{u}}\cdot \nabla T\right ) \,\text{d} V-\int _{V_i}\left (\mu^{\prime}_0 \zeta +\mu^{\prime\prime}_0 T\right )\boldsymbol{{u}}\cdot \nabla T_0 \,\text{d} V, \end{align}

and

(5.9a) \begin{align} \zeta _0&=\nabla \cdot \boldsymbol{{u}}_0=-\frac{\rho^{\prime}_0}{\rho _0}\boldsymbol{{u}}_0\cdot \nabla T_0, \end{align}

(5.9b) \begin{align} \zeta &=\nabla \cdot \boldsymbol{{u}}=- \frac{1}{\rho _0} \left [\rho^{\prime}_0 \boldsymbol{{u}}\cdot \nabla T_0+\rho^{\prime}_0 \partial _t T+\nabla \cdot (\rho^{\prime}_0 \boldsymbol{{u}}_0 T) \right ], \end{align}

(5.9c) \begin{align} \mathcal{S}_0&=\nabla \boldsymbol{{u}}_0+(\nabla \boldsymbol{{u}}_0)^{\mathrm{T}}, \end{align}

(5.9d) \begin{align} \mathcal{S}&=\nabla \boldsymbol{{u}}+(\nabla \boldsymbol{{u}})^{\mathrm{T}}. \end{align}

As before, the index $i$ indicating the phase (liquid or gas) has been suppressed for the thermophysical properties.

The terms (5.8a)–(5.8e) are the known terms for an incompressible flow and constant material properties [16]. $D_{\text{kin}}$ denotes the viscous dissipation, $I$ the kinetic energy production including effects like, for example, the lift-up process, and $M_r$ , $M_\varphi$ and $M_z$ represent the work per time done by thermocapillary forces on the interface $A_{\text{fs}}$ in the radial, azimuthal and axial direction, respectively. The well-known buoyancy production term $B$ also enters the kinetic energy budget within the Boussinesq approximation. $K_{\text{G}}$ represents the advection with the basic flow of perturbation kinetic energy through the outlet of the gas $A_{\text{out}}$ . It vanishes for a closed container holding the gas $(\left .\!\!w_0\right |_{A_{\text{out}}}=0)$ . Note we assume that no perturbation momentum is introduced by advection through the inlet of the gas. The new terms (5.8h)–(5.8j) arise due to the temperature dependence of the material parameters. $\Lambda _\rho$ and $\Lambda _\mu$ vanish, if $\rho =\text{const.}$ or $\mu =\text{const.}$ , respectively, while $\Lambda _{\rho \mu }$ vanishes if either $\rho =\text{const.}$ or $\mu =\text{const.}$

6. Discussion

The orders of magnitude of the terms arising in the energy budgets (4.8) and (5.7) depend on the physicochemical properties of the two fluids as well as on their variability. To estimate the effect of fully temperature-dependent (FTD) properties on the linear stability boundary, we compare the results computed with those obtained using the Oberbeck–Boussinesq approximation (OB) in which all material parameters are assumed constant except for the density in the buoyancy term $\rho \boldsymbol{{g}}$ of (3.1b), which is considered up to first order in $\hat T - \bar T$ .

By considering a Taylor expansion up to first order around the reference values, Gray and Giorgini [3] found that a deviation of 5% of the thermophysical parameters from the value at the reference temperature is an acceptable tolerance to use the OB approximation. Here we make the same assumption, but keep the higher-order terms of the Taylor expansion. The condition that the absolute relative deviation of any quantity $f\in \{\rho,\lambda,c_p,\mu \}$ from its value at the reference temperature is less than or equal to the threshold value $\xi/2 = 0.05$ leads to

(6.1) \begin{align} \left |\frac{f(\hat T)- f(\bar T)}{f(\bar T)} \right | = \left |\frac{f^{\prime}(\bar T)}{f(\bar T)}(\hat T - \bar T) + \frac{1}{2}\frac{f^{\prime\prime}(\bar T)}{f(\bar T)}(\hat T - \bar T)^2 + \ldots \right | \leq \frac{\xi }{2}, \end{align}

where $\hat T$ can be any temperature arising in the system, bounded by $\bar T \pm \Delta T/2$ . Assuming $f(\hat T)$ is a monotonic function and using the algebraic mean temperature $\bar T$ as the reference temperature (as in Gray and Giorgini [3]), we consider the extreme case when $\hat T -\bar T = \pm \Delta T/2$ . Then we get the restriction of the maximum tolerable relative deviation from the reference value

(6.2) \begin{align} \psi _f \,:\!= \max \left | \pm \frac{f^{\prime}(\bar T)}{f(\bar T)}\Delta T + \frac{1}{4}\frac{f^{\prime\prime}(\bar T)}{f(\bar T)}\Delta T^2 \pm \ldots \right | \leq \xi. \end{align}

In lowest order and for $\xi =0.1$ , we recover the criterion of Gray and Giorgini [3]. If the series is truncated at second order, we obtain $(\Delta T \gt 0)$

(6.3) \begin{align} \psi _f = \psi ^{\text{I}}_f + \psi ^{\text{II}}_f = \left |\frac{ f^{\prime}(\bar T) }{f(\bar T)}\right | \Delta T + \frac{1}{4}\left |\frac{ f^{\prime\prime}(\bar T)}{f(\bar T)} \right | \Delta T^2 \leq \xi. \end{align}

Therefore, if the second-order contribution $\psi ^{\text{II}}_f$ is significant, the criterion of Gray and Giorgini [3] is tightened.

If, instead of the OB approximation, a linearised model for a quantity $f$ is used, it makes sense to ensure that the relative deviation of the quantity $f$ due to its second-order variation from the linear model is sufficiently small. Assuming a threshold of $\xi/2=0.05$ as in Gray and Giorgini [3], this leads to the condition

(6.4) \begin{align} \left | \frac{f(\hat T) - \left [ f(\bar T) + f^{\prime}(\bar T) (\hat T -\bar T)\right ]}{f(\bar T) + f^{\prime}(\bar T) (\hat T -\bar T)} \right | = \frac{1}{2}\left | \frac{ f^{\prime\prime}(\bar T) (\hat T -\bar T)^2 + \ldots }{f(\bar T) + f^{\prime}(\bar T) (\hat T -\bar T)} \right | \le \frac{\xi }{2}. \end{align}

Assuming a monotonic variation with $\hat T$ , by setting $\hat T - \bar T = \pm \Delta T/2$ as above, and by neglecting cubic terms, we obtain

(6.5) \begin{align} \frac{1}{4}\left | \frac{ f^{\prime\prime}(\bar T) \Delta T^2 }{f(\bar T) \pm f^{\prime}(\bar T) \Delta T} \right | = \frac{1}{4}\left | \frac{ [f^{\prime\prime}(\bar T)/f(\bar T)] \Delta T^2 }{1 \pm [f^{\prime}(\bar T)/f(\bar T)] \Delta T} \right | = \frac{ \psi _f^{\text{II}}}{\left | 1 \pm [f^{\prime}(\bar T)/f(\bar T)] \Delta T\right |} \le \xi, \end{align}

Maximising the left-hand side, we get the condition

(6.6) \begin{align} \frac{ \psi _f^{\text{II}}}{\left | 1 - \psi _f^{\text{I}}\right |} \le \xi. \end{align}

It is well known that in experiments on thermocapillary liquid bridges even the first-order bound $\psi ^{\text{I}}\le 0.1$ provided by Gray and Giorgini [3] can be violated by the viscosity $(f=\mu )$ . Therefore, the dependence of the liquid viscosity on the temperature has already been taken into account up to first order in the stability analysis of Kozhoukharova et al. [7] $({\textit{Pr}}_{\text{L}}=\bar{\mu }_{\text{L}} \bar{c}_{p \text{L}}/ \bar{\lambda }_{\text{L}}=4)$ . In their numerical simulations for ${\textit{Pr}}_{\text{L}}\in [1,5]$ Melnikov et al. [13] found a significant impact of the linear temperature dependence of the viscosity on the linear stability boundary.

Since the functional dependence of the thermophysical properties on the temperature is not restricted in our investigation, also the effect of a higher-order temperature dependence is of interest. It is difficult, however, to quantify the effect of the FTD approach on the stability boundary without specifying the fluids, owing to the wide range of different fluids employed for liquid bridges. Therefore, we focus on two different cases: a high- and a low-Prandtl-number liquid bridge being heated from above.

6.1 High-Prandtl-number instability

Linear stability analyses have been carried out for the following setting. The length and radius of the liquid bridge are $d=1.65$ mm and $R=d/\Gamma$ , respectively, where $\Gamma =0.66$ is the aspect ratio. The liquid is 2-cSt silicone oil (KF96L-2cs, Shin-Etsu Chemical, Co., Ltd., Japan) which has a Prandtl number of ${\textit{Pr}}_{\text{L}}=28.84$ at the arithmetic mean (reference) temperature $\bar{T}=25^\circ$ C. The discrete data of $\rho _{\text{L}}$ , $\lambda _{\text{L}}$ and $c_{p\text{L}}$ for 2-cSt silicone oil provided by Shin-Etsu [24] have been fitted by least-squares to polynomials of second order. A low polynomial order is used to avoid non-physical oscillations. Since the manufacturer does not specify the temperature dependence of the surface tension, we have to stick to the linear dependence provided in Romanò et al. [20] (see Table 1). The function $\mu _{\text{L}}(\hat{T})$ is constructed from the exponential temperature dependence of the kinematic viscosity as in Ueno et al. [31] and by a quadratic fit of the density. The volume ratio of the liquid is kept constant at $\mathcal{V}=V_{\text{L}}/\pi R^2d =0.9$ . The liquid bridge is placed in a wide test chamber filled with air and confined by no-penetration ( $w_{\text{G}}\equiv 0$ ) adiabatic walls. The temperature dependence of the properties of the gas is based on explicit formulae of VDI Heat Atlas [32]. The reference values of all physical properties are given in Table 1 for both working fluids. The geometry of the test chamber (subscript tc) is defined through the radius ratio $\eta =R_{\text{tc}}/R=4$ , and the total height of the gas space is $d_{\text{tc}}=3.65$ mm within which the liquid bridge is positioned coaxially and vertically centred. Further details on the numerical methods and the explicit temperature dependence of the fluid properties will be provided in Stojanović et al. [28].

Table 1. Thermophysical reference quantities of 2-cSt silicone oil and air at $25^\circ$ C

Fixing $\bar{T}=25^\circ$ C, the condition (6.3) can be rewritten in terms of maximum allowable temperature differences for the OB approximation. Truncating (6.2) after the first and second order, we define the temperature thresholds (symmetric about $\bar T$ ), respectively, as

(6.7) \begin{equation} \Delta T_{\text{OB}}^{\text{I}} \,:\!= \xi \left | \frac{\bar f}{\bar f^{\prime}}\right |\quad \text{and}\quad \Delta T_{\text{OB}}^{\text{II}} \,:\!= 2\frac{\sqrt{|\bar{f}^{\prime}|^2 + \xi |\bar{f}^{\prime\prime}| |\bar{f}|}-|\bar{f}^{\prime}|}{|\bar{f}^{\prime\prime}|}, \end{equation}

where $\bar{f}=f(\bar{T})$ , $\bar{f}^{\prime}=f^{\prime}(\bar{T})$ and $\bar{f}^{\prime\prime}=f^{\prime\prime}(\bar{T})$ . Similarly, the temperature limit of validity for a model accounting for linearly temperature-dependent (LTD) properties can be derived by solving (6.6) for $\Delta T$ to yield

(6.8) \begin{equation} \Delta{T}_{\text{LTD}} \,:\!= \left \{ \begin{array}{ll}\displaystyle -2\xi \frac{|\bar{f}^{\prime}|}{|\bar{f}^{\prime\prime}|} + 2\sqrt{ \xi \frac{|\bar{f}|}{|\bar{f}^{\prime\prime}|} + \xi ^2 \frac{|\bar{f}^{\prime}|^2}{|\bar{f}^{\prime\prime}|^2} }, & \psi _f^{\text{I}} \lt 1, \\ \displaystyle 2\xi \frac{|\bar{f}^{\prime}|}{|\bar{f}^{\prime\prime}|} - 2\sqrt{ -\xi \frac{|\bar{f}|}{|\bar{f}^{\prime\prime}|} + \xi ^2 \frac{|\bar{f}^{\prime}|^2}{|\bar{f}^{\prime\prime}|^2} }, & \psi _f^{\text{I}} \gt 1. \end{array}\right. \end{equation}

The temperature differences $\Delta T_{\text{OB}}^{\text{I}}$ , $\Delta T_{\text{OB}}^{\text{II}}$ and $\Delta{T}_{\text{LTD}}$ are assigned to each thermophysical property of each phase, and they are given in Table 2 for $\xi =0.1$ . The most severe restriction of $\Delta T$ for the validity of the OB approximation is imposed by the condition $\psi ^{\text{I}}_{\mu \text{L}}+\psi ^{\text{II}}_{\mu \text{L}} \lt 0.1$ , not allowing $\Delta T$ to exceed $\Delta T_{\text{OB}}^{\text{II}}=4.6\,$ K. Furthermore, temperature differences greater than $20.2\,$ K would violate condition (6.6) on the viscosity of the liquid. In this case, assuming a linear dependence $\mu _{\text{L}}(\hat T)\sim (\hat T-\bar T)$ would not be sufficient to accurately describe the flow inside the liquid bridge. Besides, the criteria $\psi _f^{\text{II}}$ on $c_{p\text{L}}$ and $c_{p\text{G}}$ get violated for $\Delta T \gt 121.6\,$ K and $\Delta T \gt 578\,$ K, respectively. The latter condition is unrealistic and could only be realised by a phase change.

Table 2. Maximum allowable temperature differences $\Delta T_{\text{OB}}^{\text{I}}$ , $\Delta T_{\text{OB}}^{\text{II}}$ and $\Delta{T}_{\text{LTD}}$ based on a tolerance of $\xi =0.1$ and a reference temperature of $\bar{T} = 25^\circ$ C for different thermophysical parameters. $\Delta T_{\text{OB}}^{\text{I}}$ and $\Delta T_{\text{OB}}^{\text{II}}$ represent the validity thresholds for the applied temperature difference when using the OB approximation and assuming a first-order (up to linear) or, respectively, a second-order (up to quadratic) dependence of the thermophysical quantity on the temperature. $\Delta{T}_{\text{LTD}}$ is the validity threshold when using the linear temperature model (LTD). All temperature differences are given in Kelvin for 2-cSt silicone oil (L) and air (G)

In addition, a minimum temperature difference $\Delta T_{\text{min}}=10\chi \bar{T}$ can be obtained from the first condition of (3.3), which is required to justify the omission of the pressure work in (3.1c). For the present liquid and gas, this condition certainly holds true at $\bar{T}=25^\circ$ C, since the minimum required temperature differences are negligibly small with $\Delta T_{\text{min,L}}=2\times 10^{-6}\,$ K and $\Delta T_{\text{min,G}}= 10^{-5}\,$ K, respectively. The second condition of (3.3) does not involve $\Delta T$ , but rather turns into a condition for the length of the liquid bridge, which is $d\leq 585\,$ m in the present case. Thus, neglecting viscous dissipation in (3.1c) is also reasonable for liquid bridges, which confirms the validity of (3.1c).

In Table 3, we compare the critical temperature differences of the linear stability analyses for different approximations of the governing equations. The linear stability boundary for the onset of hydrothermal waves obtained by the present FTD approach is taken as a reference. It is compared with the result obtained using the OB approximation. To demonstrate the effect of the temperature dependence on the stability boundary of a single thermophysical property, we also combine the OB approximation with the temperature dependence of only one property at a time, keeping the remaining thermophysical properties at their reference values. For instance, within the approximation ‘ $\text{OB}+\rho (\hat T)$ ’ the temperature dependence of the fluid densities is taken into account in all the governing equations (3.1a)–(3.1c). From Table 3, it is seen that critical Reynolds number ${\textit{Re}}_c ={\textit{Ma}}_c/{\textit{Pr}}_{\text{L}} = \gamma d \Delta T_c \bar{\rho }_{\text{L}}/\bar{\mu }_{\text{L}}^2$ for the OB approximation deviates strongly (by $\epsilon _c=24.7$ %) from the reference result (FTD). The main reason is that the relatively large change of the liquid viscosity in the range $\bar{T}\pm \Delta T_c/2$ is not taken care of by the OB approximation, resulting in strongly violated conditions with $\psi ^{\text{I}}_{\mu \text{L}}=1.16$ and $\psi ^{\text{I}}_{\mu \text{L}}+\psi ^{\text{II}}_{\mu \text{L}}=1.59$ . Given the exponential behaviour of $\mu _{\text{L}}(\hat T)$ , also condition (6.6) gets violated for $\Delta T_c=55.5\,$ K with $\psi ^{\text{II}}_\mu/|1-\psi ^{\text{I}}_\mu |=2.8$ . Other than that, the OB approximation slightly fails to satisfy the conditions for $\psi ^{\text{I}}_{\lambda \text{L}}=0.14$ , $\psi ^{\text{I}}_{\lambda \text{G}}=0.16$ , $\psi ^{\text{I}}_{\rho \text{G}}=0.19$ and $\psi ^{\text{I}}_{\mu \text{G}}=0.15$ . This explains why the critical Reynolds number for the case ‘ $\text{OB}+\mu (\hat T)$ ’ is the best approximation to the reference value ${\textit{Re}}_c^{\text{FTD}}$ . The small deviation of 2.5% from FTD is due to the remaining approximations made. In contrast, the relative error in ${\textit{Re}}_c$ of $\epsilon _c\approx 22$ % with respect to the FTD model is very large if, instead, the model accounts for the full temperature dependence of only $\rho$ , $\lambda$ or $c_p$ at a time.

Table 3. Critical temperature difference $\Delta T_c$ and critical Reynolds number ${\textit{Re}}_c = \gamma \bar{\rho }_{\text{L}} \Delta T_c d/\bar{\mu }_{\text{L}}^2$ for a slender liquid bridge with $\Gamma =0.66$ and $\mathcal{V}=0.9$ made of 2-cSt silicone oil (see text). Results are given for different approximations. For all models, the critical wave number is $m_c=3$ . The relative deviation $\epsilon _c=({\textit{Re}}_c -{\textit{Re}}_c^{\text{FTD}})/{\textit{Re}}_c^{\text{FTD}}$ is given in [%]

Figure 2. Temperature and velocity distributions of the basic state for $\Delta T=44.49\,$ K along the free surface (a) and across the midplane at $z=0\,$ mm (b). Solid lines: FTD approach. Dashed lines: OB approximation. In (a), $u_{t0}=\boldsymbol{{t}}\cdot \boldsymbol{{u}}_0$ denotes the tangential velocity, where $\boldsymbol{{t}}$ is the unit vector tangent to the interface. The vertical black dashed line in (b) represents the position of the interface $h_0(z=0)$ .

The question arises as to why the critical Reynolds number using the OB approximation is larger than the one for the FTD approach. Inspecting both basic flows in Figure 2, it seems that the dimensional basic flow fields for $\Delta T = 44.49\,$ K do not differ much. The main differences concern the higher plateau temperature (full red line in Figure 2(a)) and the faster surface velocity (full blue line in Figure 2(a), in particular for $z\gt 0$ ) for the FTD model as compared to the OB approximation. These deviations are caused by a liquid viscosity $\mu [T_0(r,z)]$ which is reduced in the hotter regions with $T_0(r,z)\gt \bar T$ from the constant reference viscosity $\bar \mu _{\text{L}}$ in the OB approximation. The relative local viscosity deviation in the liquid (subscript L)

(6.9) \begin{equation} \Delta \mu _{\text{L}} = \frac{\mu _{\text{L}}[T_0(r,z)]-\bar{\mu }_{\text{L}}}{\bar{\mu }_{\text{L}}}. \end{equation}

is illustrated by colour in Figure 3(a). In the FTD model, the local viscosity is more than 60% larger than nominal near the cold wall, whereas near the hot wall and the free surface it is up to 30% smaller than nominal. The reduced viscosity near the hot wall and along the free surface provides less resistance to the flow such that the basic vortex for the FTD model is stronger than for the OB approximation. This is confirmed by the equidistant streamlines in Figure 3(a), where the full/dashed white lines correspond to the FTD/OB model obtained for the same temperature difference. From Figure 3(b), the critical mode arises in the region where the basic temperature gradients are large, and extends further into the region $\mu _{\text{L}}\lt \bar \mu _{\text{L}}$ of lower viscosity. This is confirmed by the loci of maximum thermal production (white crosses in Figure 3) in a region of slightly reduced viscosity $\mu _{\text{L}}\lt \bar \mu _{\text{L}}$ .

Figure 3. Basic state (a) and critical mode (b) for $\Delta T=44.49\,$ K using the FTD model. (a) Local deviation of the viscosity $\Delta \mu _{\text{L}}=[\mu _{\text{L}}(r,z)-\bar{\mu }_{\text{L}}]/\bar{\mu }_{\text{L}}$ (colour) and streamlines (full white lines) in the liquid. The dashed white lines show streamlines obtained with the OB approximation. (b) Critical velocity field (arrows) and critical temperature field (colour) for $m_c=3$ in the ( $r,z$ ) plane in which the local thermal production $j_1+j_2=-\rho _0 T\boldsymbol{{u}}\cdot \nabla (c_{p0}T_0)$ takes one of its maxima (white crosses in (a, b) located at $(r,\,z)=(1.73,\,0.28)\,$ mm) in the bulk. Black lines indicate isotherms of the basic state.

These properties favour the instability by two mechanisms: (a) The stronger basic vortex leads to larger internal temperature gradients in the upper half of the liquid bridge. Therefore, the hydrothermal wave can extract more energy from the basic temperature field than in the case of the OB model. (b) The perturbation vortices which created the temperature perturbations of the hydrothermal wave arise in a region of reduced viscosity and experience less resistance. For these reasons, the critical Reynolds number for the FTD model is significantly lower than for the OB approximation.

To study the instability mechanism itself, we investigate the budget of the thermal perturbation energy which is crucial for the present hydrothermal wave instability and typical for high-Prandtl-number liquids [25, 27, 33]. Figure 4 shows the main contributions to the integral thermal energy budget of the critical mode for the liquid phase (a) and for the gas phase (b). The tilde indicates that the quantities have been normalised by the dissipation term $D_{\text{th}}$ , as usual. The integral rates of change of thermal energy by the most important transfer processes are almost identical among the FTD method (red) and the OB approximation (blue). This is consistent with the integral contributions $\tilde{\Pi }_\rho$ , $\tilde{\Pi }_{c_p}$ and $\tilde{\Pi }_\lambda$ to the thermal energy budget being very small in the present FTD approximation (Table 4). They are thus negligible. Within the OB approximation, they vanish by definition. Therefore, the temperature dependence of the material parameters does not alter the general instability mechanism discussed, for example, in Stojanović et al. [27]. Note the close agreement of the energy budgets between the FTD and OB models on the stability boundary does not preclude different critical Reynolds numbers, as the terms in the energy budgets are only relative (normalised) quantities.

Table 4. Minor contributions to the thermal energy budgets of the critical mode for the FTD approach

Figure 4. Main contributions to the thermal energy budget of the critical mode in the liquid phase (a) and the gas phase (b). Results are given for the OB approximation (blue) and FTD approach (red). $J_1$ and $J_2$ are defined in (4.9b).

6.2 Low-Prandtl-number instability

For low-Prandtl-number liquids, the instability mechanism is inertial and the critical mode is stationary [33]. In that case, the kinetic energy budget of the perturbation flow is relevant for the instability. As an example, we consider a liquid bridge made of molten tin and use the reference temperature $\bar{T}=250^\circ$ C which is slightly above the melting temperature $T_m = 231.97^\circ$ C [21]. Thus, the Prandtl number is ${\textit{Pr}}_{\text{L}} = 0.0185$ . The functional dependence of the thermophysical properties of molten tin on the temperature is taken from Gancarz et al. [1] and Savchenko et al. [21], either through explicitly given correlations or by fitting quadratic polynomials to tabulated data. Since buoyancy plays a lesser role for low-Prandtl-number liquids [16], we assume weightlessness conditions. Moreover, we select $\Gamma =\mathcal{V}=1$ which allows for a comparison of the critical parameters with data from the literature. The length of the liquid bridge, the chamber geometry, the boundary conditions, and the gas are the same as for the high-Prandtl-number liquid bridge from Section 6.1.

The main contributions, normalised by $D_{\text{kin}}$ , to the kinetic energy budgets of the critical modes are shown in Figure 5 for both approximations FTD (red) and OB (blue). The tilde sign is here employed to denote the terms of the kinetic energy budget normalised by $D_{\text{kin}}$ . Both methods yield almost the same result, which is consistent with the kinetic energy budget obtained by Wanschura et al. [33]. This is consistent with Table 5 , where the obtained critical temperature differences safely fall into the validity range of the OB approximation given in Table 6. Note that increasing the reference temperature to $\bar{T}=500^\circ$ C leads to an extension of the validity range as the variability of the viscosity decreases for higher reference temperatures. Owing to the extremely small dynamic surface deformations, the radial Marangoni production terms $\tilde{M}_r,\, \tilde{M}_{r,\text{G}}\lt 10^{-4}$ are negligible. The production due to buoyancy $\tilde{B}$ vanishes by definition and, for the closed chamber considered, $\tilde{K}_{\text{G}}=0$ . It is clear from Figure 5(a) that most kinetic energy is produced by the inertial process described by $\tilde I_4$ with the work done by Marangoni forces (mainly $\tilde M_\varphi$ ) being very small. As can be seen from Figure 5(b), practically no inertial energy production takes place in the gas phase ( $\tilde{I}_{\text{1,G}}, \dots, \tilde{I}_{\text{5,G}}\lt 10^{-2}$ ). The perturbation flow in the gas is driven by axial $(\tilde M_{z,\text{G}})$ and mainly azimuthal thermocapillary forces $(\tilde M_{\varphi,\text{G}})$ , but the produced kinetic energy is readily dissipated $(\tilde D_{\text{kin,G}})$ . Thus, in the present two-phase system, the gas phase only plays a passive role for the instability mechanism. This also holds true for high-Prandl-number liquids [27]. Owing to the small critical temperature difference $\Delta T_c$ , the new contributions (5.8h)–(5.8j) remain negligibly small.

Table 5. Critical temperature differences and Reynolds numbers for the first instability in a liquid bridge made from tin at $\bar T = 250^\circ$ C with ${\textit{Pr}}_{\text{L}}=0.0185$ , $\Gamma =1$ and ${\mathcal V}=1$ . For the other parameters, see the text. Results are given for different approximations. The critical Reynolds number of Wanschura et al. [33] was obtained by linear interpolation of their data for different  ${\textit{Pr}}_{\text{L}}$ (their table 3)

Table 6. Validity ranges $\Delta T\leq \Delta T_{\text{OB}}^{\text{I}}(\bar{T})$ and $\Delta T\leq \Delta T_{\text{OB}}^{\text{II}}(\bar{T})$ of the OB approximation for each thermophysical property of molten tin at $\bar{T} = 250^\circ$ C and at $\bar{T} = 500^\circ$ C, respectively, using $\xi =0.1$ . All temperature differences are given in K

Figure 5. Main contributions to the kinetic energy budgets of the critical modes assuming constant properties (OB approximation) and fully temperature-dependent fluid properties (FTD). (a) Liquid phase. (b) Gas phase. $I_1$ to $I_5$ are defined in (5.8e).

7. Conclusions

Variable material properties are important in high-Prandtl-number liquid bridges, because the temperature difference is typically large such that the viscosity can vary over a wide range. This variation is particularly important for very small-scale liquid bridges for which the critical temperature difference $\Delta T_c \sim d^{-1}$ must be even larger. In that case, there is some ambiguity (through the reference temperature) in defining the Reynolds, Prandtl and Marangoni numbers, and the critical Reynolds numbers for different approximations of the governing equations may deviate significantly. The dependence of the critical Marangoni number on the choice of the reference temperature has already been noted by Melnikov et al. [13] who demonstrated that using the cold wall temperature as the reference temperature, $\bar T = T_{\text{cold}}$ , leads to a significant reduction of the critical Marangoni number (depending on the amount of variation of the viscosity) as compared to when the mean temperature is used as a reference. While using $\bar T = T_{\text{cold}}$ is convenient from an experimental point of view, because $\Delta T_c$ is initially unknown and the reference Prandtl number does not depend on the (varying) temperature difference, it is not so well suited to correlate the critical Marangoni numbers for different experimental realisations with different critical temperature differences.

Another aspect is the use of the OB approximation beyond its strict range of validity. Even when using the algebraic mean temperature to define the reference material parameters [7], the critical Reynolds number can still significantly depend on the approximation made. It was shown that in the high-Prandtl-number case considered, higher-order variations of the liquid’s viscosity need to be taken into account beyond a certain value for $\Delta T$ . On the other hand, it is more than sufficient to assume a linear dependence of $\rho$ and $\lambda$ on $\hat T$ for silicone oil and for air near room temperature. Moreover, the temperature dependence of $c_{p\text{L}}$ and $c_{p\text{G}}$ is negligible. Finally, we note that the free surface temperature depends on the thermal conditions in the gas phase. For instance, a weak forced axial gas flow can strongly affect the critical conditions [2, 6, 30, 35, 37].

In the future, it would be interesting to investigate the linear stability of very small-scale liquid bridges under extreme temperature gradients. In this case, the model needs to be extended by including the effects of evaporation to correctly describe the physics close to the liquid’s boiling temperature.