Onset of Marangoni Convection in a Horizontal Fluid Layer

Ind. Eng. Chem. Res. 2007, 46, 5775-5784

5775

Onset of Marangoni Convection in a Horizontal Fluid Layer Experiencing Evaporative Cooling Min C. Kim,*,† Chang K. Choi,‡ Do Y. Yoon,§ and Tae J. Chung‡ Department of Chemical Engineering, Cheju National UniVersity, Cheju 690-756, Korea, School of Chemical and Biological Engineering, Seoul National UniVersity, Seoul 151-744, Korea, and Department of Chemical Engineering, Kwangwoon UniVersity, Seoul 139-701, Korea

The onset of Marangoni convection in an initially isothermal, quiescent fluid layer experiencing evaporative cooling from above is analyzed theoretically. For this system, the stability analysis on regular cell modes is conducted, based on the linear stability theory, the energy method, and their modifications. The stability equations are derived and the critical time to mark the onset of convection (τc) is obtained, as a function of the Marangoni number (Ma) and the Prandtl number (Pr). The linear stability theory predicts more-stable conditions than the energy method. The present results show that τc decreases as Pr increases for a given Ma value. The present predictions are compared with the available experimental data of propanol. 1. Introduction The convective stability of fluids bounded by a free surface has been studied extensively, in association with Marangoni convection. Stability problems in this deformable surface are encountered in various chemical engineering situations, such as gas-liquid absorption/desorption,1,2 liquid-liquid extraction,3 and chemical reaction.4 Marangoni convection can enhance the mass-transfer rate, but nonuniformities in product quality are incurred. Therefore, it becomes an important question to predict the critical conditions to mark the onset of Marangoni convection. Recently, surface-tension-driven convection was analyzed in relation with the CO2 absorption process.5,6 When an initially quiescent liquid layer bounded by an upper evaporating free surface is cooled from above, surface tension gradients can often generate Marangoni convection. In very thin liquid films, surface-tension effects will dominate buoyant forces and surface-tension-driven cells will occur when the Marangoni number exceeds a certain critical value. The stability problem in a liquid layer with temperature-dependent surface tension was analyzed first by Pearson,7 using the linear stability theory. Later, Davis8 reconsidered this problem, using the energy method. Although their analyses are indeed appropriate for describing convection within thin layers with linear temperature profiles, in many practical systems, the Marangoni number (Ma) is quite large, and the most important parameter becomes the critical time (tc), which marks the onset of natural convection. This time-dependent instability problem, under the deep-pool condition of t f 0, is quite complex. The related instability analysis has been conducted using the quasi-static approximation,9 the propagation theory,10 the energy method,11 and the maximumtransient-Marangoni-number criterion.12 The first two models are based on linear theory and yield the critical time tc as the parameter. Except for the second model, the critical conditions are independent of the Prandtl number (Pr). The last model gives the critical condition, which is dependent on the boundary conditions on the lower boundary. The energy method yields the lower bounds on the experimental onset times, and its * To whom correspondence should be addressed. Tel.: +82-64-7543685. Fax: +82-64-755-3670. E-mail address: [email protected]. † Cheju National University. ‡ Seoul National University. § Kwangwoon University.

Figure 1. Sketch of the basic conduction state considered here.

modification has been tested by considering the growth of disturbance energy.11 We developed the propagation theory to predict the critical condition to mark buoyancy-driven convection in developing temperature fields. This theory uses the thermal penetration depth as a length scaling factor and transforms the linearized perturbation equations similarly, under the principle of exchange of stabilities. Its core concept is that a most dangerous mode of disturbances will set in, experiencing instantaneous variations in their values upon their onset. The resulting stability criteria have compared favorably with experimental data in Beńard convection,13,14 Soret-driven convection,15 and even Taylor-like vortices.16 Hence, the stability analysis based on this theory is extended to the onset of Marangoni convection in evaporating liquid layers. Also, the energy method is modified and the strong stability criterion is relaxed by introducing the relative stability criteria. For the present problem, the new stability equations will be derived based on the propagation theory, the dominant mode analysis, and the modified energy method. Here, we will concentrate on the instability problem in an initially isothermal, quiescent fluid layer. Starting from time t ) 0, the upper free boundary is cooled uniformly by evaporation. For this specific system, the surface-tension-driven instability criteria will be obtained, using the linear stability theory and the energy method and their modifications. These predictions are compared and their validity is discussed in comparison with extant experimental and theoretical results. This study may be considered to be an extension of the work by Pearson7 and Gumerman and Homsy.11 2. Theoretical Analysis 2.1. Governing Equations. The system considered here is a Newtonian fluid layer with an initial temperature Ti. For time t

10.1021/ie070006y CCC: $37.00 © 2007 American Chemical Society Published on Web 07/20/2007

5776

Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007

g 0, the horizontal layer of fluid depth d experiences evaporative cooling with heat flux qw through the upper free boundary and its lower boundary is kept at the initial temperature Ti. The evaporation rate and corresponding heat flux qw are assumed to be constant. The schematic diagram of the basic system of pure conduction is shown in Figure 1. For a high qw value, surface-tension-driven convection will set in at a certain time and the governing equations of flow and temperature fields are expressed as

∇‚U)0

(1)

∂ 1 + U ‚ ∇ U ) - ∇P + ν∇2U ∂t F

(2)

{∂t∂ + U ‚ ∇}T ) R∇ T

(3)

{

}

2

where U is the velocity vector, P the dynamic pressure, F the density, ν the kinematic viscosity, and R the thermal diffusivity. The surface temperature (Ts) at the vertical distance Z ) 0 decreases with time during the conduction period. The important parameters to describe the present system are the Prandtl number (Pr) and the Marangoni number (Ma), which are defined by

ν R

(4a)

sqwd2 kRµ

