Non-Isothermal Spontaneous Imbibition Process Including Condensation Effects and Variable Surface Tension

Due to the continued miniaturization of semiconductor devices, power electronics, biosensors and aerospace equipment, problems associated with overheating of these components have increased. Accordingly, the innovative cooling techniques are required to meet the demands of heat load removal from highly integrated electronic circuits and the electronic components of spacecraft designed for advanced long-term spacecraft missions. Those demands result in the rapid development of improved heat rejection techniques such as the interfacial vaporization and condensation heat transfer of thin liquid films. Closed two-phase devices such as heat pipes and thermosyphons have been and are being used successfully for the above application. Whatever configuration is used, the heat energy removed at the chip is transported away and rejected from the system by condensation at a remote location. Therefore, a fundamental understanding of the condensation process in minichannels and capillaries is important to optimize design considerations. In this direction, for possible use in electronic cooling applications, Begg et al. [1] developed a mathematical model of annular film condensation in a miniature tube. In this model, the liquid flow has been coupled with the vapor flow along the liquid-vapor interface through the interfacial temperature, heat flux, shear stress and pressure jump conditions. The model predicts the position of the liquid-vapor interface. The numerical results show that complete condensation of the incoming vapor is possible at comparatively low heat loads. L. P Yarin et al. [2] present a quasi-one dimensional model of laminar flow in a heated capillary. In the frame of this model, the effects of channel size, initial temperature of the working fluid, wall heat flux and gravity on two phase capillary flow are studied. It is shown that hydrodynamical and thermal characteristics of laminar flow in a heated capillary are determined by the physical properties of the liquid and its vapor, as well as the heat flux at the wall. The effect of dimensionless parameters such as the Peclet, Jakob numbers, and dimensionless heat flux or Nusselt number on the velocity, temperature and pressure within the liquid and vapor domains has been studied. In addition, the above authors conducted an experimental analysis for showing that the flow in micro-channels appear to have distinct phase domains. On the other hand, the theory of


Introduction
Due to the continued miniaturization of semiconductor devices, power electronics, biosensors and aerospace equipment, problems associated with overheating of these components have increased. Accordingly, the innovative cooling techniques are required to meet the demands of heat load removal from highly integrated electronic circuits and the electronic components of spacecraft designed for advanced long-term spacecraft missions. Those demands result in the rapid development of improved heat rejection techniques such as the interfacial vaporization and condensation heat transfer of thin liquid films. Closed two-phase devices such as heat pipes and thermosyphons have been and are being used successfully for the above application. Whatever configuration is used, the heat energy removed at the chip is transported away and rejected from the system by condensation at a remote location. Therefore, a fundamental understanding of the condensation process in minichannels and capillaries is important to optimize design considerations. In this direction, for possible use in electronic cooling applications, Begg et al. [1] developed a mathematical model of annular film condensation in a miniature tube. In this model, the liquid flow has been coupled with the vapor flow along the liquid-vapor interface through the interfacial temperature, heat flux, shear stress and pressure jump conditions. The model predicts the position of the liquid-vapor interface. The numerical results show that complete condensation of the incoming vapor is possible at comparatively low heat loads. L. P Yarin et al. [2] present a quasi-one dimensional model of laminar flow in a heated capillary. In the frame of this model, the effects of channel size, initial temperature of the working fluid, wall heat flux and gravity on two phase capillary flow are studied. It is shown that hydrodynamical and thermal characteristics of laminar flow in a heated capillary are determined by the physical properties of the liquid and its vapor, as well as the heat flux at the wall. The effect of dimensionless parameters such as the Peclet, Jakob numbers, and dimensionless heat flux or Nusselt number on the velocity, temperature and pressure within the liquid and vapor domains has been studied. In addition, the above authors conducted an experimental analysis for showing that the flow in micro-channels appear to have distinct phase domains. On the other hand, the theory of two-phase laminar flow in a heated microchannels was presented by Yarin et al. [3]. They studied the thermohydrodynamic characteristics of a two-phase capillary flow with phase change at the meniscus by using a quasi-one-dimensional model for the flow. It takes into account the principal characteristics of the phenomenon, namely, the effects of the inertia, pressure, gravity, friction forces and capillary pressure due to the curvature of the interface surface, as well as the thermal and dynamical interactions of the liquid and vapor phases. To describe the flow outside of the meniscus in the domains of the pure liquid or vapor, the one-dimensional mass, momentum and energy equations are used. The possible states of the flow are considered, and the domains of steady and unsteady states are outlined. Meanwhile, Qu et al. [4] established the physical and mathematical models to account for the formation of evaporating thin liquid film and meniscus in capillary tubes. Their results show that in regard to the capillary tubes of micron scale, the calculation results show that, the bigger the inner radius or the smaller the heat flow, the longer the evaporating interfacial region will be. There only exists meniscus near the wall, and nearby the axial center is flat interface. While as to the capillary tubes of scale about 100 μm, the evaporating interfacial region will increase with heat flux. Compared to the capillaries of micron scale, the meniscus region will extend to the center of capillary axis. Recently, the Lucas-Washburn equation [5,7], was extended by Ramon and Oron [8], describing the motion of a liquid body in a capillary tube so as to account for the effect of interfacial mass transfer due to phase change, either evaporation or condensation. They showed that the phase change affects the equilibrium height of the meniscus, the transition threshold from monotonic to oscillatory dynamics, and the frequency of oscillations. At higher mass transfer rates and/or large capillary radii, vapor recoil is found to be the dominant factor. In general, evaporation decreases the equilibrium height, increases the oscillation frequency and diminishes the transition threshold to oscillations. For condensation, two regimes are identified: at high mass transfer rates similar trends to those of evaporation are observed, whereas the opposite is found for low mass transfer rates, resulting in an increased equilibrium height, lower oscillation frequencies and a shift of the transition threshold toward monotonic dynamics. It must be noted that in the last mentioned work, a difference of temperature or interfacial temperature resistance at the interface is imposed for causing the phase change process. However, the temperature of the liquid bulk was assumed constant. We consider that a more realistic case occurs when the temperature field in the liquid is taken into account, causing a local variation of the interfacial temperature. In fact, this point was commented by Ramon and Oron [8]. Therefore, the main purpose of this paper is to solve the Lucas-Washburn equation in conjunction with the energy equation for the liquid penetrating a capillary tube, considering that the imbibition front is also controlled by a direct condensation process and considering a linear dependence on temperature of the surface tension.

