[ UPMC Univ Paris 06, Univ Paris-Sud, CNRS, UMR 7608, Lab FAST, Bat 502, Campus Univ, Orsay F-91405, France
Institut für Thermo- und Fluiddynamik, Technische Universität Ilmenau, Postfach 100565, 98684 Ilmenau, Germany
UPMC Univ Paris 06, CNRS, UMR 7190, IJLRA, 4 place Jussieu, Paris F-75005, France
?? and in revised form ??
?? and in revised form ??

The convective instability in a plane liquid layer with time-dependent temperature profile is investigated by means of a general method suitable for linear stability analysis of an unsteady basic flow. The method is based on a non-normal approach, and predicts the onset of instability, critical wavenumber and time. The method is applied to transient Rayleigh-Bénard-Marangoni convection due to cooling by evaporation. Numerical results as well as theoretical scalings for the critical parameters as function of the Biot number are presented for the limiting cases of purely buoyancy-driven and purely surface-tension-driven convection. Critical parameters from calculations are in good agreement with those from experiments on drying polymer solutions, where the surface cooling is induced by solvent evaporation.


eurm10 \checkfontmsam10 ]Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis. F. Doumenc, T. Boeck, B. Guerrier and M. Rossi] F.\nsD\lsO\lsU\lsM\lsE\lsN\lsC , \nsT.\nsB\lsO\lsE\lsC\lsK , \nsB.\nsG\lsU\lsE\lsR\lsR\lsI\lsE\lsR and M.\nsR\lsO\lsS\lsS\lsI 2009 \volume?? \pagerange??–??

1 Introduction

Thermally driven flows near liquid interfaces continue to be an active area of research since the first experimental studies by H. Bénard over 100 years ago. They are characterized by a multitude of interacting physical mechanisms and display a large variety of regular and complex flow patterns. Experimental and theoretical studies of such flows have stimulated and accompanied the development of the theory of pattern formation and nonlinear phenomena in general. Recent developments are summarized in a number of monographs and review papers (see for example Colinet et al. (2001); Nepomnyashchy et al. (2006); Bodenschatz et al. (2000)).

The driving forces in these thermally driven flows are surface-tension gradients (Marangoni forces) and buoyancy forces, which originate from the temperature dependency of surface tension and density, respectively. Both mechanisms are classically studied in convection taking place within a finite layer subjected to a steady temperature difference: for the buoyancy-driven case, this is the celebrated Rayleigh-Bénard convection formulated by Rayleigh (1916), for the surface-tension driven case this is the Bénard-Marangoni convection studied by Pearson (1958). In both cases, a certain critical temperature difference must be applied in order to sustain convective motion. The theoretical predictions for such critical temperature differences are in excellent agreement with experiments (Chandrasekhar (1961); Schatz et al. (1995)). Thermal convection can also appear when the external conditions are time-dependent, e.g. by a modulation or abrupt change of cooling or heating. In this case, the basic conductive temperature distribution is time dependent, and the stability problem becomes a non-autonomous one, i.e. the solution cannot in general be sought in the form of exponentials . For time-periodic modulation of the basic state one can resort to Floquet theory (Rosenblat & Tanaka (1971); Bhadauria & Bathia (2002)) but, in the general case, the formulation and prediction of critical conditions for the onset of convection becomes much less clear-cut than for a steady base state. The following two basic approaches are common for a general, time-dependent basic state:

  • reduction to an autonomous problem by the frozen-time assumption, whereby the basic state is supposed to evolve much more slowly than the perturbations, and

  • solution of the full non-autonomous linear perturbation problem for some initial conditions which are supposed to be representative. This is called amplification theory (denoted by AT in the following).

Both approaches go at least back to the 1960s. The first one was used by Lick (1965) and Currie (1967) and the second by Foster (1965). Later work by Homsy (1973) introduced energy stability ideas, but the resulting bounds are not necessarily useful for predicting instability thresholds and the mathematics is considerably more involved. In each of these works the focus was on buoyancy-driven convection.

The frozen-time approach was used for the Marangoni convection by Kang & Choi (1997). These authors studied the dynamics of a fluid layer subjected to a sudden change in surface temperature. However the frozen-time method was only applied for late times in the evolution. For early times they formulated the so-called propagation theory, whereby the problem is again reduced to an autonomous one by using a certain similarity solution for base state and perturbations. This approximation is actually appropriate for an infinite layer. The frozen-time model relies on the validity of the assumption of fast instability development (relative to the base state). This could be inadequate at instability thresholds, at which the instability can develop on the same time scales as the base state. The other method namely the amplification model, has a sounder mathematical foundation, but comes with the necessity of identifying representative initial conditions and amplification levels for the perturbations.