(4b)

Pr ) Ma )

where µ is the viscosity, s the negative rate of change of surface tension (with respect to temperature), and k the thermal conductivity. In case of very slow cooling or heating, the basic temperature profile finally becomes linear and time-independent, wherein the critical conditions are well-known:7,8

MaLc ) 79.61

(5a)

MaEc ) 56.77

(5b)

∞

θ0 ) z - 1 - 2 θ0 ) - x4τ

∑ n)1

1 2

cos(λnz) exp(- λn2τ)

{ ( )

∞

(-1)n ∑ n)0

ierfc

n

+

xτ

ζ

-

2

ierfc

(

n+1

-

xτ

∂θ0 ∂2θ0 ) 2 ∂τ ∂z

(6)

θ/0 ) -

θ0

x4τ

θ0 ) 0

(at τ ) 0 and z ) 1)

∂θ0 )1 ∂z

(at z ) 0)

(7a) (7b)

In the aforementioned equations, τ ) Rt/d2, z ) Z/d, and θ0 ) k(T - Ti)/(qwd). The subscript “0” denotes the basic state. Using the separation of variables technique and the Laplace transform method, the exact solutions of eqs 6 and 7 can be obtained as

2

(8b)

) ierfc

(2ζ)

(with δT ) 3.21xτ)

(9)

where δT denotes the dimensionless thermal boundary layer thickness with θ/0(ζ)/θ/0(0) ) 0.01. The aforementioned Leveque-type soution of eq 9 is in good agreement with the exact solutions of eq 8 for τ e 0.01. For τ e 0.01, eq 8b with n ) 0 yields almost the same temperature profile as eq 9. It is noted that the upper wall temperatures (θs) are -x4τ/π for small τ. The experimental results show that this profile is more adequate to model the temperature profile during the evaporation,3 although a minor discrepancy is observed at the early stage of cooling. Because we are concerned with the deep-pool case of a large Marangoni number and small τ, the aforementioned Leveque-type solution of eq 9 is used primarily for the stability analysis. 2.2. Linear Stability Theory. 2.2.1. Perturbed Equations. Under linear stability theory, infinitesimal perturbation quantities that are caused by incipient convective motion at the dimensionless critical time τc can be formulated, in dimensionless form, in terms of the temperature component θ1 and the vertical velocity component w1 by linearizing eqs 1-3:

{Pr1 ∂τ∂ - ∇ }∇ w ) 0

(10)

∂θ1 ∂θ0 + Ma(a2w1) ) ∇ 2θ 1 ∂τ ∂z

(11)

2

2

where the Laplacian is given as ∇2 ) ∂2/∂x2 + ∂2/∂y2 + ∂2/∂z2. Here, the velocity component has the scale of R/d and the temperature component has that of Rµ/(a2sd). The proper boundary conditions are given as follows:5

w1 )

∂2w1 ∂θ1 ) 0, ) θ1 ∂z ∂z2

w1 )

with the following initial and boundary conditions:

)}

ζ

where λn ) (n - 1/2)π, ζ ) z/xτ, and ierfc(η) ) exp(-η2)/ xπ - ηerfc(η). For deep-pool systems of ∆T ∝ xRt, the Leveque-type solutions are given by

1

which are independent of Pr. Here, the superscripts L and E denote the linear stability theory and the energy method, respectively. However, for rapid cooling or for a heating system with a large Ma value, the resulting transient stability problem is complicated and the critical time tc to mark the onset of convective motion becomes an important question. For the basic state of heat conduction, the dimensionless temperature profile is represented by

(8a)

µn

∂w1 ) θ1 ) 0 ∂z

(at z ) 0) (at z ) 1)

(12a) (12b)

wherein the impermeable condition with qw ) constant ((∂θ1/ ∂z) ) 0) and the force balance including surface tension effects are applied to the upper boundary and the no-slip conditions with T ) constant (θ1 ) 0) are applied to the lower boundary. For a given Pr, the critical time τc should be determined using eqs 10-12. Now, convective motion is assumed to exhibit horizontal periodicity and the normal-mode analysis is applied. The perturbed quantities then can be expressed as

Ind. Eng. Chem. Res., Vol. 46, No. 17, 2007 5777

[w1(τ,x,y,z),θ1(τ,x,y,z)] ) [w1(τ,z),θ1(τ,z)] exp[i(axx + ayy)] (13) where “i” is the imaginary number and the horizontal wavenumber (a) has the relation of a ) [ax2 + ay2]1/2. 2.2.2. Dominant Mode Method. For the limiting case of τ f 0, it may be natural to observe that θ1(τ,z) ) θ1(τ,ζ) and w1(τ,z) ) w1(τ,ζ), considering eq 9. Now, eqs 10-12 are transformed to the similarity variable of the base state. Then, for the limiting case of Pr f ∞, the perturbation equations can be expressed as

(

)

2 1 ∂2 - a2 w1 ) 0 2 τ ∂ζ

(

(14)

)

∂θ/0 ∂θ1 1 ∂2 1 ∂ 2 2 + τ θ ) a Maw ζ a (15) 1 1 ∂τ τ ∂ζ2 2 ∂ζ ∂ζ The corresponding boundary conditions are

∂θ1 ∂2w1 w1 ) ) 0, 2 ) τθ1 ∂ζ ∂ζ w1 )

∂w1 ) θ1 ) 0 ∂ζ

(at ζ ) 0) (for ζ f ∞)

(16a)

〈 〉 dθ/0

w*

dζ

σ)

∫0

∞

)

dθ/0 w* dζ dζ

∫0∞ exp(-ζ2/4) dζ

(

w* 1 )ζ exp(-axτζ) A0(τ) 2axτ

(18)

(k ) 0, 1, 2, ...) (19)

