Excitation of Surface Waves Due to Thermocapillary Effects on a Stably Stratified Fluid Layer

In chemical engineering applications, the operation of condensers and evaporators can be made more efficient by exploiting the transport properties of interfacial waves excited on the interface between a hot vapor overlying a colder liquid. Linear theory for the onset of instabilities due to heating a thin layer from above is computed for the Marangoni–Bénard problem. Symbolic computation in the long wave asymptotic limit shows three stationary, non-growing modes. Intersection of two decaying branches occurs at a crossover long wavelength; two other modes co-exist at the crossover point—propagating modes on nascent, shorter wavelength branches. The dispersion relation is then mapped numerically by Newton continuation methods. A neutral stability method is used to map the space of critical stability for a physically meaningful range of capillary, Prandtl, and Galileo numbers. The existence of a cut-off wavenumber for the long wave instability was verified. It was found that the effect of applying a no-slip lower boundary condition was to render all long waves stationary. This has the implication that any propagating modes, if they exist, must occur at finite wavelengths. The computation of 8000 different parameter sets shows that the group velocity always lies within 2 to 3 of the longwave phase velocity.


Introduction
Marangoni-Bénard convection still generates much interest in fluid and nonlinear dynamics due to its complexity.When the fluid layer is heated from below, convective instabilities can be driven by surface or buoyancy forces [1].The role of surface-tension gradients in inducing convective instability through Marangoni stresses at the air-liquid interface in a thin layer initially at rest, heated from below, was characterized in the seminal works of Sternling and Scriven [2] and Smith [3] and the role of surface deformation and surface tension gradients in the onset of patterned convection and oscillatory instability is reviewed in [4].
When the layer is heated from above, however, only overstability can be excited at sufficiently high Marangoni or Rayleigh (buoyancy) numbers.The most common chemical engineering context for a cold liquid heated from above is a condenser, which comprises a hot vapor that condenses over a colder liquid chilled from the solid support below.There are many different configurations for condensers, but all of them would condense faster if interfacial waves are excited on the interface.Thus, a stability theory that shows under what conditions self excitation occurs could better inform the design and operation of condensers.The scenario also describes an evaporator that is operated by contacting hot gas over the cold liquid, as long as there is no fluid motion imposed.Of course, imposing fluid motion adds additional complication but has been found to accelerate performance in novel contactors (evaporators, condensers and distillers) [5].
Classical experiments by Linde and coworkers [6] demonstrate a series of wave instabilities excited at high Marangoni numbers, although they could not discern whether the excited waves were surface or internal waves.It was later posited that the waves were surface manifestations of soliton solutions to a dissipative variation of the Korteweg-de Vries equation.Experiments have shown that solitary waves, excited and sustained by Marangoni stresses, undergo interactions and collisions that return the waves to the pre-collision celerities and shapes, experiencing at most a phase shift, but with either sense possible [7].Nepomnyashchy and Velarde [8] demonstrated via multiple scale perturbation methods the definitive derivation of the dissipative nonlinear evolution equation (termed the KdV-KSV equation) in a sufficiently thin layer that buoyancy effects are neglected.Their study assumed a stress-free lower boundary.The KdV-KSV theory predicts a critical Marangoni number M crit = 12, irrespective of capillary, Prandtl and gravity numbers.The interpretation of this result is simply that surface-tension gradients, if sufficiently strong, overcome viscous dissipation in the surface layer to start the fluid oscillations.These oscillations must then combine two modes-the gravity waves modified by Marangoni stresses that are normal to the surface (termed transverse waves) and elastic-like waves that are tangential to the surface (termed longitudinal waves).A study of solutocapillary Marangoni-induced interfacial waves is given in [9].
The purpose of this paper is to test the regime of validity of the KdV-KSV theory on two points-(i) the assumption of the multiple scale theory that there is a long wave cut-off for a wave-packet of unstable waves just above critical stability; (ii) the affect of the no-slip lower boundary on the coefficients of the linear terms of the KdV-KSV equation-by computing the full linear stability theory (LST) of the Bénard-Marangoni problem when heated from above.Since, in the basic state, the fluid is at rest and only the static pressure among the field variables has a vertical dependency, it is possible to formulate the analytic solution to the linearized equations and subsequently the secular equation for excitation of non-trivial solutions for the streamfunction, temperature, and surface displacement.The analysis presented is a precursor to the weakly nonlinear theory for the evolution of capillary gravity waves given in [10].The model equations are presented in Section 2. In Section 3, the linear stability of both the stress-free model and the no-slip models are formulated.In each case, the long wave asymptotic forms are computed for branches of the dispersion relation found.Results are discussed in Section 4 and the conclusions are drawn in Section 5.

