Energy Stability Analysis on the Onset of Buoyancy-Driven Convection in a Horizontal Fluid Layer Subject to Evaporative Cooling

Abstract

The onset of buoyancy-driven convection in an initially isothermal and quiescent horizontal fluid layer was analyzed theoretically. It is well-known that at the critical Rayleigh number $Ra_c=669$ convective motion sets in with a constant-heat-flux cooling through the upper boundary. Here, based on the momentary instability concept, the dimensionless critical time ${\tau}_m$ to mark the onset of convective motion for Ra > 669 was analyzed theoretically. The energy method under the momentary stability concept was used to find the critical conditions as a function of the Rayleigh number Ra and the Prandtl number Pr. The predicted critical conditions were compared with the previous theoretical and experimental results. The momentary stability criterion gives more reasonable wavenumber than the conventional energy method.

1. Introduction

Buoyancy-driven convection plays an important role in many engineering problems, such as chemical vapor deposition, solidification, electroplating, and also conventional heat and mass transfer systems. Most of these processes involve nonlinear developing temperature profile, so it becomes important to predict when the buoyancy-driven convection sets in. But a general approach to predict the critical conditions to mark the onset of buoyancy-driven convection under these circumstances is still under controversy.

When an initially quiescent, horizontal fluid layer is cooled from above or heated rapidly from below, the basic temperature profile of heat conduction develops with time and buoyancy-driven convection setting in at a critical time. In this transient system, the critical time t_c to mark the onset of convective motion becomes an important question, which may be called an extension of classical RayleighBénard problems. The related instability analyses have been conducted under linear stability theory and nonlinear energy method. Based on the linear stability theory, the frozen-time model [1], amplification theory [2] and propagation theory [3] have been derived and applied to the various systems. Homsy [4], Gummerman and Homsy [5], and Wankat and Homsy [6] applied the nonlinear energy method to analyze this kind of problem. Straughan [7] summarized the theoretical aspects of the energy method for the various systems. Also, Harfash and Straughan [8] used the energy method to study the magnetic field effect on the convective instability.

Based on the relative stability concept, Kim and colleagues [9-14] analyzed the energy stability of the various systems. Their relaxed energy method gives the critical time t_c for the whole range of Pr and Pa, but the conventional energy method yields the stability criteria independently of Pr. Here we concentrated 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 stability criteria were obtained based on the original energy method and its modification, and they were compared with available experimental and theoretical results.

2. Theoretical analysis

2-1. Governing equations and base system

The system considered here is a Newtonian fluid layer with an initial temperature. For time t ≥ 0, the horizontal layer of fluid depth, d, experiences evaporative cooling with heat flux, q, through the upper free boundary, and its lower boundary is kept at the initial temperature, T_i. A schematic diagram of the basic system of pure conduction is shown in Fig. 1. For a high q, buoyancy-driven convection will set in at a certain time, and the governing equations of flow and temperature fields are expressed by employing the Boussinesq approximation as

\(\nabla \cdot \mathbf{U}=0\),             (1)

\(\left\{\frac{\partial}{\partial t}+\mathbf{U} \cdot \nabla\right\} \mathbf{U}=-\frac{1}{\rho} \nabla P+v \nabla^{2} \mathbf{U}+\mathbf{g}\left\{1-\beta\left(T-T_{i}\right)\right\}\),             (2)

\(\left\{\frac{\partial}{\partial t}+\mathbf{U} \cdot \nabla\right\} T=\alpha \nabla^{2} T\),             (3)

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

where \(\mathbf{U}(=(U, V, W))\) is the velocity vector, ρ the density, P the dynamic pressure, ν the kinematic viscosity, T the temperature, g the gravitational acceleration, β the thermal expansion coefficient, and α the thermal diffusivity

Let’s assume the evaporation rate and corresponding heat flux q are constant. The validity of this assumption will be discussed later. Then, the basic state of heat conduction the dimensionless temperature profile is represented by [15].

\(\frac{\partial \theta_{0}}{\partial \tau}=\frac{\partial^{2} \theta_{0}}{\partial z^{2}}\),             (4)

with the following initial and boundary conditions,

\(\theta_{0}= 0 at \tau=0 and z=1\),             (5a)

\(\frac{\partial \theta_{0}}{\partial z}=-1 \text { at } z=0\),             (5b)

where τ = αt/d² , z = Z/d and θ₀ = k(T-T_i)/(qd). Here, k is the thermal conductivity of the fluid. The subscript ‘0’ denotes the basic state. The exact solution of Eqs. (4) and (5) is well-known:

\(\theta_{0}=1-z-2 \sum_{n=1}^{\infty} \frac{1}{\mu_{n}} \cos \left(\mu_{n} z\right) \exp \left(-\mu_{n}^{2} \tau\right)\),             (6a)

\(\theta_{0}^{*}=\sqrt{4 \tau} \sum_{n=0}^{\infty}(-1)^{n}\left\{\text { ierfc }\left(\frac{n}{\sqrt{\tau}}+\frac{\zeta}{2}\right)-\operatorname{ierfc}\left(\frac{n+1}{\sqrt{\tau}}-\frac{\zeta}{2}\right)\right\}\),             (6b)

where θ₀(τ,z) = θ₀(τ,ζ), μ_n = (n-1/2)π, and ζ = z/$\sqrt{\tau}$. Equation (6b) is obtained in terms of the integral of complementary error functions by using the Laplace transform.

2-2. Stability equations

Consider the following velocity, pressure and temperature perturbations: U₁ = U−U₀, P₁ = P−P₀ and T₁ = T−T₀, and let’s introduce these perturbations into Eqs. (1)-(3). Then, using α/d, ρα² /d² , and qd/k as the scaling factors of velocity, pressure and temperature, respectively, we can obtain the following dimensionless equations:

\(\nabla \cdot \mathbf{u}_{1}=0\),             (7)

\(\frac{1}{\operatorname{Pr}}\left\{\frac{\partial}{\partial \tau}+\mathbf{u}_{1} \cdot \nabla\right\} \mathbf{u}_{1}=-\nabla p_{1}+\nabla^{2} \mathbf{u}_{1}+\mathbf{k} R a \theta_{1}\),             (8)

\(\frac{\partial \theta_{1}}{\partial \tau}=\nabla^{2} \theta_{1}-w_{1} \frac{\partial \theta_{0}}{\partial z}-\mathbf{u}_{1} \cdot \nabla \theta_{1}\),             (9)

under the following boundary conditions:

\(\mathbf{u}_{1}=\frac{\partial \theta_{1}}{\partial z}=0 \text { at } z=0\),             (10a)

\(\mathbf{u}_{1}=\theta_{1}=0 \text { at } z=1\),             (10b)

where k is the unit vector of the positive z-direction, and subscripts 0 and 1 represent the base and perturbation quantities, respectively. Here, \(\operatorname{Pr}(=v / \alpha)\) and \(R a\left\{=g \beta q d^{4} /(k \alpha v)\right\}\) are the Prandtl number and the Rayleigh number, respectively

Now, multiply Eq. (8) by u₁ and Eq. (9) by θ₁ and integrate over the system volume Ω, then Eqs. (8) and (9) become,

\(\int_{\Omega} \frac{1}{2 P r} \frac{\partial \mathbf{u}_{1}^{2}}{\partial \tau} d \Omega=-\int_{\Omega} \mathbf{u}_{1} \cdot \nabla\left(p_{1}+\frac{1}{2} \mathbf{u}_{1}^{2}\right) d \Omega+\int_{\Omega} \mathbf{u}_{1} \cdot \nabla^{2} \mathbf{u}_{1} d \Omega+R a \int_{\Omega} \theta_{1} w_{1} d \Omega\),             (11)

\(\int_{\Omega} \frac{1}{2} \frac{\partial \theta_{1}^{2}}{\partial \tau} d \Omega=-\int_{\Omega} \mathbf{u}_{1} \nabla \theta_{1} d \Omega+\int_{\Omega} \theta_{1} \nabla^{2} \theta_{1} d \Omega-\int_{\Omega} w_{1} \theta_{1} \frac{\partial \theta_{0}}{\partial z} d \Omega\),             (12)

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

\(\frac{1}{2} \frac{1}{P r} \frac{\partial\left\langle\left|\mathbf{u}_{1}\right|^{2}\right\rangle}{\partial \tau}=-\left\langle\left|\nabla \mathbf{u}_{1}\right|^{2}\right\rangle+R\left\langle w_{1} \theta_{1}^{\prime}\right\rangle\),             (13)

\(\frac{1}{2} \frac{\partial\left\langle\left|\theta_{1}^{\prime}\right|^{2}\right\rangle}{\partial \tau}=-\left\langle\left|\nabla \theta_{1}^{\prime}\right|^{2}\right\rangle-R\left\langle w_{1} \theta_{1}^{\prime} \frac{\partial \theta_{0}}{\partial z}\right\rangle\),             (14)