The eigenfunctions φk of L are the Hermite polynomials Hk(ζ/2) in the semi-infinite domain with weight function exp(ζ2/4). The eigenvalues are λk ) -(k + 1)/2 for k ) 0, 1, 2, ... The most unstable eigenfunction (i.e., the dominant mode solution satisfying the aforementioned boundary conditions of eq 16) is given by

φ0 ) exp

( )

-ζ 1 , Lφ0 ) - φ0(ζ) 4 2

)

(23)

r0 )

1 dE0 E0 dτ

(24a)

r1 )

1 dE1 E1 dτ

(24b)

In the present system, the dimensionless natural energy is defined as6,9

(20)

〈 〉

dθ/0 dA0 1 2 2 2 ) στA0 ) - + a τ A0 - a Maτ w* τ dτ 2 dζ

(21a)

1 1 1 2 〈u 〉 + β〈θ2〉 2 Pr 2

(25)

where 〈(‚)〉 ) ∫Ω (‚) dΩ, where Ω is the system volume, and the coupling constant β is a positive one. For the limiting case of Pr f ∞, the natural energy of eq 25 can be reduced to

1 E(τ) ) 〈θ2〉 2

(26)

From the base temperature distribution of eq 9, r0 can be obtained as

r0 )

2

The long-wave mode solution shows that the zeroth mode of the self-similar diffusion operator decays with time, while the remainder of the spectrum decays more rapidly. Therefore, the zeroth mode is the dominant one. By substituting the dominant mode solution θ1 ) A0(τ) exp(-ζ2/4) into eq 15 and integrating across the domain, we obtain

(

(22)

where the temporal growth rate of the base energy (r0) and the perturbation energy (r1) in the (τ,z) domain are defined as

E(τ) ) -ζ2 ζ Hk 4 2

(at τ ) τr)

(17)

Ak(τ)φk(ζ)

( ) ()

(21c)

Therefore, the driving force for the instability, the integral 〈w* dθ* 0/dζ〉, is also a function of time and the wavenumber. Now, let us consider the stability criteria. Under the relative stability concept, the relative instability time τr is determined as14,19

with

Lφk ) λkφk(ζ) ) λk exp

(21b)

Here, w* ) w1/τ and σ is the growth rate in the (τ,ζ) domain. With θ1 ) A0(τ) exp(-ζ2/4), w* can be obtained analytically from eqs 14 and 16 as

r 0 ) r1

∞

∑

〈 〉

1 w* dθ0 1 dA0 )+ a2 - a2Maτ A0 dτ 2τ A0 dζ

and the temperature disturbance is expanded as

k)0

)

xπ

and

ζ ∂ ∂2 + ∂ζ2 2 ∂ζ

θ1(τ,ζ) )

)

∫0∞ w* erfc(η/2)dη

(16b)

The amplitudes of the perturbations have a time-dependent, spatial structure. By following Riaz et al.’s procedure,18 the streamwise operator of the concentration disturbance in the transformed coordinates is expressed as

L)

where

3 for τ f 0 2τ

(27)

Using the dominant mode solution of θ1 ) A0(τ) exp(-ζ2/4) and eqs 24 and 27, the following relation for the relative instability time τr, where r1 ) 2σ + 1/2τr, is obtained:

Maτr )

xπ(1 + a2τr) a2τr

∫0∞ η exp(-axτrη) erfc(η2) dη

(28)

At τ ) τr, σ ) 1/(2τr), which is quite different from Riaz et al.’s instability criteria,18 i.e., σ ) 0. For a given Ma value, the critical conditions can be determined analytically from eq 28

5778


w1 ≈ δT2 ≈ τ θ1

(32a)

∂θ0 ≈1 ∂z

Ma*a*2

(32b)

where δT ) ∆T/d ≈ xτ, Ma* ) Maτ, and a* ) axτ. The aforementioned scaling analysis is described in the work of Kang and Choi.17 From the previous discussion, there are many possible forms of dimensionless amplitude functions of disturbances, such as

[w1(τ,z),θ1(τ,z)] ) [τn+1w*(τ,z),τnθ*(τ,z)]

Figure 2. Marginal stability curve for deep-pool systems of Pr f ∞, under the propagation theory and the dominant mode method.

when

∂τr )0 ∂a

(29)

under the critical wavenumber ac. The resulting τr-value is the minimum, i.e., τc. With the present method, which is called the dominant mode method, the marginal stability curve for the present absorption system of small time is given in Figure 2 and the critical conditions are obtained as

Maτc ) 17.31

(for τ f 0)

(30a)

and

acxτc ) 0.5

(for τ f 0)

(30b)

Although the aforementioned dominant mode method gives some mathematical and physical insights into the growth of disturbances, its validity is limited by the condition of axτ f 0. Because the critical condition of eq 30 deviates from the condition of axτ f 0, the eigenfunction of the temperature disturbance under the critical conditions may be different from θ1 ) A0(τ) exp(-ζ2/4). This will be discussed later. 2.2.3. Propagation Theory. Now the present stability problem is examined using the propagation theory, which has been used to determine the onset condition of the BeńardRayleigh convection and the Taylor-Go¨rtler vortices under time-dependent basic fields.13-16 This model is based on the assumption that, in deep-pool systems, the infinitesimal temperature disturbances are propagated mainly within the thermal penetration depth at the onset of convective motion ∆T(∝ xRt) and the following scale relations are valid for perturbed quantities from eqs 2 and 3: 2

W1 a γT1 µ 2≈ 2 ∆T d

(31a)

∂T0 T1 ≈R 2 ∂Z ∆

(31b)

W1

T

from the balance between the viscous and surface-tension terms in eq 12a and also from that among the terms in eq 3. Now, based on eq 31a, the following amplitude relation is obtained in dimensionless form:

(33)

which satisfy the relation described by eqs 32. Under the relative stability concept of Chen et al.,19 the case of n ) 1/2 satisfies the condition of eq 23, i.e., r0 ) r1 at τ ) τc, which will be shown later. If the related process is still diffusion-dominant with Ma* ) constant at τ ) τc, it is probable that θ(τ,z) ) τ1/2θ*(ζ). This means that the amplitude function of the temperature disturbances follows the behavior of θ0 for small τ. This trend is different from the previous work,17 where the temperature disturbances were set to θ1(τ,z) ) θ*(ζ), to follow the behavior of θ*0(ζ). The relation of Ma* ) constant is shown, even in theoretical results from the frozen-time model,8 the energy method,11 and the maximum-transient-Marangoninumber criterion.12 By the aforementioned scaling reasoning, we set w1 ) τ3/2w*(ζ) and θ1 ) τ1/2θ*(ζ), which are quite different from those in the dominant mode method, where w1 ) τw*(ζ) and θ1 ) θ*(ζ). For a deep-pool system of δT ∝ xτ, the dimensionless time τ is related to the time for the development of thermal penetration depth and it has dual roles of time and length. Now, the self-similar stability equations are obtained from eqs 10 and 11 as

{(

)

2 d2 2 a* + dζ2 1 d2 d d3 ζ 3 - 2 - a*2ζ + 3a*2 w* ) 0 (34) 2Pr dζ dζ dζ

(

(

)}

)

dθ/0 1 d 1 d2 2 2 + θ* ) Ma*a* w* ζ a* 2 dζ dζ2 2 dζ

(35)

where a* ) axτ. The proper boundary conditions are

w* )

dθ* d2w* - θ* ) 0 ) dζ dζ2

w* )

dw* ) θ* ) 0 dζ

(at ζ ) 0)

(as ζ f ∞)

(36a) (36b)

For the limiting case of Pr f ∞, the aforementioned stability equations can be obtained from eqs 14 and 15 by setting ∂θ*/ ∂τ ) 0, axτ ) a*, and Maτ ) Ma*. Therefore, the present propagation theory is the quasi-steady-state approximation (QSSA) in the (τ,ζ) domain. For a given Pr value, Ma* and a* are treated as eigenvalues and the minimum value of Ma* should be determined from the plot of Ma* vs a*, under the principle of exchange of stabilities. In other words, the minimum value of τ (i.e., τc) and its corresponding wavenumber ac are obtained for given values of Pr and Ma. The stability equations described by eqs 34-36 are solved using the outward shooting scheme. To integrate these stability


Figure 3. Amplitude profiles of disturbances at τ ) τc for deep-pool systems under the propagation theory and the dominant mode method.

equations, the proper values of dw*/dζ, d3w*/dζ3, and θ* at ζ ) 0 are assumed for given Pr and a* values. Because of the fact that the stability equations and their boundary conditions are all homogeneous, the value of dw*/dζ at ζ ) 0 can be assigned arbitrarily and the value of the parameter Ma* is assumed. This procedure can be understood easily by taking into account the characteristics of eigenvalue problems. After all the values at ζ ) 0 are provided, this eigenvalue problem can be proceeded numerically. Integration is performed from ζ ) 0 to a fictitious upper boundary with the fourth-order RungeKutta-Gill method. If the guessed values of Ma*, d3w*/dζ3, and θ* at ζ ) 0 are correct, w*, dw*/dζ, and θ* will vanish at the lower isothermal boundary. To improve the initial guesses, the Newton-Raphson iteration is used. When convergence is achieved, the upper boundary for computation is increased by a predetermined value and the aforementioned procedure is repeated. Because the temperature disturbances decay exponentially outside the thermal penetration depth, the incremental change of Ma* also decays rapidly as the fictitious upper boundary thickness increases. This behavior enables us to extrapolate the eigenvalue to the infinite depth. A similar procedure can be applied to solve other sets of stability equations. The marginal stability curve for Pr f ∞ is shown in Figure 2. In this limiting case, the governing equations are reduced to simpler forms, because the inertia terms involving Pr in eq 34 are negligible. Based on the critical conditions predicted by the propagation theory for a deep-pool system, the critical conditions to mark the onset of convective motion are correlated as

[ (0.31 Pr ) ]

τc ) 19.18 1 +

0.65 1/0.65

Ma-1 (for τc < 0.01) (37)

within an error bound of 8%. It is believed that, for given values of Ma and Pr, a fastest-growing mode of infinitesimal disturbances would set in at t ) tc with a ) ac. The aforementioned equations show that τc decreases as Ma and Pr increase. The Pr effect becomes pronounced for Pr < 1, which means the inertia terms make the system more stable. For Pr > 100, τc is almost independent of Pr. Under the critical conditions illustrated above, the amplitude functions of w* and θ* are featured in Figure 3, wherein the eigenfunctions from the dominant mode method and the propagation theory have been normalized by the corresponding

Figure 4. Effect of Pr on the stability condition from the linear stability theory.