Formulation
A long vertical capillary tube of radius R, as shown in Figure 1, is filled with a saturated vapor at temperature T s whose density is ρ v . At time t = 0, the bottom of the capillary tube comes into contact with a large reservoir of a liquid whose density is ρ, viscosity μ and thermal diffusivity α, assumed constant. The liquid reservoir is kept at temperature T 0 , which is below the vapor temperature T s ; for t > 0, the contact originates a spontaneous non-isothermal imbibition process of the liquid into the capillary tube. It is assumed that the liquid wets the capillary inner wall completely, for which case the contact angle, θ, is set equal to zero [9]. For instance, a typical substance that satisfies this condition for the equilibrium contact angle, is the silicone oil [10]. Due the temperature difference, T s − T 0 , condensation occurs at the moving vapor-liquid interface. In addition, we assume that the temperature interface is maintained at T s . Using a one-dimensional formulation of the average conservation laws, we derive the corresponding dimensionless momentum and energy equations for the liquid penetrating the capillary. In addition, the imbibition front is characterized by a uniform capillary pressure. In this manner, the fluid flow is governed by the continuity, momentum and energy equations, given by where v (r, t) and p (r, t) represent the velocity and pressure fields, respectively. In cylindrical coordinates r = (r, θ, z), with the origin of the coordinates at the tube inlet and the z axis along its axis, the radial distance (0 ≤ r ≤ R) is measured from the z axis ( Figure 1). The velocity components v = (v r , v θ , v z ) are along the r, θ and z, respectively. In this work, unidirectional flow is assumed and admits the azimuthally symmetrical solution v =v z z, with v θ = v r = 0. Therefore, Eq. (1) yields ∂v z /∂z = 0, which means that the v z is a function of the radial coordinate and time, i.e., v z = v z (r, t). However, at the capillary inlet and in a small zone beneath the meniscus interface, the velocity deviates from its value v z (r, t) inside the column and depends on z. Therefore, with Eq. (2), the z component acquires the form We assume further that the fully developed velocity profile in the capillary tube that gives rise to the viscous force is given by the equation is the flow velocity averaged over the capillary cross section. In the specialized literature, it is accepted that the parabolic velocity profile in a pipe is well established at a distance Z from the fluid entry such that Z/d ≈ Re/30, where d is the pipe diameter and Re is the Reynolds number, defined as For typical values of capillary rise, where Re < 1, the profile of Eq. (5) is well established at a distance smaller than the capillary radius.
Because it is required to determine the temporal change of the meniscus height h(t), the azimuthally symmetric terms of Eq. (4) are averaged over the volume of the moving liquid column [11]. The time derivative of the velocity on the left-hand side of Eq. (4) is: The averaged pressure gradient term is obtained by taking into account that p (z = 0) = 0 and Supposing that the flow rate at the tube inlet v z | z−→0 −→ 0 and at the meniscus For the averaged viscosity term, With the mass balance equation, and considering the possibility of phase change where j is the interfacial mass flux due to the condensation process, given by [6].
In the above equation (12), k l and h lv are the thermal conductivity and the latent heat of vaporization of the liquid, respectively. Additionally, we need the liquid temperature field, T. Therefore, substituting Eqs. (7)-(10) into Eq. (4), we obtain the Lucas-Washburn equation that describes the rise of the liquid within the capillary [8] ρgh + 8μ with the following initial conditions where the surface tension is assumed to be dependent on temperature through the relationship: [12]. It must be noted that Eq. (15) is obtained by a balance between inertia and capillary forces, ignoring the influence of viscosity and gravity close to the contact moment between the capillary and the liquid reservoir. Quéré [13] deduced this expression, analytically and experimentally, establishing that for a short time the height of the liquid column increases linearly.
In order to solve the Eq. (13), additionally we need the energy equation of the liquid phase: with the following initial and boundary conditions In the above formulation we have assumed thermodynamic equilibrium condition at the interface, meaning that the phase change does not alter the interface temperature. In order to obtain the suitable characteristic scales for the problem, we conduct an order of magnitude analysis from the momentum and energy equations for the liquid penetrating the capillary tube. In this context, the competition between thermal and dynamics penetrations generates a non-isothermal capillary flow, which is developed inside the capillary tube. After an elapsed time t, the non-isothermal imbibition front reaches an average distance, h(t). Therefore, the thermal and imbibition effects introduce two time scales: the thermal penetration scale, t th , and the imbibition scale t i , which will be determined in the Order of magnitude analysis section.