where R = \(\sqrt{Ra}\), θ₁' = \(\sqrt{Ra}θ_1\) and <(·)> = \(\int_{0}{(\cdot) \mathrm{d} \Omega}\). In above derivation, the boundary condition of Eq. (10) and the periodicity in x- and y-direction are used.

In the present system the dimensionless natural energy can be defined as a linear combination of Eqs. (13) and (14) with the coupling constant γ > 0:

\(E(\tau)=\frac{1}{2 \operatorname{Pr}}\left\langle\mathbf{u}_{1}\right\rangle^{2}+\frac{1}{2} \gamma\left\langle\theta_{1}^{\prime}\right\rangle^{2}\),             (15)

and the following energy identity can be derived

\(\frac{\partial E}{\partial \tau}=-\gamma\left\langle\left|\nabla \theta_{1}\right|^{2}\right\rangle-\gamma R\left\langle w_{1} \frac{\partial \theta_{0}}{\partial z} \theta_{1}\right\rangle+R\left\langle w_{1} \theta_{1}\right\rangle-\left\langle\left|\nabla \mathbf{u}_{1}\right|^{2}\right\rangle\),             (16)

where w₁ is the vertical component of the velocity perturbation vector and the primes are dropped for the sake of simplicity. By setting \(\hat{\theta}_{1}=\sqrt{\gamma} \theta_{1}\), the above energy identity can be expressed as

\(\frac{\partial \hat{E}}{\partial \tau}=-\left\langle\left|\nabla \hat{\theta}_{1}\right|^{2}+\left|\nabla \mathbf{u}_{1}\right|^{2}\right\rangle+R\left\langle w_{1} \frac{\hat{\theta}_{1}}{\sqrt{\gamma}}-w_{1} \frac{\partial \theta_{0}}{\partial z} \sqrt{\gamma} \hat{\theta}_{1}\right\rangle\),             (17)

where \(\hat{E}=\frac{1}{2 P r}\left\langle\mathbf{u}_{1}\right\rangle^{2}+\frac{1}{2}\left\langle\hat{\theta}_{1}\right\rangle^{2}\). After dropping the hats, the above relation can be represented as

\(\begin{aligned} \frac{d E}{d \tau} &=R I-D \\ &=-D\left(1-\frac{I}{D} R\right) \end{aligned}\),             (18)

where

\(I=\left\langle w_{1} \frac{\theta_{1}}{\sqrt{\lambda}}-w_{1} \frac{\partial \theta_{0}}{\partial z} \sqrt{\lambda} \theta_{1}\right\rangle\),             (19)

\(D=\left\langle\left|\nabla \theta_{1}\right|^{2}+\left|\nabla \mathbf{u}_{1}\right|^{2}\right\rangle\).             (20)

Even though we can describe the temporal evolution of the perturbation energy through Eq. (18), care must be taken in defining the stability criterion for the system having time-dependent base states. Shen [16] first observed, in a study of time dependent parallel shear flow, if the kinetic energy of a perturbation decreases in time but that of the base state decreases at a faster rate, then the kinetic energy of perturbation will appear amplified in time. Conversely, if the kinetic energy of the perturbation increases in time but that of the base state increases faster still, then the kinetic energy of the perturbation will appear to decay in time. To determine the stability characteristics of perturbations of time variant base states, Shen [16] introduced the concept of “momentary stability” where the stability of the system is guaranteed if

\(\frac{d E_{R}}{d \tau}<0\),             (21)

where

\(E_{R}=\frac{E}{E_{0}}\),             (22)

is called the relative energy [17] and the momentary stability has been known as relative stability [18]. Here E and E₀ are the energy of the disturbance and that of base state, respectively. For the present system, E₀ is defined as

\(E_{0}=\frac{1}{2}\left\{\int_{0}^{1} \theta_{0}^{2} d z\right\}\),             (23)

With these definitions, the criterion for momentary stability of unsteady base state is given by

\(\frac{1}{E_{R}} \frac{d E_{R}}{d \tau}=\sigma-\sigma_{0}\),             (24)

Here σ and σ₀ are the growth rate of the disturbance and that of base energy defined as

\(\sigma=\frac{1}{E} \frac{d E}{d \tau} \text { and } \sigma_{0}=\frac{1}{E_{0}} \frac{d E_{0}}{d \tau}\),             (25)

