Transient Rayleigh-Bénard-Marangoni solutal convection Benoît Trouette, 1, 2, a) Eric Chénier, 2, b) Frédéric Doumenc, 3, c) Claudine Delcarte, 1, d) and Béatrice Guerrier 3, e) 1) Univ Paris-Sud, CNRS, Lab. LIMSI, Orsay, F-91405

2) Univ Paris-Est, Lab. Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 5 bd Descartes, Marne-la-Vallée, F-77454

3) UPMC Univ Paris 06, Univ Paris-Sud, CNRS, Lab. Fast, Orsay,

F-91405

Solutal driven flow is studied for a binary solution submitted to solvent evaporation at the upper free surface. Evaporation induces an increase in the solute concentration close to the free surface and solutal gradients may induce a convective flow driven by buoyancy and/or surface tension. This problem is studied numerically, using several assumptions deduced from previous experiments on polymer solutions.

The stability of the system is investigated as a function of the solutal Rayleigh and Marangoni numbers, the evaporative flux and the Schmidt number. The sensitivity of the thresholds to initial perturbations is analyzed. The effect of viscosity variation during drying is also investigated. At last numerical simulations are presented to study the competition between buoyancy and Marangoni effects in the nonlinear regime.

a) Electronic mail: benoit.trouette@univ-paris-est.fr b) Electronic mail: eric.chenier@univ-paris-est.fr c) Electronic mail: doumenc@fast.u-psud.fr d) Electronic mail: delcarte@limsi.fr e) Electronic mail: guerrier@fast.u-psud.fr

I. INTRODUCTION

Numerous natural and industrial processes are affected by evaporation of solutions, like drying of salty lakes, drying of paint films, crystal growth or distillation. Evaporation induces a decrease of both temperature and solvent concentration near the surface. Because the density and surface tension of the solution are usually temperature and concentration dependent, buoyancy forces and/or surface tension variations can destabilize the liquid layer leading to pattern formation [START_REF] Colinet | Nonlinear Dynamics of Surface-Tension-Driven Instabilities[END_REF][START_REF] Manneville | Dynamics of Spatio-Temporal Cellular Structures[END_REF] . When one of the components is non-volatile, the final state of the dried product may depend strongly on the convective flow induced by the evaporation, and understanding the formation of the convective patterns is an important challenge.

The thermal problem (namely Rayleigh-Bénard or Bénard-Marangoni convection induced by temperature gradient) is also encountered during the evaporation of pure liquids and has been widely studied. In the majority of theoretical and numerical works published on evaporation problems, the transition from diffusion to convection is studied for asymptotic temporal regimes. Indeed several authors have developed linear stability analysis for a steady basic state, taking into account hydrodynamics in the liquid and in the gas [START_REF] Merkt | Bénard-Marangoni convection in a strongly evaporating fluid[END_REF][START_REF] Moussy | Influence of evaporation on Bénard-Marangoni instability in a liquid-gas bilayer with a deformable interface[END_REF][START_REF] Ozen | The physics of evaporative and convective instabilities in bilayer systems[END_REF] .

Margerit et al. [START_REF] Margerit | Improved 1.5-sided model for the weakly nonlinear study of Bénard-Marangoni instabilities in an evaporating liquid layer[END_REF] also addressed dynamics in both phases, but considered only diffusion in a gas at rest. It should be emphasized that an asymptotic regime is not always present.

For instance in configurations with adiabatic solid walls, evaporation induces a transition from an initial motionless state with uniform temperature to a final new isothermal state at a lower temperature [START_REF] Toussaint | Experimental characteri-33 zation of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF] . Between these two states, the regime is time-dependent, diffusive or convective. The cooling by evaporation of a liquid layer has been investigated by Vidal and Acrivos 8 in the framework of a frozen time approach. The frozen-time assumption consists in a decoupling of the temporal growth rate of the disturbances from the characteristic time-scales of the basic solution. Other approaches that take the time dependence of the basic state explicitly into account have been developed using amplification theory [START_REF] Foster | Stability of a homogeneous fluid cooled uniformly from above[END_REF] or a more recent linear analysis based on a non-normal approach [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis[END_REF] . Results have been compared to numerical simulations [START_REF] Touazi | Simulation of Transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF][START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF] as well as experimental observations [START_REF] Toussaint | Experimental characteri-33 zation of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF] . Other experiments to characterize flow structures have also been reported in the literature. Berg et al. [START_REF] Berg | Natural convection in pools of evaporative liquids[END_REF] used the Schlieren effect to study flows in various fluids of different thicknesses. Turbulent flow has been observed by Colinet et al. [START_REF] Colinet | Interfacial turbulence in evaporating liquids: theory and preliminary results of the itel-master 9 sounding rocket experiment[END_REF] in microgravity experiments, in order to focus on Marangoni effect. Transition from hexagonal to square cells has been obtained by Mancini and Maza 15 by increasing the liquid thickness.