A specific analysis suitable for transient problems is developed in this paper, which is an extension of the previous amplification model. The underlying concept is called the non-normal approach. In such a case, the choice of initial perturbations is not somewhat arbitrary as in AT but it is guided by an optimization technique. This procedure consists in identifying, for a given time and a given wavenumber , the initial perturbation with the strongest amplification (called the optimal perturbation). Hence this method may exhibit initial conditions which are not a priori obvious configurations. As a consequence it gives the maximal value for the perturbation norm that was previously defined for the problem under study. This cannot be achieved through an AT simulation. Concerning the determination of initial conditions, we note that the non-normal approach has also been used for instability problems of steady base flow, e.g. in parallel shear flows, such as boundary layers (Schmid & Henningson (2001)) or plane Poiseuille flow (Reddy et al. (1998)). For such problems, there was a long-standing gap between the standard two-dimensional Tollmien-Schlichting modes obtained using the classical normal mode analysis – though these modes were only observed in some very controlled experiments – and the general experimental observation of streak generation. It is by resorting to non-normal mode analysis that the coherent structures observed in boundary layer transition, i.e. streamwise rolls and streaks, could be deduced from an optimal approach and their generation linked to the lift-up mechanism proposed by Landahl (Schmid & Henningson (2001)).

The non-normal approach can be extended to transient cases. This is precisely what we propose here. Such a method provides the strongest amplification, i.e. the optimal growth at a given duration after initial conditions have started their evolution. If this amplification is sufficiently large, the optimal perturbations may result in the modification of the basic flow at time and we can speak of an unstable regime in a way to be defined by a norm. To the best of our knowledge, this approach is the only one capable to provide clear-cut answers on instability problems for truly unsteady basic flows. However, as mentioned by one referee, this definition of the stability condition is blurred owing to a certain arbitrariness in choosing the norm and defining the critical amplification gain. This is inherent to any unsteady linear stability problem and does not depend on the particular method. It is then more suitable to view the results as the estimation of a transition region between a domain exhibiting strong convection and a domain where initial perturbations are damped or have no time to significantly develop during the transient regime. The non-normal approach allows one to characterize this transition domain properly. The results presented here for different norms and different critical amplification values show that this transition region is thin, so that the notion of stability threshold is still valid for this transient problem. On the contrary, the frozen-time assumption may fail, for instance if the base state evolves on the same time scale or faster then the unstable modes characterizing its frozen-time stability. In that case, the computed growth rates might not correspond to any true amplifications, at the system time scale. This might affect the determination of critical conditions 111 In the frozen-time assumption, the critical conditions are determined using the quasi-static growth rates . However the pertinent character of this method stongly depends on the rapid variation of these “quasi-static growth rates” with the basic state changes. For instance, the amplification between times and would be, at zeroth order WKB theory, equal to . So if the instantaneous growth rate strongly varies during this interval, its positive value at a given moment may not be interpreted in the right way. .

In the present work, we provide the results obtained by means of the non-normal approach as well as the frozen-time approach. In the specific flow case presented here, the quasi-static method is shown to provide similar results except for some particular quantities.

The proposed analysis is general and can be easily extended to many other unsteady problems, e.g., chemically driven hydrodynamic instabilities (Eckert et al. (2004)). In the present paper, we determine the onset of convection in a drying experiment, which leads to an unsteady Bénard-Marangoni problem. More specifically we study the sudden cooling due to evaporation of a liquid layer, where the decrease of surface temperature is induced by the vaporization latent heat. This problem has been the subject of many experimental or theoretical studies (see for example Berg et al. (1966), Vidal & Acrivos (1968), and more recently Mancini & Maza (2004) or Moussy et al. (2004)). The specific motivation for our work is the drying process studied in experiments by Toussaint et al. (2008) performed on a polymer solution (Polyisobutylene/Toluene). The solution initially at the ambient temperature is poured in a dish located in an extractive hood. When evaporation of toluene begins, convective patterns are observed at the very beginning of the experiment (quasi-instantaneous or less than 100 s after pouring the solution). They disappear well before the end of the drying. The very large Lewis number ( and denote the thermal and mass diffusivity) of the polymer solution (about ) is an indication that the thermal diffusion is faster and that convective patterns observed in the first minutes should be mainly driven by thermal effects. Two experimental observations detailed in Toussaint et al. (2008) support this thesis. First, a few experiments were conducted with deuterated solvent, whose density is higher than the polymer density. In that case, the density of the solution decreases when the polymer concentration increases, leading to a stable situation if the solutal Rayleigh-Bénard problem is considered. Since no differences were found with the experiments conducted with the standard solvent, we can exclude solutal buoyancy as a dominant mechanism. Second, free surface temperature fields measured by infrared camera showed that the end of free convection was related to the duration of the transient thermal regime. In these free convection experiments it can be inferred from the work on steady convection by Pearson (1958), that the Marangoni effect is dominant for thicknesses typically less than 1cm and the buoyancy dominant for higher thicknesses. For a more accurate description of the experiments see Toussaint et al. (2008).

