- Email: [email protected]
Dynamic Response of Geometrically Constrained Vapor Bubbles
Journal of Colloid and Interface Science 254, 346–354 (2002) doi:10.1006/jcis.2002.8631
Dynamic Response of Geometrically Constrained Vapor Bubbles Vladimir S. Ajaev,∗ G. M. Homsy,† and S. J. S. Morris‡ ∗ Department of Mathematics, Southern Methodist University, Dallas, Texas 75275; †Department of Mechanical and Environmental Engineering, University of California, Santa Barbara, California 93106; and ‡Department of Mechanical Engineering, University of California, Berkeley, California 94720 Received April 29, 2002; accepted July 22, 2002
We consider a two-dimensional model of a vapor bubble between two horizontal parallel boundaries held at different temperatures. When the temperatures are constant, a steady state can be achieved such that evaporation near the contact lines at the hot bottom plate is balanced by condensation in colder areas of the interface near the top. The dynamic response of the bubble is probed by treating the case of time-dependent wall temperatures. For periodic modulations of the wall temperature the bubble oscillates about the steady state. In order to describe such time-dependent behavior we consider the limit of small capillary number, in which the effects of heat and mass transfer are significant only near the contact lines at the bottom plate and in a small region near the top. When the bottom temperature is modulated and the top temperature is held fixed, the amplitude of forced oscillations is constant for low-frequency modulations and then rapidly decays in the high-frequency regime. When the top temperature is modulated with fixed bottom temperature, the dynamic-response curve is flat in the low-frequency regime as well, but it also flattens out when the frequency is increased. This shape of the response curve is shown to be the result of the nonmonotonic behavior of the thickness of the liquid film between the bubble interface and the top plate: when the temperature is decreased, the film thickness increases rapidly, but then slowly decays to a value which is smaller than the initial thickness. The dynamic response is also studied as a function of the forcing amplitude. C 2002 Elsevier Science (USA) Key Words: vapor bubbles; contact lines; disjoining pressure; microchannels.
1. INTRODUCTION AND PROBLEM DEFINITION
Geometrically constrained vapor and gas bubbles are used in various microelectromechnical systems (MEMS). In microactuators, mechanical parts are moved by expanding vapor bubbles created by localized heating (1). Periodically expanding and contracting vapor bubbles can also be used to pump liquid in microchannels (2). In these and other applications, transient dynamic response of the bubble to changes in temperature is of interest. An important feature of constrained bubbles created by localized heating is the presence of contact lines where the vapor– liquid interface meets the wall of a container. In order to determine the shape of the bubble one has to determine the position 0021-9797/02 $35.00
C 2002 Elsevier Science (USA)
All rights reserved.
of such contact lines and identify proper boundary conditions in their vicinity. Several local analyses of the steady version of the problem, developed by Potash and Wayner (3), Moosman and Homsy (4), and Morris (5), included the effects of heat/mass transfer, liquid flow, and disjoining pressure. Global steadystate solutions incorporating all these effects for vapor bubbles in rectangular channels were developed only recently for both two-dimensional (6) and three-dimensional (7) bubbles. The position of the contact line in (6, 7) was found by solving the freeboundary problems for the interface shape. It was shown that the contact line adjusts its position in such a way that the global mass balance is satisfied, so that evaporation near the contact line is balanced by condensation in colder areas of the interface. The main goal of this work is to generalize the approach of (6) to unsteady bubbles. Clearly, this requires a model for moving contact lines in the presence of evaporation, on which there is relatively little work. Anderson and Davis (8) considered this type of problem in the context of spreading. They assumed that corrections to the equilibrium value of the contact angle due to contact-line motion and evaporation are both relatively small and therefore their linear superposition can be used to calculate the contact angle. Contact lines with evaporation in the context of dry-out of thin films on heated surfaces were studied by Sharma (9), Oron and Bankoff (10), and Lyushnin et al. (11). These authors considered models of disjoining pressure and formulas for evaporative flux which are valid only for relatively thick films, as opposed to the adsorbed films used in (6) to describe the dry spots on heated surfaces. Wilson et al. (12) discuss the possible appearance of contact lines on hot surfaces in the context of unsteady growth of a long vapor bubble between two parallel plates. However, the focus of their mathematical model is on the case when both plates are covered with macroscopic films, so that there are no dry areas. Clearly, when dry areas appear, the motion of contact lines on hot surfaces is coupled with the deformation of the vapor– liquid interface in colder areas where condensation occurs. To our knowledge, there is no global time-dependent model that describes evolution of the bubble shape in this regime. In the present study we develop such a model. We use a framework similar to Wilson et al. (12) to model unsteady bubbles in the regime when dry areas and contact lines are present. We also use a two-dimensional model which while
346
347
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
T1 .
.
.
. .
.
vapor . .
. .
.
liquid 2. FORMULATION
.
. .
. .
.
T0 FIG. 1. peratures.
type functions of time. The dynamic response of the bubble to such forcing over a range of frequencies and forcing amplitudes is investigated.
A sketch of a bubble between two plates held at different tem-
simpler, gives insights into the much more complicated threedimensional geometries encountered in applications. This has been already shown for the steady case: the local structure of fluid flow and heat transfer near contact lines found from twodimensional models of (6) turned out to be relevant for the threedimensional case as well (7). We consider a two-dimensional model of a vapor bubble between two parallel plates held at different temperatures, see Fig. 1. The bottom plate is heated above the equilibrium saturation temperature of the liquid, TS∗ , the top plate is cooled below TS∗ . When the temperatures are constant, a steady state can be achieved such that evaporation near the contact lines at the heated bottom plate is balanced by condensation in colder areas of the interface near the top. This type of steady state can be realized experimentally when liquid in a microchannel is heated locally by a heater element of constant power. Heat conduction in the wall typically occurs on a faster time-scale than other processes and is therefore not considered in our model: we prescribe wall temperatures and assume that they can be changed instantaneously. In the present study we consider plate temperatures which are either constant or square-wave functions of time. However, by varying the frequency of the square-wave modulations over a wide range of values we gain insight into the response of the bubble to a forcing by an arbitrary time-dependent wall temperatures, which can be represented as a superposition of forcing functions at different frequencies. The paper is organized as follows. In Section 2 we discuss the geometry of a bubble between two plates held at different temperatures. The detailed formulation of the governing equations, as well as the boundary conditions, is presented. In Section 3 we outline the asymptotic method used to describe bubbles in the limit of small capillary number. The bubble shape is dominated by capillary forces everywhere except small regions near the top and bottom plates. We show that local solutions near the plates can be matched asymptotically to the capillary-statics shape. A steady-state solution is obtained such that evaporation near the contact line is balanced by condensation near the top. Finally, in Section 4 we discuss bubble oscillations about the steady state under conditions when the plate temperatures are square-wave-
Consider a two-dimensional model of a vapor bubble between two parallel plates, as shown in Fig. 1. The liquid phase has viscosity µ and density ρ; surface tension at the vapor–liquid interface is σ and is assumed to be constant. Thus, thermocapillary effects are neglected. In the nondimensional formulation we scale all lengths with d, the distance between plates. The effects of gravity are neglected since typical values of the Bond number, Bo = ρgd 2 /σ , are ∼10−4 (g is the acceleration of gravity). The equilibrium saturation temperature of the liquid, TS∗ , is chosen as the temperature scale. The energy balance at the interface gives the velocity scale (6) U=
kTS∗ , ρLd
where L is the latent heat of vaporization and k is the thermal conductivity of the liquid. In our nondimensional formulation we use d/U as the time scale which can be interpreted as the characterestic time required to evaporate a liquid film of uniform thickness d. Its typical value is on the order of 10 ms. There are also other time scales in the problem: the inverse frequency of modulations, which is typically much larger than d/U , and the time scales of heat conduction and viscous processes in the liquid, which are typically much smaller. Following our previous work, we use a capillary scale, σ/d, for the pressure. The quasistatic governing equations are −∇ p + C∇ 2 v = 0,
[1]
∇ · v = 0,
[2]
∇ 2 T = 0.
[3]
Here p and v are the nondimensional pressure and velocity in the liquid phase; C = µU/σ is the capillary number, which is typically small (10−3 –10−4 ). Both inertia and convective heat transfer are considered negligible; i.e., the Reynolds and Peclet numbers are small. Even if these numbers are small, the unsteady terms in the equations, vt and Tt , may need to be included when the frequency of modulations becomes large so that d/U is no longer the proper time scale. This is not the case in the regime that we are investigating. Contributions from the unsteady terms may become significant only when the Stokes numbers (St = ω∗ d 2 ρ/µ) are O(1), which corresponds to dimensional modulation frequencies ω∗ on the order of 1 MHz. Similar governing equations could be used to describe vapor phase. However, following our previous work, we adopt a onesided model of (13) which neglects all dynamical processes in the vapor. We also assume that the vapor pressure in the bubble
348
AJAEV, HOMSY, AND MORRIS
is constant. While inappropriate for rapidly expanding bubbles, this is generally a reasonable assumption for constrained bubbles that oscillate around a steady state. We define the nondimensional mass flux J in terms of its dimensional value J ∗ according to J = J ∗ /(ρU ). With this choice the conditions of conservation of mass and energy at the vapor–liquid interface take the nondimensional forms J = (v − v I ) · n,
[4]
J = −n · ∇T.
[5]
Here n is the unit normal vector to the interface pointing into the liquid and v I is the interfacial velocity. We write the scaled normal and shear stress balances as n · T · n = ∇ · n + − pv ,
[6]
t · T · n = 0.
[7]
Here pv is the scaled vapor pressure, and the condition [7] is satisfied for any tangent vector t. The nondimensional stress tensor, T, and disjoining pressure, , are the ratios of corresponding dimensional quantities to the capillary pressure scale, σ/d. The disjoining pressure term, , arises from unbalanced molecular interactions. It becomes significant only near the contact lines where it is inversely proportional to the cube of the liquid film thickness. Marangoni stresses are neglected in this analysis since for temperature differences of ∼10 K and length scales of a few micrometers the Marangoni number is typically small (∼0.1). Contributions from vapor recoil are also neglected. The temperature at the interface is different from the value given by equilibrium thermodynamics at the local vapor pressure since there is mass flux and a pressure jump at the interface. A modified relation, which we refer to as the nonequilibrium condition, was derived by Schrage (14). It is sometimes referred to as the Hertz–Knudsen–Langmuir equation. As shown in (6), the linearized version of this condition can be written in the following nondimensional form: K J = δ( p − pv ) + T − 1, ρU ¯ s∗ δ = σ . K = 2π RT 2ρv L Lρd
interface is at saturation temperature. This effect is sometimes called capillary condensation. Typical values of both δ and K are small (10−4 –10−2 ). Even though the nonequilibrium condition in our formulation is linearized, the response of the bubble to modulations of wall temperature can still be nonlinear due to conditions [6] and [7]. Let us now discuss the boundary conditions at the walls of the channel. Clearly all components of the velocity are zero there. The values of temperature at the walls are prescribed for the purpose of modeling the flow. In general, both the nondimensional temperature of the bottom, T0 , and the nondimensional temperature of the top, T1 , can be functions of time. In the simulations below we consider square-wave modulations of each of the temperatures when the other one is held fixed. The nondimensional frequency of modulations, scaled by U/d, is denoted by ω. 3. SMALL CAPILLARY NUMBER LIMIT
3.1. Scalings for Local Solution near Top For small values of the capillary number the curvature of the vapor–liquid interface is constant away from the plates since the interface is dominated by capillary statics there. There are regions near the contact lines and near the top plate where the interface shape is determined by the balance of viscous and capillary terms. In this section we discuss a set of scales for variables and distinguished limits for the parameters which allow one to obtain equations for these local interface shapes. The local solutions can then be matched to the capillary-statics solution away from the walls. Let us consider the local solution near the top and choose Cartesian coordinates such that the x-axis is horizontal and the y-axis is directed downwards, as shown in Fig. 2; the liquid film thickness is denoted by h. The condition of balance of viscous and capillary forces in the governing equations together with the matching condition for the curvature away from the top wall suggest that the variables have to be rescaled according to ˆ ˆ (y, h) = C 1/3 ( yˆ , h), x = C 1/6 x, t = C 2/3 tˆ.
i
[8]
Here R¯ is the gas constant per unit mass. According to [8], the departures of local temperature at the interface from the equilibrium value are characterized by two nondimensional parameters, K and δ. The kinetic parameter, K , measures the relative importance of kinetic effects, which lead to a shift in local vapor pressure from its equilibrium value; this shift changes the interfacial temperature. The value of K is found from kinetic theory by analyzing fluxes of vapor molecules towards the interface and away from it. The parameter δ shows how changes of the pressure in the liquid relative to the vapor affect the local phase-change temperature. We note that a pressure jump at the interface, typically due to capillary forces, can result in mass flux even if the
A balance of the viscous and capillary terms in the equation of motion, together with the continuity equation, determine the
x h
y FIG. 2. A sketch of the interface shape near the top. Local Cartesian coordinates are shown.
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
scales for the velocity components as
The conditions at the top wall are written in the form
ˆ v = C −1/3 v. u = C −1/2 u, ˆ
[9]
This type of balance appears in the literature in several different contexts, as discussed in (7). The governing equations at the leading order in the rescaled variables take the standard lubrication form: − pxˆ + uˆ yˆ yˆ = 0,
[10]
p yˆ = 0,
[11]
u xˆ + vˆ yˆ = 0,
[12]
Tyˆ yˆ = 0.
[13]
Let us now turn to the interfacial boundary conditions. In order to include the vapor mass-flux into the leading-order massconservation condition one has to rescale the flux according to J =C
−1/3
Jˆ, Jˆ = O(1).
With this choice and the above length and velocity scales the ˆ x) ˆ becomes kinematic condition, Eq. [4], at yˆ = h( hˆ tˆ + Jˆ + uˆ hˆ xˆ − vˆ = 0.
[14]
The rescaled condition for conservation of energy at the vapor– liquid interface, Eq. [5], at leading order becomes Jˆ = −Tyˆ .
[15]
Most of the terms in the normal and shear stress balances, Eqs. [6] and [7], drop out of the leading order in the transition region, so that the interfacial stress conditions take the relatively simple form p − pv = −hˆ xˆ xˆ −
349
ε , ˆh 3
uˆ yˆ = 0.
[16] [17]
Here ε = A/(σ d 2 C), A is the Hamaker constant. The value of ε is typically very small (10−4 –10−3 ), so the disjoining pressure is significant only for very thin films. In the region near the top the liquid film is always sufficiently thick and thus the disjoiningpressure term is negligible, so we drop it in our calculations in the next section. We rescale the kinetic constant, K , according to K = C 1/3 Kˆ , Kˆ = O(1). Then the nonequilibrium condition at the interface is written in the form ε Kˆ Jˆ = −δ hˆ xˆ xˆ + + T i − 1. [18] hˆ 3
u = v = 0, y = 0, T = T1 (tˆ), y = 0. Here T1 (tˆ) is the temperature of the top plate. 3.2. Evolution Equation near Top Let us now use the rescaled system of governing equations and boundary conditions to derive an evolution equation for film thickness near the top plate. First we note that according to [11] the pressure p is a function of xˆ only, so the momentum equation [10] can be integrated twice to give the standard lubrication-type velocity profile as described in (6). We substitute this velocity profile into the mass-conservation condition at the interface and neglect the disjoining pressure to obtain 3( Jˆ + hˆ tˆ ) = (hˆ 3 pxˆ )xˆ .
[19]
The mass flux Jˆ is then expressed in terms of the wall temperature and interface shape as described in (6). Substituting the formula for the flux, together with the expression for liquid pressure from the normal-stress balance, Eq. [16], into Eq. [19], results in the differential equation for the film thickness ˆ xˆ , tˆ): h( 1 T1 (tˆ) − 1 − δ hˆ xˆ xˆ hˆ tˆ + [hˆ 3 hˆ xˆ xˆ xˆ ]xˆ + = 0. 3 Kˆ + hˆ
[20]
This evolution equation has a clear physical interperation: the film thickness changes in time due to liquid flow driven by capillary pressure gradient, which is the second term on the left-hand side, and the condensation flux at the interface. The latter has two components: thermal and capillary condensation. The matching condition with the capillary-statics region is always written in the form hˆ xˆ xˆ → κ0 as xˆ → ∞,
[21]
where κ0 is the scaled mean curvature of the interface in the capillary-statics region. This condition together with the symmetry and no-flux conditions at xˆ = 0 allows one to solve the Eq. [20] for any time-dependence of the wall temperature, T1 (tˆ), and required value of the total mass condensing near the top, which has to be specified in order to obtain a unique solution. The choice of this value is discussed below in Section 3.4. 3.3. Evolution Equation near Contact Line Let us now discuss the region near the contact line, sketched in Fig. 3. An asymptotic method for describing steady contact lines has been developed in (6). It is relatively easy to show that the unsteady version of this method can be developed using the same approach as for the film near the top. In order to do this
350
AJAEV, HOMSY, AND MORRIS
3.4. Global Solutions evaporation
liquid flow adsorbed film
contact−line region
FIG. 3. A sketch of the region near the contact line and local coordinates showing the direction of liquid flow and evaporative mass flux. Dry areas are modeled by microscopic adsorbed film.
˜ y˜ ), shown in Fig. 3, as well as the film the local coordinates (x, ˜ have to be rescaled according to thickness h, ¯ ˜ = C 1/3 ( y¯ , h). ¯ ( y˜ , h) x˜ = C 1/6 x, Scales for the velocity components are chosen from the balance of viscous and capillary terms similar to those in Section 3.1. These scales lead to the standard lubrication-type velocity profile which allows one to reduce the problem to an evolution equation for local thickness in terms of the local coordinate x¯ and the time-dependent wall temperature T0 (tˆ), written in the form 1 T0 (tˆ) − 1 − δ h¯ x¯ x¯ − δε h¯ −3 h¯ tˆ + [h¯ 3 (h¯ x¯ x¯ + ε h¯ −3 )x¯ ]x¯ + = 0. 3 Kˆ + h¯ [22] One has to solve this equation with the matching condition ¯ for curvature at large values of x, h¯ x¯ x¯ → κ0 as x¯ → ∞,
[23]
and three conditions in the “dry” area. The latter in our approach is modeled by the isothermal adsorbed film (6). Thus h¯ decays to the local adsorbed film thickness that corresponds to T0 (tˆ). It is important to note that the steady-state solutions of Eq. [22] with the above boundary conditions are invariant under translations: there is no preference as to where the contact line should be when the plate is heated uniformly. To avoid resulting numerical difficulties without changing the physical problem we fix the interfacial position far away from the contact line where heat fluxes are very small. The physics of the problem near the contact line is illustrated in Fig. 3. Dry areas are modeled by microscopic adsorbed films and the interface rapidly changes its curvature in the region near the contact line in order to match the adsorbed film. Evaporation is suppressed in the adsorbed film, but it is strong in the region of macroscopic liquid film near the contact line.
The global solution of the capillary-statics problem for an unconstrained bubble in two dimensions is a circular interface. The actual bubble shape shown in Fig. 1 is the result of deformation of the circle in the regions near the top and bottom plates, i.e., the areas where the interfacial heat and mass transfer becomes significant. At the leading order in the capillary number the curvature of the deformed shape outside the transition regions, and thus the length of the bubble, remain unchanged. In our scaled coordinates the value of the curvature is κ0 = 2. After the value of κ0 is determined, the global problem can be solved numerically. It consists of Eq. [20] with the matching condition [21] and two conditions at the symmetry point, Eq. [22] with condition [23] and three conditions in the adsorbed film. They are coupled together by the condition of the instantaneous mass balance: we assume that processes in the vapor occur sufficiently fast so that the total amount of evaporation from contact lines at any given moment in time is equal to the absolute value of the condensation mass flux near the top. We note that at the leading order in the capillary number the vapor pressure and volume remain constant. The solution procedure is as follows. The evolution equation [22] from the previous section is solved numerically for any given bottom temperature, T0 (tˆ). The total amount of evaporation, J¯ , is found from the solution. According to the condition of instantaneous mass balance, the value of J¯ is used as the boundary condition for the local solution near the top. Equation [20] now has four boundary conditions, so it is solved numerically to obtain the shape of the liquid film near the top. We use finite-difference methods to solve the equations [20] and [22]. The length of the domain of integration is typically no larger than L = 5; for larger values of the spatial coordinate the interface can be approximated by a parabolic shape without decreasing the accuracy of the solution; the (small) correction to mass flux from that region is found analytically. In order to evaluate the nonlinear term h 3 on our finite-difference mesh we take the average of its values at two neighboring points, as described in (15). Evaluating the nonlinear terms and derivatives at the endpoints involves mesh points outside the physical domain, where the function is extrapolated using the boundary conditions at the endpoints. The total number of mesh points is varied between 102 and 103 to verify convergence of the solutions. After the spatial domain is discretized, the partial differential equation is reduced to a system of first-order ordinary differential equations in time. DVODE solver is used to solve this system (16). 4. RESULTS
4.1. Steady States Let us choose representative value of parameters Kˆ = 2.0, δ = 0.04, ε = 10−3 for all simulations in this and the subsequent two subsections. These correspond to typical dimensional parameters given in (7) corresponding to organic liquids in microchannels with modest superheating above room temperature.
351
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
a
of the inverse thickness). There is a dimple in the thickness in Fig. 4b, which is indicative of the significant capillary pressure gradient that drives the flow of the condensing liquid away from the middle of the bubble.
5.0
4.0
_ h
4.2. Modulation of the Bottom Temperature
3.0
2.0
1.0
0.0 0.0
b
1.0
0.5
2.0
1.5
_ x
4.0
a
3.0
h^
When the temperatures of the plates become time-dependent, both the position of the contact line and the shape of the interface near the top change. The description of contact-line motion is relatively simple: it slides on the bottom plate in response to evaporation while the basic flow structure, shown in Fig. 3, does not change significantly. Such motion is illustrated in Fig. 5a, where the local interface shape near the contact line is plotted for different moments in time after a rapid change in the bottom plate temperature from T0 = 1.5 to T0 = 1.4. The initial position of the interface is shown by the dotted line. We also plot the ratio of the disjoining pressure to the capillary pressure, r , as a function 0.5
2.0
^
h 0.25
1.0
0.0 0.0
1.0
2.0
3.0
x^ 0.2
FIG. 4. Steady-state interface shape for the regions near the contact line (a) and near the top of the bubble (b) for Kˆ = 2.0, δ = 0.04, ε = 10−3 , T0 = 1.5, T1 = 0.8.
First let us consider the case when wall temperatures are constant, taking T0 = 1.5, T1 = 0.8. The time evolution from a constant-curvature initial condition indicates that the solutions approach steady-state shapes for which condensation and evaporation balance each other. Local interface shapes in this steadystate regime are plotted in Fig. 4. Several features of the steady solutions were discussed in (7), but it is useful to recount them here. Figure 4a illustrates the solution near the contact line at the bottom. Clearly, when x¯ becomes large, the interface curvature becomes constant, while as x¯ goes to negative infinity, the film thickness approaches the small constant value of the isothermal adsorbed film that mimics a dry region. In the region near the top there are no dry areas: the steady film is relatively thick, as seen in Fig. 4b. We note that hˆ and h¯ in Fig. 4 are the values of film thickness rescaled by the same number, C 1/3 . Thus there is an order of magnitude difference between the values of thickness at the top and at the bottom, which shows that near the top the disjoining pressure is negligible (it being proportional to the cube
0.4
0.6
^
x b
4.0
3.0
r 2.0
1.0 0.0
1.0
0.5
1.5
^
t
FIG. 5. The motion of the contact line on the bottom plate after instantaneous change of bottom temperature T0 from 1.5 to 1.4. (a) Interface shapes are shown at tˆ = 0.3 (dashed line) and tˆ = 0.6 (solid line). The initial position of the interface is also shown (dotted line). (b) The ratio of the disjoining pressure to the capillary pressure is shown as a function of time at xˆ = 0.2.
352
AJAEV, HOMSY, AND MORRIS
of time at a fixed location near the contact line (xˆ = 0.2) in Fig. 5b. The plot shows that the relative value of the disjoiningpressure term decreases rapidly as the film thickness increases. After tˆ ≈ 1 there is essentially no change in the position of the contact line and the pressure. Thus the motion of the contact line after a change in temperature can be described as rapid monotonic approach to a new steady-state position. Changes in the interfacial shape near the top occur on a much slower time scale and are more complicated to describe, so we focus on them in our discussion of unsteady simulations in this and the following subsections. In this section we consider modulations of the bottom temperature when the top temperature is held fixed at T1 = 0.8. This implies that the bottom plate temperature T0 (tˆ) in the evolution equation [22] is now a function of time. We choose it in the form of a square wave of dimensionless amplitude and frequency ω. Since the rescaled time tˆ is the time variable, it is convenient to also introduce rescaled frequency of forcing ωˆ according to ωˆ = C
2/3
a 0.7
0.6
h0 0.5
0.4
0
1000
500
^
t
b
0.8
0.7
ω. 0.6
The square-wave forcing can then be written in the form T0 (tˆ) =
T0 (0) − , n ωˆ −1 ≤ tˆ < (n + 1/2)ωˆ −1 T0 (0),
(n + 1/2)ωˆ −1 ≤ tˆ < (n + 1)ωˆ −1 ;
h0 0.5
0.4
n = 0, 1, 2, . . . . 0.3
Let us choose T0 (0) = 1.5 and = 0.1 and study the resulting changes in the liquid film shape near the top. We note that the bubble is always symmetric with respect to a vertical line, as seen in Fig. 1, and therefore the liquid film is also symmetric with respect to the y-axis shown in Fig. 2. Although our results give the entire film profile as a function of time, it is convenient to characterize the shape of the bubble by the film thickness at ˆ the point of symmetry, h 0 = h(0). We plot h 0 as a function of time in Fig. 6 for two different values of the scaled frequency of modulations: ωˆ = 0.01 and ωˆ = 0.002. Initially, there is a transient response, as clearly seen in Fig. 6a, but then the response settles into a periodic solution with period equal to that of the forcing. The amplitude of forced oscillations of the film thickness is smaller for larger values of frequency; the system does not have time to fully relax to the steady state before the next change in forcing. This regime is termed the high-frequency regime. A different physical picture is obtained for low-frequency modulations, when the interface almost relaxes to the steady state before the next change in temperature. The response amplitude then remains constant when the frequency is decreased further. The response in the low-frequency regime is illustrated in Fig. 6b. The two regimes are clearly seen in Fig. 7a, where the normalized amplitude is plotted as a function of modulation frequency for different values of forcing. The normalized amplitude a is defined as the ratio of the amplitude of variations of h 0 to the nondimensional forcing amplitude, . Since we vary the fre-
0
1000
500
1500
^
t
FIG. 6. Liquid film thickness at the symmetry point near the top as a function of time for square-wave modulations of the bottom temperature for different values of modulation frequency: (a): ωˆ = 0.01; (b): ωˆ = 0.002.
quency over a wide range of values, it is convenient to use a linear–logarithmic plot. The transition between the two regimes occurs at the values of the scaled frequency on the order of 10−3 . The corresponding dimensional value is approximately 100 Hz. The high-frequency regime correspond to the typical dimensional values of several kHz. We note that this is well below the value of 1 MHz at which the quasistatic approximation in the governing equations can break down as discussed in Section 2. It is clear from Fig. 7a that the response is strongly nonlinear—the normalized amplitude is different for different values of forcing. In fact, Fig. 7a shows that the value of normalized amplitude can be decreasing when is increased: nonlinear effects suppress oscillations. In order to better describe the response we also plot it in loglog coordinates as shown in Fig. 7b. Clearly there is a range of frequencies between 0.005 and 0.05 where a power-law decay of the amplitude is observed. However, as frequency is increased to 0.1, the rapid processes in the contact line-region start to affect the response, so the simple power law is no longer applicable.
353
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
a
6.0
5.0
4.0
a
3.0
2.0
1.0
0.0
-3
10
b
-2
^
ω
10
-1
10
1
10
pressure gradient needed for steady state, which actually means a smaller thickness. This interpretation is also supported by the fact that initial increase occurs on a time scale of the interface motion due to phase change (t¯ ≈ 1), while the slow decay is on the scale of viscous and conductive processes in the liquid. Thus in the short run the thickness is increased by decreasing the temperature, but over a longer period, the lower temperature leads to a thinner film. The nonmonotonic behavior of the interface leads to a different dynamic response curve. This is illustrated in Fig. 9, where the normalized amplitude of oscillations is plotted as a function of modulation frequency for different values of forcing . As before, the normalized amplitude a is the amplitude of variation of h 0 divided by . We observe again the high-frequency and low-frequency regimes of oscillations: transition between them occurs for ωˆ ≈ 10−3 . Contrary to the rapid power-law decay of amplitude for the bottom-temperature modulation, the response curve here flattens out for values of frequency near 0.05. This is mostly due to the presence of the initial rapid jump in interface
a
a
0.81
0
10
0.8
0.79
0.78
T1 0.77
-1
10
-3
10
-2
^ ω
10
-1
10
0.76
FIG. 7. The dynamic response over a range of frequencies: normalized amplitude of forced oscillations of film thickness at the point of symmetry is plotted vs frequency using linear–log (a) and log–log (b) scales for two different values of forcing amplitude: = 0.1 (solid lines) and = 0.05 (dashed lines).
4.3. Modulation of the Top Temperature Let us now study the effect of modulations of the top temperature on the bubble shape. We take a fixed value of the bottom temperature, T0 = 1.5, and then apply square-wave modulation of the top temperature of amplitude = 0.05, which is shown in Fig. 8a. The measure of film thickness, which is the value at the point of symmetry, h 0 , is plotted as a function of t¯ in Fig. 8b. We are not interested in transient effects, so all data are recorded starting at the zero value of the shifted time-variable t¯ = tˆ − 10ωˆ −1 . Typical results are shown in Fig. 8b. Clearly, there is a nonmonotonic behavior. When the temperature T1 decreases, the thickness initially increases rapidly, then slowly decreases. This result can be interpreted as follows. The initial rapid decrease of the temperature causes a rapid increase in condensation; however, the pressure gradient is not sufficient to remove the condensing liquid. The interface then adjusts its shape to provide the
0.75
0.74
b
0
200
400
_ t
600
800
0
200
400
600
800
0.7
0.6
h0 0.5
0.4
_ t
FIG. 8. The top wall temperature as a function of time for square-wave modulations of frequency ωˆ = 0.005 (a) and corresponding liquid film thickness at the point of symmetry (b) for K = 2.0, δ = 0.04, T0 = 1.5.
354
AJAEV, HOMSY, AND MORRIS
12.0
10.0
8.0
a
6.0
4.0
2.0
0.0
-3
10
-2
^
ω
10
-1
10
FIG. 9. Dynamic response to the modulations of the top temperature: normalized amplitude of oscillations of film thickness near the top is plotted as a function of scaled frequency for different values of forcing amplitude: = 0.25 (solid line), = 0.5 (dashed line), and = 0.75 (dot–dashed line).
position, which is essentially independent of the modulation frequency. Thus the amplitude decreases to some finite value of this initial jump as the frequency is increased. The nonlinear character of the response is also clearly seen: the amplitude a is increased when the forcing amplitude increases. This shows that nonlinear effects lead to larger amplitude of oscillations than may be expected from a linearized response model.
5. CONCLUSIONS
We develop a mathematical model for unsteady motion of a bubble between two parallel plates subject to time-dependent temperature conditions. Global time-dependent interface shapes are found in the limit of small capillary number. The interface is dominated by capillary forces everywhere except near the top plate and the contact-line regions at the bottom plate. The dynamic response of the bubble to square-wave modulations of plate temperatures is studied. When the bottom temperature is modulated and the top temperature is held fixed, the amplitude of forced oscillations is approximately constant
when the scaled frequency is below ωˆ ≈ 10−3 and then shows power-law decay for larger frequencies. When the top temperature is modulated with fixed bottom temperature, the dynamicresponse curve is flat in the low-frequency regime as well, but it also flattens out when the frequency is increased. This shape of the response curve is a result of the nonmonotonic behavior of the liquid film between the bubble interface and the top plate. When the top-plate temperature is decreased, the film thickness increases rapidly, but then slowly decays to a value which is smaller than the initial thickness. We have also studied the dynamic reponse as a function of forcing amplitude and found it to be strongly nonlinear even for very small values of the amplitudes of temperature variations. Nonlinear effects can result in oscillation amplitudes which are much larger than can be expected from the linear analysis of the response. ACKNOWLEDGMENTS GMH is grateful to the U.S. Department of Energy, Office of Basic Energy Sciences, for partial funding of this work. VA thanks Professor Ian Gladwell for useful advice on numerical computations.
REFERENCES 1. Lin, L., Pisano, A. P., and Carey, V. P., J. Heat Transfer 120, 735 (1998). 2. Ory, E., Yuan, H., Prosperetti, A., Popinet, S., and Zaleski, S., Phys. Fluids 12, 1268 (2000). 3. Potash, M., and Wayner, P. C., Int. J. Heat and Mass Transfer 15, 1851 (1972). 4. Moosman, S., and Homsy, G. M., J. Colloid Interface Sci. 73, 212 (1980). 5. Morris, S. J. S., J. Fluid Mech. 432, 1 (2001). 6. Ajaev, V. S., and Homsy, G. M., J. Colloid Interface Sci. 240, 259 (2001). 7. Ajaev, V. S., and Homsy, G. M., J. Colloid Interface Sci. 244, 180 (2001). 8. Anderson, D. M., and Davis, S. H., Phys. Fluids 7, 248 (1995). 9. Sharma, A., Langmuir 14, 4915 (1998). 10. Oron, A., and Bankoff, S. G., J. Colloid Interface Sci. 218, 152 (1999). 11. Lyushnin, A. V., Golovin, A. A., and Pismen, L. M., Phys. Rev. E 65, 021602 (2002). 12. Wilson, S. K., Davis, S. H., and Bankoff, S. G., J. Fluid Mech. 391, 1 (1999). 13. Burelbach, J. P., Bankoff, S. G., and Davis, S. H., J. Fluid Mech. 195, 463 (1988). 14. Schrage, R. W., “A Theoretical Study of Interface Mass Transfer,” Columbia Univ. Press, New York, 1953. 15. Zhornitskaya, L., and Bertozzi, A. L., SIAM J. Numer. Anal. 37, 523 (2000). 16. Brown, P. N., Byrne, C. D., and Hindmarsh, A. C., SIAM J. Sci. Stat. Comput. 10, 1038 (1989).
Dynamic Response of Geometrically Constrained Vapor Bubbles Vladimir S. Ajaev,∗ G. M. Homsy,† and S. J. S. Morris‡ ∗ Department of Mathematics, Southern Methodist University, Dallas, Texas 75275; †Department of Mechanical and Environmental Engineering, University of California, Santa Barbara, California 93106; and ‡Department of Mechanical Engineering, University of California, Berkeley, California 94720 Received April 29, 2002; accepted July 22, 2002
We consider a two-dimensional model of a vapor bubble between two horizontal parallel boundaries held at different temperatures. When the temperatures are constant, a steady state can be achieved such that evaporation near the contact lines at the hot bottom plate is balanced by condensation in colder areas of the interface near the top. The dynamic response of the bubble is probed by treating the case of time-dependent wall temperatures. For periodic modulations of the wall temperature the bubble oscillates about the steady state. In order to describe such time-dependent behavior we consider the limit of small capillary number, in which the effects of heat and mass transfer are significant only near the contact lines at the bottom plate and in a small region near the top. When the bottom temperature is modulated and the top temperature is held fixed, the amplitude of forced oscillations is constant for low-frequency modulations and then rapidly decays in the high-frequency regime. When the top temperature is modulated with fixed bottom temperature, the dynamic-response curve is flat in the low-frequency regime as well, but it also flattens out when the frequency is increased. This shape of the response curve is shown to be the result of the nonmonotonic behavior of the thickness of the liquid film between the bubble interface and the top plate: when the temperature is decreased, the film thickness increases rapidly, but then slowly decays to a value which is smaller than the initial thickness. The dynamic response is also studied as a function of the forcing amplitude. C 2002 Elsevier Science (USA) Key Words: vapor bubbles; contact lines; disjoining pressure; microchannels.
1. INTRODUCTION AND PROBLEM DEFINITION
Geometrically constrained vapor and gas bubbles are used in various microelectromechnical systems (MEMS). In microactuators, mechanical parts are moved by expanding vapor bubbles created by localized heating (1). Periodically expanding and contracting vapor bubbles can also be used to pump liquid in microchannels (2). In these and other applications, transient dynamic response of the bubble to changes in temperature is of interest. An important feature of constrained bubbles created by localized heating is the presence of contact lines where the vapor– liquid interface meets the wall of a container. In order to determine the shape of the bubble one has to determine the position 0021-9797/02 $35.00
C 2002 Elsevier Science (USA)
All rights reserved.
of such contact lines and identify proper boundary conditions in their vicinity. Several local analyses of the steady version of the problem, developed by Potash and Wayner (3), Moosman and Homsy (4), and Morris (5), included the effects of heat/mass transfer, liquid flow, and disjoining pressure. Global steadystate solutions incorporating all these effects for vapor bubbles in rectangular channels were developed only recently for both two-dimensional (6) and three-dimensional (7) bubbles. The position of the contact line in (6, 7) was found by solving the freeboundary problems for the interface shape. It was shown that the contact line adjusts its position in such a way that the global mass balance is satisfied, so that evaporation near the contact line is balanced by condensation in colder areas of the interface. The main goal of this work is to generalize the approach of (6) to unsteady bubbles. Clearly, this requires a model for moving contact lines in the presence of evaporation, on which there is relatively little work. Anderson and Davis (8) considered this type of problem in the context of spreading. They assumed that corrections to the equilibrium value of the contact angle due to contact-line motion and evaporation are both relatively small and therefore their linear superposition can be used to calculate the contact angle. Contact lines with evaporation in the context of dry-out of thin films on heated surfaces were studied by Sharma (9), Oron and Bankoff (10), and Lyushnin et al. (11). These authors considered models of disjoining pressure and formulas for evaporative flux which are valid only for relatively thick films, as opposed to the adsorbed films used in (6) to describe the dry spots on heated surfaces. Wilson et al. (12) discuss the possible appearance of contact lines on hot surfaces in the context of unsteady growth of a long vapor bubble between two parallel plates. However, the focus of their mathematical model is on the case when both plates are covered with macroscopic films, so that there are no dry areas. Clearly, when dry areas appear, the motion of contact lines on hot surfaces is coupled with the deformation of the vapor– liquid interface in colder areas where condensation occurs. To our knowledge, there is no global time-dependent model that describes evolution of the bubble shape in this regime. In the present study we develop such a model. We use a framework similar to Wilson et al. (12) to model unsteady bubbles in the regime when dry areas and contact lines are present. We also use a two-dimensional model which while
346
347
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
T1 .
.
.
. .
.
vapor . .
. .
.
liquid 2. FORMULATION
.
. .
. .
.
T0 FIG. 1. peratures.
type functions of time. The dynamic response of the bubble to such forcing over a range of frequencies and forcing amplitudes is investigated.
A sketch of a bubble between two plates held at different tem-
simpler, gives insights into the much more complicated threedimensional geometries encountered in applications. This has been already shown for the steady case: the local structure of fluid flow and heat transfer near contact lines found from twodimensional models of (6) turned out to be relevant for the threedimensional case as well (7). We consider a two-dimensional model of a vapor bubble between two parallel plates held at different temperatures, see Fig. 1. The bottom plate is heated above the equilibrium saturation temperature of the liquid, TS∗ , the top plate is cooled below TS∗ . When the temperatures are constant, a steady state can be achieved such that evaporation near the contact lines at the heated bottom plate is balanced by condensation in colder areas of the interface near the top. This type of steady state can be realized experimentally when liquid in a microchannel is heated locally by a heater element of constant power. Heat conduction in the wall typically occurs on a faster time-scale than other processes and is therefore not considered in our model: we prescribe wall temperatures and assume that they can be changed instantaneously. In the present study we consider plate temperatures which are either constant or square-wave functions of time. However, by varying the frequency of the square-wave modulations over a wide range of values we gain insight into the response of the bubble to a forcing by an arbitrary time-dependent wall temperatures, which can be represented as a superposition of forcing functions at different frequencies. The paper is organized as follows. In Section 2 we discuss the geometry of a bubble between two plates held at different temperatures. The detailed formulation of the governing equations, as well as the boundary conditions, is presented. In Section 3 we outline the asymptotic method used to describe bubbles in the limit of small capillary number. The bubble shape is dominated by capillary forces everywhere except small regions near the top and bottom plates. We show that local solutions near the plates can be matched asymptotically to the capillary-statics shape. A steady-state solution is obtained such that evaporation near the contact line is balanced by condensation near the top. Finally, in Section 4 we discuss bubble oscillations about the steady state under conditions when the plate temperatures are square-wave-
Consider a two-dimensional model of a vapor bubble between two parallel plates, as shown in Fig. 1. The liquid phase has viscosity µ and density ρ; surface tension at the vapor–liquid interface is σ and is assumed to be constant. Thus, thermocapillary effects are neglected. In the nondimensional formulation we scale all lengths with d, the distance between plates. The effects of gravity are neglected since typical values of the Bond number, Bo = ρgd 2 /σ , are ∼10−4 (g is the acceleration of gravity). The equilibrium saturation temperature of the liquid, TS∗ , is chosen as the temperature scale. The energy balance at the interface gives the velocity scale (6) U=
kTS∗ , ρLd
where L is the latent heat of vaporization and k is the thermal conductivity of the liquid. In our nondimensional formulation we use d/U as the time scale which can be interpreted as the characterestic time required to evaporate a liquid film of uniform thickness d. Its typical value is on the order of 10 ms. There are also other time scales in the problem: the inverse frequency of modulations, which is typically much larger than d/U , and the time scales of heat conduction and viscous processes in the liquid, which are typically much smaller. Following our previous work, we use a capillary scale, σ/d, for the pressure. The quasistatic governing equations are −∇ p + C∇ 2 v = 0,
[1]
∇ · v = 0,
[2]
∇ 2 T = 0.
[3]
Here p and v are the nondimensional pressure and velocity in the liquid phase; C = µU/σ is the capillary number, which is typically small (10−3 –10−4 ). Both inertia and convective heat transfer are considered negligible; i.e., the Reynolds and Peclet numbers are small. Even if these numbers are small, the unsteady terms in the equations, vt and Tt , may need to be included when the frequency of modulations becomes large so that d/U is no longer the proper time scale. This is not the case in the regime that we are investigating. Contributions from the unsteady terms may become significant only when the Stokes numbers (St = ω∗ d 2 ρ/µ) are O(1), which corresponds to dimensional modulation frequencies ω∗ on the order of 1 MHz. Similar governing equations could be used to describe vapor phase. However, following our previous work, we adopt a onesided model of (13) which neglects all dynamical processes in the vapor. We also assume that the vapor pressure in the bubble
348
AJAEV, HOMSY, AND MORRIS
is constant. While inappropriate for rapidly expanding bubbles, this is generally a reasonable assumption for constrained bubbles that oscillate around a steady state. We define the nondimensional mass flux J in terms of its dimensional value J ∗ according to J = J ∗ /(ρU ). With this choice the conditions of conservation of mass and energy at the vapor–liquid interface take the nondimensional forms J = (v − v I ) · n,
[4]
J = −n · ∇T.
[5]
Here n is the unit normal vector to the interface pointing into the liquid and v I is the interfacial velocity. We write the scaled normal and shear stress balances as n · T · n = ∇ · n + − pv ,
[6]
t · T · n = 0.
[7]
Here pv is the scaled vapor pressure, and the condition [7] is satisfied for any tangent vector t. The nondimensional stress tensor, T, and disjoining pressure, , are the ratios of corresponding dimensional quantities to the capillary pressure scale, σ/d. The disjoining pressure term, , arises from unbalanced molecular interactions. It becomes significant only near the contact lines where it is inversely proportional to the cube of the liquid film thickness. Marangoni stresses are neglected in this analysis since for temperature differences of ∼10 K and length scales of a few micrometers the Marangoni number is typically small (∼0.1). Contributions from vapor recoil are also neglected. The temperature at the interface is different from the value given by equilibrium thermodynamics at the local vapor pressure since there is mass flux and a pressure jump at the interface. A modified relation, which we refer to as the nonequilibrium condition, was derived by Schrage (14). It is sometimes referred to as the Hertz–Knudsen–Langmuir equation. As shown in (6), the linearized version of this condition can be written in the following nondimensional form: K J = δ( p − pv ) + T − 1, ρU ¯ s∗ δ = σ . K = 2π RT 2ρv L Lρd
interface is at saturation temperature. This effect is sometimes called capillary condensation. Typical values of both δ and K are small (10−4 –10−2 ). Even though the nonequilibrium condition in our formulation is linearized, the response of the bubble to modulations of wall temperature can still be nonlinear due to conditions [6] and [7]. Let us now discuss the boundary conditions at the walls of the channel. Clearly all components of the velocity are zero there. The values of temperature at the walls are prescribed for the purpose of modeling the flow. In general, both the nondimensional temperature of the bottom, T0 , and the nondimensional temperature of the top, T1 , can be functions of time. In the simulations below we consider square-wave modulations of each of the temperatures when the other one is held fixed. The nondimensional frequency of modulations, scaled by U/d, is denoted by ω. 3. SMALL CAPILLARY NUMBER LIMIT
3.1. Scalings for Local Solution near Top For small values of the capillary number the curvature of the vapor–liquid interface is constant away from the plates since the interface is dominated by capillary statics there. There are regions near the contact lines and near the top plate where the interface shape is determined by the balance of viscous and capillary terms. In this section we discuss a set of scales for variables and distinguished limits for the parameters which allow one to obtain equations for these local interface shapes. The local solutions can then be matched to the capillary-statics solution away from the walls. Let us consider the local solution near the top and choose Cartesian coordinates such that the x-axis is horizontal and the y-axis is directed downwards, as shown in Fig. 2; the liquid film thickness is denoted by h. The condition of balance of viscous and capillary forces in the governing equations together with the matching condition for the curvature away from the top wall suggest that the variables have to be rescaled according to ˆ ˆ (y, h) = C 1/3 ( yˆ , h), x = C 1/6 x, t = C 2/3 tˆ.
i
[8]
Here R¯ is the gas constant per unit mass. According to [8], the departures of local temperature at the interface from the equilibrium value are characterized by two nondimensional parameters, K and δ. The kinetic parameter, K , measures the relative importance of kinetic effects, which lead to a shift in local vapor pressure from its equilibrium value; this shift changes the interfacial temperature. The value of K is found from kinetic theory by analyzing fluxes of vapor molecules towards the interface and away from it. The parameter δ shows how changes of the pressure in the liquid relative to the vapor affect the local phase-change temperature. We note that a pressure jump at the interface, typically due to capillary forces, can result in mass flux even if the
A balance of the viscous and capillary terms in the equation of motion, together with the continuity equation, determine the
x h
y FIG. 2. A sketch of the interface shape near the top. Local Cartesian coordinates are shown.
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
scales for the velocity components as
The conditions at the top wall are written in the form
ˆ v = C −1/3 v. u = C −1/2 u, ˆ
[9]
This type of balance appears in the literature in several different contexts, as discussed in (7). The governing equations at the leading order in the rescaled variables take the standard lubrication form: − pxˆ + uˆ yˆ yˆ = 0,
[10]
p yˆ = 0,
[11]
u xˆ + vˆ yˆ = 0,
[12]
Tyˆ yˆ = 0.
[13]
Let us now turn to the interfacial boundary conditions. In order to include the vapor mass-flux into the leading-order massconservation condition one has to rescale the flux according to J =C
−1/3
Jˆ, Jˆ = O(1).
With this choice and the above length and velocity scales the ˆ x) ˆ becomes kinematic condition, Eq. [4], at yˆ = h( hˆ tˆ + Jˆ + uˆ hˆ xˆ − vˆ = 0.
[14]
The rescaled condition for conservation of energy at the vapor– liquid interface, Eq. [5], at leading order becomes Jˆ = −Tyˆ .
[15]
Most of the terms in the normal and shear stress balances, Eqs. [6] and [7], drop out of the leading order in the transition region, so that the interfacial stress conditions take the relatively simple form p − pv = −hˆ xˆ xˆ −
349
ε , ˆh 3
uˆ yˆ = 0.
[16] [17]
Here ε = A/(σ d 2 C), A is the Hamaker constant. The value of ε is typically very small (10−4 –10−3 ), so the disjoining pressure is significant only for very thin films. In the region near the top the liquid film is always sufficiently thick and thus the disjoiningpressure term is negligible, so we drop it in our calculations in the next section. We rescale the kinetic constant, K , according to K = C 1/3 Kˆ , Kˆ = O(1). Then the nonequilibrium condition at the interface is written in the form ε Kˆ Jˆ = −δ hˆ xˆ xˆ + + T i − 1. [18] hˆ 3
u = v = 0, y = 0, T = T1 (tˆ), y = 0. Here T1 (tˆ) is the temperature of the top plate. 3.2. Evolution Equation near Top Let us now use the rescaled system of governing equations and boundary conditions to derive an evolution equation for film thickness near the top plate. First we note that according to [11] the pressure p is a function of xˆ only, so the momentum equation [10] can be integrated twice to give the standard lubrication-type velocity profile as described in (6). We substitute this velocity profile into the mass-conservation condition at the interface and neglect the disjoining pressure to obtain 3( Jˆ + hˆ tˆ ) = (hˆ 3 pxˆ )xˆ .
[19]
The mass flux Jˆ is then expressed in terms of the wall temperature and interface shape as described in (6). Substituting the formula for the flux, together with the expression for liquid pressure from the normal-stress balance, Eq. [16], into Eq. [19], results in the differential equation for the film thickness ˆ xˆ , tˆ): h( 1 T1 (tˆ) − 1 − δ hˆ xˆ xˆ hˆ tˆ + [hˆ 3 hˆ xˆ xˆ xˆ ]xˆ + = 0. 3 Kˆ + hˆ
[20]
This evolution equation has a clear physical interperation: the film thickness changes in time due to liquid flow driven by capillary pressure gradient, which is the second term on the left-hand side, and the condensation flux at the interface. The latter has two components: thermal and capillary condensation. The matching condition with the capillary-statics region is always written in the form hˆ xˆ xˆ → κ0 as xˆ → ∞,
[21]
where κ0 is the scaled mean curvature of the interface in the capillary-statics region. This condition together with the symmetry and no-flux conditions at xˆ = 0 allows one to solve the Eq. [20] for any time-dependence of the wall temperature, T1 (tˆ), and required value of the total mass condensing near the top, which has to be specified in order to obtain a unique solution. The choice of this value is discussed below in Section 3.4. 3.3. Evolution Equation near Contact Line Let us now discuss the region near the contact line, sketched in Fig. 3. An asymptotic method for describing steady contact lines has been developed in (6). It is relatively easy to show that the unsteady version of this method can be developed using the same approach as for the film near the top. In order to do this
350
AJAEV, HOMSY, AND MORRIS
3.4. Global Solutions evaporation
liquid flow adsorbed film
contact−line region
FIG. 3. A sketch of the region near the contact line and local coordinates showing the direction of liquid flow and evaporative mass flux. Dry areas are modeled by microscopic adsorbed film.
˜ y˜ ), shown in Fig. 3, as well as the film the local coordinates (x, ˜ have to be rescaled according to thickness h, ¯ ˜ = C 1/3 ( y¯ , h). ¯ ( y˜ , h) x˜ = C 1/6 x, Scales for the velocity components are chosen from the balance of viscous and capillary terms similar to those in Section 3.1. These scales lead to the standard lubrication-type velocity profile which allows one to reduce the problem to an evolution equation for local thickness in terms of the local coordinate x¯ and the time-dependent wall temperature T0 (tˆ), written in the form 1 T0 (tˆ) − 1 − δ h¯ x¯ x¯ − δε h¯ −3 h¯ tˆ + [h¯ 3 (h¯ x¯ x¯ + ε h¯ −3 )x¯ ]x¯ + = 0. 3 Kˆ + h¯ [22] One has to solve this equation with the matching condition ¯ for curvature at large values of x, h¯ x¯ x¯ → κ0 as x¯ → ∞,
[23]
and three conditions in the “dry” area. The latter in our approach is modeled by the isothermal adsorbed film (6). Thus h¯ decays to the local adsorbed film thickness that corresponds to T0 (tˆ). It is important to note that the steady-state solutions of Eq. [22] with the above boundary conditions are invariant under translations: there is no preference as to where the contact line should be when the plate is heated uniformly. To avoid resulting numerical difficulties without changing the physical problem we fix the interfacial position far away from the contact line where heat fluxes are very small. The physics of the problem near the contact line is illustrated in Fig. 3. Dry areas are modeled by microscopic adsorbed films and the interface rapidly changes its curvature in the region near the contact line in order to match the adsorbed film. Evaporation is suppressed in the adsorbed film, but it is strong in the region of macroscopic liquid film near the contact line.
The global solution of the capillary-statics problem for an unconstrained bubble in two dimensions is a circular interface. The actual bubble shape shown in Fig. 1 is the result of deformation of the circle in the regions near the top and bottom plates, i.e., the areas where the interfacial heat and mass transfer becomes significant. At the leading order in the capillary number the curvature of the deformed shape outside the transition regions, and thus the length of the bubble, remain unchanged. In our scaled coordinates the value of the curvature is κ0 = 2. After the value of κ0 is determined, the global problem can be solved numerically. It consists of Eq. [20] with the matching condition [21] and two conditions at the symmetry point, Eq. [22] with condition [23] and three conditions in the adsorbed film. They are coupled together by the condition of the instantaneous mass balance: we assume that processes in the vapor occur sufficiently fast so that the total amount of evaporation from contact lines at any given moment in time is equal to the absolute value of the condensation mass flux near the top. We note that at the leading order in the capillary number the vapor pressure and volume remain constant. The solution procedure is as follows. The evolution equation [22] from the previous section is solved numerically for any given bottom temperature, T0 (tˆ). The total amount of evaporation, J¯ , is found from the solution. According to the condition of instantaneous mass balance, the value of J¯ is used as the boundary condition for the local solution near the top. Equation [20] now has four boundary conditions, so it is solved numerically to obtain the shape of the liquid film near the top. We use finite-difference methods to solve the equations [20] and [22]. The length of the domain of integration is typically no larger than L = 5; for larger values of the spatial coordinate the interface can be approximated by a parabolic shape without decreasing the accuracy of the solution; the (small) correction to mass flux from that region is found analytically. In order to evaluate the nonlinear term h 3 on our finite-difference mesh we take the average of its values at two neighboring points, as described in (15). Evaluating the nonlinear terms and derivatives at the endpoints involves mesh points outside the physical domain, where the function is extrapolated using the boundary conditions at the endpoints. The total number of mesh points is varied between 102 and 103 to verify convergence of the solutions. After the spatial domain is discretized, the partial differential equation is reduced to a system of first-order ordinary differential equations in time. DVODE solver is used to solve this system (16). 4. RESULTS
4.1. Steady States Let us choose representative value of parameters Kˆ = 2.0, δ = 0.04, ε = 10−3 for all simulations in this and the subsequent two subsections. These correspond to typical dimensional parameters given in (7) corresponding to organic liquids in microchannels with modest superheating above room temperature.
351
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
a
of the inverse thickness). There is a dimple in the thickness in Fig. 4b, which is indicative of the significant capillary pressure gradient that drives the flow of the condensing liquid away from the middle of the bubble.
5.0
4.0
_ h
4.2. Modulation of the Bottom Temperature
3.0
2.0
1.0
0.0 0.0
b
1.0
0.5
2.0
1.5
_ x
4.0
a
3.0
h^
When the temperatures of the plates become time-dependent, both the position of the contact line and the shape of the interface near the top change. The description of contact-line motion is relatively simple: it slides on the bottom plate in response to evaporation while the basic flow structure, shown in Fig. 3, does not change significantly. Such motion is illustrated in Fig. 5a, where the local interface shape near the contact line is plotted for different moments in time after a rapid change in the bottom plate temperature from T0 = 1.5 to T0 = 1.4. The initial position of the interface is shown by the dotted line. We also plot the ratio of the disjoining pressure to the capillary pressure, r , as a function 0.5
2.0
^
h 0.25
1.0
0.0 0.0
1.0
2.0
3.0
x^ 0.2
FIG. 4. Steady-state interface shape for the regions near the contact line (a) and near the top of the bubble (b) for Kˆ = 2.0, δ = 0.04, ε = 10−3 , T0 = 1.5, T1 = 0.8.
First let us consider the case when wall temperatures are constant, taking T0 = 1.5, T1 = 0.8. The time evolution from a constant-curvature initial condition indicates that the solutions approach steady-state shapes for which condensation and evaporation balance each other. Local interface shapes in this steadystate regime are plotted in Fig. 4. Several features of the steady solutions were discussed in (7), but it is useful to recount them here. Figure 4a illustrates the solution near the contact line at the bottom. Clearly, when x¯ becomes large, the interface curvature becomes constant, while as x¯ goes to negative infinity, the film thickness approaches the small constant value of the isothermal adsorbed film that mimics a dry region. In the region near the top there are no dry areas: the steady film is relatively thick, as seen in Fig. 4b. We note that hˆ and h¯ in Fig. 4 are the values of film thickness rescaled by the same number, C 1/3 . Thus there is an order of magnitude difference between the values of thickness at the top and at the bottom, which shows that near the top the disjoining pressure is negligible (it being proportional to the cube
0.4
0.6
^
x b
4.0
3.0
r 2.0
1.0 0.0
1.0
0.5
1.5
^
t
FIG. 5. The motion of the contact line on the bottom plate after instantaneous change of bottom temperature T0 from 1.5 to 1.4. (a) Interface shapes are shown at tˆ = 0.3 (dashed line) and tˆ = 0.6 (solid line). The initial position of the interface is also shown (dotted line). (b) The ratio of the disjoining pressure to the capillary pressure is shown as a function of time at xˆ = 0.2.
352
AJAEV, HOMSY, AND MORRIS
of time at a fixed location near the contact line (xˆ = 0.2) in Fig. 5b. The plot shows that the relative value of the disjoiningpressure term decreases rapidly as the film thickness increases. After tˆ ≈ 1 there is essentially no change in the position of the contact line and the pressure. Thus the motion of the contact line after a change in temperature can be described as rapid monotonic approach to a new steady-state position. Changes in the interfacial shape near the top occur on a much slower time scale and are more complicated to describe, so we focus on them in our discussion of unsteady simulations in this and the following subsections. In this section we consider modulations of the bottom temperature when the top temperature is held fixed at T1 = 0.8. This implies that the bottom plate temperature T0 (tˆ) in the evolution equation [22] is now a function of time. We choose it in the form of a square wave of dimensionless amplitude and frequency ω. Since the rescaled time tˆ is the time variable, it is convenient to also introduce rescaled frequency of forcing ωˆ according to ωˆ = C
2/3
a 0.7
0.6
h0 0.5
0.4
0
1000
500
^
t
b
0.8
0.7
ω. 0.6
The square-wave forcing can then be written in the form T0 (tˆ) =
T0 (0) − , n ωˆ −1 ≤ tˆ < (n + 1/2)ωˆ −1 T0 (0),
(n + 1/2)ωˆ −1 ≤ tˆ < (n + 1)ωˆ −1 ;
h0 0.5
0.4
n = 0, 1, 2, . . . . 0.3
Let us choose T0 (0) = 1.5 and = 0.1 and study the resulting changes in the liquid film shape near the top. We note that the bubble is always symmetric with respect to a vertical line, as seen in Fig. 1, and therefore the liquid film is also symmetric with respect to the y-axis shown in Fig. 2. Although our results give the entire film profile as a function of time, it is convenient to characterize the shape of the bubble by the film thickness at ˆ the point of symmetry, h 0 = h(0). We plot h 0 as a function of time in Fig. 6 for two different values of the scaled frequency of modulations: ωˆ = 0.01 and ωˆ = 0.002. Initially, there is a transient response, as clearly seen in Fig. 6a, but then the response settles into a periodic solution with period equal to that of the forcing. The amplitude of forced oscillations of the film thickness is smaller for larger values of frequency; the system does not have time to fully relax to the steady state before the next change in forcing. This regime is termed the high-frequency regime. A different physical picture is obtained for low-frequency modulations, when the interface almost relaxes to the steady state before the next change in temperature. The response amplitude then remains constant when the frequency is decreased further. The response in the low-frequency regime is illustrated in Fig. 6b. The two regimes are clearly seen in Fig. 7a, where the normalized amplitude is plotted as a function of modulation frequency for different values of forcing. The normalized amplitude a is defined as the ratio of the amplitude of variations of h 0 to the nondimensional forcing amplitude, . Since we vary the fre-
0
1000
500
1500
^
t
FIG. 6. Liquid film thickness at the symmetry point near the top as a function of time for square-wave modulations of the bottom temperature for different values of modulation frequency: (a): ωˆ = 0.01; (b): ωˆ = 0.002.
quency over a wide range of values, it is convenient to use a linear–logarithmic plot. The transition between the two regimes occurs at the values of the scaled frequency on the order of 10−3 . The corresponding dimensional value is approximately 100 Hz. The high-frequency regime correspond to the typical dimensional values of several kHz. We note that this is well below the value of 1 MHz at which the quasistatic approximation in the governing equations can break down as discussed in Section 2. It is clear from Fig. 7a that the response is strongly nonlinear—the normalized amplitude is different for different values of forcing. In fact, Fig. 7a shows that the value of normalized amplitude can be decreasing when is increased: nonlinear effects suppress oscillations. In order to better describe the response we also plot it in loglog coordinates as shown in Fig. 7b. Clearly there is a range of frequencies between 0.005 and 0.05 where a power-law decay of the amplitude is observed. However, as frequency is increased to 0.1, the rapid processes in the contact line-region start to affect the response, so the simple power law is no longer applicable.
353
GEOMETRICALLY CONSTRAINED VAPOR BUBBLES
a
6.0
5.0
4.0
a
3.0
2.0
1.0
0.0
-3
10
b
-2
^
ω
10
-1
10
1
10
pressure gradient needed for steady state, which actually means a smaller thickness. This interpretation is also supported by the fact that initial increase occurs on a time scale of the interface motion due to phase change (t¯ ≈ 1), while the slow decay is on the scale of viscous and conductive processes in the liquid. Thus in the short run the thickness is increased by decreasing the temperature, but over a longer period, the lower temperature leads to a thinner film. The nonmonotonic behavior of the interface leads to a different dynamic response curve. This is illustrated in Fig. 9, where the normalized amplitude of oscillations is plotted as a function of modulation frequency for different values of forcing . As before, the normalized amplitude a is the amplitude of variation of h 0 divided by . We observe again the high-frequency and low-frequency regimes of oscillations: transition between them occurs for ωˆ ≈ 10−3 . Contrary to the rapid power-law decay of amplitude for the bottom-temperature modulation, the response curve here flattens out for values of frequency near 0.05. This is mostly due to the presence of the initial rapid jump in interface
a
a
0.81
0
10
0.8
0.79
0.78
T1 0.77
-1
10
-3
10
-2
^ ω
10
-1
10
0.76
FIG. 7. The dynamic response over a range of frequencies: normalized amplitude of forced oscillations of film thickness at the point of symmetry is plotted vs frequency using linear–log (a) and log–log (b) scales for two different values of forcing amplitude: = 0.1 (solid lines) and = 0.05 (dashed lines).
4.3. Modulation of the Top Temperature Let us now study the effect of modulations of the top temperature on the bubble shape. We take a fixed value of the bottom temperature, T0 = 1.5, and then apply square-wave modulation of the top temperature of amplitude = 0.05, which is shown in Fig. 8a. The measure of film thickness, which is the value at the point of symmetry, h 0 , is plotted as a function of t¯ in Fig. 8b. We are not interested in transient effects, so all data are recorded starting at the zero value of the shifted time-variable t¯ = tˆ − 10ωˆ −1 . Typical results are shown in Fig. 8b. Clearly, there is a nonmonotonic behavior. When the temperature T1 decreases, the thickness initially increases rapidly, then slowly decreases. This result can be interpreted as follows. The initial rapid decrease of the temperature causes a rapid increase in condensation; however, the pressure gradient is not sufficient to remove the condensing liquid. The interface then adjusts its shape to provide the
0.75
0.74
b
0
200
400
_ t
600
800
0
200
400
600
800
0.7
0.6
h0 0.5
0.4
_ t
FIG. 8. The top wall temperature as a function of time for square-wave modulations of frequency ωˆ = 0.005 (a) and corresponding liquid film thickness at the point of symmetry (b) for K = 2.0, δ = 0.04, T0 = 1.5.
354
AJAEV, HOMSY, AND MORRIS
12.0
10.0
8.0
a
6.0
4.0
2.0
0.0
-3
10
-2
^
ω
10
-1
10
FIG. 9. Dynamic response to the modulations of the top temperature: normalized amplitude of oscillations of film thickness near the top is plotted as a function of scaled frequency for different values of forcing amplitude: = 0.25 (solid line), = 0.5 (dashed line), and = 0.75 (dot–dashed line).
position, which is essentially independent of the modulation frequency. Thus the amplitude decreases to some finite value of this initial jump as the frequency is increased. The nonlinear character of the response is also clearly seen: the amplitude a is increased when the forcing amplitude increases. This shows that nonlinear effects lead to larger amplitude of oscillations than may be expected from a linearized response model.
5. CONCLUSIONS
We develop a mathematical model for unsteady motion of a bubble between two parallel plates subject to time-dependent temperature conditions. Global time-dependent interface shapes are found in the limit of small capillary number. The interface is dominated by capillary forces everywhere except near the top plate and the contact-line regions at the bottom plate. The dynamic response of the bubble to square-wave modulations of plate temperatures is studied. When the bottom temperature is modulated and the top temperature is held fixed, the amplitude of forced oscillations is approximately constant
when the scaled frequency is below ωˆ ≈ 10−3 and then shows power-law decay for larger frequencies. When the top temperature is modulated with fixed bottom temperature, the dynamicresponse curve is flat in the low-frequency regime as well, but it also flattens out when the frequency is increased. This shape of the response curve is a result of the nonmonotonic behavior of the liquid film between the bubble interface and the top plate. When the top-plate temperature is decreased, the film thickness increases rapidly, but then slowly decays to a value which is smaller than the initial thickness. We have also studied the dynamic reponse as a function of forcing amplitude and found it to be strongly nonlinear even for very small values of the amplitudes of temperature variations. Nonlinear effects can result in oscillation amplitudes which are much larger than can be expected from the linear analysis of the response. ACKNOWLEDGMENTS GMH is grateful to the U.S. Department of Energy, Office of Basic Energy Sciences, for partial funding of this work. VA thanks Professor Ian Gladwell for useful advice on numerical computations.
REFERENCES 1. Lin, L., Pisano, A. P., and Carey, V. P., J. Heat Transfer 120, 735 (1998). 2. Ory, E., Yuan, H., Prosperetti, A., Popinet, S., and Zaleski, S., Phys. Fluids 12, 1268 (2000). 3. Potash, M., and Wayner, P. C., Int. J. Heat and Mass Transfer 15, 1851 (1972). 4. Moosman, S., and Homsy, G. M., J. Colloid Interface Sci. 73, 212 (1980). 5. Morris, S. J. S., J. Fluid Mech. 432, 1 (2001). 6. Ajaev, V. S., and Homsy, G. M., J. Colloid Interface Sci. 240, 259 (2001). 7. Ajaev, V. S., and Homsy, G. M., J. Colloid Interface Sci. 244, 180 (2001). 8. Anderson, D. M., and Davis, S. H., Phys. Fluids 7, 248 (1995). 9. Sharma, A., Langmuir 14, 4915 (1998). 10. Oron, A., and Bankoff, S. G., J. Colloid Interface Sci. 218, 152 (1999). 11. Lyushnin, A. V., Golovin, A. A., and Pismen, L. M., Phys. Rev. E 65, 021602 (2002). 12. Wilson, S. K., Davis, S. H., and Bankoff, S. G., J. Fluid Mech. 391, 1 (1999). 13. Burelbach, J. P., Bankoff, S. G., and Davis, S. H., J. Fluid Mech. 195, 463 (1988). 14. Schrage, R. W., “A Theoretical Study of Interface Mass Transfer,” Columbia Univ. Press, New York, 1953. 15. Zhornitskaya, L., and Bertozzi, A. L., SIAM J. Numer. Anal. 37, 523 (2000). 16. Brown, P. N., Byrne, C. D., and Hindmarsh, A. C., SIAM J. Sci. Stat. Comput. 10, 1038 (1989).