Order of magnitude analysis and theoretical analysis
In this section we carry out an order of magnitude analysis, similar to that used by O. Bautista et al. [14] and J. P. Escandon et al. [15], in order to determine the characteristic scales of this problem. From Eq. (13) we can readily identify, by using an order of magnitude analysis, one characteristic time scale, t e , associated with the characteristic equilibrium height, h e . Similarly, from Eq. (16), we recognize two characteristic time scales associated with the thermal energy transported by the imbibed fluid: a convective scale t conv and a conductive scale t cond . For estimating the characteristic scales mentioned before, we conduct the following order of magnitude analysis. From Eq. (13), the order of magnitude for the characteristic equilibrium time t e of the imbibition process is determined by a balance between the surface tension and viscosity forces, as follows: where h e can be easily evaluated from a balance between the surface tension and viscosity forces, that is, Combining relationships (20) and (21), we obtain From this last relationship, the terms j 2 R/2σ 0 ρ v and 8μj/R 2 ρ 2 g are very small compared to the unity; in such case, the order of magnitude for the equilibrium height is h e ∼ 2σ 0 /ρgR , and the corresponding equilibrium time, from Eq. (20) is On the other hand, an energy balance between the transported thermal energy by the motion of the liquid and the accumulation energy term dictates that where, ΔT = T s − T 0 and u c , h th and t conv represent the characteristic average velocity associated with the velocity of the imbibition front, the characteristic thermal penetration and the characteristic convective time scale, respectively. Therefore, from Eq. (24), the characteristic convective time scale is given as, In a similar way, from Eq. (16), a balance between diffusive and accumulation terms, dictates that, and using the above relationship, the characteristic diffusive time scale, t cond , is given by To obtain the order of magnitude of the characteristic thermal penetration, we compare the convective and diffusive terms of the energy equation, obtaining that The characteristic average imbibition velocity, in a first approximation, is easily derived from Eq. (11), assuming that the mass flux j at the interface is very small, therefore By considering Eqs. (23), (25) and (27), we obtain The analysis of this transient heat transfer process can be characterized by adopting as characteristic time the imbibition time scale t e . The characteristic scales determined previously will be used to nondimensionalize the governing equations properly in the following section.

Model formulation and numerical solution
Introducing the following dimensionless variables the system of Eqs. (13)- (15) and (16)-(19), is transformed to where and β ∂θ ∂τ together with the following initial and boundary conditions, In the above equations, the dimensionless parameters ε, β and the Jakob number Ja are defined as and The solution of the problem (33)-(39) shall provide Y = Y (τ; ε, β, Ja, Γ) and θ = θ (η, τ; ε, β, Ja, Γ). To solve the system of Eqs. (33)-(39), we use a regular perturbation method [16], based on the fact that in direct condensation, the characteristic values of the Jakob number are usually small. This permits to consider the Jakob number, Ja, as the parameter of perturbation. We assume the following perturbation expansions, in power of Ja, for the dimensionless imbibition front and dimensionless temperature: and where Y 0 and θ 0 are the leading-order solutions, i. e., the case where no phase change occurs at the interface; Y 1 and θ 1 are higher order corrections to the leading order that include the phase change. Substituting the expansions (43) and (44) in the system of Eqs. (33)-(39), and collecting terms of the same powers of Ja, we obtain, for O Ja 0 : and The initial and boundary conditions associated with Eqs. (45) and (46) are, respectively For the first order solution, O Ja 1 , The initial and boundary conditions for Eqs. (52) and (53) are, respectively As a particular case and neglecting the inertial term, from Eq. (45), we can derive the well known case of the imbibition process without phase change, whose solution is given by, where, W (x) represents the Lambert function [17], being x the argument of this function, given by x = − exp (−1 − τ). Eq. (59) represents the zeroth-order solution for the momentum equation, without phase change and constant surface tension.
In order to obtain the zeroth-order solution for the dimensionless temperature in the liquid bulk, we apply the following finite difference scheme to the Eq. (46): where, Y 0 is obtained by integrating numerically the Eq. (47) with the classical fourth order Runge-Kutta method. In Eq. (60), i is the spatial node in the liquid bulk and n denotes the time step. Introducing the numerical solution of Eqs. (47) and (60) in Eq. (52), we obtain the next higher order equation for the imbibition front, Y 1 , which was solved by the conventional fourth order Runge-Kutta method. For simplicity, the details are omitted.