The paper is organized as follows. In section 2, we present the basic assumptions of the model and the governing equations. Thereafter the unsteady basic state is described and a specific stability analysis is introduced in section 3. In particular the non-normal method is explained and the choice of norms is discussed. Section 4 contains the main numerical findings for the two limiting situations: the pure Rayleigh-Bénard and the pure Marangoni case. Critical conditions for the optimal modes are presented. Part of these numerical results are also obtained by a scaling analysis. A comparison of this method with the frozen time approach is also performed. Finally, in section 5, the comparison with experimental results is discussed.

2 Mathematical Model

2.1 Basic Assumptions of the Model

The mathematical model of Rayleigh-Bénard-Marangoni convection used throughout this paper is based on a one-layer model in which three assumptions are made : (1) the upper surface remains planar, (2) the layer thickness remains constant, (3) the heat and mass fluxes across the upper surface are given by transfer coefficients. Moreover our analysis is restricted to fluids characterized by a Prandtl number with the kinematic viscosity. This turns out to be the case of most liquids (including water and organic solvents). In that instance, the thermal diffusion time scale is always larger than the viscous diffusion time scale.

The first hypothesis can be tested as follows. Two modes of instability are known to occur in the Bénard-Marangoni problem (Scriven & Sternling (1964); Reichenbach & Linde (1981); Goussis & Kelly (1990)). One is generated from the interaction of the velocity perturbations with the basic temperature, while the other mode, characterized by a wavelength long compared with the layer thickness, is due to the coupling of the Marangoni effect with the deflection of the free surface. Mathematically, surface deformation can be neglected on scales of order when the Laplace pressure associated with a curvature is large compared with the dynamic pressure. This condition (Davis (1987)) corresponds with the smallness of the crispation number where and respectively denote fluid density and surface tension. Moreover, the Galileo number characterizes the relative importance of gravity ( gravity constant) and diffusion. A large value of the Galileo number indicates that gravity stabilizes the long-wave mode. The free surface deformation can be neglected if and . Such conditions are shown be satisfied for the experiments considered here.

The second hypothesis can be adopted if , where the Peclet number is defined as the ratio between the interface velocity due to evaporation which is equal to minus the time derivative of , and the thermal velocity scale . Indeed when , the surface displacement remains negligible compared to the total thickness during the problem characteristic time i.e. the diffusion time . In the experiment (see section 5), the Peclet number is smaller than .

Finally let us discuss the third assumption. The boundary condition at the free surface results from the coupling between the system and its surroundings. In evaporation experiments, the solvent flux and thus the temperature gradient in the fluid depends on the heat and mass transfer with the ambient air. Several authors have developed numerical or theoretical studies taking into account this coupling, see for example Merkt & Bestehorn (2003); Colinet et al. (2003); Ozen & Narayanan (2004); Moussy et al. (2004). In this paper we adopt a simple description by global heat and mass transfer coefficients, as our main interest is directed on the transient character of the problem under study and not on the detailed description of the transfer per se.

2.2 Governing equations

We formulate the basic equations in a Cartesian coordinate system, where the bottom of the layer coincides with the plane and the upper free surface with . The fluid is characterized by a density , a kinematic viscosity , a thermal diffusivity and a thermal expansion coefficient . The Boussinesq approximation is assumed to govern the velocity field and temperature field


where denotes the acceleration due to gravity and the temperature of the ambient air. In this approximation, the density is taken to be the density at . At the bottom, the velocity satisfies the no-slip condition and the wall is assumed adiabatic since, in the experiment, the bottom of the dish is thermally insulated by an air gap.


