1 Introduction
Evolution of temperature in an incompressible liquid (Landau & Lifshitz Reference Landau and Lifshitz1987) is expressed as a balance between advection and conduction of heat as well as energy dissipation due to internal friction, termed viscous heating in this article
$$\begin{eqnarray}\unicode[STIX]{x1D70C}c_{p}(\unicode[STIX]{x1D717}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\text{grad}\unicode[STIX]{x1D717})=k_{th}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}+\unicode[STIX]{x1D70F}_{ik}^{\prime }\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}},\quad i,k=1,2,3,\end{eqnarray}$$
where summation is implied on repeated indices,
$\unicode[STIX]{x1D717}$
is the liquid temperature,
$c_{p}$
is the specific heat at constant pressure,
$\unicode[STIX]{x1D70C}$
the density of the liquid,
$k_{th}$
the thermal conductivity of the liquid and the liquid velocity
$\boldsymbol{u}$
is determined from conservation of mass and linear momentum. Here
$\unicode[STIX]{x1D70F}_{ik}^{\prime }=\unicode[STIX]{x1D702}((\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{k})+(\unicode[STIX]{x2202}u_{k}/\unicode[STIX]{x2202}x_{i}))$
,
$i,k=1,2,3,$
is the viscous stress tensor and
$\unicode[STIX]{x1D702}$
the liquid viscosity.
Consider a two-dimensional thin liquid film of thickness
$h(x,t)$
on a heated solid substrate surrounded by a colder ambient gas phase where interfacial shear stresses arise due to the variation of the surface tension
$\unicode[STIX]{x1D70E}$
with temperature
$\unicode[STIX]{x1D717}$
(see figure 1). If the surface tension decreases with growing temperature, valid for common liquids, then the static state is unstable; a surface flow from hotter to colder areas develops. Since the liquids are viscous, their bulk is dragged along; thus, interfacial temperature gradients induce bulk-liquid motion. This is the thermocapillary effect (Davis Reference Davis1987; Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002).
Thermocapillary instabilities arise because a corrugated interface in a fixed temperature gradient of the substrate and the gas phase experiences hot regions at valleys and cool regions at peaks. The thermocapillary effect drives interfacial liquid from depressions to peaks and since the liquid is viscous, bulk fluid as well, see figure 2. This leads to film thinning, rupture and substrate dewetting.
The thermocapillary effect (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002): a free interface
$z=h(x,t)$
separates a cooler gas phase at temperature
$\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}_{\infty }$
from a viscous liquid lying on a hotter solid substrate at temperature
$\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}_{0}>\unicode[STIX]{x1D717}_{\infty }$
;
$u$
and
$w$
are horizontal and vertical liquid velocities, respectively. Interfacial stresses arising from surface tension gradients, which in turn are due to temperature gradients, drive the liquid from hot depressions to cool peaks leading to depression deepening, film rupture and substrate dewetting (thermocapillary instability), see figure 2. In this article we argue that viscous heating, that is, energy dissipated by viscous stresses during the motion of the liquid, leads internal currents to inhibit heat transfer and promotes film stabilization. In doing so, viscous heating drives two subsurface vortices that effectively halve the wavelength of disturbances (see figure 6).
In this article we show that viscous heating, represented by the last term in (1.1) and neglected in earlier thin liquid-film studies, leads to qualitatively new effects. In particular, it averts the thermocapillary instability and heals the film. The effect can be explained as follows. A critical temperature difference
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}_{cr}$
or equivalently a critical Marangoni number
$M_{cr}$
needs to be exceeded for the instability to be maintained. The viscous heating term in (1.1) slows down the liquid, inhibits heat transfer and prevents the system from reaching the critical threshold of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}$
or equivalently
$M$
. The system crosses into the stability regime leading to suppression of thinning and of substrate dewetting. That viscous heating leads to liquid deceleration and inhibition of heat transfer has also been observed in other physical systems where temperature gradients are in operation, for instance, in Bénard convection (where the liquid density is temperature dependent) when heat is passed from the solid substrate to the fluid (Gebhart Reference Gebhart1962; Turcotte et al.
Reference Turcotte, Hsui, Torrance and Schubert1974).
Surface tension gradients induced by temperature gradients lead to thinning of depressions in a thin liquid film (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002). Numerical solutions of the thin film equation (2.27) with
$\unicode[STIX]{x1D6F4}_{X}=\unicode[STIX]{x1D6F4}_{X}^{(0)}$
from (3.19). (a) Snapshots of the gas–liquid interface of a liquid film towards rupture. (b) Velocity field at time
$t_{end}$
. (c) Streamlines at time
$t_{end}$
Viscous heating leads to a persistent circulating motion of two counter-rotating vortices lying diametrically opposite to a depression of the liquid–gas interface, see figure 6. These coherent structures effectively halve the wavelength of disturbances as has also been observed in the aforementioned studies of Bénard convection with viscous heating (Turcotte et al. Reference Turcotte, Hsui, Torrance and Schubert1974).
The surface tension
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D717}),\end{eqnarray}$$
is a decreasing function of temperature and in the relevant temperature regime it can be approximated by a linear law. In the classical formulation of the thermocapillary problem (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002), the temperature distribution
$\unicode[STIX]{x1D717}$
at the liquid–gas interface is determined by solving the energy equation (1.1) and neglecting its viscous heating term. Then
$\unicode[STIX]{x1D717}$
depends only on the thickness
$h$
of the film:
$\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}(h)$
. The surface tension gradient
$\unicode[STIX]{x1D70E}_{x}$
that drives the thermocapillary flow can then be easily calculated from (1.2) employing the chain rule
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}}\unicode[STIX]{x1D717}_{h}h_{x},\end{eqnarray}$$
where the leading factor is a constant proportional to the Marangoni number (in this article we place variables into subscripts to denote the corresponding partial derivatives).
On the other hand, when viscous heating is taken into account, the temperature distribution depends not only on film thickness, but also on higher-order derivatives due to gravity, the curvature of the interface and most importantly the surface tension gradient itself;
$\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}(h,h_{x},h_{xxx},\unicode[STIX]{x1D70E}_{x})$
. This functional dependence is dictated by the form the viscous heating term in (1.1) acquires in the case of a thermocapillary thin liquid film since then the temperature depends on the horizontal velocity, which in turn depends on the aforementioned fields (see equations (2.25) and (3.10) in the subsequent discussion). The surface tension gradient
$\unicode[STIX]{x1D70E}_{x}$
is again calculated from (1.2) as
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}}[\unicode[STIX]{x1D717}_{h}h_{x}+\unicode[STIX]{x1D717}_{h_{x}}h_{xx}+\unicode[STIX]{x1D717}_{h_{xxx}}h_{xxxx}+\unicode[STIX]{x1D717}_{\unicode[STIX]{x1D70E}_{x}}\unicode[STIX]{x1D70E}_{xx}],\end{eqnarray}$$
where the leading factor is still proportional to the Marangoni number. It is clear that such a problem for the determination of the surface tension gradient
$\unicode[STIX]{x1D70E}_{x}$
has to be solved self-consistently. This is the main goal of the present article.
Thermocapillary flows have been successfully integrated into engineering applications such as directing microscopic flow on patterned surfaces (Kataoka & Troian Reference Kataoka and Troian1999), moving and routing droplets in microfluidic devices by laser beam localized heating (Cordero et al. Reference Cordero, Burnham, Baroud and McGloin2008) and propelling micron-sized motors by employing incoherent wide-field illumination (Maggi et al. Reference Maggi, Saglimbeni, Dipalo, De Angelis and Di Leonardo2015).
Identifying mechanisms that would stabilize an unstable thin liquid film is important in medicine (pulmonary airway closure (Halpern & Grotberg Reference Halpern and Grotberg1992) or tear film rupture (Braun Reference Braun2012)) as well as in engineering applications, for instance, in the manufacturing of light-weight structures formed by solidifying a foam network of liquid bridges (Stewart & Davis Reference Stewart and Davis2013). Theoretical methodologies have been devised to stabilize an unstable liquid film based on mechanisms of surface shear (for instance by blowing air (Davis, Gratton & Davis Reference Davis, Gratton and Davis2010) or inducing temperature gradients (Kataoka & Troian Reference Kataoka and Troian1998)), body forces (Benney Reference Benney1966) and magnetic torques (Kirkinis Reference Kirkinis2017). These cases have been reviewed and classified by Kirkinis & Davis (Reference Kirkinis and Davis2015).
A simplified picture can be gained when viscous heating is taken to be a constant. Such modelling was employed in molecular dynamics simulations (Vo & Kim Reference Vo and Kim2016) to predict nanoscale fluid transport phenomena establishing good agreement between simulation and the continuum equations. Viscous heating in Couette flow with temperature-dependent viscosity leads to double-valued steady solutions (two different velocities associated with the same value of stress) and beyond a certain value of stress no solutions exist (Joseph Reference Joseph1965; Johns & Narayanan Reference Johns and Narayanan1997; Subrahmaniam, Johns & Narayanan Reference Subrahmaniam, Johns and Narayanan2002).
The present article is organized as follows. In § 2 we briefly review known results from the literature on the formulation of the thermocapillary problem by reducing the full Navier–Stokes equations according to the geometry displayed in figure 1. The reduction continues further in § 2.1 by resorting to the long wavelength approximation. This is possible to carry out because a thin liquid film possesses a geometrical disparity of spatial scales parallel and normal to the solid substrate. Thus, the very difficult free-boundary problem reduces to a tractable partial differential equation albeit fully nonlinear and including higher-order spatial derivatives.
In § 3 we provide a new formulation of the thermocapillary problem whereby viscous heating gives rise to a surface tension depending on the temperature distribution which now acquires the form
$\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}(h,h_{x},h_{xxx},\unicode[STIX]{x1D70E}_{x})$
.
In § 4 we perform a mean-field analysis to calculate the critical parameters giving rise to stabilization based on a corresponding analysis by Kirkinis & Davis (Reference Kirkinis and Davis2015). These estimates are in agreement with numerical simulations of the evolution equation.
In appendix A we show analytically how the viscous heating term in (1.1) leads to a commensurate heat-transfer inhibition.
In appendix B we extend the formulation developed in this article to the case of liquids endowed with rotational degrees of freedom (such as magnetic liquids). In this case new viscous heating terms arise in the energy equation, analogous to (1.1) taking into account not only spatial variations of the liquid velocity (the strain-of-rate tensor) but also the imbalance between magnetic particle angular velocity and liquid angular velocity (one half of its vorticity) as well as the spatial variation of the former. This same formulation also leads to new insights regarding the motion of magnetic liquids in channels of biological significance. In particular, we briefly discuss the existence of additional mechanisms that lead to an increase of particle temperature which has been observed to destroy malignant cells.
2 Momentum balance and liquid-film evolution equations
Consider a two-dimensional thin liquid film of thickness
$h(x,t)$
on a solid substrate (at
$z=0$
) surrounded by an ambient gas phase (see figure 1) with interfacial shear stresses arising due to the variation of the surface tension with temperature
$\unicode[STIX]{x1D717}$
. A Newtonian viscous liquid is characterized by the Cauchy stress tensor
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ik}=-p\unicode[STIX]{x1D6FF}_{ik}+\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}\right),\quad i,k=1,3,\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}$
is the liquid viscosity,
$p$
the hydrostatic pressure and the incompressibility condition
$\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{i}=0$
has been implicitly assumed to hold.
The Navier–Stokes equations become
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}[u_{t}+uu_{x}+wu_{z}]=-p_{x}+\unicode[STIX]{x1D702}(u_{xx}+u_{zz})-\unicode[STIX]{x1D719}_{x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}[w_{t}+uw_{x}+ww_{z}]=-p_{z}+\unicode[STIX]{x1D702}(w_{xx}+w_{zz})-\unicode[STIX]{x1D719}_{z}, & \displaystyle\end{eqnarray}$$
where
$u$
and
$w$
are the horizontal and vertical components of the velocity distribution (cf. figure 1),
$\unicode[STIX]{x1D70C}$
is the density of the liquid and
$\unicode[STIX]{x1D719}$
is a potential representing the effects of gravity, intermolecular forces or other external fields. In this article we tacitly assume that
$\unicode[STIX]{x1D719}$
is a gravitational potential so that
$\unicode[STIX]{x1D719}=gz$
where
$g$
is the acceleration due to gravity.
At the liquid–gas interface
$z=h(x,t)$
, the shear and normal stresses are, respectively,
$$\begin{eqnarray}\boldsymbol{t}\unicode[STIX]{x1D749}\boldsymbol{n}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x2202}s}\quad \text{and}\quad \boldsymbol{n}\unicode[STIX]{x1D749}\boldsymbol{n}=2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70E},\quad \text{at }z=h(x,t),\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}$
is the temperature dependent surface tension,
$s$
is arc length along the interface and
$\unicode[STIX]{x1D705}$
the mean curvature of the surface (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002).
No slip and no penetration boundary conditions at the liquid–solid interface are expressed as
$$\begin{eqnarray}u=0,\quad \text{and}\quad w=0,\quad \text{at }z=0,\end{eqnarray}$$
respectively.
The liquid–gas interface is described by the function
$z=h(x,t)$
or, equivalently, by the space curve
$\boldsymbol{r}(x,t)=(x,h(x,t))$
. The tangent and normal vectors to the interface are
$$\begin{eqnarray}\boldsymbol{t}=\frac{(1,h_{x})}{\sqrt{1+h_{x}^{2}}}\quad \text{and}\quad \boldsymbol{n}=\frac{(-h_{x},1)}{\sqrt{1+h_{x}^{2}}}.\end{eqnarray}$$
The shear component of the stress tensor on the liquid–gas interface is
$$\begin{eqnarray}\boldsymbol{t}\unicode[STIX]{x1D749}\boldsymbol{n}=\frac{1}{1+h_{x}^{2}}(1,h_{x})\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70F}_{11} & \unicode[STIX]{x1D70F}_{13}\\ \unicode[STIX]{x1D70F}_{31} & \unicode[STIX]{x1D70F}_{33}\end{array}\right)\left(\begin{array}{@{}c@{}}-h_{x}\\ 1\end{array}\right)=\frac{1}{1+h_{x}^{2}}[\unicode[STIX]{x1D70F}_{13}-h_{x}^{2}\unicode[STIX]{x1D70F}_{31}+h_{x}(\unicode[STIX]{x1D70F}_{33}-\unicode[STIX]{x1D70F}_{11})]=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x2202}s}.\end{eqnarray}$$
The normal component of the stress at the liquid–gas interface is
$$\begin{eqnarray}\boldsymbol{n}\unicode[STIX]{x1D749}\boldsymbol{n}=\frac{1}{1+h_{x}^{2}}(-h_{x},1)\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D70F}_{11} & \unicode[STIX]{x1D70F}_{13}\\ \unicode[STIX]{x1D70F}_{31} & \unicode[STIX]{x1D70F}_{33}\end{array}\right)\left(\begin{array}{@{}c@{}}-h_{x}\\ 1\end{array}\right)=\frac{1}{1+h_{x}^{2}}[\unicode[STIX]{x1D70F}_{33}+h_{x}^{2}\unicode[STIX]{x1D70F}_{11}-h_{x}(\unicode[STIX]{x1D70F}_{13}+\unicode[STIX]{x1D70F}_{31})]=2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70E},\end{eqnarray}$$
where
$2\unicode[STIX]{x1D705}=h_{xx}/(1+h_{x})^{3/2}$
. To the above, one may add the kinematic condition
$$\begin{eqnarray}w=h_{t}+uh_{x}\end{eqnarray}$$
at the liquid–gas interface.
In terms of the velocity field the stress conditions become (Oron et al. Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002)
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{t}\unicode[STIX]{x1D749}\boldsymbol{n}=\frac{1}{1+h_{x}^{2}}\unicode[STIX]{x1D702}[(-u_{x}+w_{z})h_{x}+(1-h_{x}^{2})(w_{x}+u_{z})]=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x2202}s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}\unicode[STIX]{x1D749}\boldsymbol{n}=\frac{1}{1+h_{x}^{2}}[-p+2\unicode[STIX]{x1D702}w_{z}-2\unicode[STIX]{x1D702}h_{x}(w_{x}+u_{z})+h_{x}^{2}(-p+2\unicode[STIX]{x1D702}u_{x})]=2\unicode[STIX]{x1D705}\unicode[STIX]{x1D70E}. & \displaystyle\end{eqnarray}$$
2.1 Lubrication approximation
In this approximation the Navier–Stokes equations are scaled with respect to a long wavelength
$\unicode[STIX]{x1D706}$
which is a natural consequence of the longitudinal and vertical disparity of scales embodied in a thin liquid film. This introduces slow dimensionless longitudinal spatial and temporal variables
$$\begin{eqnarray}X=\frac{2\unicode[STIX]{x03C0}x}{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D716}\frac{x}{h_{0}},\quad Z=\frac{z}{h_{0}}\quad \text{and}\quad T=\frac{2\unicode[STIX]{x03C0}U_{0}}{\unicode[STIX]{x1D706}}t,\quad \text{for }\unicode[STIX]{x1D716}=\frac{2\unicode[STIX]{x03C0}h_{0}}{\unicode[STIX]{x1D706}},\end{eqnarray}$$
where
$h_{0}$
and
$U_{0}$
are characteristic scales for the liquid film thickness and longitudinal velocity respectively and the partial derivatives
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}X$
and
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}Z$
are now quantities of unit order (Davis Reference Davis, Batchelor, Moffatt and Worster2002). The long wavelength
$\unicode[STIX]{x1D706}$
is defined in terms of the wavenumber that leads to instability, which in turn can be expressed with respect to the system parameters such as viscosity, thermal conductivity etc. see § 4.1. The dimensionless velocity, pressure and liquid–gas interface location are, respectively,
$$\begin{eqnarray}U=\frac{u}{U_{0}},\quad W=\frac{\unicode[STIX]{x1D706}}{2\unicode[STIX]{x03C0}h_{0}}\frac{w}{U_{0}}=\frac{1}{\unicode[STIX]{x1D716}}\frac{w}{U_{0}},\quad P=\frac{2\unicode[STIX]{x03C0}h_{0}}{\unicode[STIX]{x1D706}}\frac{h_{0}p}{\unicode[STIX]{x1D702}U_{0}},\quad H=\frac{h}{h_{0}}.\end{eqnarray}$$
We define the dimensionless Reynolds numbers, scaled capillary number, external potential and surface tension
$$\begin{eqnarray}Re=\frac{U_{0}h_{0}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D702}},\quad \bar{C}=\unicode[STIX]{x1D716}^{-3}\frac{U_{0}\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D70E}},\quad \unicode[STIX]{x1D6F7}=\frac{\unicode[STIX]{x1D716}h_{0}\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D702}U_{0}},\quad \text{and}\quad \unicode[STIX]{x1D6F4}=\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D702}U_{0}}.\end{eqnarray}$$
With the above notation, equations (2.2) and (2.3) obtain the following dimensionless form
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}Re[U_{T}+UU_{X}+WU_{Z}]=-P_{X}+(\unicode[STIX]{x1D716}^{2}U_{XX}+U_{ZZ})-\unicode[STIX]{x1D6F7}_{X}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{2}Re[W_{T}+UW_{X}+WW_{Z}]=-\frac{1}{\unicode[STIX]{x1D716}}P_{Z}+\unicode[STIX]{x1D716}(\unicode[STIX]{x1D716}^{2}W_{XX}+W_{ZZ})-\frac{1}{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6F7}_{Z}. & \displaystyle\end{eqnarray}$$
The shear and normal components of the stress tensor from relations (2.10)–(2.11) become
$$\begin{eqnarray}\displaystyle \boldsymbol{t}\unicode[STIX]{x1D749}\boldsymbol{n} & = & \displaystyle \frac{\unicode[STIX]{x1D702}U_{0}}{h_{0}}\frac{1}{1+\unicode[STIX]{x1D716}^{2}H_{X}^{2}}[\unicode[STIX]{x1D716}(-U_{X}+W_{Z})\unicode[STIX]{x1D716}H_{X}(1-\unicode[STIX]{x1D716}^{2}H_{X}^{2})(\unicode[STIX]{x1D716}^{2}W_{X}+U_{Z})],\end{eqnarray}$$
and
$$\begin{eqnarray}\boldsymbol{n}\unicode[STIX]{x1D749}\boldsymbol{n}=-p+2\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x1D716}U_{0}}{h_{0}}W_{Z}-2\frac{\unicode[STIX]{x1D716}U_{0}}{h_{0}}\unicode[STIX]{x1D702}H_{X}(U_{Z}-\unicode[STIX]{x1D716}^{2}W_{X})+\unicode[STIX]{x1D716}^{2}H_{X}^{2}\left(-p+2\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x1D716}U_{0}}{h_{0}}\unicode[STIX]{x1D716}U_{X}\right).\end{eqnarray}$$
The kinematic condition (2.9) becomes
$W=H_{T}+UH_{X}$
and leads to the following evolution equation for the profile
$H$
:
$$\begin{eqnarray}H_{T}+\left[\int _{0}^{H}U\,\text{d}Z\right]_{X}=0.\end{eqnarray}$$
This is obtained in a standard manner (Oron et al.
Reference Oron, Davis and Bankoff1997) by integrating the incompressibility condition
$U_{X}+W_{Z}=0$
from
$Z=0$
to
$Z=H(X,T)$
, using the second of (2.5) and integration by parts.
We proceed by expanding all fields in a regular perturbation series with respect to the small parameter
$\unicode[STIX]{x1D716}$
:
$(U,W,P)=(U^{0},W^{0},P^{0})+\unicode[STIX]{x1D716}(U^{1},W^{1},P^{1})+O(\unicode[STIX]{x1D716}^{2})$
. To leading order, which corresponds to the lubrication approximation, and dropping superscripts denoting order with respect to
$\unicode[STIX]{x1D716}$
, the scaled equations (2.15)–(2.16) simplify to
$$\begin{eqnarray}-P_{X}-\unicode[STIX]{x1D6F7}_{X}+U_{ZZ}=0,\quad -P_{Z}-\unicode[STIX]{x1D6F7}_{Z}=0.\end{eqnarray}$$
The boundary conditions (2.5) become
$$\begin{eqnarray}U=0,\quad W=0,\quad \text{at }Z=0,\end{eqnarray}$$
with shear stress
$$\begin{eqnarray}U_{Z}=\unicode[STIX]{x1D6F4}_{X}\quad \text{at }Z=H(X,T),\end{eqnarray}$$
and normal stress
$$\begin{eqnarray}-P=\bar{C}^{-1}H_{XX}\quad \text{at }Z=H(X,T).\end{eqnarray}$$
We may also define the reduced pressure (Oron et al. Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002)
$$\begin{eqnarray}\bar{P}=P+\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}|_{Z=H}-\bar{C}^{-1}H_{XX}.\end{eqnarray}$$
Integration of the first of (2.20) and application of the normal stress boundary condition (2.24) leads to the velocity distribution in the long wavelength approximation (cf. Oron et al. Reference Oron, Davis and Bankoff1997, equation (2.26))
$$\begin{eqnarray}U(X,Z,T)=\unicode[STIX]{x1D6F4}_{X}Z+\bar{P}_{X}({\textstyle \frac{1}{2}}Z^{2}-HZ),\end{eqnarray}$$
where
$$\begin{eqnarray}\bar{P}_{X}=(\unicode[STIX]{x1D6F7}|_{Z=H})_{X}-\bar{C}^{-1}H_{XXX}.\end{eqnarray}$$
The sought-after evolution equation for the liquid–gas interface (2.19) with the velocity profile equation (2.25) and pressure gradient from (2.26) becomes
$$\begin{eqnarray}H_{T}+[{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6F4}_{X}H^{2}-{\textstyle \frac{1}{3}}\bar{P}_{X}H^{3}]_{X}=0.\end{eqnarray}$$
3 Thermocapillary flow with viscous heating
In the system considered, there are temperature variations due to the applied temperature gradient. For such a system conservation of energy must be added to the usual Navier–Stokes equations in the form of (1.1). After some rearrangement and symmetrization, we obtain
$$\begin{eqnarray}\unicode[STIX]{x1D70C}c_{p}(\unicode[STIX]{x1D717}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\text{grad}\unicode[STIX]{x1D717})=k_{th}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}+\frac{\unicode[STIX]{x1D702}}{2}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}\right)^{2}.\end{eqnarray}$$
In the geometry of our problem (cf. figure 1), equation (3.1) reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D70C}c_{p}(\unicode[STIX]{x1D717}_{t}+u\unicode[STIX]{x1D717}_{x}+w\unicode[STIX]{x1D717}_{z})=k_{th}(\unicode[STIX]{x1D717}_{xx}+\unicode[STIX]{x1D717}_{zz})+\unicode[STIX]{x1D702}(2u_{x}^{2}+(u_{z}+w_{x})^{2}+2w_{z}^{2}),\end{eqnarray}$$
with boundary conditions
$$\begin{eqnarray}\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}_{0}\quad \text{at }z=0\quad \text{and}\quad k_{th}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}+\unicode[STIX]{x1D6FC}_{th}(\unicode[STIX]{x1D717}-\unicode[STIX]{x1D717}_{\infty })=0,\quad \text{at }z=h(x,t),\end{eqnarray}$$
where
$\unicode[STIX]{x1D717}_{0}$
is the temperature of the rigid plane,
$\unicode[STIX]{x1D717}_{\infty }$
the temperature of the gas phase and
$\unicode[STIX]{x1D6FC}_{th}$
the heat-transfer coefficient describing the rate of heat transfer from the liquid to the ambient gas phase.
We now introduce the dimensionless temperature field
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}(X,Z,T)=\frac{\unicode[STIX]{x1D717}-\unicode[STIX]{x1D717}_{\infty }}{\unicode[STIX]{x1D717}_{0}-\unicode[STIX]{x1D717}_{\infty }},\end{eqnarray}$$
and define the dimensionless Prandtl, Biot, viscous heating (Brinkman), gravity and Marangoni numbers
$$\begin{eqnarray}Pr=\frac{\unicode[STIX]{x1D702}c_{p}}{k_{th}},\quad B=\frac{\unicode[STIX]{x1D6FC}_{th}h_{0}}{k_{th}},\quad \unicode[STIX]{x1D6EC}=\frac{\unicode[STIX]{x1D702}U_{0}^{2}}{k_{th}\unicode[STIX]{x0394}\unicode[STIX]{x1D717}},\quad G=\frac{\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}gh_{0}^{2}}{\unicode[STIX]{x1D702}U_{0}}\quad \text{and}\quad M=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D702}U_{0}}\unicode[STIX]{x1D716},\end{eqnarray}$$
where
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70E}$
is the variation of surface tension over the temperature domain
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}=\unicode[STIX]{x1D717}_{0}-\unicode[STIX]{x1D717}_{\infty }$
so that
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}.\end{eqnarray}$$
Equation (3.6) is derived from the equation of state approximated by a linear law
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}-\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D717}-\unicode[STIX]{x1D717}_{0}),\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FE}$
is positive for common liquids (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002).
With these dimensionless quantities, the energy equation (3.2) and boundary conditions (3.3) become, respectively,
$$\begin{eqnarray}\unicode[STIX]{x1D716}Re\,Pr(\unicode[STIX]{x1D6E9}_{T}+U\unicode[STIX]{x1D6E9}_{X}+W\unicode[STIX]{x1D6E9}_{Z})=\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6E9}_{XX}+\unicode[STIX]{x1D6E9}_{ZZ}+\unicode[STIX]{x1D6EC}[2\unicode[STIX]{x1D716}^{2}U_{X}^{2}+(U_{Z}+\unicode[STIX]{x1D716}^{2}W_{X})^{2}+2\unicode[STIX]{x1D716}^{2}W_{Z}^{2}],\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}=1\quad \text{at }Z=0,\quad \text{and}\quad \unicode[STIX]{x1D6E9}_{Z}-\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6E9}_{X}H_{X}+B\unicode[STIX]{x1D6E9}(1+\unicode[STIX]{x1D716}^{2}H_{X}^{2})^{1/2}=0,\quad \text{at }Z=H.\end{eqnarray}$$
Employing the long wavelength approximation and considering the leading-order correction, the energy equation (3.8) and boundary conditions (3.9) become, respectively,
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{ZZ}+\unicode[STIX]{x1D6EC}(U_{Z})^{2}=0,\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}=1,\quad \text{at }Z=0,\quad \unicode[STIX]{x1D6E9}_{Z}+B\unicode[STIX]{x1D6E9}=0\quad \text{at }Z=H(X,T).\end{eqnarray}$$
Substituting the velocity distribution (2.25) into the energy equation (3.10) and integrating twice leads to the following succinct form for the temperature distribution
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}=-\,\frac{\unicode[STIX]{x1D6EC}\,(\bar{P}_{X}(Z-H)+\unicode[STIX]{x1D6F4}_{X})^{4}}{12{\bar{P}_{X}}^{2}}+C_{1}(\bar{P}_{X},\unicode[STIX]{x1D6F4}_{X})\,Z+C_{2}(\bar{P}_{X},\unicode[STIX]{x1D6F4}_{X}),\end{eqnarray}$$
for two arbitrary functions
$C_{1}$
and
$C_{2}$
. Employing the boundary conditions (3.11) they obtain the form
$$\begin{eqnarray}C_{1}=-\unicode[STIX]{x1D6EC}\,\frac{(BH^{4}{\bar{P}_{X}}^{3}-4\,BH^{3}{\bar{P}_{X}}^{2}\unicode[STIX]{x1D6F4}_{X}+6\,BH^{2}\bar{P}_{X}{\unicode[STIX]{x1D6F4}_{X}}^{2}-4\,(1+BH){\unicode[STIX]{x1D6F4}_{X}}^{3})}{12\bar{P}_{X}(BH+1)}-\frac{B}{BH+1},\end{eqnarray}$$
and
$$\begin{eqnarray}C_{2}=1+\unicode[STIX]{x1D6EC}\,\frac{(\bar{P}_{X}H-\unicode[STIX]{x1D6F4}_{X})^{4}}{12\bar{P}_{X}^{2}}.\end{eqnarray}$$
Clearly, the temperature distribution (3.12) immediately reduces to the well-known form (Oron et al. Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002)
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}(X,Z,T)=1-\frac{BZ}{1+BH},\end{eqnarray}$$
in the absence of viscous heating, that is, in the absence of the last term in (3.1) (by setting
$\unicode[STIX]{x1D6EC}=0$
in (3.8) and (3.12)).
3.1 Determination of the surface tension gradient
The value of
$\unicode[STIX]{x1D6F4}_{X}$
can be calculated by considering the variation of surface tension at the liquid–gas interface, that is, the variation of
$\unicode[STIX]{x1D6F4}=\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6E9}(H,\bar{P}_{X},\unicode[STIX]{x1D6F4}_{X}))$
. Thus,
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}=-M\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}H}H_{X}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}\bar{P}_{X}}\bar{P}_{XX}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F4}_{X}}\unicode[STIX]{x1D6F4}_{XX}\right]\end{eqnarray}$$
where
$M=-d\unicode[STIX]{x1D6F4}/d\unicode[STIX]{x1D6E9}$
is the Marangoni number defined above.
Note that in the absence of viscous heating only the first term on the right-hand side of (3.16) is present, that is
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}=-M\frac{\text{d}\unicode[STIX]{x1D6E9}}{\text{d}H}H_{X},\quad \text{where }\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D6E9}(H).\end{eqnarray}$$
Since, in the absence of viscous heating,
$\unicode[STIX]{x1D6E9}$
is known explicitly from (3.15), the surface tension gradient can be determined in closed form (cf. equation (3.19) below and Oron et al. (Reference Oron, Davis and Bankoff1997), Davis (Reference Davis, Batchelor, Moffatt and Worster2002)).
3.2 Surface tension gradient expansion with respect to viscous heating parameter
$\unicode[STIX]{x1D6EC}$
The determination of the surface tension gradient proceeds by expressing
$\unicode[STIX]{x1D6F4}_{X}$
from (3.16) in a regular perturbation series with respect to the perturbing parameter
$\unicode[STIX]{x1D6EC}$
,
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}=\unicode[STIX]{x1D6F4}_{X}^{(0)}+\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F4}_{X}^{(1)}+O(\unicode[STIX]{x1D6EC}^{2}).\end{eqnarray}$$
This leads to the following leading-order correction
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}^{(0)}=\frac{BMH_{X}}{(1+BH)^{2}},\end{eqnarray}$$
which is the standard result obtained in the absence of the nonlinear viscous term in (3.8) (setting
$\unicode[STIX]{x1D6EC}\equiv 0$
) (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002).
The first-order term in
$\unicode[STIX]{x1D6EC}$
can be most usefully expressed in divergence form
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}^{(1)}=\left[-\frac{1}{12}\,\frac{MH^{4}\bar{P}_{X}^{2}}{1+BH}+\frac{1}{3}\,\frac{M^{2}BH^{3}H_{X}\bar{P}_{X}}{(1+BH)^{3}}-\frac{1}{2}\,\frac{M^{3}B^{2}H^{2}(H_{X})^{2}}{(1+BH)^{5}}\right]_{X},\end{eqnarray}$$
and similarly for the higher-order corrections.
3.3 Surface tension gradient expansion with respect to the Biot number
An alternative derivation of the surface tension gradient can be obtained by expanding
$\unicode[STIX]{x1D6F4}_{X}$
in a perturbation series with respect to the Biot number
$B$
as
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}=B\unicode[STIX]{x1D6F4}_{X}^{(1)}+B^{2}\unicode[STIX]{x1D6F4}_{X}^{(2)}+O(B^{3}),\end{eqnarray}$$
but keeping the ratio
$\unicode[STIX]{x1D6EC}/B$
constant. The two derivations are asymptotically equivalent and as such they lead to the same bifurcation result (see the discussion of § 4).
The leading-order correction reads
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{X}^{(1)}=\left[MH-\frac{\unicode[STIX]{x1D6EC}}{B}\frac{H^{4}\bar{P}_{X}^{2}}{12}M\right]_{X}.\end{eqnarray}$$
When energy dissipation due to internal friction is neglected (
$\unicode[STIX]{x1D6EC}\equiv 0$
) equation (3.22) reduces to the well-known form of the surface tension gradient (as
$B\rightarrow 0$
) (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002).
The second-order correction reads
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6F4}_{X}^{(2)}=\left[-MH^{2}+\frac{\unicode[STIX]{x1D6EC}}{B}\left(\frac{1}{12}H^{5}M{\bar{P}_{X}}^{2}+\frac{1}{3}H^{3}\bar{P}_{X}H_{X}M^{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\left(\frac{\unicode[STIX]{x1D6EC}}{B}\right)^{2}\left(\frac{1}{18}H^{7}{\bar{P}_{X}}^{2}\bar{P}_{XX}+\frac{1}{9}H^{6}{\bar{P}_{X}}^{3}H_{X}\right)M^{2}\right]_{X},\end{eqnarray}$$
and similarly for higher-order corrections. Thus, the first-order correction
$\unicode[STIX]{x1D6F4}_{X}^{(1)}$
does not contribute to any linear stability analysis and its effects on the thermocapillary instability can only be studied nonlinearly.
Growth rate
$s=k^{2}-k^{4}$
obtained from linear stability analysis (Oron et al.
Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002) of the thin film evolution equation (2.27) showing the band of unstable wavenumbers
$0<k<k_{cr}=1$
.
4 Stability analysis
4.1 Linear stability
A standard linear stability analysis of (2.27) (Oron et al. Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002) leads to the growth rate
$$\begin{eqnarray}s=\left[\frac{MB}{2(1+B)^{2}}-\frac{G}{3}\right]k^{2}-\frac{1}{3\bar{C}}k^{4}.\end{eqnarray}$$
Thus, there is a band of unstable wavenumbers
$0<k<k_{cr}$
where the (dimensionless) critical wavenumber
$k_{cr}$
for the onset of instability is found when
$s$
becomes positive, that is, when
$$\begin{eqnarray}k_{cr}^{2}=\frac{3MB\bar{C}}{2(1+B)^{2}}-\bar{C}G.\end{eqnarray}$$
This requires that the Marangoni number
$M$
exceeds its critical value
$$\begin{eqnarray}M_{cr}=\frac{2G(1+B)^{2}}{3B}.\end{eqnarray}$$
For values
$M=4(1+B)^{2}/B,G=3,\bar{C}=1/3$
we obtain the critical wavenumber
$k_{cr}=1$
and the band of unstable wavenumbers
$0<k<1$
, see figure 3.
It is worthwhile to note here that both the velocity scale
$U_{0}$
and the long wave parameter
$\unicode[STIX]{x1D716}$
drop out when the foregoing or subsequent formulas are expressed in terms of system parameters (that is, in terms of viscosity, surface tension, thermal conductivity and coefficient of heat transfer).
4.2 Nonlinear stability
4.2.1 Expansion with respect to viscous heating parameter
$\unicode[STIX]{x1D6EC}$
We can perform a weakly nonlinear analysis in the vicinity of the onset of instability following the arguments in Kirkinis & Davis (Reference Kirkinis and Davis2015). At this threshold we can neglect the time variation of the problem. Then (2.27) can be integrated once to give
$$\begin{eqnarray}{\textstyle \frac{1}{2}}H^{2}(\unicode[STIX]{x1D6F4}_{X}^{(0)}+\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F4}_{X}^{(1)}+\unicode[STIX]{x1D6EC}^{2}\unicode[STIX]{x1D6F4}_{X}^{(2)})-{\textstyle \frac{1}{3}}H^{3}\bar{P}_{X}=K,\end{eqnarray}$$
where
$K$
is an integration constant and
$\unicode[STIX]{x1D6F4}_{X}$
is given by (3.18)–(3.20). Two orders in perturbation theory in
$\unicode[STIX]{x1D6EC}$
are necessary and sufficient to provide the sought-after bifurcation result.
We expand
$K$
and
$M$
in a perturbation series or equivalently
$K$
and the bifurcation parameter
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{M-M_{cr}}{M_{cr}}.\end{eqnarray}$$
This is effected in (4.4) by the substitution
$$\begin{eqnarray}M\rightarrow (\unicode[STIX]{x1D707}+1)M_{cr},\end{eqnarray}$$
(Kirkinis & Davis Reference Kirkinis and Davis2015, § d-v). Through this step we can derive known results as special cases to our formulation and obtain expressions which are independent of the velocity
$U_{0}$
and long wavelength
$\unicode[STIX]{x1D716}$
characteristic scales. Thus,
$$\begin{eqnarray}(K,\unicode[STIX]{x1D707})=(K_{0},\unicode[STIX]{x1D707}_{0})+\unicode[STIX]{x1D6FF}(K_{1},\unicode[STIX]{x1D707}_{1})+\unicode[STIX]{x1D6FF}^{2}(K_{2},\unicode[STIX]{x1D707}_{2})+O(\unicode[STIX]{x1D6FF}^{3}),\end{eqnarray}$$
where the positive small parameter
$\unicode[STIX]{x1D6FF}$
denotes distance from criticality. Expand
$H$
about a steady state
$$\begin{eqnarray}H=1+\unicode[STIX]{x1D6FF}(A_{c}\cos qX+A_{s}\sin qX)+\unicode[STIX]{x1D6FF}^{2}(B_{c}\cos 2qX+B_{s}\sin 2qX)+O(\unicode[STIX]{x1D6FF}^{3}),\end{eqnarray}$$
for real amplitudes
$A_{c},B_{c}$
etc. and wavenumber
$q$
associated with the primary bifurcation.
To leading order we find
$$\begin{eqnarray}\unicode[STIX]{x1D707}_{0}=\frac{q^{2}}{\bar{C}G},\end{eqnarray}$$
which is the standard result of linear stability analysis (Oron et al. Reference Oron, Davis and Bankoff1997). The first-order corrections vanish and, to second order, we obtain the amplitudes
$$\begin{eqnarray}B_{s}=A_{c}\,A_{s}\left[\,\frac{(\bar{C}\,G+q^{2})^{3}(B+1)}{36\bar{C}^{2}}\left(\frac{\unicode[STIX]{x1D6EC}}{B}\right)-\,\frac{\,(\bar{C}\,G+q^{2})(3\,B+1)}{6q^{2}(B+1)}\right],\end{eqnarray}$$
and
$$\begin{eqnarray}B_{c}=(A_{c}-As)(A_{c}+A_{s})\left[\frac{(\bar{C}\,G+q^{2})^{3}(B+1)}{72\,\bar{C}^{2}}\left(\frac{\unicode[STIX]{x1D6EC}}{B}\right)-\,\frac{(\bar{C}\,G+q^{2})(3\,B+1)}{12q^{2}(B+1)}\right],\end{eqnarray}$$
which recover equations (29) and (30) of Kirkinis & Davis (Reference Kirkinis and Davis2015) as a special case (by setting
$\unicode[STIX]{x1D6EC}=0$
in (4.10) and (4.11)). The Marangoni parameter
$\unicode[STIX]{x1D707}_{2}$
becomes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{2} & = & \displaystyle ({A_{c}}^{2}+{A_{s}}^{2})\left[\frac{(CG+q^{2})^{5}(CG+4\,q^{2})q^{2}(B+1)^{2}}{432GC^{5}}\left(\frac{\unicode[STIX]{x1D6EC}}{B}\right)^{2}\right.\nonumber\\ \displaystyle & & \displaystyle -\frac{(CG+q^{2})^{3}(3\,BCG+6\,Bq^{2}+CG-11\,q^{2})}{144C^{3}G}\left(\frac{\unicode[STIX]{x1D6EC}}{B}\right)\nonumber\\ \displaystyle & & \displaystyle \left.-\frac{(9\,B^{2}CG+45\,B^{2}q^{2}+6\,BCG+30\,Bq^{2}+CG+7\,q^{2})(CG+q^{2})}{24CG(B+1)^{2}q^{2}}\right].\end{eqnarray}$$
We recognize the last term as being the well-known correction (setting
$\unicode[STIX]{x1D70F}=0$
in (2.27) in Kirkinis & Davis (Reference Kirkinis and Davis2015)) signifying the subcritical behaviour of the bifurcation.
Second-order Marangoni number
$\unicode[STIX]{x1D707}_{2}$
versus the viscous heating (Brinkman) parameter
$\unicode[STIX]{x1D6EC}$
from (4.12). For
$\unicode[STIX]{x1D6EC}<\unicode[STIX]{x1D6EC}_{cr}$
the bifurcation is a subcritical pitchfork of revolution. It becomes supercritical for values above
$\unicode[STIX]{x1D6EC}_{cr}=8.89\times 10^{-4},1.015\times 10^{-2}$
and
$1.29\times 10^{-1}$
corresponding to
$B=0.01,0.1$
and
$1$
, respectively.
Equation (4.12) is the key result of this analysis. It provides the conditions upon which the bifurcation turns supercritical. This occurs for values
$\unicode[STIX]{x1D6EC}>\unicode[STIX]{x1D6EC}_{cr}$
, where
$\unicode[STIX]{x1D6EC}_{cr}$
is the positive critical viscous heating (Brinkman) number obtained by setting equation (4.12) equal to zero, solving the resulting quadratic equation for
$\unicode[STIX]{x1D6EC}$
and retaining the positive root. Figure 4 is the bifurcation diagram that arises from (4.12) for values
$G=3,\bar{C}=1/3$
at
$q=1$
giving rise to the threshold of instability at
$\unicode[STIX]{x1D6EC}_{cr}=8.89\times 10^{-4},1.015\times 10^{-2}$
and
$1.29\times 10^{-1}$
corresponding to values of
$B=0.01,0.1$
and
$1$
, respectively.
In a finite system of size
$L$
with periodic boundary conditions the wavenumbers become quantized,
$$\begin{eqnarray}q=(2\unicode[STIX]{x03C0}/L)n,\quad n=1,2,\ldots .\end{eqnarray}$$
When the critical threshold
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}_{cr}$
is exceeded, there is a band of unstable wavenumbers
$0<k<k_{cr}$
(see figure 3) that form and there will always be a value for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}>\unicode[STIX]{x0394}\unicode[STIX]{x1D717}_{cr}$
where the lowest of the quantized wavenumber in (4.13) will lie in the interval
$(0,k_{cr})$
and will thus lead to instability while the higher-order wavenumbers
$n\geqslant 2$
in (4.13) will not. Since the basic role of the nonlinearity is to saturate the exponential growth of the unstable mode (Cross & Hohenberg Reference Cross and Hohenberg1993, p. 870), viscous heating associated with these higher-order wavenumbers can lead to stabilization. The bifurcation criterion (4.12) leads to the following hierarchy of effectiveness of viscous heating in stabilizing an unstable liquid film:
(i) When
$\unicode[STIX]{x1D6EC}=0$
,
$\unicode[STIX]{x1D707}_{2}$
is always negative (see the last term in (4.12)) and the film is unstable; it proceeds to rupture and leads to substrate dewetting. This corresponds to a subcritical bifurcation.(ii) When
$0<\unicode[STIX]{x1D6EC}<\unicode[STIX]{x1D6EC}_{cr}$
, the film still ruptures albeit at a reduced rate. This behaviour corresponds to the dashed lines of figure 4 where
$\unicode[STIX]{x1D707}_{2}$
is negative but approaches its zero value. The bifurcation is still subcritical.(iii) When
$\unicode[STIX]{x1D6EC}>\unicode[STIX]{x1D6EC}_{cr}$
the film proceeds to a stable steady state. This behaviour corresponds to the continuous lines of figure 4 where
$\unicode[STIX]{x1D707}_{2}$
is now positive. The bifurcation is supercritical.
The inertia terms of the equations of motion and energy should be taken into account when higher-order corrections with respect to the wavelength
$\unicode[STIX]{x1D706}$
are considered. In the leading-order approximation (the lubrication approximation) that is tacitly assumed to hold in this paper, these corrections are subdominant and do not need to be considered.
Figure 5 displays the critical Brinkman number
$\unicode[STIX]{x1D6EC}_{cr}$
with respect to various Biot numbers
$B$
. It can be seen that the relation is linear. This behaviour is further verified by setting equation (4.12) equal to zero, solving for
$\unicode[STIX]{x1D6EC}$
and expanding with respect to
$B$
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{cr}(B\rightarrow 0)=\frac{3\bar{C}^{2}(\bar{C}G-11q^{2}+\sqrt{9G^{2}\bar{C}^{2}+66\bar{C}Gq^{2}+345q^{4}})}{(2\bar{C}G+8q^{2})(\bar{C}G+q^{2})^{2}q^{2}}B+O(B^{2}).\end{eqnarray}$$
Therefore, we expect the effect discussed in this article to be prominent for thin films (
$B\rightarrow 0$
) since small values of the Brinkman number
$\unicode[STIX]{x1D6EC}$
can be attained relatively easily in an experiment.
Numerical solutions of the thin film equation (2.27) with
$\unicode[STIX]{x1D6F4}_{X}=\unicode[STIX]{x1D6F4}_{X}^{(0)}+\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F4}_{X}^{(1)}+\unicode[STIX]{x1D6EC}^{2}\unicode[STIX]{x1D6F4}_{X}^{(2)}$
from (3.22) and (3.23) with
$\unicode[STIX]{x1D6EC}>\unicode[STIX]{x1D6EC}_{cr}$
determined by the bifurcation criterion equation (4.12). (a) Liquid–gas interface towards a stable steady state. Thinning of the film due to the thermocapillary effect ceases. (b) Velocity field at time
$t_{end}$
: two counter-rotating vortices lie perimetrically to the depression effectively reducing the wave-length of disturbances by one half. (c) Streamlines at time
$t_{end}$
.
We solve numerically equation (2.27)
$$\begin{eqnarray}H_{T}+[{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6F4}_{X}H^{2}-{\textstyle \frac{1}{3}}\bar{P}_{X}H^{3}]_{X}=0,\end{eqnarray}$$
with periodic boundary conditions in a spatial domain of length corresponding to the most unstable wavenumber
$k_{cr}/\sqrt{2}$
from (4.2) and initial condition
$H(X,0)=1+0.1\cos (k_{cr}X/\sqrt{2})$
. Here,
$\bar{P}_{X}=GH_{X}-\bar{C}^{-1}H_{XXX}$
and
$\unicode[STIX]{x1D6F4}_{X}$
is given by (3.18)–(3.20). We have performed standard pseudospectral Fourier computations in conjunction with one of Matlab’s variable time-stepping routines (Shampine & Reichelt Reference Shampine and Reichelt1997). Due to the higher-order derivatives,
$512$
or
$1024$
nodes were required to provide the sought-after spatial resolution. Following Kerchman & Frenkel (Reference Kerchman and Frenkel1994) we performed appropriate dealiasing. This is done by setting the Fourier coefficients of
$H$
and its derivatives corresponding to those high wavenumbers, equal to zero at each time step and prior to calculating the nonlinear terms in physical space.
For values of the dimensionless variables
$G=3,\bar{C}=1/3$
,
$B=1$
and
$M=16$
the time variation of the liquid–gas interface
$H(X,T)$
for
$\unicode[STIX]{x1D6EC}=0$
and
$\unicode[STIX]{x1D6EC}=0.5$
is shown in figures 2(a) and 6(a), respectively. The former shows the film’s evolution towards rupture and the latter its approach to a steady state. The effect of viscous dissipation is to halt the declining evolution of liquid–gas profiles and provide stabilization to the thin liquid film. This is because viscous heating causes heat-transfer inhibition requiring higher values of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D717}$
or
$M$
to bring about instability (cf. figure 6
b,c). Two subsurface vortices persist and effectively reduce the wavelength of disturbances to one half of its original value.
4.2.2 Expansion with respect to Biot number
$B$
Exactly the same bifurcation criterion (4.12) can be obtained by employing the surface tension gradient expansion with respect to the Biot number
$B$
in (3.21). Again, equation (2.27) can be integrated once to give
$$\begin{eqnarray}{\textstyle \frac{1}{2}}H^{2}(B\unicode[STIX]{x1D6F4}_{X}^{(1)}+B^{2}\unicode[STIX]{x1D6F4}_{X}^{(2)}+B^{3}\unicode[STIX]{x1D6F4}_{X}^{(3)}+B^{4}\unicode[STIX]{x1D6F4}_{X}^{(4)})-{\textstyle \frac{1}{3}}H^{3}\bar{P}_{X}=K,\end{eqnarray}$$
where
$K$
is an integration constant. In this formulation four orders in perturbation theory are required to obtain (4.12).
5 Discussion
In this article we presented a theoretical analysis of viscous heating effects, not considered previously in thin liquid-film studies. This effect is expected to find applications in microfluidic devices where thermocapillarity is used to transport small objects. With respect to the thermocapillary effect we showed that viscous heating can delay the rupture and even heal an unstable thin liquid film. Viscous heating is present during the motion of all liquids and its importance depends on the liquid properties. The bifurcation criterion we have derived in (4.12) is a complex relation linking the viscous heating (Brinkman) parameter
$\unicode[STIX]{x1D6EC}$
with all remaining dimensionless parameters of the problem. The main idea behind viscous heating can also be found for example in Landau & Lifshitz (Reference Landau and Lifshitz1987, § 55 and Problem 1 therein).
Two extensions of the present work may provide a better understanding of this effect. First, by considering the thermal conductivity of the substrate one expects that the temperature distribution in the liquid will contain a thermal resistance due to the conduction taking place in the solid (Oron et al. Reference Oron, Davis and Bankoff1997, § II.I). Second, considering the thermal properties of the gas phase, in a two layer system, see for instance the two layer experiments of the thermocapillary effect in thin silicon oil films (VanHook et al. Reference VanHook, Schatz, Swift, McCormick and Swinney1997). In both cases however a new formulation of the thermocapillary problem would be required and this is beyond the scope of the present article.
In this article we considered viscosity to be constant in the whole body of the liquid. This assumption, also adopted in the classical thermocapillarity literature (Oron et al. Reference Oron, Davis and Bankoff1997; Davis Reference Davis, Batchelor, Moffatt and Worster2002), provided a single valued and well defined velocity field (2.25) through which a certain analytical progress could be made. Although this is beyond the scope of the present article, it would be possible to also proceed with a temperature-dependent viscosity which, in certain cases, could lead to more realistic results (Johns & Narayanan Reference Johns and Narayanan1997).
When the liquid is endowed with rotational degrees of freedom (for example a magnetic liquid (Rinaldi Reference Rinaldi2002)) there are additional viscous heating terms in the energy equation that could bring about stabilization of an unstable thermocapillary film. In appendix B we therefore include such an extension of the discussion developed in this article. The mechanisms at play however in the case of a magnetic liquid, can also be applied to the motion of a magnetic liquid in a channel within a biological context. This is also further developed and highlighted in appendix B.
Acknowledgements
The authors are grateful to S. H. Davis, A. D. Gilbert and P. Yatsyshin for helpful discussions. This work was supported by the US Department of Energy Office of Science, Basic Energy Sciences under award no. DE-FG02-07ER46452.
Appendix A. Heat-transfer inhibition due to viscous heating
When
$\unicode[STIX]{x1D6EC}\equiv 0$
the energy equation (3.10) along with the boundary condition at the free surface (3.11) leads to a uniform (in
$Z$
) temperature gradient
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{Z}(X,Z)=-B\unicode[STIX]{x1D6E9}(H),\quad \unicode[STIX]{x1D6EC}\equiv 0.\end{eqnarray}$$
When
$\unicode[STIX]{x1D6EC}\neq 0$
the same two equations lead to the following temperature gradient field after some rearrangement
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{Z}(X,Z)={\textstyle \frac{1}{3}}\unicode[STIX]{x1D6EC}(H-Z)[\unicode[STIX]{x1D6F4}_{X}^{2}+\unicode[STIX]{x1D6F4}_{X}\bar{P}_{X}(H-Z)+\bar{P}_{X}^{2}(H-Z)^{2}]-B\unicode[STIX]{x1D6E9}(H,\unicode[STIX]{x1D6F4}_{X},\bar{P}_{X}).\end{eqnarray}$$
The first term on the right-hand side of (A 2) is positive. In the limit of small
$\unicode[STIX]{x1D6EC}$
the last term reduces to the right-hand side of (A 1). Thus, the temperature gradient reduces with elevation leading to a commensurate heat transfer inhibition.
Appendix B. Viscous heating in the presence of rotational degrees of freedom
The purpose of this appendix is to provide additional viscous heating mechanisms that can lead to stabilization of an unstable thermocapillary liquid endowed with rotational degrees of freedom (such as a magnetic liquid, (Rinaldi Reference Rinaldi2002)).
In the presence of rotational degrees of freedom there are dissipation processes associated with the imbalance of angular velocities and couple stresses, so (1.1) is replaced by (Aero, Bulygin & Kuvshinskii Reference Aero, Bulygin and Kuvshinskii1965)
$$\begin{eqnarray}\unicode[STIX]{x1D70C}c_{p}(\unicode[STIX]{x1D717}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\text{grad}\unicode[STIX]{x1D717})=k_{th}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}+\unicode[STIX]{x1D70F}_{ik}^{\prime }\left[\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D716}_{lki}\unicode[STIX]{x1D714}_{l}\right]+\unicode[STIX]{x1D60A}_{ik}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x2202}x_{k}},\end{eqnarray}$$
with respect to the viscous stress tensor
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ik}^{\prime }=\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}\right)+\unicode[STIX]{x1D701}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}-\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}-2\unicode[STIX]{x1D716}_{kil}\unicode[STIX]{x1D714}_{l}\right),\quad i,k,l=1,2,3,\end{eqnarray}$$
and couple stress tensor
$$\begin{eqnarray}\unicode[STIX]{x1D60A}_{ik}^{\prime }=\unicode[STIX]{x1D702}^{\prime }\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{k}}{\unicode[STIX]{x2202}x_{i}}\right)+\unicode[STIX]{x1D701}^{\prime }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{j}}{\unicode[STIX]{x2202}x_{j}}\unicode[STIX]{x1D6FF}_{ik},\quad i,k,j=1,2,3,\end{eqnarray}$$
where summation is implied on repeated indices and
$\unicode[STIX]{x1D714}$
is the suspended magnetic particles angular velocity distribution and
$\unicode[STIX]{x1D701}$
,
$\unicode[STIX]{x1D702}^{\prime }$
and
$\unicode[STIX]{x1D701}^{\prime }$
are constitutive parameters (Rinaldi Reference Rinaldi2002). After some rearrangement and symmetrization equation (B 1) becomes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}c_{p}(\unicode[STIX]{x1D717}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\text{grad}\unicode[STIX]{x1D717}) & = & \displaystyle k_{th}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}+\frac{\unicode[STIX]{x1D702}}{2}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}\right)^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D701}}{2}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}-\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{i}}\right)^{2}+\frac{\unicode[STIX]{x1D702}^{\prime }}{2}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{k}}{\unicode[STIX]{x2202}x_{i}}\right)^{2}\!.\end{eqnarray}$$
In the long wavelength approximation equation (B 4) reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{ZZ}+(\unicode[STIX]{x1D6EC}_{a}+\unicode[STIX]{x1D6EC}_{s})(U_{Z})^{2}+\unicode[STIX]{x1D6EC}_{r}(\unicode[STIX]{x1D6FA}_{Z})^{2}=0,\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FA}$
is dimensionless suspended particles angular velocity distribution, cf. Kirkinis (Reference Kirkinis2017). To (B 5) one can add the boundary conditions (3.11) for the case of a thermocapillary thin liquid film or alternative boundary conditions in the case of a channel in a biological context (Rinaldi Reference Rinaldi2002). The dimensionless coefficients (Brinkman numbers)
$\unicode[STIX]{x1D6EC}_{s}$
(
$s$
= symmetric),
$\unicode[STIX]{x1D6EC}_{a}$
(
$a$
= antisymmetric) and
$\unicode[STIX]{x1D6EC}_{r}$
(
$r$
= rotational) are each proportional to their constitutive counterparts
$\unicode[STIX]{x1D702},\unicode[STIX]{x1D701}$
and
$\unicode[STIX]{x1D702}^{\prime }$
, respectively,
$$\begin{eqnarray}(\unicode[STIX]{x1D6EC}_{s}+\unicode[STIX]{x1D6EC}_{a},\unicode[STIX]{x1D6EC}_{r})=\left(\frac{(\unicode[STIX]{x1D702}+\unicode[STIX]{x1D701})U_{0}^{2}}{k_{th}\unicode[STIX]{x0394}\unicode[STIX]{x1D717}},\frac{\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D6FA}_{0}^{2}}{k_{th}\unicode[STIX]{x0394}\unicode[STIX]{x1D717}}\right).\end{eqnarray}$$
This construction therefore, provides further corrections to the viscous heating problem considered in this article. Even further corrections are provided by the Maxwell stress tensor (Landau & Lifshitz Reference Landau and Lifshitz1960). For the thermocapillary problem, equation (B 5) replaces equation (3.10) and leads to temperature and surface tension gradient corrections emanating from the rotational degrees of freedom.
The aforementioned mechanisms can also be employed in channels of biological significance. A growing body of experimental evidence suggests that a suspension of magnetic nanoparticles in the presence of an alternating magnetic field leads to a localized temperature increase. For instance, Huang et al. (Reference Huang, Delikanli, Zeng, Ferkey and Pralle2010) grafting fluorophores to nanoparticles and others free in solution, observed that the magnetic field induced a particle surface temperature increase of
$20\,^{\circ }\text{C}$
by measuring the fluorescence intensity of the particle-bound fluorophore. The fluorophores let free in the solution did not show any temperature variation. Similar conclusions were reached by Polo-Corrales & Rinaldi (Reference Polo-Corrales and Rinaldi2012), who used iron oxide nanoparticles coated with the fluorescent thermoresponsive polymer consisting of poly(n-isopropylacrylamide) co-polymerized with a fluorescent benzofurazan monomer.
A theoretical explanation of the above observations may be related to viscous heating effects that arise due to the rotation of particles suspended in a classical liquid. We thus include here such a discussion. Consider a body immersed in such a magnetic liquid. Then, since temperature, velocity and angular velocity vary significantly over distances of order
$\ell$
(the size of the body), equation (B 5) leads to the estimate
$k_{th}\unicode[STIX]{x0394}\unicode[STIX]{x1D717}/\ell ^{2}\sim (\unicode[STIX]{x1D702}+\unicode[STIX]{x1D701})U_{0}^{2}/\ell ^{2}+\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D6FA}_{0}^{2}/\ell ^{2}$
where
$\unicode[STIX]{x1D6FA}_{0}$
is a characteristic scale for the particle angular velocity. Thus,
$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D717}\sim [(\unicode[STIX]{x1D702}+\unicode[STIX]{x1D701})U_{0}^{2}+\unicode[STIX]{x1D702}^{\prime }\unicode[STIX]{x1D6FA}_{0}^{2}]/k_{th}.\end{eqnarray}$$
This temperature increase arises because the liquid brought to rest at the surface of the body is heated by three separate mechanisms: the well-known viscous friction due to the rate-of-strain tensor (proportional to the viscosity
$\unicode[STIX]{x1D702}$
, (Landau & Lifshitz Reference Landau and Lifshitz1987, § 55)), spin friction, due to the antisymmetric part of the velocity gradient multiplied by
$\unicode[STIX]{x1D701}$
, and couple friction, due to the symmetric part of the nanoparticle angular velocity gradient, multiplied by
$\unicode[STIX]{x1D702}^{\prime }$
;
$U$
and
$\unicode[STIX]{x1D6FA}$
will depend on the applied magnetic field.
