Stability of an evaporating and condensing liquid film flowing down an inclined plane

Available online at www.sciencedirect.com

ScienceDirect

Availableonline onlineatatwww.sciencedirect.com www.sciencedirect.com Available Energy Procedia 00 (2017) 000–000

ScienceDirect ScienceDirect

www.elsevier.com/locate/procedia

Energy Procedia 142 Energy Procedia 00(2017) (2017)3944–3949 000–000 www.elsevier.com/locate/procedia

9th International Conference on Applied Energy, ICAE2017, 21-24 August 2017, Cardiff, UK

Stability of an evaporating and condensing liquid film flowing down The 15th International Symposium on District Heating and Cooling an inclined plane. Assessing the feasibility ofAbderrahmane using the heat demand-outdoor * Hamid Ait temperature function forMasdar a long-term district heat demand Khalifa University of Science and Technology, Institute, Masdar City, P.O. Box 54224, Abu Dhabi, United Arabforecast Emirates I. Andrića,b,c*, A. Pinaa, P. Ferrãoa, J. Fournierb., B. Lacarrièrec, O. Le Correc a Abstract IN+ Center for Innovation, Technology and Policy Research - Instituto Superior Técnico, Av. Rovisco Pais 1, 1049-001 Lisbon, Portugal b

Veolia Recherche & Innovation, 291 Avenue Dreyfous Daniel, 78520 Limay, France

c Département Systèmes et Environnement - IMT Atlantique, 4 rueexchangers Alfred Kastler, 44300 Nantes, France Film fluid flows, driven down an Énergétiques inclined plane by gravity is encountered in heat and cooling towers. These flows can undergo evaporation (condensation), which might affect their stability and in turn affect the heat transfers. This paper presents an investigation of the influence of evaporation (condensation) on the stability of the film down an inclined plane. Using energy integral method, to reduce the governing equations of the problem into two-equation model, we showed that Abstract (condensation) destabilizes (stabilizes) the film fluid flow. evaporation

heating networks are commonly in the literature as one of the most effective solutions for decreasing the ©District 2017 The Authors. Published by Elsevier addressed Ltd. greenhouseunder gas emissions from of thethe building sector. These of systems high investments are returned through the heat Peer-review responsibility scientific committee the 9threquire International Conferencewhich on Applied Energy. sales. Due to the changed climate conditions and building renovation policies, heat demand in the future could decrease, prolonging theexchangers; investment return period. Keywords: Heat condencers; film flow; evaporation, condensation. The main scope of this paper is to assess the feasibility of using the heat demand – outdoor temperature function for heat demand forecast. The district of Alvalade, located in Lisbon (Portugal), was used as a case study. The district is consisted of 665 that vary in both construction period and typology. Three weather scenarios (low, medium, high) and three district 1.buildings Introduction renovation scenarios were developed (shallow, intermediate, deep). To estimate the error, obtained heat demand values were compared with results from a dynamic demand model, and validated authors. Film fluid flow, driven down anheat inclined plane bypreviously gravity, isdeveloped a canonical problem by fortheseveral engineering fluid The results showed that when onlyencountered weather change is considered, the margin of errortowers. could be acceptable for some applications flow problems, including those in heat exchangers and cooling This canonical problem is the (the error in annual demand wasefforts lower than all weather scenarios considered). However, after introducing model; renovation subject of extensive modeling since20% the for fifties. The first model for this problem is a single-equation it scenarios, the error value increased to 59.5% (depending on low the weather andnumbers renovation scenarios combination was proposed by Benney [1]. Thisupmodel is valid only for Reynolds around its critical valueconsidered). at which The value of slope coefficient increased on average within the range of 3.8% up to 8% per decade, that corresponds to the the flow becomes unstable. Benney's model also suffers from finite-time blow up behavior. To extend the modeling decrease in the number of heating hours of 22-139h during the heating season (depending on the combination of weather and efforts to moderate and high Reynolds numbers an ad-hoc model, based on boundary-layer approximation (Integral renovation scenarios considered). On the other hand, function intercept increased for 7.8-12.7% per decade (depending on the Boundary Layer approximation) of the Navier–Stokes equations was developed [2]. This model consists of a twocoupled scenarios). The values suggested could be used to modify the function parameters for the scenarios considered, and equation model for the film flow height and flow rate. This model, known as Shkadov’s model, does not predict a improve the accuracy of heat demand estimations.