The upper boundary conditions are more involved. Assuming a planar surface, the local evaporation mass flux reads


Moreover the global mass balance reads (no fluid is introduced to counterbalance the evaporation mass loss)


where is the mean evaporation flux over the free surface. From relations 4 and 5 we get :


If we assume that the flux variations are much smaller than the flux itself, then , the Peclet number being defined using the evaporation velocity (see section 2.1). Since and is the velocity scale, the kinematic boundary condition so reduces to


In addition, the balance of tangential forces at the upper surface requires that the velocity field satisfies


where denotes the temperature inside the fluid layer at and denotes the surface tension which is assumed to be a linearly decreasing function of temperature,


Finally, the conservation of energy flux should be imposed at the interface. It reads :


The first l.h.s. term represents the heat conduction in the liquid, denoting the thermal conductivity related to the thermal diffusivity via with the liquid specific heat. The second l.h.s. term expresses the heat flux density in the gas using a simple phenomenological model based on the heat transfer coefficient and the difference between the air temperature far from the liquid and the surface temperature . Finally the cooling effect due to solvent vaporization is expressed in the r.h.s. term, where stands for the latent heat of vaporization and for the solvent mass flux per unit area. This latter quantity depends on the surface temperature . In the framework of the one-layer model, one imposes


where (resp. ) denotes the solvent concentration in the gas phase near the surface (resp. far from the surface) and is the phenomenological mass transfer coefficient in the gas. Assuming local thermodynamic equilibrium, directly depends on the surface temperature through the saturated vapour pressure. The variable can be linearized around which leads to the final expression:


Initially the liquid layer is isothermal i.e. . At large times, the system reaches a steady state in which the temperature in the layer is again uniform but with a temperature difference with the gas located far from the surface. This difference is imposed by the condition (12), where


To put the above equations in a non-dimensional form, scales for temperature, length, velocity and time are needed. The temperature difference provides the temperature scale and the layer thickness the relevant length scale. Two velocity scales can be introduced in this problem namely the viscous velocity scale and the thermal velocity scale . Here only the thermal scale is used. Finally the appropriate time scale is the thermal diffusion time . The non-dimensionalization leads to equations for and :


in which the Rayleigh, Marangoni, Prandtl and Biot numbers



In the following, we discriminate between two opposite cases: or . This corresponds respectively to an effective conductance in the gas much smaller or much larger than the heat conductance in the fluid.

3 Basic state and optimal linear perturbations

We study the stability of a purely conductive unsteady basic state which is initially uniform i.e. . This state evolves since the upper free surface is cooled down by the latent heat released through evaporation : The velocity field remains always zero but the unsteady field satisfies a pure heat equation


The basic state only depends on the Biot number . Note that the term characterizes the heat transfer in the gas phase and the term represents the cooling effect due to evaporation. The unsteady temperature field is shown in figure 1 for different Biot numbers and times. A cooled layer develops from the upper surface whose thickness is similar at a given time for different Biot numbers and increases until it fills the whole layer. Conversely, if denotes the characteristic temperature difference within the fluid, the maximum reached over the whole time evolution by this quantity increases with Biot number (figure 1d). Actually, two different regimes are observed according to the value of 222Details are contained in Appendix A.. For small Biot numbers, typically , the cooled layer reaches the bottom while the jump is less than one. Afterwards, decreases. Such time evolution can be summarized by the following scalings


For , the temperature field relaxes towards the steady uniform temperature and all the energy needed by evaporation is carried on by convection in the gas phase.

For large Biot numbers, typically , the temperature jump reaches a maximum before the cooled layer reaches the bottom (figure 1a). In that case, the maximal jump is equal to one and at that time. The time evolution can be summarized by the following scalings :


For , the temperature decreases in the whole layer thickness to reach the steady state regime .

In order to study the stability of the unsteady conductive state we split the fields into a basic flow and three-dimensional perturbations


and linearize in the perturbations to get a set of linear equations. Since this problem has no preferential direction in the plane, the perturbation Fourier modes in such directions decouple in the linear regime. Without loss of generality, we thus consider a nondimensional wavenumber in the direction and no dependence in the direction reducing the flow to be two-dimensional


These infinitesimal perturbations are governed by the linear system


