- Email: [email protected]
Improved 1.5-sided model for the weakly nonlinear study of Bénard–Marangoni instabilities in an evaporating liquid layer
Journal of Colloid and Interface Science 290 (2005) 220–230 www.elsevier.com/locate/jcis
Improved 1.5-sided model for the weakly nonlinear study of Bénard–Marangoni instabilities in an evaporating liquid layer J. Margerit ∗ , M. Dondlinger, P.C. Dauby Institut de Physique B5a, Université de Liège, Allée du 6 Août 17, B-4000 Liège, Belgium Received 18 November 2004; accepted 9 April 2005 Available online 1 June 2005
Abstract Bénard–Marangoni instability, with coupled gravity and surface tension effects, in an evaporating liquid layer surmounted by its vapor and an inert gas is investigated theoretically. We show that this system can be described by a model that consists in the liquid layer equations plus the diffusion equation for the vapor in the gas (the so-called 1.5-sided model) and that this model is equivalent to a one-sided model when the vapor mass fraction field can be considered as quasi-stationary, provided that the equivalent Biot number is a nonlocal function of the interface temperature. A comparison of weakly nonlinear results for the 1.5-sided model with a previous one-sided model [M. Dondlinger, J. Margerit, P.C. Dauby, J. Colloid Interface Sci. 283 (2) (2005) 522–532] that considered a Biot number depending on the wavenumber evaluated at the threshold is performed. Very good agreement is found between both models. For this reason, the present analysis can also be considered as a detailed theoretical justification for the use of a one-sided model in the study of evaporative thermoconvection. 2005 Elsevier Inc. All rights reserved. Keywords: Marangoni convection; Evaporation; Pattern formation; Amplitude equations
1. Introduction Thermoconvection is a standard example of self-organization phenomena. Bénard–Marangoni instability appears in horizontal fluid layers when a sufficiently large vertical temperature gradient is imposed. Arrangements of rolls, hexagons, or square convective patterns are then observed [1–6]. The temperature gradient can be produced by external heating of the bottom plate on which the liquid is lying, but if the liquid is volatile enough, the latent heat consumption caused by evaporation at the interface can lead to a vertical temperature gradient that is sufficient to induce instability. This is equivalent to considering a liquid layer cooled from above. The thermoconvective motion that appears gives rise to an increase of the evaporation rate. Evaporative convection plays an important role in heat exchangers, distillation columns, and drying technologies. * Corresponding author. Fax: +3243 66 2335.
E-mail address: [email protected] (J. Margerit). 0021-9797/$ – see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2005.04.031
Experimental results can be found in [7–9]. The qualitative results at Berg et al. [7] have revealed a variety of convection patterns for different fluid depths and Mancini and Maza [9] observed the successive transitions from square to hexagonal patterns and to conductive state for hexamethyldisiloxane (C6 H18 OSi2 ) when the liquid depth decreased progressively due to evaporation. Unfortunately properties of this fluid are not well enough tabulated to permit a quantitative comparison with our theoretical model. Colinet et al. describe in [8] the ITEL experiment, which made it possible to visualize in a microgravity environment evaporative chaotic patterns in ethanol. However, due to the short microgravity time of the experiment, only transient effects have been investigated. From a theoretical point of view, the instability of a liquid layer that evaporates in an inert gas has been investigated for linear stability by Ha and Lai [10] and by Colinet et al. [11], Colinet et al. [8], and Haut and Colinet [12]. Moussy et al. [13] have studied the role of interfacial deformation in the same problem. Under some assumptions, the system can be described by a one-sided model that introduces an equivalent Biot number to take into account the gas–vapor mixture
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
phase and the evaporation at the interface. In the different papers by Colinet and co-workers [8,11,12] a Biot number depending on the wavenumber is derived in the linear analysis. Its use has been extended by Dondlinger et al. [14] in the Galerkin–Eckhaus weakly nonlinear analysis of the system, with a Biot number depending on both the wavenumber evaluated at the threshold and the supercritical bottom plate temperature. The nonlinear study of such a model with a constant Biot number has also been carried out by Merkt and Bestehorn [15] thanks to a 3D numerical simulation. The purpose of the present work is to analyze the validity of the assumption of a Biot number, either considered as a constant or depending on the critical wavenumber. This will be done by deriving an improved model, the so-called 1.5-sided model, that consists of the liquid layer equations plus the diffusion equation for the vapor in the gas. A linear and a weakly nonlinear analysis of this 1.5-sided model will be presented and the results will be shown to be in very good agreement with the results of the simplified Biot approximation model presented in [14], for both the linear and nonlinear approaches. For this reason, the present analysis must be considered as a strong theoretical justification of the Biot approximation. The structure of the paper is the following. In Section 2, the physical system, the set of governing equations, and some dimensionless physical parameters are introduced. Section 3 is dedicated to the derivation of the 1.5-sided model in a configuration where thermal transport in the gas can be neglected, and to its equivalence to a one-sided model with an equivalent Biot number that is a nonlocal function of the interface temperature when the vapor mass fraction field can be considered as quasi-stationary. In the context of a linear analysis, this 1.5-sided model is equivalent to the model with a Biot number depending on the wavenumber derived in [11] and the corresponding results are compared in Section 4 with those of the complete two-sided problem and with those obtained with a constant Biot number [15]. Section 5 deals with the weakly nonlinear study of this 1.5-sided model. The influence of gravity is considered too. Finally, conclusions are drawn in Section 6.
2. Problem formulation We consider an evaporating liquid layer of infinite horizontal extent and depth dl , surmounted by a gas layer of thickness dg at a time t0 (Fig. 1). The gas layer is a mixture of the vapor of the liquid and an inert gas. Subscripts l, g, v, and i describe respectively the liquid, gas, vapor, and interface. The lower, rigid, and perfectly heat-conducting plate is maintained at the temperature Te . The gas layer is bounded by a perfectly heat-conducting upper plate maintained at the same temperature Te and at constant pressure and concentration. We assume that the Boussinesq approximation is valid and that the gas may be taken as perfect. The interface is assumed to be nondeformable so that the moving horizontal
221
Fig. 1. System under study.
liquid–gas interface is described by the equation z = h(t). Contrary to what was assumed in [14], we do not neglect gravity and the mass density ρ is assumed to be a linear function of temperature T , ρ = ρ0 1 − αT (T − Te ) , where αT is the constant coefficient of volumic expansion. The surface tension σ decreases linearly with the temperature, i.e., σ = σ0 − σT (T − Te ), where σT is a constant coefficient. All other physical properties of the fluids are considered constant and are evaluated at temperature Te . The properties of the fluids considered in the numerical calculations are taken from [16,17]. The physical properties of the gas mixture are also calculated at the reference temperature Te and at the reference mass fraction Yv,t of the vapor at the top plate by assuming an ideal gas mixture. As shown in [14], it is possible to assume that fluctuations in the system evolve on a much faster time scale than the displacement of the interface due to evaporation. Therefore the liquid layer thickness will be frozen during the stability analysis. The variables are expressed in dimensionless form: lengths (coordinates (x, y, z)) are scaled by the liquid layer thickness dl , time t by dl2 /κl , with κl the heat diffusivity of the liquid, velocity V = (U, V , W ) by κl /dl (vectors are denoted by bold letters), pressure P by ρ0,l κl2 /dl2 , the mass flux J at the interface by ρ0,l κl /dl , temperature T by L/cp,l where L is the latent heat of vaporization and cp,l the liquid heat capacity per unit mass, and the mass fraction of vapor in the gas Yv by 1 (Yv is already dimensionless). The dimensionless equations and boundary conditions governing the system are standard [14] and given as follows: in the liquid layer, ∇.Vl = 0, ∂t Vl + Vl .∇Vl + ∇Pl − Pr ∇ 2 Vl − Pr Ra Tl ez = 0, ∂t Tl + Vl .∇Tl − ∇ 2 Tl = 0;
(1) (2) (3)
in the gas layer, ∇.Vg = 0,
(4)
1 µ ∂t Vg + Vg .∇Vg + ∇Pg − Pr ∇ 2 Vg ρ ρ − ραT Pr Ra Tg ez = 0,
(5)
222
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
∂t Tg + Vg .∇Tg − κ∇ 2 Tg = 0, κ ∂t Yv + Vg .∇Yv − ∇ 2 Yv = 0; Le
(6)
3. Derivation of the improved 1.5-sided model
(7)
In this section we derive the 1.5-sided model, which is a simplification of the two-layer model, and we use this new model to derive an improved one-sided model, which is characterized by an equivalent Biot number that is a nonlocal function of the interfacial temperature. The one-sided model used in [15] (constant Biot number) will be recovered as a local approximation of the present improved one-sided model. Note also that, for the linear analysis, the 1.5-sided model is equivalent to the one used by Colinet et al. [11] and Dondlinger et al. [14] (Biot number depending on the wavenumber), but it improves the nonlinear analysis in the sense that it makes it possible to get rid of the assumption of a Biot number evaluated at the critical wavenumber.
at the bottom plate (z = 0), Ul = Vl = Wl = 0,
Tl = Te ;
(8)
at the top (z = H ), Vg = (0, 0, C),
Tg = T e ,
Yv = Yv,t ,
Pg = Pt ;
(9)
and at the interface (z = h(t)), J = Wl − Wi ,
(10)
Wl − Wi = ρ(Wg − Wi ),
(11)
Ul = Ug ,
(12)
3.1. The 1.5-sided model
Vl = Vg ,
−µ(∂x Wg + ∂z Ug ) + ∂x Wl + ∂z Ul + Ma ∂x Tl = 0, −µ(∂y Wg + ∂z Vg ) + ∂y Wl + ∂z Vl + Ma ∂y Tl = 0,
(13)
Tl = T g = T i ,
(14)
J − λ∂z Tg + ∂z Tl = 0,
(15)
psat (Tl ) −
Pt 1+
1−Yv rw Yv
= 0,
(16)
J κ ∂z Yv + (1 − Yv ) = 0, Le ρ
(17)
where Pg and Pl are respectively the gas and liquid dynamical pressures, C is a constant with respect to space, Wi = ∂t h is the interface velocity, and psat (T ) is the dimensionless liquid–vapor saturation curve. Yv,t and Pt are respectively the vapor mass fraction and the total thermodynamic pressure imposed at time t0 at the top of the system, and rw = Mg /Ml is the ratio of the inert gas and liquid molecular weights. The dimensionless parameters that appear in the above equations are defined in Table 1 with ν the kinematic viscosity, g the acceleration due to gravity, and Dg the mass diffusivity of the liquid vapor in the gas. These equations are the same as in [14], except that gravity is added. Moreover, we use Eq. (16) instead of the Hertz– Knudsen relation, because it was shown in [12,14,18] that the interface can be considered to be at thermodynamic equilibrium.
The set of equations in the liquid phase are coupled with the equations in the gas phase via the interfacial boundary conditions (10)–(17). The gas viscous stress in Eq. (13) can be neglected with respect to the liquid viscous stresses, as the viscosity ratio µ is small. As a consequence, Eq. (13) does not imply coupling between the two layers. We also assume that the gas layer is passive, i.e., that convection in the gas is negligible. This is true provided the convective terms of relations (6) and (7) are small as regards corresponding diffusive terms. With gas velocity Wg order of magnitude given by J /ρ, according to interfacial relations (10)–(11) where we have considered that ρ is small, the criteria for the gas layer to be passive, deduced from Eqs. (6) and (7), are J dg /dl κρ and J dg /dl κρ/ Le. Moreover, as the gas conductivity is small with respect to the liquid conductivity, the gas phase heat flux contribution in Eq. (15) is small compared to the other two terms, provided that the depth of the gas layer is not too small with respect to the liquid layer depth (dg /dl λ). We neglect this term in the following. The different assumptions are listed in Table 2 for water and ethanol and we have checked that for the situations analyzed in the following, these conditions are well satisfied. However, Eq. (15), expressing the interfacial energy conservation, still implies coupling between the liquid and gas phases through the term describing the heat consumed by evaporation. Indeed, this term is a function of the evaporation rate J that is given by the nonsolubility property of the
Table 1 Dimensionless parameters ν
Prandtl number
Pr = κl l
Rayleigh number
Ra = ν Tκ,lc l l l p,l σ dL Ma = κ Tµ cl
Marangoni number Lewis number Depth ratio
gα
d3L
l l p,l
κ Le = Dg g d +d H = ld g l
α
Volumetric expansion ratio
αT = αT ,g T ,l
Density ratio Diffusivity ratio
ρ = ρg l κ κ = κg
Viscosity ratio
µ = µg l
Conductivity ratio
λ = λg l
ρ
l
µ
λ
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
3.2. The one-sided model with a nonlocal Biot number
Table 2 Physical properties
µ κρ κρ/ Le λ
Water at 320 K
Ethanol at 250 K
Condition
0.033 0.18 0.21 0.042
0.0051 0.014 0.09 0.011
1 J dg /dl J dg /dl dg /dl
inert gas in the liquid (17), so that to find the evaporation rate, the mass fraction field Yv must be known. Relation (10) gives the liquid velocity Wl as a function of the evaporation rate J too. Fortunately, provided that the gas phase advection is negligible, as we have assumed above, the diffusion equation for the vapor in the gas (7) and the corresponding boundary conditions (9) and (16) are written as the linear equations Le ∂t Yv , in h(t) < z < H, κ Yv = Yv,t , at z = H,
(18b)
Yv = Yi (Tl ),
(18c)
∇ 2 Yv =
at z = h(t),
where
Yi (Tl ) = 1 − rw + rw
Pt psat (Tl )
(18a)
−1 .
(19)
After this Eq. (18) is solved, the evaporation rate is easily deduced from relation (17). Finally, our system is described by the above set of Eqs. (18), (19) in the gas layer and by the following equations in the liquid layer and the associated boundary conditions: in the liquid layer (0 < z < h(t)), ∇.Vl = 0,
(20)
∂t Vl + Vl .∇Vl + ∇Pl − Pr ∇ Vl − Pr Ra Tl ez = 0,
(21)
∂t Tl + Vl .∇Tl − ∇ 2 Tl = 0;
(22)
2
at the bottom plate (z = 0), Ul = Vl = Wl = 0,
223
Tl = Te ;
(23)
at the interface (z = h(t)), Wl = J + Wi ,
We now demonstrate that the model presented in the previous section is equivalent to a one-sided model with a Biot number that is a nonlocal function of the interfacial temperature. In order to avoid too complicated calculations that would introduce an additional memory effect in the Biot number, we assume that the vapor mass fraction field is quasi-stationary, which amounts to neglecting the time derivative in Eq. (18a). This hypothesis is valid provided the time scale dg2 /Dg of vapor diffusion in the gas is very small with respect to the time scale dl2 /κl of thermal diffusion in the liquid,√i.e., provided that the gas layer is not too thick (dg /dl κ/ Le). First, let us show that the evaporation rate can take the form of an expression depending on a convolution product of the interfacial temperature field. As the problem is unbounded in the horizontal direction, we take a Fourier transform of system (18) in the horizontal direction and we get 2 D − k 2 Yˆv = 0, for h(t) < z < H, (28a) at z = H, Yˆv = Yv,t δ(k), (28b) ˆ ˆ Yv = Yi Ti (x, t) , at z = h(t), (28c) where D stands for d/dz, ˆ· denotes the Fourier transform of function ·, k = kx ex + ky ey is the wavenumber, δ is the Dirac function, x = xex + yey , and t is the time. The solution of Eq. (28a) has the structure Yˆv (z) = A sinh (z − H )k + B cosh (z − H )k , (29) where k is the norm of vector k and the two constants A and B are deduced from the two boundary conditions (28b) and (28c): B = Yv,t δ(k), Yˆi (Ti (x, t)) − cosh[(H − h)k]Yv,t δ(k) . A=− sinh[(H − h)k]
(24)
∂y Wl + ∂z Vl + Ma ∂y Tl = 0,
(25)
∂z Tl = −J, ρκ ∂z Yv J =− . Le 1 − Yv
(26) (27)
Note that the mass conservation Eq. (20) and relation (24) lead to the compatibility relation Wi = −J , where · is a horizontal average, which makes it possible to determine Wi at any time. We call this simplified model the 1.5-sided model since only part of the fundamental conservation equations in the gas are taken into account.
(30b)
To determine the mass flux at the interface, we can calculate the derivative of the known function Yˆv (z) at z = h: −D Yˆv (h) =
∂x Wl + ∂z Ul + Ma ∂x Tl = 0,
(30a)
k{cosh[(H − h)k]Yˆi (Ti (x, t)) − Yv,t δ(k)} . sinh[(H − h)k]
(31)
Define now the two functions ha and hb as the inverse (H −h)k Fourier transforms of respectively sinh[(H −h)k] and cosh[(H − h)k]. The expressions of the functions ha and hb are π/2 1 , ha = (32a) (H − h) 1 + cosh[π x ] (H −h)
hb =
1 (H − h)
+∞ n=0
(−1)n 2n!
δ(x)(2n) ,
(32b)
where ·(2n) stands for the derivative of order 2n. The function ha is plotted in Fig. 2. It is worth noting that Yv,t is
224
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
Consequently the reference evaporation rate (36) is a function of the interfacial temperature only, given by Jref (Ti ) = Biref (Ti )[Ti − Tv,t ],
Fig. 2. Function ha .
+∞ independent of x and −∞ ha (s) ds = 1, so that ha ∗ Yv,t = +∞ Yv,t −∞ ha (s) ds = Yv,t , where ∗ is the convolution product. The inverse Fourier transform of Eq. (31) can be rewritten as
−(H − h)DYv (h) = ha ∗ hb ∗ Yi Ti (x, t) − Yv,t . (33) The expression of the evaporation rate in terms of a convolution product of the interfacial temperature field is then obtained from the nonsolubility property of the inert gas in the liquid Eq. (27), J = K Yi Ti (x, t) − Yv,t , (34) where the mass transfer coefficient K is defined by ha ∗ hb ∗ Yi (Ti (x, t)) − Yv,t ρκ . (35) Le(H − h) [1 − Yi (Ti (x, t))][Yi (Ti (x, t)) − Yv,t ] In order to introduce a Biot number, we can now define the temperature Tv,t by Yv,t = Yi (Tv,t ). The evaporation rate takes the form J = Bi Ti (x, t) − Tv,t , (36) K=
where the Biot coefficient Bi is defined by K[Yi (Ti (x, t)) − Yv,t ] (37) Ti (x, t) − Tv,t and can therefore be seen as a nonlocal function of the interfacial temperature. The above definition of the temperature Tv,t is judicious because its value is inferior to the bottom temperature Te when the vapor mass fraction Yv,t at the top plate is small enough to have evaporation and Tv,t = Te in the nonevaporating case. In the stability analysis that is carried out in the next sections, one needs to introduce a reference solution, which is characterized by the fact that no motion exists in the system. In that situation, all unknown fields are function of the vertical coordinate z only and it is interesting to reformulate Eqs. (18)–(27) in this case. Since the temperature distribution is uniform along the interface, ha ∗ hb ∗ Yi (Ti (x, t)) = Yi (Ti ) and the Biot number is a function of the interfacial temperature only: Bi =
Biref (Ti ) =
Yi (Ti ) − Yv,t ρκ . Le(H − h) [1 − Yi (Ti )][Ti − Tv,t ]
(38)
(39)
where the subscript ref design the reference state. This relation is based on the Clausius–Clapeyron relation and the diffusion of the vapor into the inert gas. When the field of temperature is no longer uniform along the interface, the Biot number at a point of the interface is function not only of the interfacial temperature at this point but of the interfacial temperature distribution so that the model that we have derived in this section is a nonlocal one-sided model. Nevertheless, the shape of the function ha (cf. Fig. 2) shows that the temperature field close to this point is the more important. If we neglect the nonlocal contribution due to the function ha , then the evaporation rate becomes a local function of the interfacial temperature and of its spatial derivatives (which appear due to hb in Eq. (35)). If Eq. (32b) is truncated at the first term, the evaporation rate becomes a function of the interface temperature only, and is given explicitly by relations (38) and (39), which are now also valid in the nonreference state. We thus recover the local one-sided model with the assumption, used for example by Ha and Lai [10] and by Merkt and Bestehorn [15], that the evaporation rate is a function of the interface temperature only. However, relations (38) and (39) provide a better description since they give explicitly the evaporation rate and the Biot number dependencies with respect to the physical parameters of the problem. To conclude this section, it is worth noting that, for the linear and nonlinear thermoconvective instability studies that are carried out below, the differential form (18) of the equations is more easy to use than relations (36) and (37). The latter have, however, the advantage of leading by a truncation procedure to a natural hierarchy of simplified models, the first of which is the local one-sided model used in [10] and [15]. Finally, note that an additional reason to prefer the 1.5-sided model is that it does not require the quasistationarity of the mass fraction field that was assumed to deduce (36) and (37).
4. Linear stability analysis We will now compare the local one-sided and 1.5-sided models. As explained above, the two models have the same reference state, but we will see that the linear stability results differ because the 1.5-sided model leads to a Biot number that is a function of the wavenumber and the growth rate, whereas the local one-sided model does not. We work under the large latent heat of vaporization assumption [14,19], which leads to a time scale of the interface displacement that is much larger than the time scales of heat diffusion in the liquid layer and of vapor diffusion in the gas
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
layer and to a null vertical velocity of the liquid at the interface instead of the interfacial relation (24): (24 )
Wl = 0.
It was verified in [14,18] that the latent heats of vaporization of water and ethanol are large enough for this approximation to be valid. Within this assumption, the above system (18)– (27) admits a quasi-steady reference solution with liquid at rest where time derivatives can be neglected in Eqs. (18a) and (22). The velocity, temperature, and vapor mass fraction fields are given by Vref = 0,
(40)
Tref = (Tref,i − Te )z + Te , (41) Yv,ref,i − Yv,t Yv,ref = (42) (z − H ) + Yv,t , 1−H where Tref,i and Yv,ref,i are solutions of the three interfacial relations (18c), (26), and (27). To determine the stability of this quasi-steady reference state at time t0 , we introduce small perturbations (v = (u , v , w ), π , θ , yv , j ) for the velocity, pressure, temperature, vapor mass fraction, and evaporation rate, respectively. From Eqs. (18)–(27), it is easily deduced that for 0 z 1: ∇.v = 0,
(43)
2
∂t v + v .∇v + ∇π − Pr ∇ v − Pr Ra θ ez = 0,
2
(44)
∂t θ + v .∇θ + (Tref,i − Te )w − ∇ θ = 0,
(45)
for 1 z H : κ ∂t yv − ∇ 2 yv = 0, Le at z = 0:
(46)
v = θ = 0,
(47)
225
this reason, we make only a brief comment on the Biot number and we give results that make it possible to compare the local one-sided and 1.5-sided models. As shown in [11,12, 14] by combining the linearization of Eqs. (46), (48), and (51)–(53), it is possible to introduce an equivalent Biot condition, ∂z θ + Biev θ = 0,
at z = 1,
(54)
where the evaporative Biot number Biev depends on the wavenumber k and the growth rate σ , Biev =
i Jref dY dT |Tref,i 1 − Yv,ref,i
dg 1 − Yv,ref,i Le 2 × 1+ f k + σ , Yv,ref,i − Yv,t dl κ
(55)
with f (x) = x/ tanh(x). In contrast, the local one-sided model approximation introduced in Section 3.2 leads to the same relation (54) where the evaporative Biot number is independent of the wavenumber and given by dJref Biev = (56) , dTi Tref,i with Jref (Ti ) defined by relation (39). Relation (56) is the expression used by Merkt and Bestehorn [15] and is also the limit for k → 0 and σ → 0 of the Biot number given by relation (55). Fig. 3 represents the marginal stability curve obtained in the case where the liquid is water and the gas is air. The discrepancy between the local one-sided and 1.5-sided formulations is shown. The results corresponding to the full two-layer system are also presented to show the much better approximation given by the 1.5-sided model. If one considers a configuration where the depth of the gas layer is too small, the gas phase heat flux contribution should be
at z = H : yv = 0,
(48)
and at z = 1: w = 0,
(49)
∂z u + Ma ∂x θ = 0, ∂z v + Ma ∂y θ = 0,
−∂z θ = j , yv
(50) (51)
= Yi (Tref,i + θ ) − Yi (Tref,i ),
Jref yv − ρκ Le ∂z yv j = . 1 − Yv,ref,i − yv
(52) (53)
Except for the gravity effect, the linear model considered here for the linear stability analysis is the same as the model used in [14], with the additional assumption that the gas conduction can be neglected and that local thermodynamic equilibrium is achieved at the liquid–gas interface (because these assumptions are valid in the configuration that we study). For
Fig. 3. Marginal stability curves for water (g = 9.8 m s−2 , Yv,t = 0, Pt = 1 atm, Te = 320 K, dg = 0.01 m). (—) exact two-layer problem; (- - -) 1.5-sided model; (- · -) local one-sided model.
226
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
incorporated into the previous 1.5-sided model derivation (cf. Section 3.1). This can easily be done under the passive gas assumption by considering the thermal equation in the gas layer in addition to the diffusion equation for the vapor mass fraction in the gas layer. This contribution has already been taken into account in the linear regime by [12]. To check the role of these contribution, we have computed in the case of Fig. 3 the relative error, calculated at the nondimensional wavenumber k = 2, of the liquid depth marginal stability curves, compared to the exact two-layer one. It is 17% for the local one-sided model, 4.9% for the 1.5-sided model, and 3.4% for the 1.5-sided model with nonneglected thermal transport in the gas (not represented in Fig. 3), so that it was appropriate to neglect the gas phase heat flux contribution. Note finally that we checked that the exchange of stability (Re(σ ) = 0 ⇒ Im(σ ) = 0) is verified in the cases studied here.
of the set of nonlinear equations (57) in terms of the eigenfunctions fpk defined in Appendix A, namely f=
∞
Akp (t)fpk (x, y, z),
(58)
p=1 k
where Akp is the amplitude of the mode fpk , and after projection of the nonlinear set of Eqs. (57) onto the eigenfunctions of the adjoint problem defined in Appendix A and using the slaving principle explained in [20], one obtains the amplitude equations, whose general form is the same as that given in [14, Eq. (75)] i.e., ∂t Akp = σpk Akp k +ε Mq,p Akq − q Akq ∈Kc
k1 ,k2 ,k k1 k2 Nq,l,p Aq Al
q,l,k1 ,k2 k
k
Aq 1 ,Al 2 ∈Kc
−
k1 ,k3 ,k4 ,k k1 k3 k4 Zq,m,n,p Aq Am An .
(59)
q,m,n,k1 ,k3 ,k4
5. Nonlinear stability analysis for the improved one-sided model
k
k
k
Aq 1 ,Am3 ,An4 ∈Kc
We have performed a nonlinear analysis of the 1.5-sided model thanks to the same Galerkin–Eckhaus method that was used in [14,20]. Let us recall that the nonlinear approach in [14] was based on a one-sided model with Eqs. (46), (48), and (51)–(53) replaced by relations (54), (55) with the Biot number evaluated under the critical wavenumber and growth rate conditions, and at the temperature Te . In contrast, we use here the improved 1.5-sided model to deduce the amplitude equations governing the nonlinear convection. For clarity, the different models that we have considered and their equivalence are summarized in Table 3. To study the 1.5-sided model nonlinear regime, it is appropriate to reformulate the set of nonlinear equations (43)– (53) in term of linear and nonlinear operators. At order 3 of the small perturbations, it is written as Lc f = −L f + M∂t f + N2 (f, f ) + N3 (f, f, f ),
(57a)
Bf = 0,
(57b)
where the operators L, M, N2 , N3 , and B are defined in Appendix A, L = L − Lc , and subscript c designates the critical conditions. Let ε = (Te − Te,c )/Te,c be the relative distance to the linear stability threshold. L is a function of ε. We only retain the order ε dependence of L and we express the solution f
The expression of the coefficients is, however, different due to the use of the improved model. They are written as k k k k −1 ∂(Tref,i − Te ) Mq,p = −Te,c ζp wp , θ q ∂Te Te,c ∂ Pr 2 k k ∂ Pr Ra k k − ∇ vp , vq − θ ,w ∂Te Te,c ∂Te Te,c p q κ ∂ Le ∂Ma k k ∇ 2 y k , y k ∇ + θ , ϒ − h p q z=1 v,p v,q ∂Te Te,c ∂Te Te,c ∂Jref − y k , Θ k ∂Te Te,c v,p q z=1 ∂Yv,ref,i k k ∂z θp , Θq z=1 + ∂Te Te,c k ∂ ρκ Le ∂z yv,p , Θqk z=1 + ∂Te Te,c
i k k ∂ dY dT |Tref,i , θ ,Ψ + (60) ∂Te Te,c p v,q z=1 k k1 ,k2 ,k k2 k 1 Nq,l,p = vkq 1 .∇vkl 2 , vk p + vq .∇θl , θp k 1 ∂ θ k2 , Θ k + yv,q z l p z=1
−1 k k2 k 1 d2 Yi 1 + θ θ , Ψv,p z=1 ζpk , 2 dT 2 Tref,i,c q l (61)
Table 3 Different one-sided models Linear
Nonlinear
Bi(Ti )
Local approximation
Choice Biev = constant
Biev (k, σ, Te ) Bi nonlocal function of Ti 1.5-sided model
Equivalent when the vapor mass fraction field can be considered as quasi-stationary
Local choice Biev = (kc , 0, Te ) Equivalent when the vapor mass fraction field can be considered as quasi-stationary
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
227
k1 ,k5 ,k k5 ,k1 ,k k3 ,k4 ,k5 (Nq,r,p + Nr,q,p )Nm,n,r
k1 ,k3 ,k4 ,k Zq,m,n,p =−
k
r,k5
k
i Im(σm3 + σnk4 ) − σr 5
k
Ar 5 ∈Ks
+
k k θqk1 θm3 θnk4 , Ψv,p 1 d3 Yi z=1 , 3 k 6 dT Tref,i,c ζp
where k k k k k ζpk = vk p , vp + θp , θp + yv,p , yv,p ,
(62)
(63)
k ) are the eigenfunctions of the linear eigenvalue (vkp , θpk , yv,p problem, Im is the imaginary part, and Kc and Ks are respectively the sets of critical and slaved modes. It is interesting to stress that, due to the choice of the scale factors used to nondimensionalize the equations, which are different from those used, for instance, in [20], the Rayleigh and Marangoni effects appear through the first right-hand ∂(Tref,i −Te ) |Te,c wpk , θqk ; the side term of relation (60), i.e., ∂Te Pr Ra k k |Te,c θpk , wqk and ∂Ma terms ∂ ∂T ∂Te |Te,c ∇h θp , ϒ q z=1 ape pear only because of the variation of the fluid properties with Te . Contrary to the study of Bénard–Marangoni instabilities without evaporation [20], we have to take into account the variation of the physical properties of the fluid with Te . This is due to the fact that the temperature gradient is not obtained by external heating but only through evaporation; i.e., to obtain a temperature difference of, say, some kelvin between the bottom and top of the liquid layer, it is necessary to increase the temperature Te by about 10 K. Due to these large variations of Te , the physical properties can no longer be considered as constant. We will briefly discuss the meaning and the relative importance of the different terms of relation (60) for water. The ∂(T −Te ) first term ( ref,i |Te,c wpk , θqk ) is the “driving force” ∂Te term and the most important one. The following three terms appear because of the variation of the physical properties of the fluid with the temperature Te and are respectively about 10−1 , 10−6 , and 4 × 10−1 times less important than the first term. The following term appears because we take into account the time variations of the vapor mass fraction field (and the variation of the physical properties with the temperature) and is about 10−2 times less important than the first term. The last four terms mainly take into account the evaporation process (and also the variation of the physical properties with the temperature) and are respectively about 10−4 , 10−2 , 10−3 , and 10−1 times less important than the first term. In order to describe the interactions of roll, square, and hexagonal patterns, we choose 12 critical wave vectors with an angle of 30◦ between them, as in [14] (Kc = ±k6 ±k2 1 {A±k 1 , A1 , . . . , A1 }). A stability analysis of the solutions of the resulting amplitude equations gives rise to the stability diagram of Fig. 4 for a water–air system. In this figure, the characters C, R, H, and S mean respectively that the conductive state, roll, hexagon, and square patterns are stable. At a given liquid depth, when the ambient temperature
Fig. 4. Weakly nonlinear stability diagram for water. (—) one-sided model with Biot number; (- - -) 1.5-sided model (g = 9.8 m s−2 , Yv,t = 0, Pt = 1 atm, dg = 0.01 m).
Fig. 5. Weakly nonlinear stability diagram for ethanol. (—) without gravity; (- - -) with gravity g = 9.8 m s−2 (Yv,t = 0, Pt = 1 atm, dg = 0.01 m).
is progressively increased above the critical value, the first patterns to appear in this region are hexagons. If the liquid layer depth is great enough, then, by a further increase of the ambient temperature, the roll patterns become stable too, and over a particular temperature, only roll patterns are stable. When the liquid layer is small enough, square patterns can be stable instead of the roll patterns. To study the influence of gravity, the results corresponding to zero gravity can be compared with those described above. For the system described in Fig. 4, no difference was observed due to the low liquid layer thicknesses. In Fig. 5 we study the case of an ethanol–air system. We see that for liquid layer thicknesses great than ≈ 0.002 m,
228
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
the presence of gravity has a destabilizing influence on the system. In Fig. 4, we compare the results of the present 1.5-sided model with those of the one-sided model used in [14] with the acceleration due to gravity g = 9.8 m s−2 . We can see that we obtain the same results in both cases. To understand this very good agreement between both models, we compute the discrepancy terms that we derive in Appendix B. The new nonlinear terms that appear compared to [14] and the last term of ζpk (i.e., time variations of the vapor mass fraction field) are of about 10−6 times less important than the previous ones so that the difference between the two nonlinear weakly stability diagrams should be negligible. This demonstrates a posteriori that the previously nonjustified approximation of an evaluation of the Biot number at the critical wavenumber and at zero growth rate and of neglecting nonlinear terms due to evaporation as done in [14], is correct for the derivation of weakly nonlinear stability diagram.
Appendix A. Eigenvalue and adjoint eigenvalue problems We define the eigenvalue problem as Lc fpk = σpk Mfpk , Bfpk = 0,
(64a)
p = 1, 2, 3, . . . ,
f T = (v, π, θ, yv , u|z=1 , v|z=1 , θ |z=1 , ∂z yv |z=1 ),
Pr ∇ 2 ∇· (Te − Tref,i )(·) · ez ··· L= (·) · ex |z=1 ∂ z ∂z (·) · ey |z=1 ···
The authors thank Professor G. Lebon and Th. Desaive (Liège University) as well as P. Colinet (Brussels University) for fruitful discussions. This research has been supported by the ESA MAP under Contract AO-99-110 (CIMEX) and by the ESA PRODEX under Contract 90160. M.D. acknowledges financial support from the F.R.I.A. (Fonds pour la Formation à la Recherche dans l’Industrie et dans l’Agriculture, Belgium).
1 ··· ··· ··· M = ··· ··· ··· ··· and
··· ··· ··· ··· ··· ··· ··· ···
··· ··· 1 ··· ··· ··· ··· ···
(·)|z=0 ··· B = ··· (·) · ez |z=1
−∇ ··· ··· ··· ··· ···
··· ··· ··· 1 ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ···
(65)
Pr Ra ez ··· ∇2 ··· Ma ∂x (·)|z=1 Ma ∂y (·)|z=1
· · · (1 − Yv,ref,i )∂z (·)|z=1 ···
··· ··· ··· κ 2 Le ∇ ··· ···
Jref (·)|z=1 − ρκ Le ∂z (·)|z=1 (·)|z=1
6. Conclusions
Acknowledgments
(64b)
wherein the following notation has been used:
···
Our objective was to derive and study an improved model, the so-called 1.5-sided model, that describes the evaporation of a liquid layer surmounted by an inert gas–vapor mixture. Due to its linearity, we found that it is convenient to solve the diffusion equation for the vapor mass fraction in the gas layer in addition to the equations in the liquid layer. Moreover we have shown that this improved model is characterized by a Biot number that is a nonlocal function of the interfacial temperature, when the vapor mass fraction field can be considered as quasi-stationary. We compared the linear stability analysis of this 1.5-sided model with the local one-sided one used by [15] and found that the 1.5-sided model improves the results. We performed a weakly nonlinear stability analysis of this improved model and found that there is very good agreement with the simplified Biot approximation model results presented in [14]. This analysis provides thus a strong theoretical justification of the Biot approximation. Finally, we analyzed the influence of gravity on the results.
k ∈ R2 ,
··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ···
i − dY dT |Tref,i (·)|z=1 ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· , ··· ··· ··· ··· ··· ···
··· ··· ··· ··· , ··· ··· ··· ···
··· ··· ··· · · · (·)|z=0 · · · · · · · · · (·)|z=H ··· ··· ···
··· ··· ··· ···
(66)
(67)
··· ··· ··· ···
··· ··· ··· ···
··· ··· . ··· ··· (68)
An upper index T means transposition, subscript c designates the critical conditions, L is the operator of the linear problem associate to the set of Eqs. (43)–(53) σpk are the eigenvalues, and fpk the eigenfunctions. It should be observed that the linearizations of the boundary conditions (50)–(52) have been introduced in operator L. By doing so, the functions f do only have to verify a priori the boundary conditions (47)–(49) and the general form (65). The adjoint
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
operator of L, denoted L , is defined by Lf, f = f, L f ,
∀f, f ,
N2T (f, f ) =
0, ∂z θ
f T = v , π , θ , yv , ϒ |z=1 , Θ |z=1 , Ψv |z=1
ϒ = Pr(v · ex )ex + Pr(v · ey )ey , 1 Θ = − θ , 1 − Yv,ref,i Jref κ Ψv = θ − ∂z yv , 1 − Yv,ref,i Le
H +
L L dx dy
(71)
N3T (f, f , f
)
1 d3 Yi
= 0, 0, 0, 0, 0, 0, 0, θ |z=1 θ |z=1 θ |z=1 . 6 dT 3 Tref,i,c (78)
g4 h¯ 4 dz +
(73)
gi h¯ i dz
i=1 0
−L −L 8
gi h¯ i |z=1
(74)
(gi and h¯ i are the ith components of respectively g and the complex conjugate of h), L =
Pr ∇ 2 ∇· Pr Ra(·) · ez ··· − Pr ∂z (·) · ex |z=1 − Pr ∂z (·) · ey |z=1
Pr Ma ∂x (·) · ex |z=1 +∂y (·) · ey |z=1 ··· ··· ··· ··· κ 2 Le ∇ ··· ··· dYi κ dT |Tref,i Le ∂z (·)|z=1 κ (·)|z=1 − Le
−∇ ··· ··· ··· ··· ···
(Te − Tref,i )ez ··· ∇2 ··· ··· ··· dYi dT |Tref,i
··· ··· ··· ··· ··· ···
· · · − 1−Yv,ref,i Jref (·)|z=1 +∂z (·)|z=1 ρκ ··· Le(1−Yv,ref,i ) (·)|z=1
··· ··· ··· ··· ··· . ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ···
To understand the difference between the formulations of the 1.5-sided model and the one-sided model used in [14] in the weakly nonlinear regime, we can first introduce small perturbations of the quasi-steady reference state into nonlocal relations (36), (37). We are then led at order θ 3 to
(75)
··· ··· ··· ···
Lc fpk = σpk Mfpk , p = 1, 2, 3, . . . ,
at z = 1,
κ
k , y k = 0 for critical modes (k = k and ∂TLee |Te,c ∇ 2 yv,p c v,q and p = q = 1) so that the amplitude Eqs. (59)–(63) and [14, Eq. (75)] differ only by the noncritical modes function of the Biot number evaluated at (k, σ ) instead of (kc , σ = 0), by additional nonlinear terms (the last two terms of (61) and the last term of (62)), and by the last term of ζpk (relation (63)), which takes into account the time variations of the vapor mass fraction field.
(76a) k ∈ R2 .
(76b)
Finally, the nonlinear operators that appear in Eqs. (57) are defined by
(79)
where Biev is given by relation (55) with σ = 0 and N2 and N3 are respectively quadratic and cubic nonlinear terms, the expressions of which are not given here for reasons of brevity, so that one can say that, in the weakly nonlinear regime, the improved 1.5-sided model differs from the previous one-sided one by evaluating the Biot number at wavenumber k instead of kc , by additional nonlinear terms, and by the consideration of the time variations of the vapor mass fraction field. The difference between the two descriptions may be important a priori. By considering now the linear mode solutions of the eigenvalue and adjoint eigenvalue problems (64) and (76), k , and y k and Ψ k can be expressed as functhe modes yv,p v,p v,p tions of modes θpk and θpk respectively so that, by algebraic calculus, one can show that the last four terms of relation k k (60) can be rewritten αev,q,p θpk , θqk z=1 with αev,q,p a function of Biev and Te , the expression of which are not given here for reasons of brevity. The important point to note is k c ,σ =0,Te ) that algebraic calculus leads to αev,q,p = ∂Biev (k∂T |Te,c e ∂
We define the adjoint eigenvalue problem as
Bfpk = 0,
Appendix B. Difference between the formulations of the 1.5-sided model and the one-sided model
∂z θ + Biev θ = N2 + N3 ,
i=5
1
1 d2 Yi
θ | θ | z=1 z=1 2 dT 2 Tref,i,c (77)
and
(72)
3 1
|z=1 yv |z=1 ,
(70)
with
1 g, h = lim L→∞ 4L2
v · ∇v , 0, v · ∇θ , 0, 0,
(69)
where
229
References [1] H. Bénard, Ann. Chem. Phys. 23 (1901) 62. [2] D.A. Nield, J. Fluid Mech. 19 (1964) 341.
230
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
[3] J.R.A. Pearson, J. Fluid Mech. 4 (1958) 489. [4] E.L. Koschmieder, Bénard Cells and Taylor Vortices, Cambridge Univ. Press, Cambridge, UK, 1993. [5] C. Normand, Y. Pomeau, M. Velarde, Rev. Mod. Phys. 49 (1977) 581. [6] P. Colinet, J.C. Legros, M.G. Velarde, Nonlinear Dynamics of Surface Tension Driven Instabilities, Wiley–VCH, Berlin, 2001. [7] J.C. Berg, M. Boudart, A. Acrivos, J. Fluid Mech. 24 (1966) 721. [8] P. Colinet, L. Joannes, C.S. Iorio, B. Haut, M. Bestehorn, G. Lebon, J.C. Legros, Adv. Space Res. 32 (2003) 119. [9] H. Mancini, D. Maza, Europhys. Lett. 66 (6) (2004) 812–818. [10] V. Ha, C. Lai, J. Chin. Inst. Eng. 21 (1998) 547. [11] P. Colinet, B. Haut, J. Margerit, F. Peeters, J.C. Legros, G. Lebon, V. Halloin, in: Proc. of CHISA 2002, 15th International Congress of Chemical and Process Engineering, CD-ROM, Paper 858. [12] B. Haut, P. Colinet, J. Colloid Interface Sci. 285 (2005) 296–305.
[13] C. Moussy, G. Lebon, J. Margerit, Eur. Phys. J. B 40 (2004) 327–335. [14] M. Dondlinger, J. Margerit, P.C. Dauby, J. Colloid Interface Sci. 283 (2005) 522–532. [15] D. Merkt, M. Bestehorn, Physica D 185 (2003) 196. [16] R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, McGraw–Hill, New York, 1987. [17] V.P. Carey, Liquid–Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation in Heat Transfer Equipment, Hemisphere Publishing Corporation, Washington, DC, 1992. [18] J. Margerit, G. Lebon, P. Colinet, C.S. Iorio, J.C. Legros, Phys. Rev. E 68 (2003) 041601. [19] J.P. Burelbach, S.G.B.S.H. Davis, J. Fluid Mech. 195 (1988) 463. [20] P.M. Parmentier, V.C. Regnier, G. Lebon, J.C. Legros, Phys. Rev. E 54 (1996) 411.
Improved 1.5-sided model for the weakly nonlinear study of Bénard–Marangoni instabilities in an evaporating liquid layer J. Margerit ∗ , M. Dondlinger, P.C. Dauby Institut de Physique B5a, Université de Liège, Allée du 6 Août 17, B-4000 Liège, Belgium Received 18 November 2004; accepted 9 April 2005 Available online 1 June 2005
Abstract Bénard–Marangoni instability, with coupled gravity and surface tension effects, in an evaporating liquid layer surmounted by its vapor and an inert gas is investigated theoretically. We show that this system can be described by a model that consists in the liquid layer equations plus the diffusion equation for the vapor in the gas (the so-called 1.5-sided model) and that this model is equivalent to a one-sided model when the vapor mass fraction field can be considered as quasi-stationary, provided that the equivalent Biot number is a nonlocal function of the interface temperature. A comparison of weakly nonlinear results for the 1.5-sided model with a previous one-sided model [M. Dondlinger, J. Margerit, P.C. Dauby, J. Colloid Interface Sci. 283 (2) (2005) 522–532] that considered a Biot number depending on the wavenumber evaluated at the threshold is performed. Very good agreement is found between both models. For this reason, the present analysis can also be considered as a detailed theoretical justification for the use of a one-sided model in the study of evaporative thermoconvection. 2005 Elsevier Inc. All rights reserved. Keywords: Marangoni convection; Evaporation; Pattern formation; Amplitude equations
1. Introduction Thermoconvection is a standard example of self-organization phenomena. Bénard–Marangoni instability appears in horizontal fluid layers when a sufficiently large vertical temperature gradient is imposed. Arrangements of rolls, hexagons, or square convective patterns are then observed [1–6]. The temperature gradient can be produced by external heating of the bottom plate on which the liquid is lying, but if the liquid is volatile enough, the latent heat consumption caused by evaporation at the interface can lead to a vertical temperature gradient that is sufficient to induce instability. This is equivalent to considering a liquid layer cooled from above. The thermoconvective motion that appears gives rise to an increase of the evaporation rate. Evaporative convection plays an important role in heat exchangers, distillation columns, and drying technologies. * Corresponding author. Fax: +3243 66 2335.
E-mail address: [email protected] (J. Margerit). 0021-9797/$ – see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2005.04.031
Experimental results can be found in [7–9]. The qualitative results at Berg et al. [7] have revealed a variety of convection patterns for different fluid depths and Mancini and Maza [9] observed the successive transitions from square to hexagonal patterns and to conductive state for hexamethyldisiloxane (C6 H18 OSi2 ) when the liquid depth decreased progressively due to evaporation. Unfortunately properties of this fluid are not well enough tabulated to permit a quantitative comparison with our theoretical model. Colinet et al. describe in [8] the ITEL experiment, which made it possible to visualize in a microgravity environment evaporative chaotic patterns in ethanol. However, due to the short microgravity time of the experiment, only transient effects have been investigated. From a theoretical point of view, the instability of a liquid layer that evaporates in an inert gas has been investigated for linear stability by Ha and Lai [10] and by Colinet et al. [11], Colinet et al. [8], and Haut and Colinet [12]. Moussy et al. [13] have studied the role of interfacial deformation in the same problem. Under some assumptions, the system can be described by a one-sided model that introduces an equivalent Biot number to take into account the gas–vapor mixture
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
phase and the evaporation at the interface. In the different papers by Colinet and co-workers [8,11,12] a Biot number depending on the wavenumber is derived in the linear analysis. Its use has been extended by Dondlinger et al. [14] in the Galerkin–Eckhaus weakly nonlinear analysis of the system, with a Biot number depending on both the wavenumber evaluated at the threshold and the supercritical bottom plate temperature. The nonlinear study of such a model with a constant Biot number has also been carried out by Merkt and Bestehorn [15] thanks to a 3D numerical simulation. The purpose of the present work is to analyze the validity of the assumption of a Biot number, either considered as a constant or depending on the critical wavenumber. This will be done by deriving an improved model, the so-called 1.5-sided model, that consists of the liquid layer equations plus the diffusion equation for the vapor in the gas. A linear and a weakly nonlinear analysis of this 1.5-sided model will be presented and the results will be shown to be in very good agreement with the results of the simplified Biot approximation model presented in [14], for both the linear and nonlinear approaches. For this reason, the present analysis must be considered as a strong theoretical justification of the Biot approximation. The structure of the paper is the following. In Section 2, the physical system, the set of governing equations, and some dimensionless physical parameters are introduced. Section 3 is dedicated to the derivation of the 1.5-sided model in a configuration where thermal transport in the gas can be neglected, and to its equivalence to a one-sided model with an equivalent Biot number that is a nonlocal function of the interface temperature when the vapor mass fraction field can be considered as quasi-stationary. In the context of a linear analysis, this 1.5-sided model is equivalent to the model with a Biot number depending on the wavenumber derived in [11] and the corresponding results are compared in Section 4 with those of the complete two-sided problem and with those obtained with a constant Biot number [15]. Section 5 deals with the weakly nonlinear study of this 1.5-sided model. The influence of gravity is considered too. Finally, conclusions are drawn in Section 6.
2. Problem formulation We consider an evaporating liquid layer of infinite horizontal extent and depth dl , surmounted by a gas layer of thickness dg at a time t0 (Fig. 1). The gas layer is a mixture of the vapor of the liquid and an inert gas. Subscripts l, g, v, and i describe respectively the liquid, gas, vapor, and interface. The lower, rigid, and perfectly heat-conducting plate is maintained at the temperature Te . The gas layer is bounded by a perfectly heat-conducting upper plate maintained at the same temperature Te and at constant pressure and concentration. We assume that the Boussinesq approximation is valid and that the gas may be taken as perfect. The interface is assumed to be nondeformable so that the moving horizontal
221
Fig. 1. System under study.
liquid–gas interface is described by the equation z = h(t). Contrary to what was assumed in [14], we do not neglect gravity and the mass density ρ is assumed to be a linear function of temperature T , ρ = ρ0 1 − αT (T − Te ) , where αT is the constant coefficient of volumic expansion. The surface tension σ decreases linearly with the temperature, i.e., σ = σ0 − σT (T − Te ), where σT is a constant coefficient. All other physical properties of the fluids are considered constant and are evaluated at temperature Te . The properties of the fluids considered in the numerical calculations are taken from [16,17]. The physical properties of the gas mixture are also calculated at the reference temperature Te and at the reference mass fraction Yv,t of the vapor at the top plate by assuming an ideal gas mixture. As shown in [14], it is possible to assume that fluctuations in the system evolve on a much faster time scale than the displacement of the interface due to evaporation. Therefore the liquid layer thickness will be frozen during the stability analysis. The variables are expressed in dimensionless form: lengths (coordinates (x, y, z)) are scaled by the liquid layer thickness dl , time t by dl2 /κl , with κl the heat diffusivity of the liquid, velocity V = (U, V , W ) by κl /dl (vectors are denoted by bold letters), pressure P by ρ0,l κl2 /dl2 , the mass flux J at the interface by ρ0,l κl /dl , temperature T by L/cp,l where L is the latent heat of vaporization and cp,l the liquid heat capacity per unit mass, and the mass fraction of vapor in the gas Yv by 1 (Yv is already dimensionless). The dimensionless equations and boundary conditions governing the system are standard [14] and given as follows: in the liquid layer, ∇.Vl = 0, ∂t Vl + Vl .∇Vl + ∇Pl − Pr ∇ 2 Vl − Pr Ra Tl ez = 0, ∂t Tl + Vl .∇Tl − ∇ 2 Tl = 0;
(1) (2) (3)
in the gas layer, ∇.Vg = 0,
(4)
1 µ ∂t Vg + Vg .∇Vg + ∇Pg − Pr ∇ 2 Vg ρ ρ − ραT Pr Ra Tg ez = 0,
(5)
222
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
∂t Tg + Vg .∇Tg − κ∇ 2 Tg = 0, κ ∂t Yv + Vg .∇Yv − ∇ 2 Yv = 0; Le
(6)
3. Derivation of the improved 1.5-sided model
(7)
In this section we derive the 1.5-sided model, which is a simplification of the two-layer model, and we use this new model to derive an improved one-sided model, which is characterized by an equivalent Biot number that is a nonlocal function of the interfacial temperature. The one-sided model used in [15] (constant Biot number) will be recovered as a local approximation of the present improved one-sided model. Note also that, for the linear analysis, the 1.5-sided model is equivalent to the one used by Colinet et al. [11] and Dondlinger et al. [14] (Biot number depending on the wavenumber), but it improves the nonlinear analysis in the sense that it makes it possible to get rid of the assumption of a Biot number evaluated at the critical wavenumber.
at the bottom plate (z = 0), Ul = Vl = Wl = 0,
Tl = Te ;
(8)
at the top (z = H ), Vg = (0, 0, C),
Tg = T e ,
Yv = Yv,t ,
Pg = Pt ;
(9)
and at the interface (z = h(t)), J = Wl − Wi ,
(10)
Wl − Wi = ρ(Wg − Wi ),
(11)
Ul = Ug ,
(12)
3.1. The 1.5-sided model
Vl = Vg ,
−µ(∂x Wg + ∂z Ug ) + ∂x Wl + ∂z Ul + Ma ∂x Tl = 0, −µ(∂y Wg + ∂z Vg ) + ∂y Wl + ∂z Vl + Ma ∂y Tl = 0,
(13)
Tl = T g = T i ,
(14)
J − λ∂z Tg + ∂z Tl = 0,
(15)
psat (Tl ) −
Pt 1+
1−Yv rw Yv
= 0,
(16)
J κ ∂z Yv + (1 − Yv ) = 0, Le ρ
(17)
where Pg and Pl are respectively the gas and liquid dynamical pressures, C is a constant with respect to space, Wi = ∂t h is the interface velocity, and psat (T ) is the dimensionless liquid–vapor saturation curve. Yv,t and Pt are respectively the vapor mass fraction and the total thermodynamic pressure imposed at time t0 at the top of the system, and rw = Mg /Ml is the ratio of the inert gas and liquid molecular weights. The dimensionless parameters that appear in the above equations are defined in Table 1 with ν the kinematic viscosity, g the acceleration due to gravity, and Dg the mass diffusivity of the liquid vapor in the gas. These equations are the same as in [14], except that gravity is added. Moreover, we use Eq. (16) instead of the Hertz– Knudsen relation, because it was shown in [12,14,18] that the interface can be considered to be at thermodynamic equilibrium.
The set of equations in the liquid phase are coupled with the equations in the gas phase via the interfacial boundary conditions (10)–(17). The gas viscous stress in Eq. (13) can be neglected with respect to the liquid viscous stresses, as the viscosity ratio µ is small. As a consequence, Eq. (13) does not imply coupling between the two layers. We also assume that the gas layer is passive, i.e., that convection in the gas is negligible. This is true provided the convective terms of relations (6) and (7) are small as regards corresponding diffusive terms. With gas velocity Wg order of magnitude given by J /ρ, according to interfacial relations (10)–(11) where we have considered that ρ is small, the criteria for the gas layer to be passive, deduced from Eqs. (6) and (7), are J dg /dl κρ and J dg /dl κρ/ Le. Moreover, as the gas conductivity is small with respect to the liquid conductivity, the gas phase heat flux contribution in Eq. (15) is small compared to the other two terms, provided that the depth of the gas layer is not too small with respect to the liquid layer depth (dg /dl λ). We neglect this term in the following. The different assumptions are listed in Table 2 for water and ethanol and we have checked that for the situations analyzed in the following, these conditions are well satisfied. However, Eq. (15), expressing the interfacial energy conservation, still implies coupling between the liquid and gas phases through the term describing the heat consumed by evaporation. Indeed, this term is a function of the evaporation rate J that is given by the nonsolubility property of the
Table 1 Dimensionless parameters ν
Prandtl number
Pr = κl l
Rayleigh number
Ra = ν Tκ,lc l l l p,l σ dL Ma = κ Tµ cl
Marangoni number Lewis number Depth ratio
gα
d3L
l l p,l
κ Le = Dg g d +d H = ld g l
α
Volumetric expansion ratio
αT = αT ,g T ,l
Density ratio Diffusivity ratio
ρ = ρg l κ κ = κg
Viscosity ratio
µ = µg l
Conductivity ratio
λ = λg l
ρ
l
µ
λ
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
3.2. The one-sided model with a nonlocal Biot number
Table 2 Physical properties
µ κρ κρ/ Le λ
Water at 320 K
Ethanol at 250 K
Condition
0.033 0.18 0.21 0.042
0.0051 0.014 0.09 0.011
1 J dg /dl J dg /dl dg /dl
inert gas in the liquid (17), so that to find the evaporation rate, the mass fraction field Yv must be known. Relation (10) gives the liquid velocity Wl as a function of the evaporation rate J too. Fortunately, provided that the gas phase advection is negligible, as we have assumed above, the diffusion equation for the vapor in the gas (7) and the corresponding boundary conditions (9) and (16) are written as the linear equations Le ∂t Yv , in h(t) < z < H, κ Yv = Yv,t , at z = H,
(18b)
Yv = Yi (Tl ),
(18c)
∇ 2 Yv =
at z = h(t),
where
Yi (Tl ) = 1 − rw + rw
Pt psat (Tl )
(18a)
−1 .
(19)
After this Eq. (18) is solved, the evaporation rate is easily deduced from relation (17). Finally, our system is described by the above set of Eqs. (18), (19) in the gas layer and by the following equations in the liquid layer and the associated boundary conditions: in the liquid layer (0 < z < h(t)), ∇.Vl = 0,
(20)
∂t Vl + Vl .∇Vl + ∇Pl − Pr ∇ Vl − Pr Ra Tl ez = 0,
(21)
∂t Tl + Vl .∇Tl − ∇ 2 Tl = 0;
(22)
2
at the bottom plate (z = 0), Ul = Vl = Wl = 0,
223
Tl = Te ;
(23)
at the interface (z = h(t)), Wl = J + Wi ,
We now demonstrate that the model presented in the previous section is equivalent to a one-sided model with a Biot number that is a nonlocal function of the interfacial temperature. In order to avoid too complicated calculations that would introduce an additional memory effect in the Biot number, we assume that the vapor mass fraction field is quasi-stationary, which amounts to neglecting the time derivative in Eq. (18a). This hypothesis is valid provided the time scale dg2 /Dg of vapor diffusion in the gas is very small with respect to the time scale dl2 /κl of thermal diffusion in the liquid,√i.e., provided that the gas layer is not too thick (dg /dl κ/ Le). First, let us show that the evaporation rate can take the form of an expression depending on a convolution product of the interfacial temperature field. As the problem is unbounded in the horizontal direction, we take a Fourier transform of system (18) in the horizontal direction and we get 2 D − k 2 Yˆv = 0, for h(t) < z < H, (28a) at z = H, Yˆv = Yv,t δ(k), (28b) ˆ ˆ Yv = Yi Ti (x, t) , at z = h(t), (28c) where D stands for d/dz, ˆ· denotes the Fourier transform of function ·, k = kx ex + ky ey is the wavenumber, δ is the Dirac function, x = xex + yey , and t is the time. The solution of Eq. (28a) has the structure Yˆv (z) = A sinh (z − H )k + B cosh (z − H )k , (29) where k is the norm of vector k and the two constants A and B are deduced from the two boundary conditions (28b) and (28c): B = Yv,t δ(k), Yˆi (Ti (x, t)) − cosh[(H − h)k]Yv,t δ(k) . A=− sinh[(H − h)k]
(24)
∂y Wl + ∂z Vl + Ma ∂y Tl = 0,
(25)
∂z Tl = −J, ρκ ∂z Yv J =− . Le 1 − Yv
(26) (27)
Note that the mass conservation Eq. (20) and relation (24) lead to the compatibility relation Wi = −J , where · is a horizontal average, which makes it possible to determine Wi at any time. We call this simplified model the 1.5-sided model since only part of the fundamental conservation equations in the gas are taken into account.
(30b)
To determine the mass flux at the interface, we can calculate the derivative of the known function Yˆv (z) at z = h: −D Yˆv (h) =
∂x Wl + ∂z Ul + Ma ∂x Tl = 0,
(30a)
k{cosh[(H − h)k]Yˆi (Ti (x, t)) − Yv,t δ(k)} . sinh[(H − h)k]
(31)
Define now the two functions ha and hb as the inverse (H −h)k Fourier transforms of respectively sinh[(H −h)k] and cosh[(H − h)k]. The expressions of the functions ha and hb are π/2 1 , ha = (32a) (H − h) 1 + cosh[π x ] (H −h)
hb =
1 (H − h)
+∞ n=0
(−1)n 2n!
δ(x)(2n) ,
(32b)
where ·(2n) stands for the derivative of order 2n. The function ha is plotted in Fig. 2. It is worth noting that Yv,t is
224
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
Consequently the reference evaporation rate (36) is a function of the interfacial temperature only, given by Jref (Ti ) = Biref (Ti )[Ti − Tv,t ],
Fig. 2. Function ha .
+∞ independent of x and −∞ ha (s) ds = 1, so that ha ∗ Yv,t = +∞ Yv,t −∞ ha (s) ds = Yv,t , where ∗ is the convolution product. The inverse Fourier transform of Eq. (31) can be rewritten as
−(H − h)DYv (h) = ha ∗ hb ∗ Yi Ti (x, t) − Yv,t . (33) The expression of the evaporation rate in terms of a convolution product of the interfacial temperature field is then obtained from the nonsolubility property of the inert gas in the liquid Eq. (27), J = K Yi Ti (x, t) − Yv,t , (34) where the mass transfer coefficient K is defined by ha ∗ hb ∗ Yi (Ti (x, t)) − Yv,t ρκ . (35) Le(H − h) [1 − Yi (Ti (x, t))][Yi (Ti (x, t)) − Yv,t ] In order to introduce a Biot number, we can now define the temperature Tv,t by Yv,t = Yi (Tv,t ). The evaporation rate takes the form J = Bi Ti (x, t) − Tv,t , (36) K=
where the Biot coefficient Bi is defined by K[Yi (Ti (x, t)) − Yv,t ] (37) Ti (x, t) − Tv,t and can therefore be seen as a nonlocal function of the interfacial temperature. The above definition of the temperature Tv,t is judicious because its value is inferior to the bottom temperature Te when the vapor mass fraction Yv,t at the top plate is small enough to have evaporation and Tv,t = Te in the nonevaporating case. In the stability analysis that is carried out in the next sections, one needs to introduce a reference solution, which is characterized by the fact that no motion exists in the system. In that situation, all unknown fields are function of the vertical coordinate z only and it is interesting to reformulate Eqs. (18)–(27) in this case. Since the temperature distribution is uniform along the interface, ha ∗ hb ∗ Yi (Ti (x, t)) = Yi (Ti ) and the Biot number is a function of the interfacial temperature only: Bi =
Biref (Ti ) =
Yi (Ti ) − Yv,t ρκ . Le(H − h) [1 − Yi (Ti )][Ti − Tv,t ]
(38)
(39)
where the subscript ref design the reference state. This relation is based on the Clausius–Clapeyron relation and the diffusion of the vapor into the inert gas. When the field of temperature is no longer uniform along the interface, the Biot number at a point of the interface is function not only of the interfacial temperature at this point but of the interfacial temperature distribution so that the model that we have derived in this section is a nonlocal one-sided model. Nevertheless, the shape of the function ha (cf. Fig. 2) shows that the temperature field close to this point is the more important. If we neglect the nonlocal contribution due to the function ha , then the evaporation rate becomes a local function of the interfacial temperature and of its spatial derivatives (which appear due to hb in Eq. (35)). If Eq. (32b) is truncated at the first term, the evaporation rate becomes a function of the interface temperature only, and is given explicitly by relations (38) and (39), which are now also valid in the nonreference state. We thus recover the local one-sided model with the assumption, used for example by Ha and Lai [10] and by Merkt and Bestehorn [15], that the evaporation rate is a function of the interface temperature only. However, relations (38) and (39) provide a better description since they give explicitly the evaporation rate and the Biot number dependencies with respect to the physical parameters of the problem. To conclude this section, it is worth noting that, for the linear and nonlinear thermoconvective instability studies that are carried out below, the differential form (18) of the equations is more easy to use than relations (36) and (37). The latter have, however, the advantage of leading by a truncation procedure to a natural hierarchy of simplified models, the first of which is the local one-sided model used in [10] and [15]. Finally, note that an additional reason to prefer the 1.5-sided model is that it does not require the quasistationarity of the mass fraction field that was assumed to deduce (36) and (37).
4. Linear stability analysis We will now compare the local one-sided and 1.5-sided models. As explained above, the two models have the same reference state, but we will see that the linear stability results differ because the 1.5-sided model leads to a Biot number that is a function of the wavenumber and the growth rate, whereas the local one-sided model does not. We work under the large latent heat of vaporization assumption [14,19], which leads to a time scale of the interface displacement that is much larger than the time scales of heat diffusion in the liquid layer and of vapor diffusion in the gas
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
layer and to a null vertical velocity of the liquid at the interface instead of the interfacial relation (24): (24 )
Wl = 0.
It was verified in [14,18] that the latent heats of vaporization of water and ethanol are large enough for this approximation to be valid. Within this assumption, the above system (18)– (27) admits a quasi-steady reference solution with liquid at rest where time derivatives can be neglected in Eqs. (18a) and (22). The velocity, temperature, and vapor mass fraction fields are given by Vref = 0,
(40)
Tref = (Tref,i − Te )z + Te , (41) Yv,ref,i − Yv,t Yv,ref = (42) (z − H ) + Yv,t , 1−H where Tref,i and Yv,ref,i are solutions of the three interfacial relations (18c), (26), and (27). To determine the stability of this quasi-steady reference state at time t0 , we introduce small perturbations (v = (u , v , w ), π , θ , yv , j ) for the velocity, pressure, temperature, vapor mass fraction, and evaporation rate, respectively. From Eqs. (18)–(27), it is easily deduced that for 0 z 1: ∇.v = 0,
(43)
2
∂t v + v .∇v + ∇π − Pr ∇ v − Pr Ra θ ez = 0,
2
(44)
∂t θ + v .∇θ + (Tref,i − Te )w − ∇ θ = 0,
(45)
for 1 z H : κ ∂t yv − ∇ 2 yv = 0, Le at z = 0:
(46)
v = θ = 0,
(47)
225
this reason, we make only a brief comment on the Biot number and we give results that make it possible to compare the local one-sided and 1.5-sided models. As shown in [11,12, 14] by combining the linearization of Eqs. (46), (48), and (51)–(53), it is possible to introduce an equivalent Biot condition, ∂z θ + Biev θ = 0,
at z = 1,
(54)
where the evaporative Biot number Biev depends on the wavenumber k and the growth rate σ , Biev =
i Jref dY dT |Tref,i 1 − Yv,ref,i
dg 1 − Yv,ref,i Le 2 × 1+ f k + σ , Yv,ref,i − Yv,t dl κ
(55)
with f (x) = x/ tanh(x). In contrast, the local one-sided model approximation introduced in Section 3.2 leads to the same relation (54) where the evaporative Biot number is independent of the wavenumber and given by dJref Biev = (56) , dTi Tref,i with Jref (Ti ) defined by relation (39). Relation (56) is the expression used by Merkt and Bestehorn [15] and is also the limit for k → 0 and σ → 0 of the Biot number given by relation (55). Fig. 3 represents the marginal stability curve obtained in the case where the liquid is water and the gas is air. The discrepancy between the local one-sided and 1.5-sided formulations is shown. The results corresponding to the full two-layer system are also presented to show the much better approximation given by the 1.5-sided model. If one considers a configuration where the depth of the gas layer is too small, the gas phase heat flux contribution should be
at z = H : yv = 0,
(48)
and at z = 1: w = 0,
(49)
∂z u + Ma ∂x θ = 0, ∂z v + Ma ∂y θ = 0,
−∂z θ = j , yv
(50) (51)
= Yi (Tref,i + θ ) − Yi (Tref,i ),
Jref yv − ρκ Le ∂z yv j = . 1 − Yv,ref,i − yv
(52) (53)
Except for the gravity effect, the linear model considered here for the linear stability analysis is the same as the model used in [14], with the additional assumption that the gas conduction can be neglected and that local thermodynamic equilibrium is achieved at the liquid–gas interface (because these assumptions are valid in the configuration that we study). For
Fig. 3. Marginal stability curves for water (g = 9.8 m s−2 , Yv,t = 0, Pt = 1 atm, Te = 320 K, dg = 0.01 m). (—) exact two-layer problem; (- - -) 1.5-sided model; (- · -) local one-sided model.
226
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
incorporated into the previous 1.5-sided model derivation (cf. Section 3.1). This can easily be done under the passive gas assumption by considering the thermal equation in the gas layer in addition to the diffusion equation for the vapor mass fraction in the gas layer. This contribution has already been taken into account in the linear regime by [12]. To check the role of these contribution, we have computed in the case of Fig. 3 the relative error, calculated at the nondimensional wavenumber k = 2, of the liquid depth marginal stability curves, compared to the exact two-layer one. It is 17% for the local one-sided model, 4.9% for the 1.5-sided model, and 3.4% for the 1.5-sided model with nonneglected thermal transport in the gas (not represented in Fig. 3), so that it was appropriate to neglect the gas phase heat flux contribution. Note finally that we checked that the exchange of stability (Re(σ ) = 0 ⇒ Im(σ ) = 0) is verified in the cases studied here.
of the set of nonlinear equations (57) in terms of the eigenfunctions fpk defined in Appendix A, namely f=
∞
Akp (t)fpk (x, y, z),
(58)
p=1 k
where Akp is the amplitude of the mode fpk , and after projection of the nonlinear set of Eqs. (57) onto the eigenfunctions of the adjoint problem defined in Appendix A and using the slaving principle explained in [20], one obtains the amplitude equations, whose general form is the same as that given in [14, Eq. (75)] i.e., ∂t Akp = σpk Akp k +ε Mq,p Akq − q Akq ∈Kc
k1 ,k2 ,k k1 k2 Nq,l,p Aq Al
q,l,k1 ,k2 k
k
Aq 1 ,Al 2 ∈Kc
−
k1 ,k3 ,k4 ,k k1 k3 k4 Zq,m,n,p Aq Am An .
(59)
q,m,n,k1 ,k3 ,k4
5. Nonlinear stability analysis for the improved one-sided model
k
k
k
Aq 1 ,Am3 ,An4 ∈Kc
We have performed a nonlinear analysis of the 1.5-sided model thanks to the same Galerkin–Eckhaus method that was used in [14,20]. Let us recall that the nonlinear approach in [14] was based on a one-sided model with Eqs. (46), (48), and (51)–(53) replaced by relations (54), (55) with the Biot number evaluated under the critical wavenumber and growth rate conditions, and at the temperature Te . In contrast, we use here the improved 1.5-sided model to deduce the amplitude equations governing the nonlinear convection. For clarity, the different models that we have considered and their equivalence are summarized in Table 3. To study the 1.5-sided model nonlinear regime, it is appropriate to reformulate the set of nonlinear equations (43)– (53) in term of linear and nonlinear operators. At order 3 of the small perturbations, it is written as Lc f = −L f + M∂t f + N2 (f, f ) + N3 (f, f, f ),
(57a)
Bf = 0,
(57b)
where the operators L, M, N2 , N3 , and B are defined in Appendix A, L = L − Lc , and subscript c designates the critical conditions. Let ε = (Te − Te,c )/Te,c be the relative distance to the linear stability threshold. L is a function of ε. We only retain the order ε dependence of L and we express the solution f
The expression of the coefficients is, however, different due to the use of the improved model. They are written as k k k k −1 ∂(Tref,i − Te ) Mq,p = −Te,c ζp wp , θ q ∂Te Te,c ∂ Pr 2 k k ∂ Pr Ra k k − ∇ vp , vq − θ ,w ∂Te Te,c ∂Te Te,c p q κ ∂ Le ∂Ma k k ∇ 2 y k , y k ∇ + θ , ϒ − h p q z=1 v,p v,q ∂Te Te,c ∂Te Te,c ∂Jref − y k , Θ k ∂Te Te,c v,p q z=1 ∂Yv,ref,i k k ∂z θp , Θq z=1 + ∂Te Te,c k ∂ ρκ Le ∂z yv,p , Θqk z=1 + ∂Te Te,c
i k k ∂ dY dT |Tref,i , θ ,Ψ + (60) ∂Te Te,c p v,q z=1 k k1 ,k2 ,k k2 k 1 Nq,l,p = vkq 1 .∇vkl 2 , vk p + vq .∇θl , θp k 1 ∂ θ k2 , Θ k + yv,q z l p z=1
−1 k k2 k 1 d2 Yi 1 + θ θ , Ψv,p z=1 ζpk , 2 dT 2 Tref,i,c q l (61)
Table 3 Different one-sided models Linear
Nonlinear
Bi(Ti )
Local approximation
Choice Biev = constant
Biev (k, σ, Te ) Bi nonlocal function of Ti 1.5-sided model
Equivalent when the vapor mass fraction field can be considered as quasi-stationary
Local choice Biev = (kc , 0, Te ) Equivalent when the vapor mass fraction field can be considered as quasi-stationary
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
227
k1 ,k5 ,k k5 ,k1 ,k k3 ,k4 ,k5 (Nq,r,p + Nr,q,p )Nm,n,r
k1 ,k3 ,k4 ,k Zq,m,n,p =−
k
r,k5
k
i Im(σm3 + σnk4 ) − σr 5
k
Ar 5 ∈Ks
+
k k θqk1 θm3 θnk4 , Ψv,p 1 d3 Yi z=1 , 3 k 6 dT Tref,i,c ζp
where k k k k k ζpk = vk p , vp + θp , θp + yv,p , yv,p ,
(62)
(63)
k ) are the eigenfunctions of the linear eigenvalue (vkp , θpk , yv,p problem, Im is the imaginary part, and Kc and Ks are respectively the sets of critical and slaved modes. It is interesting to stress that, due to the choice of the scale factors used to nondimensionalize the equations, which are different from those used, for instance, in [20], the Rayleigh and Marangoni effects appear through the first right-hand ∂(Tref,i −Te ) |Te,c wpk , θqk ; the side term of relation (60), i.e., ∂Te Pr Ra k k |Te,c θpk , wqk and ∂Ma terms ∂ ∂T ∂Te |Te,c ∇h θp , ϒ q z=1 ape pear only because of the variation of the fluid properties with Te . Contrary to the study of Bénard–Marangoni instabilities without evaporation [20], we have to take into account the variation of the physical properties of the fluid with Te . This is due to the fact that the temperature gradient is not obtained by external heating but only through evaporation; i.e., to obtain a temperature difference of, say, some kelvin between the bottom and top of the liquid layer, it is necessary to increase the temperature Te by about 10 K. Due to these large variations of Te , the physical properties can no longer be considered as constant. We will briefly discuss the meaning and the relative importance of the different terms of relation (60) for water. The ∂(T −Te ) first term ( ref,i |Te,c wpk , θqk ) is the “driving force” ∂Te term and the most important one. The following three terms appear because of the variation of the physical properties of the fluid with the temperature Te and are respectively about 10−1 , 10−6 , and 4 × 10−1 times less important than the first term. The following term appears because we take into account the time variations of the vapor mass fraction field (and the variation of the physical properties with the temperature) and is about 10−2 times less important than the first term. The last four terms mainly take into account the evaporation process (and also the variation of the physical properties with the temperature) and are respectively about 10−4 , 10−2 , 10−3 , and 10−1 times less important than the first term. In order to describe the interactions of roll, square, and hexagonal patterns, we choose 12 critical wave vectors with an angle of 30◦ between them, as in [14] (Kc = ±k6 ±k2 1 {A±k 1 , A1 , . . . , A1 }). A stability analysis of the solutions of the resulting amplitude equations gives rise to the stability diagram of Fig. 4 for a water–air system. In this figure, the characters C, R, H, and S mean respectively that the conductive state, roll, hexagon, and square patterns are stable. At a given liquid depth, when the ambient temperature
Fig. 4. Weakly nonlinear stability diagram for water. (—) one-sided model with Biot number; (- - -) 1.5-sided model (g = 9.8 m s−2 , Yv,t = 0, Pt = 1 atm, dg = 0.01 m).
Fig. 5. Weakly nonlinear stability diagram for ethanol. (—) without gravity; (- - -) with gravity g = 9.8 m s−2 (Yv,t = 0, Pt = 1 atm, dg = 0.01 m).
is progressively increased above the critical value, the first patterns to appear in this region are hexagons. If the liquid layer depth is great enough, then, by a further increase of the ambient temperature, the roll patterns become stable too, and over a particular temperature, only roll patterns are stable. When the liquid layer is small enough, square patterns can be stable instead of the roll patterns. To study the influence of gravity, the results corresponding to zero gravity can be compared with those described above. For the system described in Fig. 4, no difference was observed due to the low liquid layer thicknesses. In Fig. 5 we study the case of an ethanol–air system. We see that for liquid layer thicknesses great than ≈ 0.002 m,
228
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
the presence of gravity has a destabilizing influence on the system. In Fig. 4, we compare the results of the present 1.5-sided model with those of the one-sided model used in [14] with the acceleration due to gravity g = 9.8 m s−2 . We can see that we obtain the same results in both cases. To understand this very good agreement between both models, we compute the discrepancy terms that we derive in Appendix B. The new nonlinear terms that appear compared to [14] and the last term of ζpk (i.e., time variations of the vapor mass fraction field) are of about 10−6 times less important than the previous ones so that the difference between the two nonlinear weakly stability diagrams should be negligible. This demonstrates a posteriori that the previously nonjustified approximation of an evaluation of the Biot number at the critical wavenumber and at zero growth rate and of neglecting nonlinear terms due to evaporation as done in [14], is correct for the derivation of weakly nonlinear stability diagram.
Appendix A. Eigenvalue and adjoint eigenvalue problems We define the eigenvalue problem as Lc fpk = σpk Mfpk , Bfpk = 0,
(64a)
p = 1, 2, 3, . . . ,
f T = (v, π, θ, yv , u|z=1 , v|z=1 , θ |z=1 , ∂z yv |z=1 ),
Pr ∇ 2 ∇· (Te − Tref,i )(·) · ez ··· L= (·) · ex |z=1 ∂ z ∂z (·) · ey |z=1 ···
The authors thank Professor G. Lebon and Th. Desaive (Liège University) as well as P. Colinet (Brussels University) for fruitful discussions. This research has been supported by the ESA MAP under Contract AO-99-110 (CIMEX) and by the ESA PRODEX under Contract 90160. M.D. acknowledges financial support from the F.R.I.A. (Fonds pour la Formation à la Recherche dans l’Industrie et dans l’Agriculture, Belgium).
1 ··· ··· ··· M = ··· ··· ··· ··· and
··· ··· ··· ··· ··· ··· ··· ···
··· ··· 1 ··· ··· ··· ··· ···
(·)|z=0 ··· B = ··· (·) · ez |z=1
−∇ ··· ··· ··· ··· ···
··· ··· ··· 1 ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ···
(65)
Pr Ra ez ··· ∇2 ··· Ma ∂x (·)|z=1 Ma ∂y (·)|z=1
· · · (1 − Yv,ref,i )∂z (·)|z=1 ···
··· ··· ··· κ 2 Le ∇ ··· ···
Jref (·)|z=1 − ρκ Le ∂z (·)|z=1 (·)|z=1
6. Conclusions
Acknowledgments
(64b)
wherein the following notation has been used:
···
Our objective was to derive and study an improved model, the so-called 1.5-sided model, that describes the evaporation of a liquid layer surmounted by an inert gas–vapor mixture. Due to its linearity, we found that it is convenient to solve the diffusion equation for the vapor mass fraction in the gas layer in addition to the equations in the liquid layer. Moreover we have shown that this improved model is characterized by a Biot number that is a nonlocal function of the interfacial temperature, when the vapor mass fraction field can be considered as quasi-stationary. We compared the linear stability analysis of this 1.5-sided model with the local one-sided one used by [15] and found that the 1.5-sided model improves the results. We performed a weakly nonlinear stability analysis of this improved model and found that there is very good agreement with the simplified Biot approximation model results presented in [14]. This analysis provides thus a strong theoretical justification of the Biot approximation. Finally, we analyzed the influence of gravity on the results.
k ∈ R2 ,
··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ··· ··· ···
i − dY dT |Tref,i (·)|z=1 ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· , ··· ··· ··· ··· ··· ···
··· ··· ··· ··· , ··· ··· ··· ···
··· ··· ··· · · · (·)|z=0 · · · · · · · · · (·)|z=H ··· ··· ···
··· ··· ··· ···
(66)
(67)
··· ··· ··· ···
··· ··· ··· ···
··· ··· . ··· ··· (68)
An upper index T means transposition, subscript c designates the critical conditions, L is the operator of the linear problem associate to the set of Eqs. (43)–(53) σpk are the eigenvalues, and fpk the eigenfunctions. It should be observed that the linearizations of the boundary conditions (50)–(52) have been introduced in operator L. By doing so, the functions f do only have to verify a priori the boundary conditions (47)–(49) and the general form (65). The adjoint
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
operator of L, denoted L , is defined by Lf, f = f, L f ,
∀f, f ,
N2T (f, f ) =
0, ∂z θ
f T = v , π , θ , yv , ϒ |z=1 , Θ |z=1 , Ψv |z=1
ϒ = Pr(v · ex )ex + Pr(v · ey )ey , 1 Θ = − θ , 1 − Yv,ref,i Jref κ Ψv = θ − ∂z yv , 1 − Yv,ref,i Le
H +
L L dx dy
(71)
N3T (f, f , f
)
1 d3 Yi
= 0, 0, 0, 0, 0, 0, 0, θ |z=1 θ |z=1 θ |z=1 . 6 dT 3 Tref,i,c (78)
g4 h¯ 4 dz +
(73)
gi h¯ i dz
i=1 0
−L −L 8
gi h¯ i |z=1
(74)
(gi and h¯ i are the ith components of respectively g and the complex conjugate of h), L =
Pr ∇ 2 ∇· Pr Ra(·) · ez ··· − Pr ∂z (·) · ex |z=1 − Pr ∂z (·) · ey |z=1
Pr Ma ∂x (·) · ex |z=1 +∂y (·) · ey |z=1 ··· ··· ··· ··· κ 2 Le ∇ ··· ··· dYi κ dT |Tref,i Le ∂z (·)|z=1 κ (·)|z=1 − Le
−∇ ··· ··· ··· ··· ···
(Te − Tref,i )ez ··· ∇2 ··· ··· ··· dYi dT |Tref,i
··· ··· ··· ··· ··· ···
· · · − 1−Yv,ref,i Jref (·)|z=1 +∂z (·)|z=1 ρκ ··· Le(1−Yv,ref,i ) (·)|z=1
··· ··· ··· ··· ··· . ··· ··· ··· ··· ··· ··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ···
To understand the difference between the formulations of the 1.5-sided model and the one-sided model used in [14] in the weakly nonlinear regime, we can first introduce small perturbations of the quasi-steady reference state into nonlocal relations (36), (37). We are then led at order θ 3 to
(75)
··· ··· ··· ···
Lc fpk = σpk Mfpk , p = 1, 2, 3, . . . ,
at z = 1,
κ
k , y k = 0 for critical modes (k = k and ∂TLee |Te,c ∇ 2 yv,p c v,q and p = q = 1) so that the amplitude Eqs. (59)–(63) and [14, Eq. (75)] differ only by the noncritical modes function of the Biot number evaluated at (k, σ ) instead of (kc , σ = 0), by additional nonlinear terms (the last two terms of (61) and the last term of (62)), and by the last term of ζpk (relation (63)), which takes into account the time variations of the vapor mass fraction field.
(76a) k ∈ R2 .
(76b)
Finally, the nonlinear operators that appear in Eqs. (57) are defined by
(79)
where Biev is given by relation (55) with σ = 0 and N2 and N3 are respectively quadratic and cubic nonlinear terms, the expressions of which are not given here for reasons of brevity, so that one can say that, in the weakly nonlinear regime, the improved 1.5-sided model differs from the previous one-sided one by evaluating the Biot number at wavenumber k instead of kc , by additional nonlinear terms, and by the consideration of the time variations of the vapor mass fraction field. The difference between the two descriptions may be important a priori. By considering now the linear mode solutions of the eigenvalue and adjoint eigenvalue problems (64) and (76), k , and y k and Ψ k can be expressed as functhe modes yv,p v,p v,p tions of modes θpk and θpk respectively so that, by algebraic calculus, one can show that the last four terms of relation k k (60) can be rewritten αev,q,p θpk , θqk z=1 with αev,q,p a function of Biev and Te , the expression of which are not given here for reasons of brevity. The important point to note is k c ,σ =0,Te ) that algebraic calculus leads to αev,q,p = ∂Biev (k∂T |Te,c e ∂
We define the adjoint eigenvalue problem as
Bfpk = 0,
Appendix B. Difference between the formulations of the 1.5-sided model and the one-sided model
∂z θ + Biev θ = N2 + N3 ,
i=5
1
1 d2 Yi
θ | θ | z=1 z=1 2 dT 2 Tref,i,c (77)
and
(72)
3 1
|z=1 yv |z=1 ,
(70)
with
1 g, h = lim L→∞ 4L2
v · ∇v , 0, v · ∇θ , 0, 0,
(69)
where
229
References [1] H. Bénard, Ann. Chem. Phys. 23 (1901) 62. [2] D.A. Nield, J. Fluid Mech. 19 (1964) 341.
230
J. Margerit et al. / Journal of Colloid and Interface Science 290 (2005) 220–230
[3] J.R.A. Pearson, J. Fluid Mech. 4 (1958) 489. [4] E.L. Koschmieder, Bénard Cells and Taylor Vortices, Cambridge Univ. Press, Cambridge, UK, 1993. [5] C. Normand, Y. Pomeau, M. Velarde, Rev. Mod. Phys. 49 (1977) 581. [6] P. Colinet, J.C. Legros, M.G. Velarde, Nonlinear Dynamics of Surface Tension Driven Instabilities, Wiley–VCH, Berlin, 2001. [7] J.C. Berg, M. Boudart, A. Acrivos, J. Fluid Mech. 24 (1966) 721. [8] P. Colinet, L. Joannes, C.S. Iorio, B. Haut, M. Bestehorn, G. Lebon, J.C. Legros, Adv. Space Res. 32 (2003) 119. [9] H. Mancini, D. Maza, Europhys. Lett. 66 (6) (2004) 812–818. [10] V. Ha, C. Lai, J. Chin. Inst. Eng. 21 (1998) 547. [11] P. Colinet, B. Haut, J. Margerit, F. Peeters, J.C. Legros, G. Lebon, V. Halloin, in: Proc. of CHISA 2002, 15th International Congress of Chemical and Process Engineering, CD-ROM, Paper 858. [12] B. Haut, P. Colinet, J. Colloid Interface Sci. 285 (2005) 296–305.
[13] C. Moussy, G. Lebon, J. Margerit, Eur. Phys. J. B 40 (2004) 327–335. [14] M. Dondlinger, J. Margerit, P.C. Dauby, J. Colloid Interface Sci. 283 (2005) 522–532. [15] D. Merkt, M. Bestehorn, Physica D 185 (2003) 196. [16] R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, McGraw–Hill, New York, 1987. [17] V.P. Carey, Liquid–Vapor Phase Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation in Heat Transfer Equipment, Hemisphere Publishing Corporation, Washington, DC, 1992. [18] J. Margerit, G. Lebon, P. Colinet, C.S. Iorio, J.C. Legros, Phys. Rev. E 68 (2003) 041601. [19] J.P. Burelbach, S.G.B.S.H. Davis, J. Fluid Mech. 195 (1988) 463. [20] P.M. Parmentier, V.C. Regnier, G. Lebon, J.C. Legros, Phys. Rev. E 54 (1996) 411.