maximum magnitude θ*max. As shown in this figure, the eigenfunction distributions of the dominant mode method and the propagation theory are quite different. Based on the distribution of temperature disturbances of the propagation theory, its temporal growth rate can be obtained from eq 33: r0 ) r1 ) 3/(2τ). This means that, for large Ma, the growth rate at τ ) τc is inversely proportional to τc and propagation theory bounds the relative stability concept. Also, based on θ1 ) A0(τ) exp(-ζ2/4) and the relative stability condition of σ ) (1/A0)(dA0/dτ) ) 1/(2τ), the condition of eq 33 (i.e., r0 ) r1 ) 3/(2τ)) is satisfied. Now, the domain of time is extended to τc > 0.01 by keeping eqs 34 and 35 and using eq 8b. In eq 36b, the upper boundary ζ f ∞ is replaced with z ) 1 (i.e., ζ ) 1/xτc) and in eqs 34 and 35, Ma* and a* are replaced by Maτc and axτc. Also, in eq 8b, τ is replaced by τc but ζ is maintained. Because τc is the fixed parameter, the resulting stability equations are functions of ζ only and the physics of eqs 8 and 28 is still alive. For given values of Pr and τc, the minimum value of Ma and its corresponding wavenumber ac are obtained. The solution procedure is almost the same as that in the previous section. The results are summarized in Figure 4, wherein those obtained from the conventional frozen-time model are also shown. For τc < 0.01, the former results are the same as those of the deeppool system (eq 36). For large τc, they approach the well-known critical condition of Mac ) 79.61 (see eq 5), because the basic temperature profile becomes linear. It is known that, for small τ, the frozen-time model yields the lower bounds of τc and the terms involving ∂(‚)/∂τ in eqs 10 and 11 stabilize the system. It is interesting that the propagation theory smoothly yields the stability criteria over the entire domain of time. The conventional frozen-time model neglects the terms involving ∂(‚)/∂τ in eqs 10 and 11 (i.e., the frozen-time model is the QSSA in the (τ,z) domain). This results in (d2/dζ2 a*2)2w* ) 0 and (d2/dζ2 - a*2)θ* ) Ma*a*2w*(dθ* 0/dζ), instead of eqs 34 and 35. The resulting stability criteria become independent of Pr and τc is obtained for a given Ma. For the limiting case of τc f 0, the following results can be derived under the base tempearature profiles of eq 9 and the boundary conditions of eq 36:

w* ) -ζ exp(-a*ζ)

(38a)

5780

θ* )


{

Ma*a* exp(-a*ζ) 2

∫0ζ exp(a*η)w* erfc(η2) dη +

∫

∞

2 1 ∂〈|u1| 〉 ) - 〈|∇u1|2〉 - Ma 2Pr ∂τ

η exp(a*ζ) ζ exp(-a*η)w* erfc dη + 2 ∞ η exp(-a*ζ) 0 exp(-a*η)w* erfc dη (38b) 2

()

∫

∂w

∫0 θ1 ∂z1

〈

(46)

〉

2 ∂θ0 1 ∂〈|θ1| 〉 ) -〈|∇θ1|2〉 - w1θ1 2 ∂τ ∂z

() }

(47)

where and

Ma* )

∫0 (‚) ) ∫x ∫y (‚)z ) 0 dy dx

4

(38c)

∫∞

{1 + exp(2a*)} 0 exp(-a*η) η erfc {η exp(-a*η)} dη 2

()

and

〈(‚)〉 )

Based on these relations, the following asymptotic relations can be obtained for the domain of small τ:

Maτc ) 2

(for τc f 0)

(39a)

acxτc ) 0

(for τc f 0)

(39b)

In the aforementioned derivation, the boundary condition of eq 43 and the periodicity in the x- and y-direction are used. In the present system, the dimensionless energy identity is defined as a linear combination of eqs 36 and 47 with the coupling constant β > 0:

E(τ) )

Acrivos.9

which were already suggested by Vidal and 2.3. Energy Method. Based on eqs 1-3, the disturbances of the base state of eqs 6 and 7 satisfy the following nonlinear system:

∇ ‚ u1 ) 0

{

}

(41)

∂θ0 ∂θ1 ) ∇2θ1 - w1 - u1 ‚ ∇θ1 ∂τ ∂z

(42)

where the subscripts 0 and 1 represent the base and disturbance quantities, respectively. The velocity vector u ) (u,V,w), the temperature θ, the time τ, and the pressure p have the scales of R/d, qwd/k, d2/R, and FRν/d3, respectively. The proper boundary conditions are

∂θ1 )0 ∂z

∫

2

1 ∂θ1 dΩ ) Ω 2 ∂τ

∫

∫

〈

〉

∂θ0 ∂E ) -β〈|∇θ1|2〉 - β w1 θ1 ∂τ ∂z 〈|∇u1|2〉 - Ma

(49)

∂E′ ) -〈|∇θ′1|2 + |∇u1|2〉 + ∂τ

{ 〈

〉

∂θ0 1 θ′ ∂z 1 xγ

xMa xγ w1

∂w

}

∫0 θ′1 ∂z1

(50)

where

1 1 〈u 〉2 + 〈θ′1〉2 2Pr 1 2

The aforementioned relation can be represented as

(at z ) 1)

dE′ ) xMaI - D dτ

(43b)

∫Ω u1 ‚ ∇(p1 + 21u12) dΩ + ∫Ω u1 ‚ ∇2u1 dΩ

u ∇θ1 dΩ + Ω 1

∂w

∫0 θ1 ∂z1

By setting θ′1 ) xβθ1 and β ) Maγ, the aforementioned energy identity is expressed as

(at z ) 0) (43a)

The first two conditions of eq 43a are the kinematic conditions. Now, multiplying eq 41 by u1 and eq 42 by θ1 and integrating them over the system volume Ω, eqs 41 and 42 then become 2 1 ∂u1 dΩ ) Ω 2Pr ∂τ

(48)

This leads to

E′ )

∂θ1 ∂V1 ∂θ1 ∂u1 ) Ma , ) Ma , ∂z ∂x ∂z ∂y

u1 ) θ 1 ) 0

1 1 〈u 〉2 + β〈θ1〉2 2Pr 1 2

(40)

1 ∂ + u1 ‚ ∇ u1 ) - ∇p1 + ∇2u1 Pr ∂τ

w1 ) -

∫Ω (‚) dΩ

where

{ 〈

〉

)

∂w

(51)

}

∫1 θ′1 ∂z1

(52)

D ) 〈|∇θ′1|2 + |∇u1|2〉

(44)

(53)