To quantify the amplification gain at time , it is customary to define a norm which is generally based on the kinetic energy of perturbations. This norm can be orthogonally decomposed in a Fourier basis so that the individual contributions of each wavenumber  can be studied independently. In the present case, the temperature field is playing a major role as well. We thus define two different norms corresponding to two different situations. The first one is based on the kinetic energy of perturbations


where superscript denotes complex conjugation. The integration is performed over the entire layer width and perturbations are obtained after integrating the above linear system over the time period . The second one


is based on the temperature field and not the velocity field.

For finite Prandtl number two extreme cases are considered, with the initial perturbation concerning either the velocity field, or the temperature field. In the following we will consider the amplification factors or to characterize the stability. Then, when the initial perturbations only concern the velocity field, we only take into account the amplification since . Conversely, when the initial perturbations only concern the temperature field we use the amplification factor . For the infinite Prandtl number case, the velocity perturbations are not dynamical quantities since the time derivative of the velocity drops out from equations (28)-(29). In other words, the velocity perturbations are slaved to the temperature perturbations. In this case, we use only and hence the norm .

For a given initial disturbance, one evaluates the amplification gain at time by computing where or . It is the purpose of a non-normal analysis to compute the quantity which is the upper bound for the energy amplification that a disturbance of wavenumber can reach at time . This approach solves a sort of the finite time stability problem : it investigates the transient evolution and defines for any time , an optimal perturbation mode which actually reaches the upper bound (Schmid & Henningson (2001)). This mode (i.e. the optimal initial profile) is found numerically by solving an optimization problem (Farrell & Ioannou (1996); Andersson et al. (1999); Luchini (2000); Schmid & Henningson (2001)). The optimum of is determined taking into account several constraints: (i) the disturbance energy at time is equal to unity; (ii) the disturbance satisfies the linear governing equation as well as the boundary conditions during the complete time interval . This problem is best solved with the help of a Lagrangian formalism and Lagrangian multipliers which are introduced to precisely enforce the above constraints. In the present case, these multipliers are adjoint fields . Following a standard derivation, these quantities satisfy a set of adjoint equations. It is


in which . These adjoint equations have to be solved backwards in time. Let us denote by the symbol the vector field . One obtains the optimal perturbation for time by an iterative scheme which propagates a given initial condition forward in time using the direct problem (here denoted by ), the result of which serves as an “initial” condition for the backward propagation by the adjoint equations (here denoted by ). More specifically 333Details are contained in Appendix B., a relation between the adjoint and is imposed. After one forward-backward integration, quantity is obtained and a relation between and is also imposed. An updated initial condition for the next iterative step is then available. This process should be self-consistent : one uses an iteration procedure which is schematically illustrated by a diagram


Convergence is reached when the initial condition for the forward problem does not change appreciably – up to a normalization constant – by an appropriately chosen criterion from one iterative step to the next. The converged mode is precisely the initial optimal perturbation for time . The maximum energy amplification is computed by propagating the converged initial condition forward in time and by forming the ratio of the disturbance energy at the end of the time interval to the energy at the beginning. The direct and adjoint equations have been discretized using a pseudospectral scheme based on Chebyshev polynomials and a streamfunction-based formulation to account for incompressibility444Details are contained in Appendix C..

4 Numerical Results

4.1 Quantities provided by the non-modal analysis

For the unsteady basic state , a given set , and a given type of initial perturbation (temperature or velocity), the non-modal analysis determines at a given time , the maximum energy amplification over all possible perturbations of wavenumber (see figure 2). In a way, the value is analogous to a growth rate for classical stability analysis. More generally, it appears possible to extend the usual concepts of classical stability analysis to unsteady flows. For instance, can be maximized over wavenumber and time providing the global maximum amplification (see figure 3). This value is reached at time , for an optimal wavenumber denoted by and for a specific perturbation structure in . These latter two quantities play the role of the most amplified wavenumber and of the most amplified mode for the standard stability analysis. One also obtains a ”stability” diagram, by determining the region of the space where the amplification gets above a threshold . It is a way of separating the region where amplification or attenuation occurs. The value of , for instance , is somewhat arbitrary : as already said in the introduction, an unsteady problem is indeed characterized more by a transition domain than a well defined threshold. It is demonstrated below that choosing or does not result in major differences, so that the transition region is thin compared with the absolute values of the Marangoni and Rayleigh numbers. For fixed and , one can determine the curve such that if (resp. ) then (resp. ). Similarly one may define for fixed and , the curve . These curves play a role very much similar to marginal stability curves. Each point of the critical curve is associated to a critical wavenumber and critical optimal time . Note that until this point, most of this procedure can be extended to other unsteady problems in a straightforward way. A comparison between these critical curves and the experimental diagram that separates the domains where convection is observed or not observed, is made in section 5.