right critical Reynolds. This drawback was corrected using energy and weighted residuals integral methods [3,4,5]. © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the Scientific Committee of The 15th International Symposium on District Heating and * Corresponding author. Tel.: 97128109306; Cooling. E-mail address: [email protected]

Keywords: Heat demand; Forecast; Climate change 1876-6102 © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the scientific committee of the 9th International Conference on Applied Energy.

1876-6102 © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the Scientific Committee of The 15th International Symposium on District Heating and Cooling.

1876-6102 © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the scientific committee of the 9th International Conference on Applied Energy . 10.1016/j.egypro.2017.12.301

2

Hamid Ait Abderrahmane et al. / Energy Procedia 142 (2017) 3944–3949 Author name / Energy Procedia 00 (2017) 000–000

3945

When a temperature gradient is imposed across the film, the flow might become unstable owing to thermocapillary, buoyancy forces and Marongoni forces [6]. The dynamics of a film falling down a uniformly heated wall was conducted by Kalliadasis et al. (2003) [7]. These authors used an Integral Boundary Layer approximation to reduce the problem into three-equation model for the film height, flow rate and interfacial temperature. This model was improved using high-order weighted residuals approach with polynomial expansions for both velocity and temperature fields [8,9]. When, the inclined plane surface is heated (cooled) at constant temperature that is higher (lower) than the saturation temperature of the fluid, evaporation (condensation) occurs at the liquid-vapour interface. Evaporation (condensation) disturbs the liquid-vapour interface and can affect the film flow instability, mass and heat transfers at the free surface [6]. In this paper we investigate the effects of evaporation (condensation) on the stability of a viscous Newtonian thin liquid film flowing down an inclined plane. We expand the two-equation model proposed by Ait Abderrahmane and Vatistas [3] by including the effect of phase changes at the free surface. Using the energy integral method with self-similar velocity profile, we reduced the governing equations in a two-equation model. This describes the spatial-temporal evolution of the film height and the flow rate. The influence of phase change at the the film free surface is investigated using the obtained reduced model. Nomenclature b=(a/2-a)(ρvr2/kTs)(M/2πRgTs)1/2 m-1. M: molecular weight, kg/mol. Rg: universal gas constant, J/(molK). Ts: temperature at the free surface. R: latent heat, J/kg. ρv: vapour density ϵ: shallow water parameter =h0/L h0: initial film depth L: characteristic length σ: surface tension coefficient We: Weber number. Re: Reynolds number ΔT: temperature difference= Tw-Ts ρ: fluid density cp: thermal capacity at constant pressure ρ0= ρv/ ρ Nd: ζ2/ρ0 Pr. ζ= cp ΔT/r. Tw: wall temperature. µ: fluid dynamics viscosity k: film flow thermal conductivity. 2. Formulation of the problem The problem consists of an evaporating (condensing) two-dimensional flow of a viscous incompressible liquid flowing down an inclined solid surface with angle, θ, and put at a temperature Tw. The schematic of the physical model, in Cartesian coordinate system, is shown in Fig. 1. The simplified problem is described by the continuity, boundary layer and heat equations: u x  vy  0

(1)

Hamid Ait Abderrahmane et al. / Energy Procedia 142 (2017) 3944–3949 Author name / Energy Procedia 00 (2017) 000–000

3946

u t  uu x  vu y 

1 p x  g sin    u yy  0 

3

(2) (3)

Tyy  0

The boundary conditions associated with the problem above are the non-slip condition at the solid surface, the normal and tangential stresses boundary conditions at the vapour-liquid interface.

u  v  0 , T  Tw at y=0, 

 J2  n  T.n.n  2K v

and

  T.n. t  . t

T  Ts

at y=h

(4)

P and  are the pressure and kinematics viscosity of liquid.  v is vapour density. J is mass flux at the interface

due to the evaporating (condensing) fluid. K is the mean curvature of the interface. T is the stress tensor of liquid.   The unit vectors n and t are the outward normal and tangent vectors, respectively. Ts is the free surface temperature.

Ti

Tw

Fig. 1: Schematic of the problem.

One should add the kinematic condition at the free surface to the governing equations above

ht  u hx  v  0

at y=h

(5)

The linearized constitutive equation of mass flux, derived from the kinetic theory, is given as:

a r  M J v  Ts  2R g Ts

1

2  Ti  Ts    Ti is the interior temperature in the film fluid. a is the accommodation coefficient. M , R g Rg and r are the molecular weight, universal gas constant and latent heat, respectively. We assume that there is no accumulation of heat in the liquid and all the heat exchange between the wall and the free surface contributes to the phase change at