When the fluid is a binary solution, Rayleigh-Bénard-Marangoni convection can arise not only due to convection induced by temperature variations but also because of concentration gradients. In the latter case, the flow is then known as solutal Rayleigh-Bénard-Marangoni convection. The literature dealing with concentration gradient-induced convection is more sparse than for thermal convection. Ha and Lai [START_REF] Ha | Onset of Marangony instability of a two-component evaporating droplet[END_REF] have investigated theoretically the onset of Marangoni instability due to the evaporation of a two-component evaporating droplet using a quasi-steady approximation. Machrafi et al. [START_REF] Machrafi | Bénard instabilities in a binary-liquid layer evaporating into an inert gas[END_REF] have also performed a stability analysis using a quasi-stationary reference state, taking into account thermal and solutal buoyancy and the Marangoni effect, as well as the Soret effect (water-ethanol solution). De Gennes 18 used scaling arguments to estimate the critical thickness for the onset of convection in a dilute polymer solution. An additional point may occur in complex fluids when one of the components is non volatile (for instance in polymer solution or colloidal suspension). Indeed, physical properties are often strongly dependent on solvent concentration, so that the ratio between destabilizing and stabilizing forces changes continuously during drying. Practical consequences of this can be significant and several experimental studies mention a possible correlation between solutal convection and wrinkles or surface corrugations observed on dried films [START_REF] Sakurai | Control of mesh pattern of surface corrugation via rate of solvent evaporation in solution casting of polymer film in the presence of convection[END_REF][START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF] . This paper presents the results of a numerical study of the drying of complex fluids. It focuses on the onset of convection for the transient solutal regime, and aims to analyze the influence of the concentration-dependent viscosity on the stability thresholds. We develop a 2D nonlinear numerical simulation of the Rayleigh-Bénard-Marangoni solutal convection for films with flat interfaces, taking into account the decrease of liquid thickness and the variation of viscosity during drying. Simplifying assumptions and the assumed magnitudes of physical parameters are inferred from previous drying experiments performed upon polymer solutions, as detailed in the next section. The equations and dimensionless parameters characterizing the problem are then introduced. In this transient problem, the formulation and prediction of critical conditions for the onset of convection are much less clear-cut than when the instability concerns a steady basic state [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis[END_REF] . Consequently it is necessary to first establish a criterion to define the transition from situations mainly driven by diffusion to cases with significant flow velocity.

The transient behavior of this problem makes the solution sensitive to initial conditions. Thus, we analyze the effect of random perturbation on dynamics and we characterize the transition zone between the diffusive and convective configurations.

Critical Marangoni numbers and critical Rayleigh numbers are determined as a function of the evaporative flux and initial thickness. Systematic comparisons between results with constant and variable viscosity are presented in order to highlight the part played by the increase of viscosity during drying.

At last numerical simulations are presented to study the competition between buoyancy and Marangoni effects in the nonlinear regime. For that purpose the two mechanisms are considered first separately, then the coupled problem is adressed.

II. MATHEMATICAL MODEL

In the following the superscript/subscript s, p, i, int indicate solvent, solute, initial and interface value, respectively. A symbol without superscript/subscript refers to the fluid solution. Vectors are indicated in boldface.

A. Assumptions

The mathematical model of solutal Rayleigh-Bénard-Marangoni convection used in this paper is based on a one-layer model (i.e. the vapor dynamic is neglected), with the following geometrical and physical assumptions :

• constant and uniform evaporative flux.

• Newtonian and incompressible fluid, with a constant mass diffusion coefficient and constant or variable viscosity

• constant density except in the buoyancy term,

• solutal convection,

• 2D rectangular geometry with a mobile but flat interface,

• local thermodynamic equilibrium at the interface.

The numerical approach presented here follows an experimental study dealing with drying of polymer solutions 7 and several of the above assumptions are inferred from typical characteristics of polymer solutions. Figure 1 been previously analyzed both experimentally and theoretically [START_REF] Toussaint | Experimental characteri-33 zation of buoyancy-and surface tension-driven convection during the drying of a polymer solution[END_REF][START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis[END_REF][START_REF] Touazi | Simulation of Transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF][START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF] , so that it will not be considered here.

The present paper focuses on the second regime (domain 2 in Figure 1(a)), where the temperature and the evaporative flux are nearly constant. This regime results from polymer solution properties. In the framework of the one-layer model, the local solvent evaporative flux Φ ev (kg/(m 2 .s) ) is given by

Φ ev = h m (c g s,int -c g s,∞ ) (1)
where h m is the mass transfer coefficient (m/s), c g s,int and c g s,∞ are the solvent concentration in the gas phase just above the interface and far from the interface, respectively (kg/m 3 ).

The latter is zero for experiments performed in the open atmosphere with organic solvents.

Using the ideal gas law, we get:

Φ ev = h m M S P V S0 (T ) RT a s (T, ϕ s ) (2)
where M S is the solvent molar mass, a s is the solvent activity, P V S0 is the saturated solvent vapor pressure, T is the temperature, ϕ s is the solvent volume fraction in the liquid phase at the interface and R is the ideal gas constant.

For polymer solutions, the activity is close to one over a large range of solvent concentration, as illustrated in Figure 1(b) for the polyisobutylène(PIB)/toluene solution [START_REF] Doumenc | Mutual diffusion coefficient and vapor-liquid equilibrium data for the system PIB/Toluene[END_REF] . It is usually expressed by the well-known Flory-Huggins model:

a s = (1 -ϕ p ) exp(ϕ p + χ ϕ 2 p ) (3)
where ϕ p = 1 -ϕ s is the polymer volume fraction in the solution at the interface, and χ the interaction parameter which characterizes the affinity between the solvent and the polymer (χ = 0.496 + 0.261ϕ p for PIB/toluene at 25 0 C 21 ).

In the following we assume that the temperature is constant and we neglect the variation of the activity for small polymer volume fractions: a s ≃ 1 for ϕ p,int < 0.6 (cf Figure 1(b), a s > 0.9 in this concentration domain). Then, from equation 2, we get a known and constant evaporative flux. The relative error induced by setting a s = 1 everywhere is less than 10%

for the local variations of the evaporative flux as well as for its mean value. The boundary condition at the free surface used in this paper is then a uniform, constant and a priori known evaporative flux. The validity domain of our simulations is restricted to ϕ p,int < 0.6.

In the last stage, not considered here, the evaporative flux drops rapidly when the saturated vapor pressure of the solution becomes strongly dependent on the solvent concentration (domain 3 in Figure 1(a)).

In the regime studied in this paper (domain 2), we consider a binary solution with a mutual diffusion coefficient D, which can be assumed constant in the concentration domain of our simulations (ϕ p,int < 0.6) [START_REF] Doumenc | Mutual diffusion coefficient and vapor-liquid equilibrium data for the system PIB/Toluene[END_REF] . The solution density is also assumed constant (constant and same partial specific volumes for the solute and the solvent), except for the buoyancy term in the Navier-Stokes equation. The following linear relations are used to express the solution density, ρ and the surface tension, σ:

ρ = ρ 0 (1 + β ϕ p ) (4)
σ = σ 0 + γ ϕ p (5)
where ρ 0 and σ 0 are the density and surface tension of the pure solvent, ϕ p = 1 -ϕ s is the polymer volume fraction and β and γ are positive constants.

As a first approximation the fluid is assumed Newtonian and incompressible. The description of the solution rheology is limited to the use of a variable viscosity µ that depends strongly on the solute concentration. Moreover, since one of the main issues of this paper is to analyze the role of viscosity increase during drying, systematic comparisons are made for fluids with constant or variable viscosities. For simulations corresponding to variable viscosity we have used an empirical law interpolated from measurements performed on a PIB/toluene solution (cf. Figure 2 and22 ):

log 10 (µ) = a 0 + a 1 X + a 2 X 2 + a 3 X 3 + a 4 X 4 (6)
with X = log 10 (ϕ p ), (a 0 , a 1 , a 2 , a 3 , a 4 ) = (8.235, 14.02, 6.575, 1.392, 0.1114).

The cavity is chosen to be two-dimensional for the sake of simplicity and to reduce computational cost. Previous studies performed in the thermal regime showed very similar values for the thresholds when using 2D or 3D simulations [START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF] . The domain is therefore rectangular, with an initial aspect ratio A i = L/e i , where L is the length and e i the initial thickness of the solution layer. A i = 10, except in the last section dedicated to nonlinear simulations where A i = 5 in order to reduce computational time. In what follows, e x (e z )

is the unit vector in the horizontal (vertical) direction. The evaporation takes place at the upper surface. The thickness e(t) is time dependent but the free surface is assumed nondeformable. This assumption is fully justified as long as the surface tension and the gravity are large enough to balance the free surface deformation induced by the fluid flow. Close to the thresholds, and using the diffusion velocity as a reference velocity, this corresponds to the crispation number Cr ≡ (µD)/(σe) ≪ 1 and the Galileo number Ga ≡ (ge 3 ρ)/(µD) ≫ 1 [START_REF] Regnier | Linear and nonlinear Rayleigh-Bénard-Marangoni instability with surface deformations[END_REF] , where g is the acceleration due to gravity. For higher velocities or variable viscosity we have used an a posteriori procedure to estimate the domain of validity for the flat interface assumption (cf. section IV C).

B. Governing equations

Let us define v = ue x + we z to be the mass-averaged velocity of the solution, given by :

v = c s v s + c p v p ρ = ω s v s + ω p v p (7)
where c s , ω s and v s are the local solvent concentration, its mass fraction, and velocity, respectively. Subscript p applies to the solute. With the above assumptions, the velocity and the solvent volume fraction are governed by the following equations:

div v = 0 (8) ∂v ∂t + (v • ∇) v = (1/ρ)[-∇p + div µ∇v + µ∇ T v ] -gβ ′ (ϕ si -ϕ s )e z (9) ∂ϕ s ∂t + (v • ∇) ϕ s = D∇ 2 ϕ s (10)
with ρ being the initial density, p the deviation from the initial hydrostatic pressure and

β ′ = β/(1 + β × ϕ pi )
, with β the solutal expansion coefficient (from equation ( 4)), and ϕ pi the initial polymer volume fraction. We impose no-slip boundary conditions along the rigid bottom (z = 0) and lateral walls (x = 0 and x = L). Across these boundaries the normal mass flux is zero (non-permeable surfaces).

z = 0 and x = 0, L    v = 0 , ∂ n ϕ s = 0 (11)
The upper boundary conditions at the free surface (z = e(t)) are defined in the following way. With Eq. ( 5), the stress equilibrium gives the boundary condition for the horizontal velocity u:

µ ∂u ∂z = -γ ∂ϕ s ∂x . (12)
The global and solvent mass flux conservation relations at the free surface give the two other boundary conditions, for the vertical component of the velocity and the concentration. We first express these two relations in the general case of different partial specific volumes for the solute and the solvent. They are then simplified using the assumption of constant solution density.

Following Fick's law, the local flux of solvent is given by:

c s v s = c s v -Dρ∇ω s . (13)
The solvent flux across the moving liquid/gas interface reads:

c s (v s -v int ) • e z = (c s (v -v int ) -Dρ∇ω s ) • e z ( 14
)
where v int = v int e z is the interface velocity.

Since the solvent is the only volatile fluid, the volumic evaporative flux is v evap = Φ ev . Vs , Vs being the solvent specific volume in liquid phase. Moreover, the bottom being impermeable, the interface velocity v int is equal to -v evap . In the framework of this model the interface velocity is then constant and a priori known. The solvent mass conservation across the interface reads:

c s (w s -v int ) = c s (w -v int ) -Dρ ∂ω s ∂z = Φ ev = v evap Vs ( 15
)
Taking into account that the polymer does not evaporate, the same analysis gives a similar equation for the polymer:

c p (w -v int ) -Dρ ∂ω p ∂z = 0 (16)
Adding Eqs. 15 and 16, with c s + c p = ρ and ω s + ω p = 1:

w = 1 ρ Vs -1 v evap . (17)
Neglecting the density variation as already done in the diffusion equation, we obtain

w = 0 (18)
Finally Eq. ( 15) reads

-D ∂ϕ s ∂z = v evap (1 -ϕ s ) . ( 19
)
Boundary conditions at the upper interface are then described by equations ( 12), ( 18)

and (19).

C. Dimensionless equations

To transform to dimensionless relations, coordinates x and z, the velocity v = ue x +we z , the dynamic pressure p and the time t are scaled by e i , D/e i , ρD 2 /e 2 i and e 2 i /D, respectively. The solvent volume fraction is scaled in the following way:

φ s = ϕ s -ϕ si ∆ϕ with ∆ϕ = - ∂ϕ s (t = 0) ∂z int × e i = P e int (1 -ϕ si ) (20)
P e int = v evap e i /D is the non dimensional velocity of the interface and is then called in the following the Peclet number of the interface.

For the sake of clarity, no extra notation is used to highlight the new variables. The dimensionless form of Eqs. ( 8)-( 10) is:

div v = 0 (21) ∂v ∂t + (v • ∇) v = -∇p + Sc div µ µ i ∇v + ∇ T v +Ra Sc φ s e z ( 22
) ∂φ s ∂t + (v • ∇) φ s = ∇ 2 φ s (23)
where Sc = µ i /ρD is the Schmidt number and Ra = gρβ ′ ∆ϕe 3 i /(µ i D) is the solutal Rayleigh number.

The dimensionless expression for the boundary conditions is:

z = 0 or x = 0, A i      v = 0, ∂φ s ∂n = 0 (24) z = e(t) e i              ∂u ∂z = - µ i µ M a ∂φ s ∂x , w = 0, - ∂φ s ∂z = 1 -P e int φ s , (25)
where M a = γe i ∆ϕ µ i D is the solutal Marangoni number. Each run is therefore characterized by the quintuplet Q = (ϕ si , P e int , M a, Ra, Sc). The initial volume fraction ϕ si does not appear explicitly. It is only involved in the concentration-dependent viscosity µ/µ i calculations (see the empirical law ( 6)), and in determining the validity domain of the model: t lim is the time such that ϕ s,int > 0.4 that is φ s,int > (0.4 -ϕ si )/(P e int (1 -ϕ si )) (see Eq. 20 and section II A for details).

Finally the above model is valid for t ≤ min(t lim , t lim2 ), where t lim and t lim2 are the limiting times over which the evaporative flux can be assumed constant, and the interface flat. This latter condition is discussed in details in section (IV C).

D. Numerical method

Eqs. ( 21)-( 23) are expressed using a finite volume method on moving staggered and structured grids. An arbitrary Lagrangian Eulerian method [START_REF] Hirt | An Arbitrary Lagrangian-Eulerian Computing Method for All Flow Speeds[END_REF] is used in order to follow the displacement of the interface versus time. We consider moving control volumes V u (x, t) and V w (x, t) for the components u and w of the velocity, and V p,φ (x, t) for the scalar variables, namely the pressure and the volume fraction. We also let δV u (x, t), δV w (x, t) and δV p,φ (x, t) be the boundary surfaces of the corresponding control volumes, n the unit normal vector pointing outward and v mesh (x, t) the mesh velocity. The equations governing fluid flow are then:

∂V p,φ ∂t + δV p,φ c • n dS = 0 (26) ∂ ∂t Vu u dV + δVu uc • n dS = - δVu pn • e x dS + Sc δVu µ µ i ∇u • n + ∂ ∂x (v • n) dS (27a) ∂ ∂t Vw w dV + δVw wc • n dS = - δVw pn • e z dS + Sc δVw µ µ i ∇w • n + ∂ ∂z (v • n) dS +Ra Sc Vw φ s dV (27b) ∂ ∂t V p,φ φ s dV + δV p,φ φ s c.n dS = δV p,φ ∇φ s • n dS (28) with c = v(x, t) -v mesh (x, t).
It is worth noticing that Eq. ( 26) reduces to

δV p,φ v • n dS = 0 (29)
because the volume flow rate due to the moving mesh, that is δV p,φ v mesh • n dS, is exactly balanced by (∂V p,φ )/(∂t), the change rate of the control volume V p,φ .

The equations are solved sequentially, by starting with the solute Eq. (28). The spatial and time derivatives are discretized with second order schemes. A backward-Euler scheme is used for the time derivatives. The first term of the viscous contribution in the momentum equations as well as the right-hand side of the volume fraction equation and φ s in the buoyancy force are implicitly expressed in time whereas the other terms are extrapolated with a second order Adams-Bashforth scheme. Note that the normal velocity of the mixture to the interface is zero (see Eq. ( 18) and the proof just above) as for impervious boundaries. Note too that the mass Eq. ( 29) is identical for both a moving and a frozen mesh. Consequently, the velocity and pressure fields can be uncoupled with a classical fractional-step method [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF] developed for incompressible Navier-Stokes equations on fixed grids. The discretization of the time derivative together with implicit diffusion terms give rise to Helmholtz problems.

They can be factorized into x-and z-components before being spatially discretized, leading to tridiagonal matrices. Finally, the Poisson equation involved by the pressure correction Φ = p n+1 -p n is generally solved with a direct Crout method, using a checkerboard renumbering for the cells. If the cell numbers N x and N z are too large, for example when N x = 1000 and N z = 100, an iterative BiConjugate Gradient Stabilized (BCGS) method with an incomplete LU preconditioning is preferred because of superior efficiency. Because the evolution of the interface is very slow between two successive time steps, the incomplete LU factorization does not need to be updated each time the Poisson equation is solved. This leads to substantial computational time savings. For some simulations, this decomposition is only necessary every 10, 000 time steps.

Computations have been performed with various grid sizes and time steps, essentially depending on the value of the interface Péclet number. For simulations near the thresholds, with A i = 10, and for P e int 1, the typical meshing was (N x , N z ) = (500, 50) and ∆t = 10 -4 . A finer grid, i.e. (N x , N z ) = (1000, 100), was preferred for larger interface Péclet numbers, with ∆t = 10 -4 for constant viscosity up to ∆t = 10 -6 for variable viscosity.

For simulations in the nonlinear domain with A i = 5, the meshing was (N x , N z ) = (750, 150)

and a geometric progression (with ratio 0.95) was used to allocate more nodes near the free surface. When the flow is driven by surface tension, the time step was ∆t ≃ 10 -9 . For buoyancy-driven flows we used ∆t ≤ 10 -7 . In all cases the results were verified to be relatively insensitive to any extra refinement of the grid or the time step.

III. DIFFUSIVE PROBLEM

Let us first consider the purely diffusive problem (M a = Ra = 0). Starting from a uniform solvent volume fraction ϕ si , the concentration field evolves due to solvent evaporation. With the boundary conditions and assumptions of the problem, the diffusive problem is 1D and the time-varying field φ s (z, t) satisfies the dimensionless equation:

∂φ s ∂t = ∂ 2 φ s ∂z 2 (30)
with boundary conditions ∂φs ∂z = 0 at z = 0 and -∂φs ∂z = 1 -P e int φ s at z = e(t)/e i . There is no stationary state in this evaporation problem, since the thickness and the solvent concentration evolve all through the drying process. Figure 3 then becomes nearly constant (see Appendix A for details and Figure 4 with δ such that φs(0,t)-φs[(e(t)/e i -δ(t)),t] φs(0,t)-φ s,int (t) = 0.1%).

IV. TRANSIENT CONVECTION

A. Criterion

As stated in the introduction, the formulation and the prediction of critical conditions for the onset of convection are not so clear-cut in a problem where the basic state is transient

(cf section III).
Indeed the perturbation must be amplified enough to modify the basic state significantly during the transient regime. The detection of a convective flow then depends on the sensor sensitivity and on the amplitude of perturbations. As a consequence it is not possible to determine a well defined threshold between stable or unstable configurations but rather a transition zone. This issue has been discussed in previous papers [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis[END_REF][START_REF] Touazi | Simulation of Transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF] for the thermal regime observed in domain 1 of Figure 1(a). Several methods have been compared: linear analysis on the one hand, nonlinear direct simulations on the other hand. The former was based on the frozen time approach or on a non-normal method that takes into account the time dependance of the basic state explicitly. Both non normal approach and direct simulations provide the width of the transition zone and give similar results. The same conclusion would be expected here. In this paper we use nonlinear direct simulations in order to get also the evolution of concentration and velocity fields in the nonlinear regime.

To characterize the occurrence of convection, a criterion based on the velocity has been used. For a given set of parameters Q = (ϕ si , P e int , M a, Ra, Sc), the convection will be considered as significant if, when the system is submitted to an initial perturbation of the concentration field, the mean flow velocity increases and reaches the diffusion velocity D/e i .

In dimensionless form, this condition can be written:

∃t d (Q) such that      ∀t < t d (Q), ||v(t)|| 2 < 1 t = t d (Q), ||v(t)|| 2 = 1 (31)
with ||v(t)|| 2 the L2 norm of the velocity v(t) and t d the diffusive/convective transition time.

The thresholds and critical times are defined in the following way. The critical Marangoni number M a c is taken to correspond, for a given initial perturbation with all other parameters frozen, to the smallest value of M a that fulfills condition (31). The critical time is taken to be t c = t d (M a c ) (same definition for Ra c ).

Although the choice of criterion may seem somewhat arbitrary, we have checked that for several configurations the thresholds were not very sensitive to the criterion definition. For instance threshold variations were less than 10% for criteria based, either on the L2 norm of the velocity (i.e. the criterion outlined here), or on the deviation of the concentration profile from the pure diffusive case (cf [START_REF] Trouette | Instabilités de Rayleigh-Bénard-Marangoni induites par évaporation en régime transitoire -Application aux solutions polymères[END_REF] for more details).

B. Influence of random perturbations on convection

Since the transient problem is sensitive to initial perturbations, a preliminary study was performed on two very different configurations to analyze the influence of the initial conditions on the evolution of the velocity and concentration fields. In both cases, the viscosity was assumed to depend on the solvent concentration. The patterns of real experimental perturbations not being known a priori, small perturbations were applied randomly on the uniform concentration field φ s (t = 0) as follows:

φ s (x i , t = 0) = a ∆ϕ Rand(x i ), ( 32
)
where a is the amplitude of the perturbation and "Rand(x i )" is a random draw from a uniform distribution between -0.5 and +0.5. Other choices for the spatial structure of the initial perturbation were discussed for the thermal regime (see domain 1 of Figure 1(a))

in previous papers [START_REF] Trouette | Numerical study of convection induced by evaporation in cylindrical geometry[END_REF][START_REF] Chénier | Sensitivity of diffusive-convective transition to the initial conditions in a transient Bénard-Marangoni problem[END_REF] . It has been shown that the use of a simple random perturbation, as described by Eq. ( 32), gives similar results to those obtained with more sophisticated methods which assumed spatial correlations for the disturbances.

Impact of random draws

In the first run, the amplitude a of the perturbation was fixed to 10 -6 . Simulations were performed with variable viscosity, for the following parameters: ϕ pi = 0.01, M a = 4950, P e int = 10, Ra = 0, Sc = 10. In this configuration, the Marangoni number was larger than the critical one, so that convection was observed.

For a set of 20 random draws, Figure 5 gives typical illustrations, on the one hand, of the evolution of the L2 norm of the velocity, and on the other hand, of the evolution of the difference between the spatial mean values of the solvent concentrations at the bottom and at the free surface, ∆φ s =< φ s (x, z = 0, t) > -< φ s (x, z = e(t)/e i , t) >.

Due to the evolution of the concentration field and of the film thickness, the velocity first increases up to a maximal value before decreasing and tending to zero again (Fig. 5(a)).

Convection can also be highlighted by considering the mean concentration gap between the pure diffusive reference case and the simulations (Fig. 5(b)). As can be seen, the curves deviate from each other for 0.06 t 0.085, which roughly corresponds to ||v(t)|| 2 > 1.

Moreover, it can be seen that the velocity norm and the diffusive/convective transition time do not depend much on the specific draw used for the simulation. We then use only one draw in the following. 0 0,01 0,02 0,03 0,04 0,05 0,06 0,07 0,08 0,09 time

Impact of the perturbation amplitude

The influence of the perturbation amplitude a (see relation (32)) was analyzed in a second run on another configuration defined by M a = 9.9 × 10 5 , Ra = 9.2 × 10 4 , Sc = 1.1 × 10 4 , P e int = 3 and ϕ pi = 0.006. The random draw is the same for all the simulations performed in this section. Figure 6 shows the evolution of the mean velocity for different amplitudes of the disturbance. It is evident that, for very short times t 5 × 10 -4 , the evolution depends on the perturbation amplitude. This corresponds to the linear regime in the sense that velocities scaled with respect to "a" are superposed. For larger times, that is in the nonlinear regime, the evolution of the velocity is qualitatively the same whatever the perturbation amplitude.

It can also be emphasized that the diffusive-convective transition time t d (see relation (31))

is not very sensitive to the initial perturbation amplitude. Indeed, the reduction of the perturbation amplitude by four orders of magnitude delays the transition time t d by about half a decade. This is the uncertainty range due to the lack of knowledge about initial perturbations.

C. Model validity

Some attention is now paid to determine the validity domain of our model, using the configuration of the previous section (M a = 9.9 × 10 5 , Ra = 9.2 × 10 4 , Sc = 1.1 × 10 4 , P e int = 3 and ϕ pi = 0.006). First, the evaporative flux was found to be constant as long as the solutal volume fraction was larger than 0.4. Results of the simulation therefore imply that the simulation time must not exceed t lim = 0.3 for this configuration. The second assumption of the model is related to the non-deformable interface. Assuming small deflections, the normal component of the normal stress

σ zz = -p + 2 µ µ i ∂w ∂z . ( 33
)
is balanced either by the Laplace pressure induced by the surface tension or by gravity [START_REF] Wehausen | Surface waves[END_REF] .

The Laplace pressure is given by the pressure jump across the interface

∆p L = 1 We 1 R , (34)
where We = ρD 2 /(σe i ) is the Weber number with σ being the surface tension and R the dimensionless radius of curvature. With the small slope approximation, the curvature be-

comes 1 R ∼ ∆h L λ 2 , ( 35
)
with λ and ∆h L being the dimensionless wavelength and the dimensionless interface deflection, respectively.

Thus, the deflection given by the Laplace pressure is finally expressible as:

∆h L ∼ We σ zz λ 2 . ( 36
)
If the normal component of the normal stress is now balanced by the fluid weight per unit surface, the resulting deflection is:

∆h B ∼ Fr σ zz (37)
with Fr = D 2 /(e 3 i g) being the Froude number and g the acceleration due to gravity. If we assume that the interface remains flat when ∆h L and ∆h B are smaller than 1% of the current thickness e(t)/e i , the conditions (36) and (37) provide the time t lim2 . Using parameters corresponding to the experimental configuration (cf section V E) and e i = 1mm to get the Weber and Froude numbers, and estimating σ zz and λ from numerical simulations, we get t lim2 ≃ 0.15, about half t lim . Notice however that the transition time is t d ≃ 10 -4 ≪ min(t lim , t lim2 ) which confirms the validity of our model at the onset of convection.

To sum up the main results of this section, we can say that the amplitude of the initial disturbance affects the transition time t d , but not the general trends of the flow behavior in the nonlinear regime. Furthermore, close to the transition (i.e. for dimensionless velocity ∼ 1), the validity of the assumption of non deformable interface has been checked a posteriori for all configurations. Thus it is meaningful to perform a systematic analysis of the critical parameters, as defined by relations (31). This is the purpose of the next section.

V. THRESHOLDS

A. Bénard-Marangoni convection, Ra = 0

We first consider pure Bénard-Marangoni convection with Ra = 0. In this transient problem, the perturbation growth depends in a non trivial way on the coupling between the time evolutions of the concentration field and layer thickness, which leads to the critical Marangoni and critical time behaviors depicted in Figure 7. However the influence of viscosity can be studied on its own, since simulations can be performed with constant or concentration-dependent viscosities. Figure 7 The results of Figure 7(a) highlight the significant part played by the viscosity at large values of P e int . Indeed the use of a concentration-dependent viscosity can increase the critical Marangoni number by more than an order of magnitude. In agreement with the previous results, it can be seen that changing the amplitude of the perturbation by several orders of magnitude induces only small modifications of the threshold values. The dispersion around the critical value achieved with a = 10 -6 is about 15%. This "blurring" effect, inherent to any transient problem, allows us to define a transition region between stable and unstable diffusive configurations. This transition zone is rather narrow compared to the extent of the Marangoni and Rayleigh variations covered by the different simulations when varying the initial thickness, the evaporative flux (both included in P e int ) and the initial solvent concentration. Note that similar conclusions were already obtained for the thermal transient regime that occurs at the beginning of the drying [START_REF] Doumenc | Transient Rayleigh-Bénard-Marangoni Convection due to Evaporation : a Linear Non-normal Stability Analysis[END_REF][START_REF] Touazi | Simulation of Transient Rayleigh-Bénard-Marangoni convection induced by evaporation[END_REF] .

Figure 7(b) provides an evaluation of the critical time t c as well as the limit time t lim (such that ϕ s (t lim ) = 0.4, cf II C). It appears clearly that the critical time decreases when P e int increases. The maximal value of P e int that can be explored for a constant viscosity is The circles correspond to simulations performed with a constant viscosity and the triangles with a variable viscosity. Parameters: ϕ pi = 0.047, Ra = 0, Sc = 10. about 10. Indeed for P e int ≃ 10, t c is close to t lim , which means that convection occurs just before ϕ s,int = 0.4.

B. Rayleigh-Bénard convection, M a = 0

Next, a study was performed for pure Rayleigh-Bénard convection (M a = 0) at the same initial concentration and Schmidt number as in the previous section (with ϕ pi = 0.047 and Sc = 10). Figure 8 gives the critical Rayleigh number Ra c and critical time t c as a function of P e int . The circles correspond to simulations performed with a constant viscosity and the triangles with a variable viscosity. Parameters: ϕ pi = 0.047, M a = 0, Sc = 10.

The conclusions are globally similar for the two configurations (Ra = 0 or M a = 0).

Taking into account the variable viscosity leads to comparable relative effect on the thresholds: at large P e, Ra c or M a c are about 1.5 to 2 orders of magnitude higher for variable viscosity than for constant one. As shown and explained in section V D here below, this strong increase of M a c (resp. Ra c ) when P e int > 1 and for concentration-dependent viscosities is due to the high value ϕ pi = 0.047. It impacts a little more the critical time in the Bénard-Marangoni configuration than in the Rayleigh-Bénard one. Indeed Marangoni convection is driven by surface phenomena. It can no more occur if a very viscous skin has formed on the surface, which explains that it should start earlier if the increase of viscosity with polymer concentration is involved in the equations.

C. Influence of the Schmidt number, Ra = 0 Simulations were performed for different values of the Schmidt number, with a concentrationdependent viscosity, ϕ pi = 0.047, Ra = 0, P e int = 1 and a = 10 -6 . As illustrated in Figure 9, the critical Marangoni number decreases from about 30% between Sc = 1 and Sc = 10, and then it maintains a nearly constant value. This result is in accord with previous work done on the transient thermal regime (for fluids with P r > 1) where the Prandtl number plays the part of the Schmidt number 11 .

D. Influence of the solute initial concentration, Ra = 0

As mentioned in section II, the problem depends on the initial concentration through Eq. 6, when the variation of the viscosity with concentration is taken into account.

In Figure 10, thresholds and critical times are compared for an initial polymer volume fraction varying from 10 -4 to 10 -1 with Ra = 0, Sc = 10 and a = 10 -6 . At high P e int the thresholds depend strongly on the initial polymer volume fraction through the function µ(ϕ p ) (Eq. 6). Indeed the increase in viscosity is especially felt when the polymer volume fraction has significantly changed. Let us first consider the three lower initial polymer volume fractions. The critical times are very close, which corresponds to the same scaled volume fraction φ p at the onset of convection (with φ p = (ϕ p -ϕ pi )/∆ϕ = -φ s ). Indeed, for t < t c , the problem is mainly driven by diffusion and φ p is then described by Eq. 30, whatever the value of ϕ pi . The solution viscosity depends on ϕ p , with ϕ p -ϕ pi = φ p P e int ϕ pi .

For the same φ p , variations of ϕ p and then of the viscosity are larger for high P e int and ϕ pi , which explains the increase of M a c with the two parameters. For the highest initial polymer volume fraction (ϕ pi = 0.1), simulations show that the only way to get convection is to have a much larger M a c which induces a smaller critical time t c . Parameters: Ra = 0, Sc = 10, a = 10 -6 , variable viscosity.

E. Simulation and experiments

In this section, we turn to a consideration of the experimental parameter plane (initial thickness -initial concentration). For the PIB/toluene drying experiment, the evaporation rate was v evap = 3 × 10 -7 m/s. The viscosity of the solution is given by relation ( 6 The critical thickness for different initial polymer volume fractions in the Bénard-Marangoni ("BM") and Rayleigh-Bénard ("RB") problems is presented in Figure 11.

Unlike the convection induced by thermal gradients, each marginal curve obtained in the solutal regime exhibits a minimum. Indeed when ϕ p goes to zero (pure fluid), solutal convection can no longer be initiated. For polymer volume fraction larger than 0.1 the viscosity strongly increases (Eq.2) and so does the critical thickness. The minimum is encountered for a ϕ p value of about 0.01 for the two mechanisms (buoyancy or surface tension driven convection), but for very different thicknesses. The thickness values were approximately 10 µm and 1 mm for BM and RB convection, respectively. Surface tension instabilities are clearly dominant and can be observed for a thickness as small as 10 µm.

Considering the experimental evaporation rate, P e int varies from 0.03 to 3 for thicknesses 10 µm and 1 mm respectively. For the small thicknesses, corresponding to small P e int and Bénard-Marangoni convection, the thresholds are not very sensitive to the viscosity as shown in Figure 7(a). Consequently, the results obtained with constant or variable viscosities are similar. The effect is more important for larger thicknesses corresponding to Rayleigh-Bénard convection and for larger P e int (cf. Figure 8(a)).

For small polymer volume fraction (ϕ pi ≤ 0.001), asymptotic behaviors can be inferred from numerical results. For this purpose, viscosity µ i can be assumed to be constant and not to depend significantly on ϕ pi . From Figure 10(a) we notice that M a c (P e int ) grows by less than one order of magnitude whereas P e int increases by three orders of magnitude. Thus, we can consider the critical Marangoni to be roughly constant, with a value of M a c ≃ 100.

We can then deduce an estimation of the critical thickness for dilute solutions:

e 2 i ≃ 10 2 × µ i D 2 ϕ pi γ v evap (38)
It is worthwhile comparing this expression with the scaling law obtained by de Gennes [START_REF] De Gennes | Instabilities during the evaporation of a film. Non-glassy polymer + volatile solvent[END_REF] which, with our notations, reduces to e 2 i ∼ µ i D 2 γ vevap where the initial polymer volume fraction ϕ pi does not appear. This discrepancy can be explained by the boundary condition used in de Gennes scaling analysis, -D ∂ϕs ∂z = v evap , which does not apply to evaporation of dilute solutions. Indeed, using boundary condition 19 in de Gennes analysis leads to Equation 38, which successfully predicts that the critical thickness goes to infinity in the limit of pure solvent.

For larger thicknesses and Rayleigh-Bénard instability, the same approach can be developed, even if the approximation Ra c ∼ cste is somewhat questionable. Using Ra c ≃ 10 3 (Figure 8(a)), we get:

e 4 i ≃ 10 3 × µ i D 2 ϕ pi g ρ β ′ v evap (39)
Both scalings for low solute concentrations are reported in Figure 11. Bassou and Rharbi 20 recently performed experiments in a configuration rather close to ours (drying of PS/Toluene solutions with v evap = 3.5 × 10 -7 m/s). They observed solutal Bénard-Marangoni convection with a critical thickness of the order of 100 µm, a value which is consistent with the above results. New simulations with the PS/Toluene viscosity law would be required in order to perform more quantitative comparisons, but this work is beyond the scope of the present paper.

VI. NON LINEAR SIMULATIONS

This last section is dedicated to the analysis of the nonlinear regime and addresses the question of the competition between buoyancy and Marangoni effects in this regime. The configuration is a 8 mm thick layer of PIB/Toluene, with ϕ pi = 0.047 and v evap = 3 × 10 -7 m/s (red star in Figure 11). The scalings e 2 i /D and D/e i used for non dimensional time and velocity are then 64 × 10 4 s and 0.125 × 10 -7 m/s, respectively. The corresponding Ra and M a numbers are 1.5 × 10 8 and 2.5 × 10 7 , respectively. This configuration is unstable with respect to both mechanisms, as stated in Figure 11. To decrease the computational time, simulations were performed with an initial aspect ratio A i = 5. Moreover, the Schmidt number used in the simulations is set to 10 3 while the real one for this configuration is 2.3 × 10 5 . Indeed we have checked that the results do not depend on the Schmidt number for values larger than 10 3 .

The two mechanisms are studied separately (M a = 0 or Ra = 0) or together (M a = 0 and Ra = 0). The main difference between the configurations is due to the presence, or absence, of the Marangoni effect. Indeed the results corresponding to the RBM configuration (M a = 0 and Ra = 0) and to the pure BM configuration (Ra = 0) are similar, while the pure RB configuration (M a = 0) exhibits a very different behavior (cf Figure 12 and13). Convection occurs in the first seconds when surface tension effects are present, while pure buoyancy convection appears only about 45 minutes after the beginning of the drying (t ≃ 0.004), as shown for instance in Figure 13. This is related to the distance of the studied configuration to the thresholds. Indeed it is not very far from the RB threshold, so that the time of occurrence of convection is close to the critical time (estimated by extending Figure 8(b) to P e int = 24). On the contrary the configuration is far away from the BM threshold, which explains that convection occurs very early when surface tension is taken into account. For The differences between the configurations is further illustrated by the observation of the convective structures evolution during the drying. For the two configurations with surface tension contribution (M a = 0) convection occurs at t ∼ 10 -5 , which corresponds to a penetration depth of the perturbation of about 1% and to only a small increase in the surface polymer concentration (cf Figures 3 and4). The surface viscosity is then close to the initial one and small convective structures first appear near the surface. These plumes then penetrate into the bulk, improving the mixing between the solute and the solvent and the whole layer is progressively overrun, as shown in Figure 14 (and corresponding animations enhanced on line). Once again it can be seen that the additional effect of buoyancy is weak, and results mainly in the reduction of the time needed by the perturbation to reach the bottom of the layer. Indeed, as can be seen in Figure 14, plumes have almost reached the bottom of the layer at t = 7.5 × 10 -4 (480 s) for the RBM configuration, while only about one half of the layer is concerned for the BM configuration. The differences in the structure morphologies observed at t = 1.1 × 10 -3 are not really significant as a factor of two in wavelengths can easily be obtained by changing a little the initial perturbations (see for instance [START_REF] Chénier | Sensitivity of diffusive-convective transition to the initial conditions in a transient Bénard-Marangoni problem[END_REF] ).

But when surface tension is not active (M a = 0), the behavior is different. Convection occurs much later, as already seen in Figure 12. At t = 3.5 × 10 -3 (about 37 min) the perturbation has overrun about 20% of the thickness and a strong concentration gradient, corresponding to a viscous skin, has formed at the surface (cf Figures 3 and4). Indeed the viscosity at the surface is 500 times larger than the one in the bulk (cf Figure 15 and corresponding animations enhanced on line). Small plumes appear under this viscous skin.

Then convective structures grow up and carry viscous fluid from the surface to the bulk and vice versa.

It is worthwhile to zoom on the upper surface to compare the evolution of the viscosity for the different configurations. For the RBM configuration, the viscosity of the surface increases slowly during the drying. For instance, at t = 0.01, i.e. about 107 min after the beginning of the drying, the maximal viscosity is about three times higher than the initial one (see spatio-temporal diagram in Figure 16(a)). Quite the reverse can be observed for the RB configuration: at t = 0.01 the upper surface morphology is very different with strong concentration gradients. Indeed part of the surface is very viscous (viscosity about five orders of magnitude larger than its initial value, the local velocity is very small in this part of the surface), while convection has succeeded to form holes in the crust and to carry

  (a) gives an illustration of the solvent evaporative flux obtained for a 4 mm thick polymer solution. As indicated on the Figure, three different successive regimes are observed. First, the solvent evaporative flux decreases due to the thermal transient regime (domain 1 in Figure 1(a)). During this first regime, buoyancy or surface-tension thermal-driven convection can be observed, depending on the initial thickness and viscosity. As stated in the introduction, this transient thermal regime has

Figure 1 .

 1 Figure 1. (a) Solvent evaporative flux for a 4 mm thick layer solution of polyisobutylene/toluene, at room temperature in an extractor hood. The initial polymer volume fraction is ϕ pi = 0.05. 1=thermal regime, 2=solutal regime with constant evaporative flux, 3=end of the drying. (b) Solvent activity as a function of polymer volume fraction ϕ p . Straight line: Flory-Huggins model (Eq 3); diamonds: experimental data from gravimetric experiments 21 .

Figure 2 .

 2 Figure 2. Dynamic viscosity µ as a function of polymer volume fraction ϕ p at room temperature. Symbols: experimental points 22 , dashed line: interpolation curve, Eq. 6.

Figure 3 .Figure 4 .

 34 Figure 3. Profile of the scaled solvent volume fraction in the layer at different times for the diffusive problem, with P e int = 24.

Figure 5 .

 5 Figure 5. Initial perturbation -Influence of the random draw on (a) L2 norm of the velocity (mean value ± 0.5 standard deviation for 20 simulations), (b)∆φ s =< φ s (z = 0, t) > -< φ s (z = e(t)/e i , t) >. Parameters: ϕ pi = 0.01, M a = 4950, P e int = 10, Ra = 0, Sc = 10 , a = 10 -6 , variable viscosity.

10 Figure 6 .

 106 Figure 6. Evolution of ||v|| 2 over time for different amplitudes of the intial perturbation. Parameters: M a = 9.9 × 10 5 , Ra = 9.2 × 10 4 , Sc = 1.1 × 10 4 , P e int = 3 and ϕ pi = 0.006, variable viscosity.

  (a) and 7(b) gives the critical Marangoni number M a c and the associated critical time t c as a function of P e int , for ϕ pi = 0.047 and Sc = 10, for constant and concentration-dependent viscosities and for three amplitudes of the solutal perturbation applied at initial time (see relation (32)).

Figure 7 .

 7 Figure 7. M a c (a) and t c (b) as a function of P e int for three values of the perturbation amplitude.

Figure 8 .

 8 Figure 8. (a) Ra c and (b) t c as a function of P e int for three values of the perturbation amplitude.

Figure 9 .

 9 Figure 9. M a c as a function of Sc. Parameters: ϕ pi = 0.047, Ra = 0, P e int = 1, a = 10 -6 , variable viscosity.

Figure 10 .

 10 Figure 10. (a)M a c and (b) t c as a function of P e int for different initial solvent concentrations.

  ) and the following values have been used for the other properties: ρ 0 = 0.867 × 10 3 kg/m 3 , σ 0 = 28 × 10 -3 N/m, β = 5.82 × 10 -2 , γ = 5.4 × 10 -3 N/m, D = 10 -10 m 2 /s.

Figure 11 .

 11 Figure 11. Critical thickness as a function of the initial polymer volume fraction for PIB/toluene drying experiments. The evaporation rate is v evap = 3 × 10 -7 m/s.

Figure 12 .

 12 Figure 12. Evolution of ||v|| 2 over time at the free surface ("fs", z = e(t)/e i ) and in the middle of the layer ("ml", z = e(t)/(2e i )) for BM, RB and BM+RB configurations. Parameters: M a = 2.5×10 7 , Ra = 1.5 × 10 8 , Sc = 10 3 , P e int = 24 and ϕ pi = 0.047, variable viscosity.

Figure 13 .

 13 Figure 13. Mean concentration difference between the surface and the bottom, ∆ϕ s =< ϕ s (z = 0, t) > -< ϕ s (z = e(t)/e i , t) >, for BM, RB and BM+RB configurations. Parameters: M a = 2.5 × 10 7 , Ra = 1.5 × 10 8 , Sc = 10 3 , P e int = 24 and ϕ pi = 0.047, variable viscosity.

Figure 16 .

 16 Figure 16. Viscosity field at the upper surface for the RBM (top) and RB (bottom) configurations.Each horizontal line corresponds to the profile µ(x, t)/µ i at a given time. Parameters: M a = 0 or M a = 2.5 × 10 7 , Ra = 1.5 × 10 8 , Sc = 10 3 , P e int = 24 and ϕ pi = 0.047, variable viscosity.

ACKNOWLEDGMENTS

The authors gratefully acknowledge the financial support of the French Fédération "Transferts de Masse et de Chaleur". F.D and B.G also gratefully acknowledge the financial support from the MULTIFLOW Marie Curie FP7-ITN Network and have benefited from fruitful discussions with P. Colinet and A. Rednikov from the MULTIFLOW network and M.Rossi from IJLRA (UPMC).

up less viscous bulk fluid in some places where the velocity is a few µm/s. Note that the evaporation rate is imposed a priori and is the same for the two configurations. In both cases the thickness has decreased of about 2mm at t = 0.01. The mean concentration in the whole layer is then the same, while the surface behavior is very different.

As a conclusion, for the polymer solution studied in this paper, Marangoni effects are clearly dominant even for rather thick layer (several millimeters). Only a strong decrease of the surface tension dependence with concentration would modify this behavior. It would be interesting, but beyond the scope of the present paper, to shift to reverse configurations when surface tension becomes a stabilizing effect and thus may compete with buoyancy.

VII. CONCLUSION

The onset of solutal convection for a transient drying problem has been studied. A model was developed taking into account the variation of the viscosity with the solute concentration and the decrease of the film thickness due to evaporation. The model assumed a constant evaporative flux and a non-deformable interface, and the validity domain of these assumptions was numerically investigated.

In order to describe the onset of convection, a criterion based on the mean velocity was defined and discussed. The impact of the amplitude of a random initial perturbation was studied and the sensitivity to initial conditions found to lead to a small transition region where the fluid flow changes from a diffusive to a convective regime. The dependence of critical parameters on the Péclet number has been estimated for the limiting cases of purely

Appendix A: Appendix A -Scaling analysis of the basic state

As long as the concentration at the bottom of the layer (z = 0) is close to the initial concentration, it is possible to consider the layer as a semi-infinite medium. By using the distance from the interface, i.e. y = e(t) -z, as the new space variable, we obtain the following set of equations (in dimensionless form):

Initial condition:

Boundary conditions:

∂φ s ∂y = 1 -P e int φ s at y = 0 (A3)

Let δ(t) denote the distance from the interface where the concentration varies. Using this quantity in a scaling approach, Eq. (A1) then states the balance between the three following terms :

• P e int 1

If P e int < 1, the second term is always smaller than the third (δ being always ≤ 1).

The mobility of the interface can then be neglected and the problem boils down to a classical diffusion problem with an evolution of the perturbation thickness as the square root of time: δ ∼ √ t, up to δ ∼ 1.

• P e int 1

At the beginning, the diffusion (third term) is balanced by the temporal derivative (first term) and as before one obtains a classical diffusion problem with an evolution of the perturbation thickness as the square root of time. Then the second term due to the moving interface becomes dominant and the thickness saturates to a value δ max .

We get the following orders of magnitude: