Hopf instability of a Rayleigh–Taylor unstable thin film heated from the gas side

A thin liquid film located on the underside of a horizontal solid substrate can be stabilized by the Marangoni effect if the liquid is heated at its free surface. Applying long-wave approximation and projecting the velocity and temperature fields onto a basis of low-order polynomials, we derive a dimension-reduced set of three coupled evolution equations where nonlinearities of both the Navier–Stokes and the heat equation are included. We find that in a certain range of fluid parameters and layer depth, the first bifurcation from the motionless state is oscillatory which sets in with a finite but small wave number. The oscillatory branch is determined using a linear stability analysis of the long-wave model, but also by solving the linearized original hydrodynamic equations. Finally, numerical solutions of the reduced nonlinear model equations in three spatial dimensions are presented.


Introduction
The full mathematical description of the dynamics of nonisothermal incompressible fluids is well known for a long time. The fluid velocity field is described by the Navier-Stokes equations, the temperature field is governed by the heat equation, and the location of the interface and its spatio-temporal evolution are determined by the kinematic boundary condition. All these equations are amended with suitable boundary conditions. As a result, a complex nonlinear dynamic boundary-value problem emerges and even nowadays, at the age of supercomputers, its further treatment, especially in three spatial dimensions, presents a great challenge.
On the other hand, a direct numerical solution of the set of governing equations mentioned above, may be considered merely as another experiment. To achieve a deeper insight into the physics behind pattern formation, other methods have been devised, for an overview, see [1]. Very often one of the three spatial dimensions is singled out due to the geometry of the system. A good example for this case provide thin viscous films on a solid substrate. Here, the behavior of the solutions IMA10 -Interfacial Fluid Dynamics and Processes. Guest editors: Rodica Borcia, Sebastian Popescu, Ion Dan Borcia. a e-mail: bestehorn@b-tu.de (corresponding author) b e-mail: meroron@technion.ac.il in the vertical direction (z ) is very different from that in the other two horizontal ones, and the 3D fields of dependent variables can be projected as a good approximation onto rather simple low-order polynomials in z . This approach is successfully applied in the longwave theory, where upon this simplification, 2D thinfilm equations describing 3D fluid flow in the system arise [2]. Normally in these theories, vertical velocities are small and inertia is neglected, referred to as the zero Reynolds number (Re) limit. A similar approximation for the heat equation (zero Péclet number (Pe) limit) results in a linear temperature profile and in a linear coupling between the local film thickness and the local temperature gradient. Using these approximations, both long-wave Marangoni and Rayleigh-Taylor instabilities as well as their combination were successfully described [3].
In the present paper, we shall develop a reduced description where the nonlinear terms in both the Navier-Stokes and the heat equation are kept, i.e., where Re and Pe are considered to be both non-zero. The reduction is accomplished by projection onto distinct 2nd-order polynomials that satisfy the boundary conditions exactly, in the tradition of the Galerkin method. This procedure is also known as the weightedresidual integral boundary layer (WRIBL) technique or the Karman-Pohlhausen method for both isothermal [4][5][6] and nonisothermal problems [7,8]. It has been recently found that the Faraday instability can be described in this way also in its nonlinear regime [9,10]. Projecting the temperature field onto polynomials, the method was then recently extended to the nonisothermal case [8] to study the impact of the Marangoni effect on the Rayleigh-Taylor instability.
It is a purpose of the present paper to show that the extended model with non-zero Pe is able to describe an oscillatory first instability in the case of a fluid located on the underside of a solid planar substrate stabilized by the Marangoni effect originating from heating at the gas side. This instability has a finite band of growing wave numbers that is bounded also from below away from zero, contrary to the monotonic long-wave instability found for Re = P e = 0. The reduction to a low-dimensional system has the advantage that both the Hopf frequency and the critical Marangoni number can be derived analytically. A codimension-two point is determined where the oscillatory and the monotonic instabilities reverse their dominance. The results of the linearized reduced system are compared and checked with computations based on the linearized long-wave hydrodynamic basic equations and a good agreement is achieved. Finally, we present nonlinear simulations of the reduced system in 3D, showing persistent oscillating waves of the free surface over a long time terminated by rupture. To the best of our knowledge, this oscillatory instability was not reported before.

