EPJ Web of Conferences 67, 0 20 44 (2014) DOI: 10.1051/epjconf / 2014 6702044 C Owned by the authors, published by EDP Sciences, 2014 Concept of CFD model of natural draft wet-cooling tower ﬂow T. Hyhlı́k1,a CTU in Prague, FME, Department of Fluid Dynamics and Thermodynamics, Technická 4, 166 07 Prague Abstract. The article deals with the development of CFD model of natural draft wet-cooling tower ﬂow. The physical phenomena taking place within a natural draft wet cooling tower are described by the system of conservation law equations along with additional equations. The heat and mass transfer in the counterﬂow wet-cooling tower ﬁll are described by model [1] which is based on the system of ordinary diﬀerential equations. Utilization of model [1] of the ﬁll allows us to apply commonly measured ﬁll characteristics as shown by [2].The boundary value problem resulting from the ﬁll model is solved separately. The system of conservation law equations is interlinked with the system of ordinary diﬀerential equations describing the phenomena occurring in the counterﬂow wet-cooling tower ﬁll via heat and mass sources and via boundary conditions. The concept of numerical solution is presented for the quasi one dimensional model of natural draft wet-cooling tower ﬂow. The simulation results are shown. 1 Introduction Natural draft cooling towers are widely used in industry especially in large power plants and in chemical industry. Low eﬃciency of the cooling process of these towers will cause a loss of power generation or product quantity. Careful and accurate analyses of cooling towers can result in improvement of eﬃciency of cooling process. The aim of cooling tower is to transfer waste heat from cooled water into the atmosphere. The air which is ﬂowing through the cooling tower is warmed up and its humidity is increasing which leads to cooling of ﬂowing water. The density of warmed air is decreasing unlike the surrounding air and this density diﬀerence produce the natural draft. Frequently encountered type of natural draft cooling tower is concrete structure of rotational hyperboloid shape. Typical cooling tower of this kind can be larger then 100 m in height and diameter. Cooled water is sprayed above the ﬁll using the set of sprayers. In the channels of counterﬂow wet-cooling tower ﬁll of ﬁlm type water ﬂows vertically down through the ﬁll as a liquid ﬁlm. Air is driven by a tower draft and ﬂows vertically in the opposite direction. Heat and mass transfer occurs at the water and air interface. Evaporation and convective heat transfer cool the water, what leads to increase of air humidity and temperature. The water is leaving the ﬁll and falling down in the form of rain. There is a pond in the bottom part of the cooling tower where water is collected. The most intense heat and mass transfer occurs in the ﬁll zone. The simplest model of heat and mass transfer in the ﬁll is Merkel model [3] which is industry standard now. The more advanced models are e.g.. Poppe model [4] and model [1]. As an example of CFD models of natural draft cooling tower ﬂow can be mentioned e.g. [5] and commercial code based [6]. a e-mail: [email protected] 2 Balance laws The moist air ﬂowing through the cooling tower will be considered as homogeneous mixture. The total value of extensive quantity ϕ is expressed as ρα wα ϕα , wα = , wα = 1, (1) ϕ= ρ α α where wα is the mass fraction of component α. Local form of the balance law is [7] k k ∂(ϕα ρα ) ∂(ϕα ρα vk ) ∂ ϕα jDα − jα + − σα = 0. (2) + ∂t ∂xk ∂xk Diﬀusive ﬂux density jkDα is expressed through the center of mass velocity wα vαk , (3) vk = α jkDα where = ρα (vαk − vk ), jkDα = 0. (4) (5) α Diﬀusive ﬂux jkDα can be written using Fick’s law [8] as ∂ρα , (6) ∂xk where Dα,m is diﬀusion coeﬃcient. Density of ﬂux of quantity ϕ we express as jkα (ϕα ). (7) jk (ϕ) = jkDα = −Dα,m α Density of production of ϕ is σα (ϕα ). σ(ϕ) = (8) α Local form of balance law for quantity ϕ can be written for the mixture as ∂(ρϕ) ∂(ρϕvk ) ∂ jk (ϕ) − − σ(ϕ) = 0. (9) + ∂t ∂xk ∂xk This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Article available at http://www.epj-conferences.org or http://dx.doi.org/10.1051/epjconf/20146702044 EPJ Web of Conferences 3 Moist air The moist air is mixture of dry air and superheated or saturated steam. For a description of the cooling tower ﬂow, let us consider the humid air as a mixture of dry air and water vapour. The liquid phase will be ignored in the model of moist air ﬂow. The density of the moist air is ρ = ρa + ρv , (10) where the densities can be expressed using equation of state as pa Ma pv Mv ρa = , ρv = . (11) RT RT Mass fraction of dry air and mass fraction of water vapour we can express as wa = where ρa ρv , wv = , ρ ρ wa + wv = 1. (13) (14) Speciﬁc gas constant of mixture can be expressed by means of speciﬁc gas constants of both components r = wa ra + wv rv . ρv pv Mv Vρv = = . Vρa ρa pa Ma 1 x , wv = . 1+x 1+x (17) (18) Equation of state of a mixture of dry air and water vapour is expressed using Dalton’s law, the deﬁnition of speciﬁc humidity and the deﬁnition of the overall density as ρ= p Ma (1 + x) , RT 1 + Ma x Mv (19) where Ma is molecular weight of dry air, Mv is molecular weight of water vapour, R is universal gas constant and T is temperature. Pressure can be expressed using the equation of state as a ρRT 1 + M Mv − wv (Mv − Ma ) Mv x (20) = ρRT p= Ma (1 + x) Ma Mv The total energy of moist air is expressed as the sum of internal energy and kinetic energy e=u+ vk vk . 2 (23) This deﬁnition of enthalpy is standard and is used when working with moist air, but due to the fact that the dry air and water vapour can be considered under atmospheric conditions as an ideal gases we can deﬁne ad hoc internal energy of unsaturated humid air as u = (wa cva + wv cvv )T = cv T. (24) The total energy of unsaturated moist air is expressed in the form r vk vk vk vk =ρ , (25) ρe = ρ cv T + T+ 2 κ−1 2 where κ is Poisson coeﬃcient κ= cp . cv (26) (16) Mass fractions can be expressed using speciﬁc moisture wa = h = wa c pa ta + wv (l0 + c pv ta ). 4 Governing equations Speciﬁc moisture indicates the weight of steam per kilogram of dry air x= (22) where ta = T − 273.15 is temperature is degrees of Celsius. In the case where we are working with enthalpy per kilogram of substance we can introduce (15) Speciﬁc heat capacity of the mixture is expressed as cv = wa cva + wv cvv . h1+x = c pa ta + x(l0 + c pv ta ), (12) Pressure can be expressed using Dalton’s law p = pa + pv . In the case of moist air we are usually working with enthalpy relative to one kilogram of dry air and x kilograms of water vapour (21) Continuity equation for dry air can be written as ∂ρa ∂ ∂(ρa ) ∂(ρa vk ) Da,v = 0, − + ∂t ∂xk ∂xk ∂xk (27) where Da,v is diﬀusion coeﬃcient of dry air in the water vapour. We can express the continuity equation of water vapour in the form ∂ρv ∂ ∂(ρv ) ∂(ρv vk ) Dv,a = σv (xk , t), − (28) + ∂t ∂xk ∂xk ∂xk where σv (xk , t) represents the density of source of water vapour and diﬀusion coeﬃcient of water vapour in the air is Dv,a = Da,v . The system of momentum equations can be written in the form ∂(ρvi ) ∂ |v|2 ρvi vk − σik = −ρgδi3 −ζ(xk , t)ρ δi3 , (29) + ∂t ∂xk 2 where −ρgδi3 represents the gravitational force acting on the ﬂowing moist air, σik is the stress tensor in a ﬂowing ﬂuid and ζ represents a loss coeﬃcient per meter. The energy equation is written in the form ∂ ∂(ρe) ρvk e − σk j v j + q˙k = σq (xk , t) − gρv3 , (30) + ∂t ∂xk where σq (xk , t) is density of heat source where only sensible part of heat source is considered because of the deﬁnition of total energy which is based on equation (25). 02044-p.2 EFM 2013 5 Quasi one-dimensional ﬂow equations The modiﬁcation of governing equations for the case of quasi one-dimensional ﬂow is performed for control volume A(z)dz. Dry air continuity equation has form ∂(ρa A(z)) ∂(ρa vA(z)) + = 0. ∂t ∂z (31) We include force acting on side wall −pdA(z) in the case of momentum equation. Finally we have the system of equations for quasi one-dimensional ﬂow of the form ∂(WA) ∂(FA) + = Q, ∂t ∂z where ⎡ ⎡ ⎤ ⎤ ⎢⎢⎢ ρa v ⎥⎥⎥ ⎢⎢⎢ ρa ⎥⎥⎥ ⎢⎢⎢ ρv v ⎥⎥⎥ ⎢⎢⎢ ρv ⎥⎥⎥ ⎥⎥ W = ⎢⎢⎢⎢ ⎥⎥⎥⎥ , F = ⎢⎢⎢⎢ 2 ⎢⎣ ρv + p ⎥⎥⎥⎦ ⎢⎣ ρv ⎥⎦ ρe v(ρe + p) and vector of sources is ⎡ ⎤ 0 ⎢⎢⎢ ⎥⎥⎥ ⎢⎢⎢ ⎥⎥⎥ A(z)σ (z) v ⎥. Q = ⎢⎢⎢⎢ dA v2 ⎥ ⎢⎢⎣ p dz − A(z)(ρg + ζρ 2 ) ⎥⎥⎥⎥⎦ A(z)(σq (z) + ρgv) (32) and modiﬁed vector of sources is ⎤ ⎡ A(z)σv (z) ⎥⎥ ⎢⎢⎢ 2 ⎢⎢⎢⎢ p dA − A(z)(ρg + ζρ v ) ⎥⎥⎥⎥⎥ 2 ⎥ Q = ⎢⎢⎢ dz ⎥⎥⎥ . ⎢⎢⎣ A(z)σv (z) ⎥⎦ A(z)(σq (z) + ρgv) where ṁa is mass ﬂow rate of dry air. The change in water mass ﬂow rate can be expressed using mass transfer coefﬁcient αm as dṁw = αm (x (tw ) − x)dA, (33) (34) (35) (36) The system of governing equations is necessary to close. There are two unknown sources in the source vector (36), i.e. density of water vapour source σv (xk , t) and density of heat source σq (xk , t). The sources in the ﬁll zone can be calculated using one dimensional model [1]. ṁa dh1+x = α(tw − ta )dA + hv (tw )dṁw , (40) where α is heat transfer coeﬃcient and hv (tw ) is enthalpy of water vapour. 6.2 Model of Klimanek & Białecky [1] To derive the system of ordinary diﬀerential equations we have to choose independent variable. The model [1] is based on the selection of spatial coordinate z as independent variable contrary to Poppe’s model, see e.g. [4] which is based on the choice of water temperature tw . The interface area can be expressed using variable z as (41) where aq is the transfer area per unit volume and A f r is the cross sectional area of the ﬁll. We can derive equation for the change of the water mass ﬂow rate from equation (38) dṁw = αm aq A f r (x (tw ) − x). dz (42) To obtain the equation for the change of speciﬁc humidity in the ﬁll we can substitute equation (42) into equation (37) dx αm aq A f r (x (tw ) − x) . = dz ṁa 6.1 Balance laws Due to the complexity of two phase ﬂow occurring in the wet-cooling tower ﬁll the one dimensional models of heat and mass transfer are used. These models are based on few assumptions which allow to create simpliﬁed one dimensional models. The ﬁrst assumption talks about neglecting of the eﬀects of horizontal temperature gradient in the liquid ﬁlm, horizontal temperature gradient in air temperature (39) where h1+x is enthalpy of air water vapour mixture and hw is enthalpy of water. The change in total enthalpy can be evaluated using interface parameters similarly like in the case of mass balance dA = aq A f r dz, 6 Computation of sources in ﬁll zone (38) where x (tw ) is saturated speciﬁc humidity at tw and dA is inﬁnitesimal contact area. The energy balance can be written in the form ṁa dh1+x = m˙w dhw + hw dṁw , The system of governing equations can be modiﬁed by summing of individual continuity equations to get mixture continuity equation. The water vapour continuity equation is modiﬁed to get water vapour mass fraction equation. ⎡ ⎤ ⎤ ⎡ ρv ⎢⎢⎢ ρ ⎥⎥⎥ ⎥ ⎢⎢⎢ ⎢⎢⎢ ρv ⎥⎥⎥ ⎢⎢⎢ ρv2 + p ⎥⎥⎥⎥⎥ ⎥ ⎥⎥ , F = ⎢⎢⎢ W = ⎢⎢⎢⎢ ⎢⎢⎣ ρvwv ⎥⎥⎥⎥⎦ ⎢⎣ ρwv ⎥⎥⎥⎦ ρe v(ρe + p) and humidity, see e.g. [4]. The second one states that temperatures and humidity are represented only by averaged value for each vertical position. We are also assuming that at the interface of two phases, there is a thin vapour ﬁlm of saturated air at the water temperature, see e.g. [4]. The derivation of every one dimensional model of heat and mass transfer in the ﬁll is based on balance laws. We have four variables in this problem: ta temperature of air, tw temperature of water, x speciﬁc humidity and ṁw water mass ﬂow rate. Mass balance of the incremental step of the ﬁll is given by (37) dṁw + ṁa dx = 0, (43) Diﬀerentiation of equation (22) leads to dta dx dh1+x = (c pa + xc pv ) + (l0 + c pv ta ) . dz dz dz (44) The left hand side of equation (44) can be substituted from equation (40), the last term can be expressed using (43) and 02044-p.3 EPJ Web of Conferences 35 after application of the deﬁnition of Lewis factor, which is equal to the ratio of heat transfer Stanton number S t to the mass transfer Stanton number S tm St α = , S tm c p αm (45) t / ◦C Le f = 30 25 20 where c p is constant pressure speciﬁc heat capacity of moist air (46) c p = c pa + xc pv , 15 ta tw twb we can get αm aq A f r dta Le f (tw − ta ) c pa (ta ) + = dz ṁa (c pa (ta ) + x c pv (ta )) + x c pv (ta ) + (c pv (tw )tw − c pv (ta )ta )(x (tw ) − x) . (47) 10 The application of the Lewis factor in the previous equations simpliﬁed the problem to ﬁnd experimentally only mass transfer coeﬃcient αm and calculate heat transfer coeﬃcient α using known value of Lewis factor. The most commonly used formula for the calculation of Lewis factor is Bosnjakovic formula [9]. The driving force of evaporation process is (x (tw ) − x (ta ))) in the case of supersaturation. The system of ordinary diﬀerential equations has to be changed [1]. The enthalpy of supersaturated air is 0.4 0.6 z / m 0.8 1 40 x / g/kg 35 30 25 20 15 10 5 x x (ta ) x (tw ) 0 0.2 0.4 0.6 z / m 0.8 1 Fig. 2. Distribution of speciﬁc humidity x, saturation speciﬁc humidity at air temperature x (ta ) and saturation speciﬁc humidity at water temperature x (tw ) 235 σv A / kg/s/m h1+x = c pa ta + x (ta )(l0 + c pv ta ) + (x − x (ta ))cw (ta )ta . (49) Equations for the air and water temperature can be derived similarly like in the case of under-saturated moist air, see [1]. The system of four ordinary diﬀerential equations mentioned in this section represents the boundary value problem. We know air temperature tai and speciﬁc humidity xi at air inlet and water temperature twi and water mass ﬂow rate ṁi on the opposite site of the ﬁll because air and water ﬂows in the opposite direction. There is additional unknown set of parameters in the system of equations, i.e. αm aq A f r . These parameters have to be solved by using experimentally obtained characteristics of the ﬁll [1,2]. 0.2 Fig. 1. Distribution of air temperature ta , water temperature tw and adiabatic saturation temperature twb By using equation (39) after substitution of equation (42) and equation (40) we derive equation for the change of water temperature αm aq A f r dtw Le f (c pa (ta ) + x c pv (ta ))(tw − ta ) + = dz m˙w cw (tw ) (48) + (x (tw ) − x)(c pv (tw )tw − c pw (tw )tw + l0 ) . 0 230 225 220 215 210 205 200 0 0.2 0.4 0.6 z / m 0.8 1 Fig. 3. Distribution of density of mass source 6.3 Fill zone simulation example The calculation was performed for a ﬁll of height H = 1 m. The air inlet mass ﬂow rate is equal to ṁa = 14333 kg/s and water inlet mass ﬂow rate is ṁw = 17200 kg/s. Inlet water temperature is twi = 34.9◦ C. Air inlet temperature is tai = 15.7◦ C and speciﬁc humidity at air inlet is xi = 7.622 g/kg. The atmospheric pressure is p = 98100 Pa. Outlet water temperature is two = 25.7◦ C. The numerical solution of the ﬁll zone model was performed using Runge-Kutta method combined with shooting method [2]. The calculated distributions of temperature and humidity are depicted in ﬁgures 1 and 2. From the mentioned ﬁgures is possible to observe supersaturation in the upper part of of the ﬁll, it means the air temperature is equal to adiabatic saturation temperature and speciﬁc humidity is slightly higher then saturation humidity. The distribution of the density of mass source is shown in the ﬁgure 3 and the density of sensible part of heat source is in the ﬁgure 4, where density of mass source of water vapour can be expressed as 02044-p.4 σv = ṁa dx A(z) dz (50) EFM 2013 1.65e+08 Total pressure p0 and total temperature T 0 are prescribed at the inlet boundary condition, where vector of conservative variables is calculated using extrapolated Mach number, see e.g. [11]. The outlet boundary condition is based on the prescription of static pressure and extrapolation of three components of the vector of conservative variables. σqS A / W/m 1.64e+08 1.63e+08 1.62e+08 1.61e+08 1.6e+08 1.59e+08 7.2 Calculation example 1.58e+08 1.57e+08 Natural draft cooling tower of rotational hyperboloid shape with 2 π 0.006977z2 − 1.2764z + 131.61 (59) A(z) = 4 1.56e+08 1.55e+08 0 0.2 0.4 0.6 z / m 0.8 1 Fig. 4. Density of sensible part of heat source distribution and sensible part of heat source is σ qS = ṁa (dh1+x − l0 dx) . A(z) dz (51) 7 Cooling tower ﬂow computation Numerical solution is based on the (32), where because of the deﬁnition of total energy (25) the density of heat source should be calculated as σq = σqS + 273.15c pv (ta )σv (52) 7.1 Numerical method The numerical solution utilizes cell-centered Finite Volume Method in conservative form. Explicit MacCormack scheme is applied to solve system (32). Predictor step is expressed as Win+1/2 = Win − Δt n Fi+1 Ai+1/2 − Fin Ai−1/2 + μi Δt + Qi Δz μi (53) and corrector step is Win+1 = 1 Δt n+1/2 1 n Wi + Win+1/2 − F Ai+1/2 − 2 2 μi i 1 Δt n+1/2 Ai−1/2 + Qi Δz, +Fi−1 2 μi (54) is selected [12]. Tower height is Ht = 150m and ﬁll zone is placed at the height of 11.5m. Calculation was performed for ﬁll height H = 1m. Water inlet mass ﬂow rate is ṁw = 17200 kg/s. Inlet water temperature is twi = 34.9◦ C. Air inlet temperature is tai = 15.7◦ C and speciﬁc humidity at air inlet is xi = 7.622 g/kg. The atmospheric pressure is p = 98100 Pa. Outlet water temperature is two = 25.7◦ C. Prescribed parameters allows to quantify ﬁll requirement to cool water from twi = 34.9◦ C to two = 25.7◦ C. Outlet static pressure in the height Ht is calculated using [13] 3.5(1+xi )(1−xi /(xi +0.62198)) Ht . pout = p0 1 − 0.00975 T0 (60) Total pressure at the tower inlet was recalculated using previous equation to cooling tower inlet height 10m. In ﬁgures 5, 6, 7, 8 and 9 the distribution of chosen parameters are shown over the tower height. The pressure decay in the ﬁgure 5 is aﬀected not only by gravity but also by losses in the rain and ﬁll zones. There is a drop in the density distribution in the ﬁgure 6 which is connected with heat addition in the ﬁll zone. The heat addition is affecting also the distribution of temperature in the ﬁgure 7 where the sharp increase of temperature in the ﬁll zone is observed. The mass source in the ﬁll zone leads to jump in the distribution of speciﬁc humidity as shown in the ﬁgure 8. The decrease of temperature and humidity in the upper part of the cooling tower is natural and is connected with gravitational force acting on ﬂowing moist air. The eﬀect of heat addition in the ﬁll zone is perceptible also in the distribution of velocity. There is slight increase of velocity in the ﬁll zone, see ﬁgure 9. The increase of velocity in ﬁll zone have to compensate the mentioned drop in the density distribution. where Ai+1/2 + Ai−1/2 . 2 Jameson’s artiﬁcial diﬀusion AD(Win ) [10] is n n AD(Win ) = Cγ Wi−1 , − 2Win + Wi+1 μi = Δz where γ= |pi−1 − 2pi + pi+1 | . |pi−1 | + 2|pi | + |pi+1 | (55) (56) (57) Vector of conservative variables is computed in new time level as (58) Win+1 = Win+1 + AD(Win ). 8 Conclusions The one dimensional model of natural draft cooling tower ﬂow is developed and numerically solved. The key part of the model is calculation of density of heat and mass sources which is performed using [1] and recalculated for the utilization of developed model of the ﬂow of homogeneous mixture of dry air and water vapour. Numerically obtained results are relevant and they correspond to general expectations of the behavior inside natural draft cooling tower. Author is aware of imperfections of proposed model especially omission of liquid droplet in the spray 02044-p.5 EPJ Web of Conferences 98000 1.18 97800 ρ / kg/m3 1.17 p / Pa 97600 97400 97200 97000 1.16 1.15 1.14 96800 1.13 96600 1.12 96400 96200 0 20 40 60 80 z / m 100 120 1.11 140 Fig. 5. Static pressure distribution 0 20 40 60 80 z / m 100 120 140 Fig. 6. Density of mixture distribution 25 24 23 t / ◦C and rain zone, where additional heat and mass transfer occurs. Last but not least it is necessary to mention the limit of one dimensional approach. Author is going to overcome imperfections and limits of this model in future work. 22 21 20 19 9 Acknowledgement 18 17 This work has been supported by Technology Agency of the Czech Republic under the project Advanced Technologies for Heat and Electricity Production - TE01020036. 16 15 0 20 40 60 80 z / m 100 120 140 Fig. 7. Distribution of temperature References 22 x / g/kg 20 18 16 14 12 10 8 6 0 20 40 60 80 z / m 100 120 140 Fig. 8. Distribution of speciﬁc humidity 3.5 3 v / m/s 1. A. Klimanek, R. A. Białecky, International Communications in Heat and Mass Transfer 36, (2009), 547-553 2. T. Hyhlı́k, Engineering Mechanics 2013, 19th International conference, Svratka, (2013) 3. F. Merkel, Verdunstungskühlung, VDI Forschungsarbeiten , no. 275, VDI - Verlag, Berlin, 1925 4. N. J. Williamson, Numerical Modelling of Heat and Mass Transfer and Optimisation of a Natural Draft Wet Cooling Tower. Ph.D. thesis, University of Sydney, 2008 5. A. K. Majumdar, A. K. Singhal, D. B. Spalding, ASME J. Heat Transfer, 105(4), (1983) 6. N. Williamson, M. Behnia, S. Armﬁeld, Int. Journal of Heat and Mass Transfer, 51 (2008) 7. F. Maršı́k, Flow of Heterogeneous Mixtures with Relaxation Processes (in czech), Institute of Thermomechanics, Academy of Sciences of the Czech Republic, Prague, 1994 8. A. Bejan, Convection Heat Transfer, John Wiley & Sons, 2004 9. F. Bosnjakovic, P. L. Blackshear Jr. (Ed.), Technical Thermodynamics, Holt, Rinehart and Winston, New York, 1965 10. P. Punčochářová-Pořı́zková, K. Kozel, J. Horáček, J. Fürst, Engineering Mechanics, 17(2), (2010) 11. J. Halama, F. Benkhaldoun, J. Fořt, Int. Journal for Numerical Methods in Fluids, 65(8), (2011) 12. I. Zúñiga-González, Modelling Heat and Mass Transfer in Cooling Towers with Natural Convection. Ph.D. thesis, Czech Technical University in Prague, 2005 13. D. G. Kröger, Air-Cooled Heat Exchangers and Cooling Towers. Penn Well Corporation, Tulsa, 2004 02044-p.6 2.5 2 1.5 1 0 20 40 60 80 z / m 100 120 Fig. 9. Velocity magnitude distribution 140

© Copyright 2021 Paperzz