4.2 Infinite Prandtl number: the pure Marangoni case and .

For infinite Prandtl number, the velocity is slaved to the temperature field. As a consequence, only perturbations in temperature field are pertinent. In the plane , the critical curve , critical wavenumber and critical optimal time are presented for the pure Marangoni case and two thresholds and (figure 4). The critical Marangoni number slightly depends on the value of the threshold and seems to be consistent, within the numerical uncertainties, to the following laws (here )


Moreover the wavenumber (figure 4b) is an increasing function of the Biot number while optimal time is a decreasing function of the same parameter.

Using heuristic arguments 555Details are contained in Appendix D., the following scaling laws can be deduced :
For the critical conditions for instability can be expressed as


For , the critical conditions become


The spatial structure in of the optimal perturbation at and is shown on figure 5 for two Biot numbers and two different times : time and time when the perturbation reaches its maximum amplification. The spatial structure of this optimal perturbation is shown to change slightly during the time evolution. In this respect, this optimal mode does not differ much from the classical most amplified mode of steady problems.

4.3 Infinite Prandtl number : the pure Rayleigh case and .

In the plane , the critical curves , and are presented for two thresholds and (figure 6). The critical Rayleigh slightly depends on the value of the threshold and seem to be consistent, within numerical uncertainties, to the following laws (here ) :


Moreover the wavenumber (resp. time ) is an increasing (resp. decreasing) function of the Biot number for small Biot and reaches a plateau for larger Biot number.

One can deduce using heuristic arguments 666Details are contained in Appendix D., the following scaling laws:


4.4 Comparison with classical steady results and with transient frozen-time method

It is worth comparing the results presented here with the well-known results obtained by Pearson (1958) 777In this case, we computed some additional results to cover a larger range of Biot numbers than the one given by Pearson in his article. (resp. Sparrow et al. (1963)) in the framework of the steady Marangoni (resp. Rayleigh) problem. This comparison is pertinent since the boundary conditions at the top and the bottom of the layer (equations 31-32) are similar in the present work and in these classical analyses. However, the Marangoni and Rayleigh numbers defined by these authors are based on the steady temperature difference between the top and the bottom of the layer. This steady temperature difference is missing in the transient problem under study. It is then not possible to make a direct comparison between the thresholds values obtained in our paper and the previous ones from Pearson’s or Sparrow’s publications, and a preliminary transformation is needed. Indeed, at each time , one might define the equivalent temperature difference between the top and the bottom of the layer i.e. . Using such a temperature difference, it is then easy to define a time-dependent Marangoni or a time-dependent Rayleigh numbers .

Let us compute the maximum of obtained during the time evolution for the Marangoni (resp. Rayleigh) case. This maximum is reached at time and leads to a new Marangoni number (resp. Rayleigh number ) which can be compared to the critical values (resp. ) obtained by Pearson (1958) (resp. Sparrow et al. (1963)) from the steady case. With the scalings used here, the critical Marangoni predicted from these steady results reads


Such estimations have been plotted on figure 4 and 6. They are found to be close to our results for . The same comment applies for the critical wavenumbers except for the Bénard-Marangoni case at high Biot numbers. For critical times, however, the steady-state approximation differs from our results.

When using the results by Pearson (1958) or Sparrow et al. (1963), we are clearly using the normal mode results obtained for a linear temperature field in on the whole thickness, which is a rough approximation especially at high Biot number. We can go even further and compare the non-normal mode results within the frozen-time approximation. In this latter approximation, a stability analysis in terms of normal modes is performed at the each time . The ”steady” base flow is assumed to be the temperature field computed at this specific time . When the flow is stable within this frozen-time approximation for each time , it will be assumed stable. When the control parameter reaches a critical value (here or ), there exists a unique critical time for which the frozen time state possesses a marginal eigenvector with a critical wavenumber. In the present problem we have determined the critical parameters for neutral conditions in the frozen-time case by the Lapack routine GGEV for generalized eigenvalue problems in combination with a Chebyshev collocation method. As can be seen in figures 4 and 6, results for the thresholds and the critical wavenumbers using the non-normal approach () and the frozen-time approximation are close. However the prediction of critical times still differs 888In the frozen-time approximation, corresponds to a time when the quasi-static growth rate becomes zero. This critical time is of a different nature than the critical time for the non-normal approach. The latter characterizes the perturbation evolution over the time interval from up to . In zeroth order WKB theory, the critical time for the non-normal approach would be determined by . Suitable modifications of the frozen-time analysis based on this observation should therefore lead to closer agreement regarding the critical times. We do not pursue this issue further in the present work.. The frozen-time approximation apparently provides a bound of order unity for because the available temperature difference attains its maximum as soon as the thermal boundary layer of the basic temperature distribution has reached the bottom of the layer. For small , remains fairly constant for larger times, and the instability can develop on this quasi-steady background over fairly long times. This observation can explain the apparently unbounded growth of with for in the non-normal analysis.