Under the strong instability concept, the strong stability limit (τs) is defined as the time when dE/dτ ) 0. In other words, for a given τ, the strong stability limit (Mas) is determined as

∫

∂θ

(

∂θ0 1 I ) xγ w1 θ′1 ∂z xγ

θ ∇2θ1 dΩ Ω 1

∫Ω w1θ1 ∂z0 dΩ

I ) -D 1 - xMa D

(45)

Using the divergence theorem, the following relations can be obtained:

1

xMas

( 〈

under the condition of

〉

∂θ0 1 θ′1 ∂z xγ

) max xγ w1

∫0 θ′1

)

∂w1 ∂z

(54)


〈|∇θ1|2 + |∇u1|2〉 ) 1

(55)

This maximum problem can be solved by the variational technique, and, under the normal-mode analysis, the following Euler-Lagrange equations can be obtained:11

( (

) )

2

2

∂θ0 d 1 - a2 w1 ) xγ a2θ1 2 2 ∂z dz 2

(56)

1 ∂θ0 d - a2 θ1 ) - xγ w1 2 2 ∂z dz

For a greater time of τ > 0.01, similar approach can be applied and the terms containing 3/(4τ) should be slightly modified. According to the base temperature profile of eq 9, dE/dτ can be written as

dE dτ

) (r0)E ) ∞

(16/λn2)exp(- λn2τ){1 - exp(-λn2τ)} ∑ n)1

(57) (1/12) -

with the following boundary conditions: 2 1 Ma ∂w1 ∂ w1 1 2Ma w1 ) 2 θ )0 ) 2 - a 2 xγ ∂z 2 xγ 1 ∂z ∂z (at z ) 0) (58a)

∂w1 ) θ1 ) 0 ∂z

(at z ) 1)

(58b)

where a is the horizontal wavenumber. The last two boundary conditions of eq 58a are “natural” boundary conditions20 of the dynamic boundary conditions. As shown in these equations, the strong stability conditions are independent of Pr. The strong stability limit Ma(τs) is given by

min xMa xMas ) max γ a

(at τ ) τr)

(for τ f 0)

(61)

and the relative stability limit can be obtained as

1

xMar

) max

[

I D + 3E/(2τ)

]

(for τ f 0)

(62)

Under the normal-mode analysis, the Euler-Lagrange equations to represent the relative stability limit are obtained for the limiting case of τ f 0:

(

)

(

)

2 1 ∂θ0 2 3 1 ∂2 ∂2 2 x a w ) γ θ + - a2 w1 a 1 1 2 ∂z 4τ Pr ∂z2 ∂z2

(

)

∂2 1 ∂θ0 3 - a2 θ1 ) xγ w1 + θ1 2 2 ∂z 4τ ∂z

(63) (64)

under the boundary conditions of eq 46. The relative stability limit Ma(τr) is given by

min xMa xMar ) max γ a

E(τm) ) E(0)

2

(65)

(at τ ) τm)

(67)

For given values of Ma and τ, they found the fastest-growing disturbances by solving the following maximum problem:

$(τ) ) max

[(

I(τ) D xMa -1 E D

)]

(68)

where $(τ) is the maximum growth rate at τ. At the strong stability limit, the aforementioned equation degenerates to

0 ) max

(60)

where E0 is the basic energy (i.e., E0 ) 〈θ02〉/2, because u0 ) 0). For the limiting case of τ f 0, it is known that (1/E0)(dE0/ dτ) ) 3/(2τ), based on eq 9, and, therefore, the relaxed energy identity for the relative stability becomes

3E ) xMaI - D 2τ

2

From the aforementioned equation, we obtain r0 ) 3/(2τ) for the limiting case of τ f 0 but r0 ) 0 for another limiting case of τ f ∞. The aforementioned stability equations (eqs 63 and 64) degenerate into the strong stability formulation of eqs 56 and 57, respectively. Therefore, the present stability equations can cover the entire time region without any contradictions and assumptions. Gumerman and Homsy11 suggested the marginal stability limit Ma(τm) when the energy identity recovers its initial value; i.e.,

(59)

Now, let us relax the strong stability limits. Under the relative stability concept of Chen et al.,19 the relative instability time τr is determined as

1 dE 1 dE0 ) E dτ E0 dτ

∑(8/λn )exp(-λn τ){2 - exp(-λn τ)} 4

n)1

∂2θ1

w1 )

E (66)

∞

{(

I(τs) D xMa -1 E D

)}

(69)

which is identical to eq 54. The time τm, which bounds the marginal stability limit, can be determined from the following conditions:

∫0τ

m

$(τ′) dτ′ ) 0

(70)

As mentioned by Gumerman and Homsy,11 the computational burden is tremendous, so we do not try to revisit the marginal stability limit. However, for the limiting case of Ma f MaEc , it may be expected that $ e 0 and, therefore, τm f ∞, i.e., the marginal stability limit degenerates to the strong stability limit. The stability limits based on the energy method are summarized in Figure 5. As shown in this figure, the relative stability limit suggests a higher stability bound than the strong and the marginal stability limits. The relative stability and the marginal stability limits seem to converge into the strong stability limit for the limiting case of τ f ∞, as expected. From these observations, the relation of τs e τm e τr seems to hold for the entire time region. The present relative stability results show that the relative stability time τr decreases as Ma and Pr each increase. The Pr effect becomes pronounced for Pr < 1, which means the inertia term of 3/(4τ)(1/Pr)((∂2/∂z2) - a2)w1, which corresponds to 1/(2Pr)((ζd3/dζ3) - (d2/dζ2) - a*2ζ(d/dζ) + 3a*2)w* in eq 34, makes the system more stable. 3. Discussion The stability limits based on the linear stability theory and the energy method are compared in Figure 6. As shown in this

5782