T  Ti a v r 2  the vapour-liquid interface, only. q  k w h Ts

 M   2R g Ts 

1

2  Ti  Ts   

k is the liquid thermal conductivity. The vapor pressure due to evaporation and condensation can be expressed as:

4

Hamid Ait Abderrahmane et al. / Energy Procedia 142 (2017) 3944–3949 Author name / Energy Procedia 00 (2017) 000–000

p v   v u 2v   1n

1  b k T     v  r1  b h  

3947

2

1

 M 2   T  Ts  and u v is the velocity of vapour leaving from or entering into the interface.  2R g Ts  i   When n  0 , the film is evaporating and the vapour is leaving from the interface. When n  1 , the film is condensing. When ΔT=0, the film flow is in isothermal state and the vapour pressure is zero. The averaged velocity of the unperturbed flow a r 2 Where b  v Ts

h0

u

 0

g sin   2 y y 2  2   h 0 h 02

2  dy  g sin h 0  2 

We assume that the velocity profile of the perturbed flow should slightly depart from the parabolic velocity profile of the base flow solution. u  a(x, t ) y  b(x, t ) y 2

The coefficients a(x, t) and b(x, t ) can be expressed in term of the physical quantities. a ( x, t ) 

3q h

2



3q 3 1  s and b( x, t )   2  s 2 h h

(6)

Where q , h and  s are the flow rate, the film fluid thickness and the shear stress at the interface fluid-vapour interface. It can be easily shown that the shear stress  s can be expressed in terms of the two fundamental variables

h0 , L is the L characteristic length x  direction. Integrating the continuity (1) and the kinematic boundary condition equation (5) with the aid of the Leibniz rule we obtain.

q and h [3]. The shear stress is of order two with respect to the shallow water parameter  

ht  qx  0

(7)

The second equation is obtained using the energy integral method, which consist of multiplying the x-momentum equation, by the velocity component u and integrating along the film fluid depth, as follows: h



 u  u

t

 uu x  vu y 

o

 1 p x  g sin    u yy dy  0  

(8)

The governing equations and the boundary conditions are nondimensionalized using scales based on of the initial fluid depth, h0, and the averaged velocity of the unperturbed flow. The v component of the velocity is obtained with the aid of the continuity equation (1). The expression of the pressure P , accurate up to the first order with respect to the shallow water parameter, is obtained by integrating the y-momentum with the boundary condition (4). Introducing the velocity components and the expression of the pressure in (8) one can obtain the second equation of the reduced model. The equations (7) and (8) form the two-equation model, which describes the spatial- temporal evolution of the locals flow rate q( x, t ) and film depth h(x, t). 3. Linear stability A linear stability analysis of the problem is performed through the obtained model. We investigated the

Hamid Ait Abderrahmane et al. / Energy Procedia 142 (2017) 3944–3949 Author name / Energy Procedia 00 (2017) 000–000

3948

5

evolution of infinitesimal disturbances on the primary flow. These perturbations may decay, remain at the same magnitude, or grow exponentially and lead to turbulent flow regime. The first type of disturbance is called stable, the second neutrally stable and the latter unstable. The neutral stability curves separate the instability and stability region in the flow parameters plane. We are interested in the impacts of phase change on the stability at different Reynolds numbers and inclination angles. The linear stability analysis involves the investigation of phase change on the neutral stability curves (marginal stability conditions) and perturbation celerity. We consider infinitesimal disturbances from the uniform flow h=1+H, and q=1+Q. The flow rate perturbation Q can be eliminated in favour of the perturbed film thickness H described by a single linear equation is obtained as:  27 H xxt 36 H xxx 2   5 Re  5 Re

3H xxt 2 1n Nd H xx 54 6  3Ht 9 Hx   102 H xt  cot   0  H xx  h tt     2 Re 5  Re Re 35 Re   35

(9)

2

1  c p T  cp hu   is the phase change parameter and Pr  is a Reynolds number, Nd  is the Prandtl  0 Pr  r   k number . Considering a wave like disturbance in a Fourier expansion: Re 

H(x, t ) 

1 i kxct  e

(10)