Scaling
Consider a layer of viscous fluid with dynamic viscosity µ, density ρ, and thermal diffusivity κ, which is heated from above.The background no flow and no deformation state is perturbed, giving rise to gravity, capillary, and Marangoni forces.The relevant physical parameters are summarized below: , where σ 0 is the surface tension at the reference temperature T 0 .

•
Surface height z = d + h (x, y, t) , where d is the nominal height of the surface in the absence of deformation, and x and y are horizontal coordinates.t is time.
As the focus here is on surface forces, the buoyancy in the bulk will be neglected by assuming density constant.
We take diffusive scales as follows: , where σ ij is the stress tensor.The asterisks refer to dimensionless variables.Henceforth, the asterisks will be dropped and all calculations are given in dimensionless variables.Dimensional analysis yields four dimensionless groups: Marangoni number M = − γβd 2 µκ .As the layer is heated from above, β will be negative for simple fluids and the Marangoni number as defined above is an intrinsically positive quantity for the target situation.

Model Equations
The full governing equations are adapted from Davis and Homsy [11] by neglecting buoyancy.In the bulk, the velocity, temperature and pressure are constrained as: 1 Pr where The comma-subscript represents the index convention for partial differentiation with respect to x i .The boundary conditions are a rigid lower planar surface held at constant temperature at z = 0: The upper free surface is open and deformable at z = 1 + η (x, y, t): The differential geometry of the surface S(t) : z = 1 + η is given in terms of the element of arc length N(η), curvature K(η), and the normal n i and the tangent t i vectors: ( The base state is motionless with only hydrostatic pressure: η = 0.

Linearized System
The linearized equations about the motionless base state (4) for a two-dimensional disturbance can be written in terms of a streamfunction ψ, defined in the usual way for incompressible flow.In dimensionless form, the streamfunction-vorticity equation is given by: Subscripting by coordinates refers to partial differentiation by the respective coordinate.The pressure, which is needed for the normal stress boundary condition, can be found from The mathematical analysis leading to Equations ( 5) and ( 6) is outlined in Davis and Homsy [11].The linear stability analysis is a standard mathematical approach having linearized the equations about the base state.Heat transport couples temperature convection and diffusion.The linearized version is We apply seven boundary conditions at the upper and lower surfaces.
At z = 0, material surface (no penetration): no slip: fixed temperature: At z = 1, normal stress: tangential stress: material surface (kinematic condition): no heat flux: We presume separation of variables with a factor that is a wave of real wavenumber k and (complex) phase velocity c by normal mode expansion of the field variables: The wave system is thus reduced to a two-point boundary value problem of ordinary differential equations (ODEs) in P (z) , Ψ (z) , Θ (z) and H, with algebraic constraints: Bulk equations, Boundary conditions, at z = 0, The algebraic manipulation package Mathematica 5.2 (Wolfram, Champaign, IL, USA) was used to manipulate the governing equations for this system.As the bulk equations are linear ODEs with constant coefficients, the general solutions can be represented by linear combinations of complex exponentials given by the characteristic roots, ±k ±λ = k 2 − ikc Pr and ±Λ = The seven boundary conditions (19) and (20) reduce to a matrix equation in the seven unknowns α = [A, B, C, D, F, J, H] T , Zα = 0, where the matrix Z = Z stick is given by where .
Alternatively, with slip boundary conditions, where boundary condition ( 9) is replaced by the stress-free condition on the lower surface ψ zz = 0, we have the matrix Z = Z slip :

Slip Boundary Condition
Either we have α = 0, in which case no wave solution exists, or there is a non-trivial solution with a dispersion relation for the phase velocity as a function of wavenumber c (k) implied by the singularity of the matrix, ∆ = det (Z) = 0. Expanding the determinant ∆ algebraically using the permutation rule gives rise to 720 complex exponential terms upon eliminating H through application of the kinematic condition on the upper surface.Thus, the determination of the dispersion relation through solving ∆ = 0 symbolically is not a trivial undertaking.In contrast, the computation is readily tractable numerically, apart from a particular difficulty in finding the zeros of ∆.Numerical computations with fixed precision arithmetic are unlikely to return exactly zero.Therefore, a small number below a given threshold is typically taken to be zero.When the determinant of a matrix approaches zero, its condition number becomes large and thus the matrix becomes ill-conditioned, meaning that computations could involve high levels of numerical error.Consequently, the threshold level required to identify zeros of ∆ can be difficult to determine a priori.This motivates our desire to seek closed form symbolic approximations to the dispersion relation.
Such approximations usually start with the long wave limiting cases.Presuming that c is an analytic function of k allows us to write the Taylor's series expansion of ∆ (c, k): Truncating the series at n = N yields an approximate dispersion relation implicitly by requiring the approximate determinant to vanish.We note that, with the matrix Z slip as written, the first non-trivial contribution comes at n = 3/2, with having a quadruple root at c = 0 and real roots These are the long wave asymptotic limits that demonstrate that the Marangoni effect is additive to the gravity wave effect in contributing to the phase velocity.The capillary effect, however, is classically weaker by O(k 2 ) [12].
To isolate these modes, it is convenient to translate the phase velocity to the frame of reference of the leftward (or rightward) moving wave and consider ∆ δc + (G + M) Pr, k = 0. Since, by construction, ∆ is analytic in both k and c, it is useful to develop the k-power series of the c-truncation of ∆, i.e., with which is readily interpreted as a long wave instability occurring for a Marangoni number greater than the critical value M crit = 12 for long waves k << 1.As k = 0 is a critical wavenumber with neutral stability, it follows that for the problem to be well posed in the sense of Joseph and Saut [13], Im{c(k)} < 0 as k → ∞.Thus, there must be a cut-off wavenumber for the long wave instability, where Im{c(k cuto f f )} = 0 for k > 0 and higher wavenumbers decay, presumably due to the viscous and heat dissipation modes.Conveniently, the exponential growth rate with slightly supercritical M grows quadratically in k 2 (Marangoni pumping) and decays as k 4 , since the growth constant is −kIm{c}.As the odd order contributions in k are purely imaginary, it is possible that the long wave instability is limited to a wave packet, with higher wavenumbers than k cuto f f being damped.Figure 1 demonstrates this occurrence for slightly supercritical M = 12.1.It can be seen that the growth rate is positive for wavenumbers less than approximately 0.007, but that for larger wavenumbers, the growth rate is negative, which is indicative of decaying modes.An analytic form for the cutoff wavenumber can be computed from Equation (28), which also shows the supercritical nature of the long wave packet of instability.Ref. [8] assumed the occurrence of just this form of dispersion relation as the basis for their multiple scale perturbation theory.They assumed k cuto f f << 1 as the formal perturbation parameter.Given that k cuto f f = 0.00724 for the conditions in Figure 1, their intuition has been proved correct here.

No-Slip Boundary Conditions
Introducing the no-slip lower boundary makes a major structural change to the physics of the problem.It is well known that there are three fluid dynamical dissipative mechanisms that lead to the attenuation of surface waves.Bottom friction is the dominant mechanism wherever the layer depth is substantially less than the wavelength, so that the wave induces large horizontal motions near the bottom (see [12], Figure 55 for an illustration of this point).
Deep fluid waves, however, do not induce movement near the bottom, and thus no frictional dissipation.Internal dissipation by viscous stresses acting throughout the wave cause attentuation.
Surface dissipation is associated with departures of the surface from its equilibrium value, described, for example, by Lucassen (1968) in the case of doping of the surface with a monolayer of surfactant leading to an elastic dissipative mechanism.
In addition to these mechanisms, the configuration under study has internal dissipation from the thermal conduction and possible surface dissipation from the Marangoni effect, which can play the role of a Lucassen-like tangential surface stress.
Lighthill [12] estimated the proportional energy loss per period, due solely to the bottom friction, as 2π d where Ω = c/k is the wave frequency.The square bracket factor is the ratio of the thickness of the bottom viscous boundary layer induced by the wave to the depth.The curly bracket factor corrects for finite wavenumber-it is unity for infinitely long waves and zero for infinitely short ripples.
A complementary analysis for internal dissipation leads to the opposite preferences.Infinitely long waves are unaffected by internal dissipation; infinitely short waves are massively damped by internal dissipation.No doubt that consideration of internal dissipation effects only influenced the search by [8] for long wave excitation of surface solitary waves.
A priori, this estimate would suggest that excitation of solitary disturbances cannot occur for either long waves or short waves due to the dominance of the attenuation by bottom friction and internal dissipative mechanisms, respectively.Candidate wavenumbers for solitary wave excitation should be intermediate wavenumbers where the Marangoni effect induces sufficient disturbance energy to overcome internal dissipation and bottom friction.
To investigate this hypothesis, the modal structure of the no-slip problem should be clarified.The principal reason is that the major branches can be described with analytic approximations, leading to better understanding of the structure, and highlighting changes that are possible at intermediate wavenumbers.

No-Slip Boundary Condition: Modal Structure
It is now convenient to express separation of variables for the no-slip problem for each normal mode in terms of exponential factors appropriate for growing standing waves (ω real): The general solution can be found as: where The modal structure can be examined by considering the long wave limit lim k→0 ∆ (ω, k; Pr, M, K, G) = 0 (32) implicitly defines ω (k → 0; Pr, M, K, G).This limit can be determined symbolically from: }}.We find that the neutrally stable root ω → 0 has multiplicity of at least seven.The root ω → −π 2 4 is associated with exponential decay.Since the scaling for time is diffusive, the exponential decay constant is a factor of the thermal time scale.Thus, this mode is termed the thermal mode.As the root ω → −π 2 Pr 4 decays according to the viscous time scale, we label this the viscous mode.
A key property of all of these modes is their stationarity.The long wave slip modes comprise a pair of propagating simple waves.The no-slip bottom boundary condition has the effect of rendering all long waves stationary.Thus, if any propagating modes exist, they must occur at finite wavelengths.This suggests investigating the k-dependence for long waves branching from the infinitely long wave modes identified above.