based on the frozen-time model and the energy method based on the strong energy concept. Tan and Thorpe12 suggested a simple instability analysis, assuming that, at the onset of convection, the following relations are maintained, based on eq 5:

maximum of

{ ( )} γZ2 ∂T0 kµ ∂Z

) 79.6

(71a)

ac × 7.26 × Zmax ) 1.993 2π

Figure 5. Effect of Pr on the stability condition from the energy method.

Figure 6. Comparison of critical Marangoni numbers based on the various models. Results from the propagation theory and the modified energy method are given for the case of Pr ) 30; the numbers shown in parentheses indicate the Prandtl number (Pr).

figure, the energy method gives more unstable stability bounds than the linear stability theory, because the energy method introduces the finite disturbance rather than infinitisimal disturbance, as in the linear stability approach. In other words, the relative stability bounds based on the energy method is lower than the stability bounds from the propagation theory, based on the relative stability criterion under the linear stability theory. Therefore, the energy method provides the necessary condition and minimum bound for the onset of instability. It is interesting that, for the limiting case of small τ, the strong stability limit of energy method and the frozen time give similar stability criteria. A similar trend can be found in the work of Ha and Lai.21,22 For this region of τ f 0, the propagation theory and the dominant mode analysis under the linear stability theory give higher stability bounds and predict the experimental data closer than the energy method, based on the relative stability concept. The propagation theory and the dominant mode analysis have their limitations (i.e., both methods are derived under the assumption of Ma . 1 and τ , 1). However, the stability bounds based on the energy method (i.e., the strong, relative, and marginal stability limits) do not use this assumption (i.e., the energy methods derived here can be applied for the entire time domain without any limitations). It is interesting that the stability bounds of the propagation theory and the energy method based on the relative stability concept are higher than those

(71b)

which are satisfied by ∂T/∂Z ) 0.070qw/k at Zmax ) 1.684xRt from eq 9. This results in τo ) 120Ma-1 and ac ) 0.0935Ma1/2, independent of Pr. Therefore, eq 71 seems to correspond to the upper bound of the critical time, wherein the temperature profile is assumed to be linear within Z ) Zmaz. It is well-known that, for deep-pool systems, the boundary conditions on the lower boundary do not affect the critical conditions but the critical condition from this model is dependent on to the conditions of the lower boundary. However, common relations are involved in the aforementioned results: Ma* ≈ constant and a* ≈ constant. In the present study, the time for the amplification of disturbances is not considered. Foster23 commented that, for the case of δT e 0.2, the thermal penetration depth ∆T is a proper length scale and the relation of to = 4tc would be maintained for the case of the buoyancy-driven convection in horizontal fluid layers heated isothermally from below. This means that the growth period for the infinitesimal disturbances to grow is required until they are detected experimentally. Therefore, it seems evident that the predicted onset time tc is smaller than the detection time to. This means that the motion of instabilities, which set in at t ) tc, will grow with time until manifest motion is first detected experimentally. Following Foster’s comments, eq 37, which is obtained from the linear stability theory, using the thermal penetration depth as a length scale, is chosen as the onset time of a fastest-growing mode of instability and it is rearranged as a function of the Marangoni number based on the temperature difference, M ( ) Max4τ/π):

[ ( ) ] ( x)

τc ) 367.9 1 +

0.31 0.65 Pr

2/0.65

π 4

M

-2

(for τc e 0.01) (72)

From the relation of τo = 4τc, the following relation is obtained:

[ ( ) ] ( x)

τo = 4τc ) 1471.6 1 +

0.31 0.65 Pr

2/0.65

-2

π 4 (for τc e 0.01) (73)

M

Because the temperature difference increases continuously during the growth period (τc e τ e τo), the detection time τo may be suggested as

[ (0.31 Pr ) ]

τo ) 38.36 1 +

0.65 1/0.65

Ma-1

(for τc e 0.01) (74)

where Ma ) Mxπ/(4τo). A similar trend is expected to be extended to τc > 0.01. This scenario is supported by the previous experimental and theoretical results for propanol (Pr ) 30), as shown in Figure 7. Although a limited set of experimental data, which are the only reliable experimental data, to our knowledge, is used in the present comparison, the relation of τo = 4τc holds for the various systems such as Rayleigh-Beńard convection,13


Greek Letters R ) thermal diffusivity (m2/s) ∆T ) thermal penetration depth (m) δT ) dimensionless penetration depth; δT ) ∆T/d ζ ) similarity variable; ζ ) z/xτ µ ) viscosity (Pa s) ν ) kinematic viscosity (m2/s) F ) density (kg/m3) τ ) dimensionless time; τ ) Rt/d2 θ1 ) dimensionless perturbed temperature; θ1 ) T1a2sd/(Rµ) or θ1 ) kT1/(qwd) θ0 ) dimensionless basic temperature; θ0 ) k(T - Ti)/(qwd) $ ) maximum growth rate Subscripts Figure 7. Comparison of the experimental results with the present predictions; the numbers shown in parentheses indicate the Prandtl number (Pr).