Governing equations
We consider a three-dimensional thin film of an incompressible fluid, see Fig. 1, with the deformable interface separating between the film and the ambient gas phase deposited on the underside of a horizontal, solid substrate in the gravity field g. The reference frame used in what follows is Cartesian with the x and y axes located in the plane of the substrate, whereas the zaxis is normal to the latter and points into the film. The location of the interface is a priori unknown and given by z = h(x, y, t).
In the long-wave approximation (LWA), the governing equations are written as [9] Re Here, υ H = (v x , v y ) denotes the horizontal component of the velocity field vector υ, T is the temperature field, the subscripts to ∂ denote the variables with respect to which partial differentiation is taken, and ∇ 2 = (∂ x , ∂ y ). All variables were made dimensionless by applying the following scaling: where τ = d/U 0 , U 0 stands for a certain characteristic velocity, d for the mean depth of the film, andT e is the constant environment temperature (variables with tilde bear a dimension). The temperature difference represents the difference between the uniformly cooled substrate and the warmer environment. Reynolds and Péclet numbers are defined via with the Prandtl number

Boundary conditions
In LWA, the boundary conditions for the velocity field are Here, we have introduced the new 2D variable as the temperature along the free surface. The Marangoni number Ma is defined as where γ T = −dσ/dT and σ is the surface tension between the liquid and the ambient gas. In what follows, the Marangoni number Ma is negative.
For the temperature, we find with the Biot number B where q is the rate of heat transfer from the liquid to the ambient gas at the interface and α is the thermal conductivity of the liquid. The location of the free surface is determined by the kinematic boundary condition at z = h. Note that in convection problems, the Marangoni number is sometimes defined as M = Ma P e B/(1 + B).

Hydrostatic pressure
In LWA, the pressure consists of components, namely, the hydrostatic and the capillary ones [9], so its gradient reads with the dimensionless groups Galileo and inverse Crispation number respectively. Thus, the six groups Re, G, Γ, Ma, Pe, B define our problem. Since U 0 can be chosen arbitrarily, one of the first five parameters can be put to one, in what follows we shall use viscous scaling U 0 = ν/d, leading to Re = 1 and P e = P r. g.

Polynomials and 2D variables
In addition to the surface temperature (5), we introduce next the 2D vector field which represents the horizontal flow rate. We now use the following representation for the horizontal velocity field with the polynomials Note that Thus, the representation (13) satisfies the boundary conditions (4) for arbitrary q and θ. For the temperature, we take also a quadratic dependence on z and write Since (15) fulfills (5) and the boundary conditions (7) for arbitrary θ. The boundary condition (9) can be reformulated with (12) as

The dimension-reduced model
To derive a set of evolution equations for h, q and θ, we insert Eqs. (13,15) into Eqs. (1a,1b), multiply Eqs. (1a) and (1b) by the weight functions W 1 (z;h) and W 2 (z;h), respectively, and integrate them over 0 ≤ z ≤ h. After some algebra, a set of 2D evolution equations is obtained with the numerical values of the coefficients depending on the choice of W i . Choosing we obtain the following set complemented by Eq. (17): which can be compared directly to the expansions given by Trevelyan et al. [11] and Sterman-Cohen and Oron [8]. As mentioned above, the coefficients depend on the choice of the weight function W 2 and are slightly different in all the three works. Taking, for instance, W 2 = z, we find for the case B 1 Note that in both limits Eqs. (19) and (20), we retain B in the right-hand side of the equations. Then it is possible to approximate the solution if the convective terms are absent in the heat equation what formally corresponds to κ → ∞ or P e → 0 with the well-known result, e.g., [2,8] θ = 1 1 + Bh .

Linear stability analysis
It is easy to identify the base state for Eqs. (18) in the form of the quiescent conductive state with a flat interface that reads Due to isotropy of the laterally unbounded system in the x-y plane, the linear analysis may be performed in just one horizontal spatial dimension, say x . Substituting into Eq. (18), results after linearization in the realvalued linear set of algebraic equations For the sake of simplicity, we assume in this section from here on a physically relevant case of a small Biot number B 1, according to the limits of Eq. (19). The solvability condition for Eqs. (22) represents a cubic polynomial equation At the threshold of a monotonic instability λ = 0, leading thus to a 0 = 0, and the system is monotonically unstable if a 0 < 0, leading to taking place at k = 0. For Marangoni numbers larger than Ma M c.min the stabilizing Marangoni effect becomes too small, and the Rayleigh-Taylor instability emerges. This is the well-known condition for the long-wave Marangoni instability in a thin film that is not affected by inertia and thermal convection (Re = P e = 0).
For the onset of an oscillatory (Hopf) instability, λ = iω c with a real ω c , we obtain ω c = a 0 /a 2 and, since a 2 > 0, one must have a 0 > 0. The latter is possible for Ma < Ma M c , i.e. where the monotonic instability does not emerge, and, therefore, occurs then as a first instability increasing Ma if Pe exceeds a certain critical value. The critical Marangoni number Ma os c for the Hopf instability is found from Eq. (23) equating a 0 Re = a 1 a 2 , ω 2 c = a 0 /a 2 and turns out as  The critical value M os c is a convex function of k and has a minimum at (Fig. 2 which is always smaller than the expression in Eq. (25), since G is negative. However for this case, k c → 0 and the wave length of the critical modes would become very long and the Hopf frequency (27b) tends to zero.
For P e > P e CD2 , the instability of the motionless state is oscillatory. It is interesting to note that qualitatively similar results Eqs. (27-30) can be derived if inertia is neglected, i.e., Re 1, at least for not too small Pe. In the case of Re = 0, the eigenvalue problem (22) reduces to a second-order problem, which is easily solved and the critical thresholds for both monotonic and oscillatory instabilities are obtained by the same equations (26-29), but now with Re = 0. This is probably due to the fact that the fluid velocities for this case stay rather small and inertia effects play no significant role.
To demonstrate that our results are not an artifact of the projection method used to obtain the set of Eqs. (18), we next solve the linearized long-wave set of Navier-Stokes and heat equations (1) by a standard finite-difference method for the same parameters and found a phase diagram presented in the right frame of Fig. 2. It clearly shows that for P e = 25, the first instability of the quiescent state is oscillatory.

Numerical results
Finally, we present our numerical results based on the reduced system (18). An FTCS (forward-time centralspace) [12] finite-difference method has been applied For Ma = −300, we are clearly in the monotonic region of Fig. 2 and after a certain time of coarsening, a few single stalactites survive and the code diverges, as a probable sign of film rupture, see Fig. 3, left frame.
The situation is completely different in the oscillatory regime for Ma = −400, Fig. 3, right frame. Now, the amplitudes of the emerging pattern remain small for a long time before more or less suddenly film rupture takes place, see Fig. 4. During this phase, the film surface has the form of mixed traveling and standing waves, nonlinear selection processes do not take place probably due to the rather small amplitudes. The emergence of the Hopf bifurcation is obvious from very early stages of the evolution.

Conclusions
We have shown that an oscillatory (wave) instability emerges via the Marangoni effect for a Rayleigh-Taylor unstable thin layer heated from the gas side. Increasing the (negative) Marangoni number, this instability may bifurcate first from the motionless flat film state for a certain range of fluid parameters and mean layer depth. For rather thin films and/or high fluid viscosity, this instability is obliterated by the monotonic long-wave instability. Inspection of the codimension-two condition (31) yields an upper bound for the kinematic viscosity of the liquid ν < 9 1024 which enables the emergence of the oscillatory instability. It is presented for a typical set of fluid parameters σ = 0.02 N/m, κ = 0.15 · 10 −6 m 2 /s, ρ = 920 kg/m 3 as a function of the mean depth in Fig. 5. It demonstrates that the layer must be rather thick to exhibit an oscillatory instability with an increase in the kinematic viscosity of the fluid. A similar instability behavior has been found based on the original hydrodynamic equations which justifies the projection method developed here. Nonlinear simulations of the set of reduced model equations terminate always with rupture (or dripping of some singular stalactites) in both the oscillatory and the monotonic flow regimes, probably due to the rather large mean film thickness of 0.5 mm used for our computations. Finally, the physical situation described in this paper in which the dynamics of a Rayleigh-Taylor unstable thin liquid film is heated at the gas side has been studied both theoretically and experimentally. In the limit of the inertialess long-wave theory, heating the film at its interface provides an extra-stabilization mechanism along with the surface tension which leads to the emergence of the monotonic instability and its subsequent saturation [13,14] in both 2D and 3D. The same was found experimentally in double-layered systems [15] and then reproduced theoretically [3].
The emergence of the oscillatory instability in this physical setting has been found for the first time, and a possibility of its saturation, which was not found here, may be explored in more detail in the future.
Funding Information Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Data availability statement
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.