k-Space Continuation and Modal Character
Identifying the critical surface in {M, K, Pr, G} parameter space for the co-existence of the pair of intermediate wavenumber propagating modes and the long wave stationary modes (long wave thermal mode and long wave neutral temperature mode) is essential in mapping the parameter space.Continuation methods based on Newton iteration must search for either purely real ω (long wave stationary modes) or genuinely complex ω (intermediate wavenumber propagating modes).Figure 2 shows this critical curve k − M and ω − M for four mode co-existence at fixed parametric values {K = 0.0001, Pr = 1.001,G = 1} found by numerically solving for ∆ = 0 and ∆ ω = 0 simultaneously by Newton's method.For M ∈ {500, 1500}, it can be seen that the wavenumber of the crossover point is a decreasing function, whilst the growth rate is an increasing function.The four modes that co-exist along this critical curve are the two long wave stationary modes and two intermediate modes propagating from opposing directions.Figure 3 shows an analogous k − Pr curve for the crossover point from stationary long waves to propagating steady intermediate waves.A Prandtl number of 1 is common for gases and a Prandtl number of 10 is common for liquids.The plots in Figure 3 show how the behavior of the fluids changes from a near gaseous state to a near liquid state.
Once the critical surface dividing the wavenumber space into the two different regions of wave character has been identified, it is a simple matter to use parameter space continuation methods in k to find the remainder of the dispersion relation.Figure 4 summarizes the dispersion relation by plotting maximum growth rate and the wavenumber at which it is achieved for a given Marangoni number M at the values {K = 0.0001, Pr = 1.001,G = 1}.There is jaggedness in the graph in Figure 4a as the wavenumber of maximum growth can switch modes.The normalized group velocity is also plotted as a function of the Marangoni number.It is clear that neutral stability is found by ramping up the Marangoni number to M ∼ 1100 for the other parameters fixed at these values.26) for which ∆ = 0 and Re {ω max } > Re {ω} for all k at the parametric values {K = 0.0001, Pr = 1.001,G = 1} and M ∈ {500, 1250}.These curves identify the most dangerous mode (wavelength with most rapid growth or slowest decay).Clearly, there is neutral stability at M ∼ 1100.The jaggedness of (a) reflects the granularity of the discretization in k-M-continuation.