Soret-driven convection,15 and the Taylor-Go¨rtler vortex.16 It seems evident that convective motion is very weak during tc e t e to, because the related heat transport is well-represented by the conduction state. 4. Conclusions The critical condition to mark the onset of surface-tensiondriven motion in an initially quiescent, horizontal fluid layer experiencing evaporative cooling has been analyzed theoretically. A new set of stability equations is derived under the linear stability theory, the energy method, and their modification. The resulting stability criteria compare reasonably well with extant experimental data and it seems that manifest convection seems to be observable in to for the present systems. For τ eτo, the velocity disturbances seem to be too weak to be observable experimentally. It is very interesting that our modification to the linear stability theory and energy method covers the entire domain of time and the limiting stability criteria of τc g 100 approach the well-known values of eq 5. Nomenclature a ) dimensionless horizontal wave number, xax2+ay2 a* ) modified wave number, axτ d ) depth of liquid layer (m) E ) dimensionless natural energy defined in eq 25 k ) thermal conductivity (W m-1 K-1) M ) Marangoni number based on the temperature difference; M ) Max4τ/π Ma ) Marangoni number; Ma ) sqwd2/(kRµ) Ma* ) modified Marangoni number; Ma* ) Maτ P ) pressure (Pa) Pr ) Prandtl number; Pr ) ν/R q ) heat flux (W/m2) s ) negative rate of change of surface tension, with respect to temperature (K-1) T ) temperature (K) t ) time (s) U ) velocity vector (m/s) u ) dimensionless velocity vector W ) vertical velocity (m/s) w ) dimensionless vertical velocity x, y, z ) dimensionless Cartesian coordinates

c ) critical state i ) initial state m ) marginal instability condition o ) observable condition r ) relative instability condition s ) strong instability condition w ) surface condition x ) x-direction y ) y-direction 0 ) basic state 1 ) perturbed state Superscripts * ) variables in the (τ,ζ) domain L ) conditions under the linear stability theory E ) conditions under the energy method Acknowledgment This work was supported by a grant from the Chuongbong Academic Research Fund of the Cheju National University Development Foundation. Literature Cited (1) Tan, K. K.; Thorpe, R. B. The onset of convection induced by buoyancy during gas diffusion in deep fluids. Chem. Eng. Sci. 1999, 54, 4179. (2) Brian, P. L. T. Effects of Gibbs adsorption on Marangoni instability. AIChE J. 1971, 7, 765. (3) Patzer, J. F.; Homsy, G. M. Global stability of transient drop extraction to Marangoni instabilities. Phys. Fluids 1981, 24, 567. (4) Mendes-Tatsis, M. A.; Perez de Ortiz, E. S. Marangoni instabilities in systems with an interfacial chemical reaction. Chem. Eng. Sci. 1996, 51, 3755. (5) Sun, Z. F.; Yu, K. T.; Wang, S. Y.; Miao, Y. Z. Absorption and desorption of carbon dioxide into and from organic solvents: Effects of Rayleigh and Marangoni instability. Ind. Eng. Chem. Res. 2002, 41, 1905. (6) Sun, Z. F.; Fahmy, M. Onset of Rayleigh-Benard-Marangoni convection in gas-liquid mass transfer with two-phase flow: Theory. Ind. Eng. Chem. Res. 2006, 45, 3293. (7) Pearson, J. R. A. On convection cells induced by surface tension. J. Fluid Mech. 1958, 4, 489. (8) Davis, S. H. Buoyancy-surface tension instability by the method of energy. J. Fluid Mech. 1969, 39, 347. (9) Vidal, A.; Acrivos, A. Effect of nonlinear temperature profiles on the onset of convection driven by surface tension gradients. Ind. Eng. Chem. Fundam. 1968, 7, 53. (10) Kang, K. H.; Choi, C. K.; Hwang, I. G. Onset of solutal Marangoni convection in a suddenly desorbing liquid layer. AIChE J. 2000, 46, 15. (11) Gumerman, R. J.; Homsy, G. M. The stability of uniformly accelerated flow with application to convection driven by surface tension. J. Fluid Mech. 1975, 68, 191. (12) Tan, K. K.; Thorpe, R. B. On convection driven by surface tension caused by transient heat conduction. Chem. Eng. Sci. 1999, 54, 775.

5784


(13) Kim, M. C.; Yoon, D. Y.; Choi, C. K. Onset of buoyancy-driven instability in gas diffusion systems. Ind. Eng. Chem. Res. 2006, 45, 7321. (14) Park, J. H.; Chung, T. J.; Choi, C. K.; Kim, M. C. The onset of mixed convection in a porous layer heated with constant heat flux. AIChE J. 2006, 53, 2677. (15) Kim, M. C.; Hong, J. S.; Choi, C. K. The analysis of onset of Soretdriven convection in nanoparticles suspension. AIChE J. 2006, 53, 2333. (16) Kim, M. C.; Choi, C. K. The onset of Taylor-Go¨rtler vortices during impulsive spin-down to rest. Chem. Eng. Sci. 2006, 61, 6478. (17) Kang, K. H.; Choi, C. K. A theoretical analysis of the onset of surface-tension-driven convection in a horizontal liquid layer cooled suddenly from above. Phys. Fluids 1997, 9, 7. (18) Riaz, A.; Hesse, M.; Tchelepi, H. A.; Orr, F. M., Jr. Onset of convection in a gravitationally unstable, diffusive boundary layer in porous media. J. Fluid Mech. 2006, 548, 87. (19) Chen, J.-C.; Neitzel, G. P.; Jankowski, D. F. The influence of initial

condition on the linear stability of time dependent circular Couette flow. Phys. Fluids 1985, 28, 749. (20) Courant, R.; Hilbert, D. Methods of Mathematical Physics, Vol. 1; Interscience: New York, 1937. (21) Ha, V.-M.; Lai, C.-L. The onset of stationary Marangoni instability of an evaporating droplet. Proc. R. Soc. London Ser. A 2001, 457, 885. (22) Ha, V.-M.; Lai, C.-L. Theoretical analysis of Marangoni instability of an evaporating droplet by energy method. Int. J. Heat Mass Transfer 2004, 47, 3811. (23) Foster, T. D. Onset of manifest convection in a layer of fluid with a time-dependent surface temperature. Phys. Fluids 1969, 12, 2482.

ReceiVed for reView January 2, 2007 ReVised manuscript receiVed April 3, 2007 Accepted June 8, 2007 IE070006Y

Onset of Marangoni Convection in a Horizontal Fluid Layer

Recommend Documents