For the present system, based on Eq. (6a) the growth rate of base energy is

\(\sigma_{0}=\frac{\sum_{n=1}^{\infty}\left(8 / \mu_{n}^{4}\right) \exp \left(-\mu_{n}^{2} \tau\right)\left\{1-\exp \left(-\mu_{n}^{2} \tau\right)\right\}}{1 / 24-\sum_{n=1}^{\infty}\left(16 / \mu_{n}^{4}\right) \exp \left(-\mu_{n}^{2} \tau\right)\left\{2-\exp \left(-\mu_{n}^{2} \tau\right)\right\}}\),             (26)

For the limiting case of τ→∞, σ₀→0 is obtained.

It is well-known that the present system with Ra > 669 is asymptotically unstable [19]. Therefore, our primary concern is the instantaneous instability, which is defined as

σ > σ₀,             (27)

under the momentary instability concept [16]. The neutral stability condition under the momentary instability can be determined from

\(\sigma_{0} E=R I-D\),             (28)

And, therefore the momentary stability limit can be obtained as

\(\frac{1}{R}=\max \left[\frac{I}{D+\sigma_{0} E}\right]\),             (29)

under the condition of

\(D=\left\langle\left|\nabla \theta_{1}\right|^{2}+\left|\nabla \mathbf{u}_{1}\right|^{2}\right\rangle=1\),             (30)

This maximum problem can be solved by the variational technique. And, under the normal mode analysis, the following EulerLagrange equations can be obtained:

\(\left(\frac{d^{2}}{d z^{2}}-a^{2}\right)^{2} w_{1}=-\frac{1}{2} R\left(\frac{1}{\sqrt{\gamma}}-\sqrt{\gamma} \frac{\partial \theta_{0}}{\partial z}\right) a^{2} \theta_{1}+\frac{\sigma_{0}}{2 P r}\left(\frac{d^{2}}{d z^{2}}-a^{2}\right) w_{1}\),             (31)

\(\left(\frac{d^{2}}{d z^{2}}-a^{2}\right) \theta_{1}=\frac{1}{2} R\left(\frac{1}{\sqrt{\gamma}}-\sqrt{\gamma} \frac{\partial \theta_{0}}{\partial z}\right) w_{1}+\frac{\sigma_{0}}{2} \theta_{1}\),             (32)

under the following boundary conditions:

\(w_{1}=\frac{d^{2} w_{1}}{d z^{2}}=\frac{d \theta_{1}}{d z}=0 \text { at } z=0\),             (33a)

\(w_{1}=\frac{d w_{1}}{d z}=\theta_{1}=0 \text { at } z=1\),             (33b)

The momentary stability limit Ra is given by

\(\sqrt{R a}=\max _{γ} \min _{a} R\),             (34)

Since σ₀→0 as τ→∞, for the limiting case of large τ, the above stability equations degenerate into the conventional strong stability equations.

3. Solution Method

The stability equations (31)-(33) were solved by employing the outward shooting scheme [20]. To integrate them, trial values of the eigenvalue R and the boundary conditions d³ w₁/dz³ and θ₁ at z = 0 are assumed properly for a given a and γ. Since the boundary condition, Eq. (33) are all homogeneous, the value of dw₁/dz at z = 0 can be assigned arbitrarily. This procedure is based on the outward shooting method in which the boundary value problem is transformed into the initial value problem. The trial values, together with the three known conditions at the lower boundary, give all the information to make numerical integration smooth.

The integration based on the 4th-order Runge-Kutta method is performed from z = 0 to z = 1. By using the Newton-Raphson iteration the trial values of R, d³ w₁/dz³ , and θ₁ are corrected until the stability equations satisfy the upper boundary conditions within the relative tolerance of 10^-10. For the strong stability limits, the solution procedure is almost the same as above.

4. Results and Discussion

Since, for large τ, the base temperature field becomes linear and therefore σ₀→0, the present momentary stability degenerates to the conventional energy method. For this case, the critical condition is Ra = 669 [19]. The present stability limits given in Fig. 2 reconstruct this condition. By employing the momentary stability concept, we tried to relax the conventional energy method and to reanalyze the well-known transient Rayleigh-Bénard problem. The present relaxation can show the Prandtl number effect on the stability conditions, which has been ignored in the original energy method based on the strong stability criterion. The present relaxation shows that the critical time τ_m based on the momentary stability concept decreases with an increase in Ra and also Pr. The Pr-effect becomes pronounced for Pr < 1, which means the inertia term \(\frac{\sigma_{0}}{2 P r}\left(\frac{d^{2}}{d z^{2}}-a^{2}\right) w_{1}\) in Eq. (31) makes the system more stable.