Critical Parameters via the Neutral Stability Method
Identifying the critical Marangoni number and associated wavenumber and frequency by continuation in k and M for fixed {K, Pr, G} is computationally expensive if parameter space is to be mapped.Once a single critical parameter set is known, continuation along the neutral stability surface should be possible.Ref. [14] recognized that the Marangoni number only arises linearly in Equation (20), so that it is possible to solve the tangential stress boundary condition for M and eliminate all quantities α, thus finding the neutral stability curve by the condition that at arbitrary k and with Re {ω} = 0, Im {M} = 0 is imposed by adjusting s = Im {ω}.Takashima [15] also used this technique.Neither study, however, identified the long wave stationary branches put forth here.The critical parameters {M crit , k crit , s crit } are then found by minimizing M over k.In this paper, the relation ∆ = 0 was solved for M, and neutral stability was found by using Newton's method to adjust s to achieve Im {M} = 0. Use of numerical root finding methods for s such that Im {M} = 0 are faster than those that seek values of s such that ∆ = 0.
An algorithm for parameter space continuation in {K, Pr, G} for {M crit , k crit , s crit } by traversing the neutral stability surface was developed.The essential steps are to find the M-k and s-k neutral curves by computing a neighborhood in k surrounding an estimated critical point {M, k, s} from adjacent values of {K, Pr, G}.Subsequently, Newton's method is used to find the minimum of the M-k neutral curve.This is illustrated in Figure 5 for the parameters Pr = 1.001,K = 0.0001 and G = 1.The crossover wavenumber M − k neutral stability curve, i.e., that for which (ω) = 0, is plotted, along with a curve showing the functional dependence of the frequency s = Im(ω) on the wavenumber k.A space of twenty values of each of {K, Pr, G} was mapped by this method.Table 1 exhibits a sampling of the database built up from Newton continuation in {K, Pr, G}.The range of parameters was selected to represent some physically realizable liquids (small capillary number and O(1) Prandtl number) and layer depths (microgravity to terrestrial gravity).From Table 1, it is clear that values of M crit ∼ 10 4 -10 5 are achievable under combinations of small capillary number and thin layers in liquids.Critical wavenumbers k crit are typically of intermediate scale, neither long (<0.1, say) nor short (>10, say), but usually about k crit ∼ 0.25.Figure 4 is suggestive that the group velocity s/k of the critical mode might actually be well predicted by the long wave formula (20) for the phase velocity.The scatter plot in Figure 6 tests this hypothesis by crudely applying it to the entire mapped set of 8000 critical points.The critical group velocity, c 0 = (G + M)Pr, is plotted against the longwave phase velocity, s crit /k crit , for each of the 8000 critical parameter sets considered.The plot suggests that c 0 is the correct order of magnitude for the group velocity of the critical mode regardless of the parameter set in the range mapped.Furthermore, it is observed that the magnitude of the group velocity is constrained to be within 1  2 to 2 3 of the longwave phase velocity.This suggests that c 0 is the correct order of magnitude and a fair approximation for the group velocity of the critical mode regardless of the parameter set in the range mapped.Dotted lines have slope 2 3 and 1 2 .