4.5 Results for finite Prandtl numbers

In this section, we focus on finite Prandtl numbers and more specifically on the role of initial perturbations on the transition zone estimation. In this case, the velocity field is not slaved to the temperature field so that perturbations in velocity can be considered as well as perturbations in temperature. We only discuss the curves for the pure Marangoni case () (Figure 7) but similar results apply to the pure Rayleigh case ().

Since the frozen time approximation seems to be valid for this unsteady problem and an exchange of stability in a normal mode is not affected by the Prandtl number, the effect of this latter number should not be significant. Indeed, for temperature perturbations, the critical Marangoni is not affected when infinite number case is compared to . For the optimal time and wavenumber , the same conclusions apply. For velocity perturbations, the curves differ but not in a drastic way taking into account that perturbations are of a completely different nature compared with the infinite-Prandtl-number problem. Perturbing the temperature is more efficient than perturbing the velocity. Indeed critical Marangoni numbers are smaller for temperature perturbations, but the difference is of the same order as the one obtained by changing the threshold value from 1 to 100. Actually, the “blurriness” of this transient problem induced by the choice of the threshold values and perturbation types is not very broad and does not modify the order of magnitude of the critical numbers if one excepts the time . Finally, let us note that non-linear direct simulations have also been made to solve this problem(Touazi et al. (2009)), showing very good agreement with the linear results presented here.

5 Comparison with experiments

The present section is devoted to the comparison between results obtained from the optimal mode calculations and from an experimental work described in Toussaint et al. (2008), where transient Rayleigh-Bénard-Marangoni convection is generated by drying a polymer solution of PolyIsobutylene-Toluene at ambient temperature. In the experiments, buoyancy and Marangoni effects are equally present. It is the evaporation of the solvent i.e. toluene which cools the upper surface by latent heat. During the experiments, the following parameters are kept constant:


Different thicknesses and dynamic viscosities are considered. is varied from to while dynamic viscosity is set to a value in the range by monitoring the initial polymer concentration. These data allow us to estimate the crispation, the Galileo and the Peclet numbers for each experiment. We get the following bounds for these numbers:


As a consequence, the assumption of planar free surface and constant thickness layer are justified. The other relevant non-dimensional numbers vary in the following range :


The comparison is displayed using the critical dynamic viscosity as a function of thickness (figure 8). In the plane , each point corresponds to a given parameter set . The four critical curves correspond to two thresholds ( and ) and two types of initial pertubation (temperature and velocity). The theoretical critical curves divide the experimental points corresponding to regions of convection or pure conduction in a satisfactory manner. The temperature perturbation critical curve is above the velocity one. This analysis in the thickness/viscosity plane shows again that the bandwith of uncertainty due to the choice of threshold and perturbation types is not very broad and does not modify the order of magnitude of the critical thickness.

6 Conclusion.

This paper presents a novel linear stability analysis of an unsteady base state within the general conceptual framework of amplification theory. The non-normal approach is used, which possesses the advantage over more classical methods to solve the transient problem in a mathematically rigorous way. In turn, this allows one to test other approximations. Here we have applied this approach for the first time to characterize the stability of a transient Rayleigh-Bénard-Marangoni problem in an horizontal fluid layer suddenly cooled from above. It provides the upper limit of the energy amplification that a disturbance of wavenumber can reach at time . This quantity reaches a maximum at time , for a specific optimal wavenumber and a specific perturbation structure in . These latter two quantities play the role of the most amplified wavenumber and most amplified mode for the standard steady analysis.