The characteristic equation for the wave velocity is obtained by substituting (10) into Eq. (9). The characteristic equation relates the wave number, k , and the wave celerity c . Since we are interested in temporal stability, we assume the wave celerity as complex number c  c r  i c i . c r is the wave celerity with which the perturbation propagates along the free surface, while c i is the growth rate of the perturbation. This equation is solved for c i and k . The threshold for linear stability is obtained by substituting c  c r  i c i and then setting the wave growth rate zero (𝑐𝑐𝑖𝑖 = 0). The neutral wave celerity is: c r  3  k cot  O( 2 ) . Note that the speed of neutral long waves

during evaporation or condensation is the same as in isothermal condition. The neutral stability curve equation is thus:

cot   We 

2  1n 6 1 9211 N d  Re  Re We k 2  Re k 2  0 3 Re 5 3 2240 

h 02 u

(11)

is the Weber Number. The neutral wave speed at this condition is equal to 3, and the critical condition

for instability is: Re 

5 cot  cot   25  80 tan 2  1n Nd . 12 12

For the isothermal flow, the critical Reynold

5 number reduces to the one found in the literature [3,4,5], Re  cot  6 4. Results and discussion Evaporation (condensation) induces a thermal non-equilibrium at the interface. We show in Fig. 2 (a,b) the effect of phase change on film stability. Fig. 2a shows that the neutral stability curves against Reynolds number. This figure shows that condensation decreases the extent of the region where the flow is unstable. However, evaporation has

6

Hamid Ait Abderrahmane et al. / Energy Procedia 142 (2017) 3944–3949 Author name / Energy Procedia 00 (2017) 000–000

3949

opposite effect; it increases the extent of the instability region. Figure 2a also shows that the effect of phase change on film flow instability is significant at low Reynolds numbers. At higher Reynolds numbers, the film flow stability become less dependent of phase change; the inertia controls the flow stability. Figures 2b shows the growth rate of the instability, 𝑐𝑐𝑖𝑖 , versus of the wave number, k. The maximum growth corresponds to the most unstable mode. Figures 2b shows that evaporation (condensation) increases (decreases) the wave number of the instable modes. Evaporation (condensation) slightly moves the most unstable mode to higher (lower) wave numbers. More importantly evaporation (condensation) destabilise (stabilise) the fluid film flow by reducing the growth of the instability.

a)

b)

Fig. 2. a) Neutral stability curves and b) growth rates of the instability

5. Conclusion The effects of evaporation and condensation on the stability of a film fluid falling down an inclined plane are investigated. The governing equations are reduced in two-equation model, which describes the spatiotemporal evolution of the film height and fluid flow rate. The reduced model was derived using Energy Integral Method. Similar to previous works, we have shown that evaporation destabilizes the flow while condensation stabilizes the flow. Results indicate that effects of phase change on neutral stability are only noticeable at low Reynolds numbers. References [1] Benny, D. J., 1966, “Long Waves on Liquid Films,” J. Math. Phys., 45, pp. 150-155. [2]Shkadov, V. 1967 Wave flow regimes of a thin layer of viscous fluid subject to gravity. Izv. Ak. Nauk SSSR, Mekh. Zhi. Gaz 2, 43, 29-34. [3] Ait Abderrahmane., H., and Vatistas., G. H., 2007, “Improved two-equation model for thin layer fluid flowing down an inclined plane problem”, Physic of Fluids, 19, 098106 pp 1-4 [4] Amaouche, M., Mehidi, N. & Amatousse, N. 2005 An accurate modeling of thin film flows down an incline for inertia dominated regimes. Eur.J.Mech. B 24, 49-70. [5] Ruyer-Quil, C. & Manneville, P. 2000 Improved modeling of flows down inclined planes. Eur.Phys. J. B 15, 357-369. [6] Oron A., Stephen H. D., Bankoff S. G., 1997, “Long scale evolution of thin liquid films”, Reviews of Modern Physics, 69, pp.931-980 [7] Kalliadasis, S., Kiyashko, A. and Demekhin, E. A. 2003a Marangoni instability of a thin liquid film heated from below by a local heat source. J. Fluid Mech. 475, 377–408. [8] Ruyer-Quil, C. Scheid, B., Kalliadasis, S., Velarde, M.G. & Zeytounian, R.Kh. 2005,Thermocapillary long waves in a liquid film flow. J. Fluid Mech. 538, 199-222. [9] Scheid, B., Ruyer-Quil, C., Thiele, U., Kabov, O., Legros, J. & Colinet, P. 2005b Validity domain of the Benney equation including Marangoni effect for closed and open flows. J. Fluid Mech. 527, 303-335.