Conclusions
An analytical study based on linear theory has been presented that considers the development of instabilities that arise when a thin layer of a stably stratified fluid is heated from above for the Marangoni-Bénard problem.Particular focus was directed towards the value of the critical Marangoni number and to the influence of the lower boundary condition.
The stress-free model modifies the predicted critical Marangoni number to depend on the Galileo number (gravity effects), but reduces to M crit = 12 under microgravity.It was not possible to find the cut-off wavenumber analytically as it requires solution of the dispersion relation to O(k 2 ) in small wavenumber k, but numerical solutions of the secular equation indicates, for all parameters tested, that k cuto f f < ∼0.05.This quantitative verification of the assumptions of KdV-KSV theory [8] under the stress-free boundary conditions leads one to conclude that any deficiencies of the theory must lie in the lower boundary condition.The analytic study of the long wave solutions to the secular equation for the no-slip model led to the surprising discovery that there are no infinitely long propagating modes.The three modes identified as k → 0 are a neutral pure temperature mode, a damped viscous (pure velocity) mode, and a damped thermal (pure temperature) mode.The temperature and viscous modes both decay with increasing k, but the thermal mode increases in growth rate with small k.The prediction of an intersection point is computed by Newton continuation of the nonlinear solution of the secular equation.At the intersection of the two modes, the crossover point, there is co-existence of two propagating modes with the two stationary long modes.At higher wavenumbers, only the propagating modes exist.These modes are mapped by parameter space continuation in wavenumber, and a parametric study by continuation in the physical parameter space is made.Further developments in this research area could be facilitated by full numerical simulations.

2 and λ 6 =
−λ 5 are the characteristic values of the general solution.

Figure 5 .
Figure 5. Magnitude of velocity field on z = 0 mm plane: (a) crossover wavenumber k for the M-k neutral stability curve (Re {ω} = 0) for the parameters Pr = 1.001,K = 0.0001, G = 1.M crit is found from the global minimum of the neutral curve and k crit is the associated wavenumber; (b) frequency s = Im {ω}.

Table 1 .
Critical values {M crit , k crit , s crit } for fixed {Pr, G, K} found by the neutral stability/minimization algorithm and parameter space continuation within the database.Samples were tested against the dispersion relation method and were found to be in agreement and global maximum of growth rate versus k.The table shows only selected values; the database is {Pr ∈ [1.,10.]}× K ∈ 10 −6 , 10 −2 × G ∈ 10 0 , 10 5 .There are twenty values for each parameter, spaced exponentially.