A ”stability” diagram in the space has been determined for the pure Marangoni and the pure Rayleigh problem by the non-normal approach. Note that the marginal conditions used to determine the stability curve was obtained by setting the optimal amplification equal to 1 or 100. The critical time and critical wawenumber were evaluated for this marginal conditions. Critical Marangoni and Rayleigh numbers exhibit a strong dependency on the Biot number and a weak sensitivity to Prandtl number variations in the range . Scaling exponents for critical Rayleigh or critical Marangoni versus Biot numbers have been found numerically and confirmed by scaling analysis in the limit of very small and very large Biot numbers. Comparison of the non-normal approach with the frozen-time approximation (a classical quasi-static approach) shows similar results for the critical Marangoni or Rayleigh numbers and the critical wave numbers. Moreover, the “blurriness” inherent in any transient problem was analyzed as a function of the amplification threshold values and perturbation types. It has been shown that the transition region is thin compared to the large domain of Rayleigh or Marangoni numbers covered by the analysis.

Finally, a comparison with experimental results has been performed, where convection is induced by solvent evaporation during the drying of polymer solution. A good agreement was indeed found between the present theoretical study and experimental observations. The method was developed in this paper for the cooling of a fluid induced by solvent evaporation but could easily be extended to other transient problems.

TB and MR acknowledge financial support from the Deutsche Forschungsgemeinschaft in the framework of the Emmy–Noether Program (grant Bo 1668/2).

Appendix A Basic State

We obtain here the scaling laws given in section 3 which provide the evolution of a purely conductive basic state. It is recalled that the basic temperature field is initialy uniform i.e. . and satisfies a heat equation


with the boundary conditions


Due to the evaporation, a cooled layer of characteristic thickness develops from the upper surface. Equation (51) provides the standard estimate


In the following, we determine the scaling laws for the two extreme cases and .

A) Case

Since is initially zero, it remains small during a period of time  (the time is determined below and shown to be much larger than ). Diffusion in the liquid and evaporation terms thus dominate in the free-surface boundary condition (52). Such a balance can be expressed in terms of order of magnitude as follows


One thus obtains using equation (53)


During this time interval, the cooling layer has not reached the bottom at so that temperature field reads in terms of orders of magnitude


Thereafter the thickness remains constant and the following condition holds


During this latter time interval, the heat equation (51) can be used with scaling (57) to get the temperature field in terms of orders of magnitude


These equations are valid if evaporation dominates heat transfer in the gas phase. This requires that and determines the value . For , the temperature field relaxes towards the steady uniform temperature equal to . Since the temperature gradient is equal to zero in that regime, all the energy due to evaporation is transfered by convection in the gas phase.

B) Case

For small times, the condition holds and an analysis similar to the one performed for the case is valid leading to


The value of is obtained by determining the time when . At that time, the heat flux in the gas phase becomes of the same order of the evaporation. During this time interval, the cooling layer has not reached the bottom at so that temperature field reads in terms of orders of magnitude


For , a new regime begins where the surface temperature remains constant while the cooled layer thickness keeps increasing


This regime ends when at time . Thereafter the temperature decreases in the whole layer thickness to reach the steady state regime .

The case is the limiting case of the two previous ones. The equation (55) is valid until , when the cooled layer thickness and the surface temperature both reach their extremum. Thereafter the temperature decreases in the whole layer thickness, to reach the steady state.

Appendix B Obtaining the Adjoint Equations

Let us denote by , the components of the vector field

To find the maximum amplification at a given time , we maximize the perturbation norm


at time with respect to the set of all possible initial perturbations such that . It is recalled that the integration is performed over the entire layer depth and the superscript denotes complex conjugation. Coefficients are weight coefficients chosen to put emphasis on temperature or velocity acoording to the case considered. To analyse the initial perturbation in velocity, the kinetic energy norm is used and one takes and . To analyse the initial perturbation in temperature, the temperature norm is used and and .

The variation with respect to a variation of the initial perturbation is to be evaluated. This computation cannot be performed in a straightforward manner since the energy can be explicitly written in terms of but only implicitely in terms of . It is known via several constraints: normalization of and time integration over the interval of equations (3.8)–(3.10). These dynamic equations relating to , are formally written here as . This optimization with constraints necessitates the introduction of Lagrangian multipliers, the so-called adjoint fields .

More specifically, a Lagrangian function is defined, which depends on direct and adjoint variables over the interval , and a normalization scalar :


where means complex conjugate and stands for the scalar product


When satisfies the constraints (direct problem plus normalization at ), all terms but the first one on the r.h.s. of equation (63) are zero and, by consequence, and