Results and discussion
In this chapter we conduct a semi-analytical approach to describe the imbibition process of a fluid into a capillary tube including the direct condensation at the imbibition front, considering that the surface tension at the liquid-vapor interface can be a function of the temperature. The present analysis serves to emphasize that the inertial terms in the momentum equation, the phase change at the imbibition front together with temperature-dependent surface tension effect, are very important for determining the temporal evolution of the imbibition front. This is reflected as a consequence of the main involved dimensionless parameters in the analysis: ε, the Jakob number, Ja and the parameter Γ. In all numerical calculations presented in this section, we have selected representative values for aforementioned parameters. Because the range of these parameters can be very large, the assumed values in our numerical results are representative of typical cases. The numerical results are plotted in Figs. 2-6. In these, we have emphasized the importance of the aforementioned dimensionless parameters on the imbibition process, and should be noted that the zeroth-order solution for the momentum equation is independent of the temperature field (see Eq. (46)). In Fig. 2 we evidence the importance of the parameter Γ on the imbibition front as a function of the dimensionless time. In this figure, the dimensionless height is plotted as a function of the dimensionless time for ε = 0.5, β = 600, Ja = 0, 0.1, 0.3. Here, it is very clear that the sensitivity of the imbibition front to the assumed values of the parameter Γ; for increasing values of above parameter, for a given dimensionless time, the dimensionless height of the imbibition front tends to decrease. The influence of the parameter ε, is shown in Fig. 4. Here we plotted the imbibition height as a function of the dimensionless time, for β = 600, three values of ε and two values of the parameter Γ. Black and hollows symbols represent the cases of Γ = 1 × 10 −6 and 1 × 10 −5 , respectively. As can be seen, the parameter ε results in the occurrence of oscillations, which appear for a value of ε = 0.5. For higher values of this parameter, greater imbibition heights are reached. As can be appreciated, for values of the parameter ε, it takes a larger dimensionless time for reaching the equilibrium height, given by the Jurin's law. Therefore, the first term (inertial term) at the right hand of Eq. (33), has a great influence on the hydrodynamic solution of the present problem. From the physical point of view, this behavior occurs for liquids of very low viscosity or for increasing values of the capillary radius (see relationship (40)) for a given liquid.
In Fig. 5, the dimensionless interfacial velocity of the liquid as a function of the dimensionless time is plotted, for β = 600, Ja = 0.1, ε = 0.1 and four vales of the parameter Γ. This figure  shows the case where no oscillations are present, and for the value of τ ∼ 3, the motion of the imbibition front is practically decelerated. From this figure, and according to the Fig. 2, increasing values of the parameter Γ tends to decrease the velocity of the imbibition front.
On the other hand, in Fig. 6, negative values for dY/dτ are obtained. In this figure, we show that for the case where oscillations are present, the equilibrium of the imbibition front occurs for larger dimensionless times in comparison with the case of ε = 0.  The effects of the Jakob number, the parameter ε and the dimensionless parameter Γ on the maximum dimensionless imbibition height are shown in Fig. 7. For the case of Ja = 0, the effect of the dimensionless parameter Γ has no influence on the imbibition front; however, for a fixed value of the Jakob number (Ja = 0), for increasing values of the parameter Γ, the maximum imbibition height decreases. This behavior is according to those shown in Figs. 2-4. From the same figure, Ja = 0 corresponds to that of isothermal case, i. e., there is no phase change process at the imbibition front.
The dimensionless temperature profiles in the liquid are shown in Fig. 8, as a function of the dimensionless coordinate η, for different values of the dimensionless time, two values of the parameter β(= 600, 1000), with ε = 0. In this figure it is evident that the dimensionless temperature profile in the liquid shows a linear behavior for large values of the dimensionless time, whereas for short dimensionless time, large temperature gradients exist near the imbibition front.