Fig. 2. Effect of Pr on the stability condition.

For the isothermally heated system, Neitzel [21] reported the global stability limits under the conventional energy stability method. The global limits are Ra = 1699 (at τ≈0.14) and 1013 (at τ≈0.08 ) for the rigid-rigid boundaries system and free-rigid boundaries one, respectively. This global stability limits are lower than the asymptotic stability limits, which are 1708 and 1101, respectively. However, this global minimum cannot be shown in the rigid-free boundaries system. The free-rigid boundaries system corresponds to rigid-free boundaries system cooled from free, upper boundary, which is similar situation to the present system. Kim et al. [13] reconsidered this problem using the relative stability criterion. According to their results, the global minimum was not observed for the various boundary situations. And, as shown in Fig. 2, the present system does not show the global minimum. Therefore, the global minimum phenomena seem to be dependent on the boundary conditions, heating or cooling history and stability criteria.

Wankat and Homsy [6] analyzed the stability condition and suggested the minimum bound of stability of the system similar to the present one. However, they simulated the evaporative cooling as the ramp cooling rather than the present constant flux cooling and used free-free boundary conditions for the upper and lower boundaries. Furthermore, they employed the strong instability criterion σ > 0 rather than present momentary instability σ > σ₀ as an instability criterion. They compared their results with Foster’s [22] experiments for water layer, where the top boundary was cooled by evaporation. Manifest convection was detected first at t = t₀ by visual observation of the motion of a thin layer of ink near the bottom layer. The typical surface temperature record, Figure 2 of Foster’s, might be represented well by the constant flux cooling model, \(\left(T_{s}-T_{i}\right) \sim \sqrt{t}\), rather than the ramp cooling one, \(\left(T_{s}-T_{i}\right) \sim t\). In the present study the experimental data are converted,

\(R a_{\phi} \tau=R a \sqrt{4 \tau} \operatorname{ierfc}(0)\),             (35)

where \(R a_{\phi}=g \beta \phi d^{5} /\left(\alpha^{2} v\right)\) is the Rayleigh number defined by Foster based on the cooling rate φ. Due to the differences described above, direct comparison with Wankat and Homsy’s [6] work is not possible.

Foster [22] analyzed the stability limits using the amplification theory based on the ramp cooling model. In comparing his prediction with his experimental data, he argued that the amplification of the initial disturbances of somewhere between 10 and 100 is necessary for the detection of manifest convection, as shown in Fig. 3. He defined the amplification factor $\overline{\mathcal{W}}$ as the ratio of the root-meansquare quantity of velocity disturbances at t = t_o to that of the assumed white-noise ones at t = 0. But we do not know what initial conditions exist in nature. In Fig. 3, the present predictions are compared with Foster’s theoretical and experimental work. The present predictions are quite different from the experimental data. However, the present momentary stability criterion is much closer to the experimental data than the conventional strong stability criterion. The quantitative discrepancy between the experimental data and the predictions based on the energy method is of foreknowledge, since the energy methods need not satisfy the dynamical equations which describe the actual experimental process. Therefore, the energy methods have been employed as lower bounds of stability. Since no experimental data lie to the left of the energy stability limits, the present predictions do not commit the theoretical base point. It is well-known that the energy method cannot give the information on the crtical wave number. However, the present modification gives a reasonable wavenumber at the onset of convection. As shown in Fig. 4, the present momentary stability criterion gives a more reasonable wavenumber than the original energy method, where the critical wavenumer is a_c = 2.08 for the whole range of Ra.

Fig. 3. Comparison of critical Rayleigh numbers with previous results.

Fig. 4. Comparison of the critical waver number with the experimental results.

5. Conclusions

The critical condition to mark the onset of convective motion driven by buoyancy forces in an initially quiescent, horizontal layer cooled from above was analyzed based on the energy method. By considering the growth rate of the relative energy, we modified the conventional energy method. Based on the present modification, we defined the momentary stability time τ_m from which the growth rate of the perturbation energy exceeds that of the base energy. The present modification predicts experimental trends which cannot be explained by the original energy method. Since the energy methods need not satisfy the dynamical equations such as Navier-Stokes equation and the heat transport equation, the growth of disturbance should be studied by solving dynamical governing equations fully.