Concept of CFD model of natural draft wet-cooling tower flow

The article deals with the development of CFD model of natural draft wet-cooling tower flow. 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 counterflow wet-cooling tower fill are described by model [1] which is based on the system of ordinary differential equations. Utilization of model [1] of the fill allows us to apply commonly measured fill characteristics as shown by [2].The boundary value problem resulting from the fill model is solved separately. The system of conservation law equations is interlinked with the system of ordinary differential equations describing the phenomena occurring in the counterflow wet-cooling tower fill 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 flow. The simulation results are shown.


Introduction
Natural draft cooling towers are widely used in industry especially in large power plants and in chemical industry.Low efficiency 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 efficiency of cooling process.
The aim of cooling tower is to transfer waste heat from cooled water into the atmosphere.The air which is flowing through the cooling tower is warmed up and its humidity is increasing which leads to cooling of flowing water.The density of warmed air is decreasing unlike the surrounding air and this density difference 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 fill using the set of sprayers.In the channels of counterflow wet-cooling tower fill of film type water flows vertically down through the fill as a liquid film.Air is driven by a tower draft and flows 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 fill 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 fill zone.The simplest model of heat and mass transfer in the fill 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 flow can be mentioned e.g.[5] and commercial code based [6].a e-mail: tomas.hyhlik@fs.cvut.cz

Balance laws
The moist air flowing through the cooling tower will be considered as homogeneous mixture.The total value of extensive quantity ϕ is expressed as where w α is the mass fraction of component α.Local form of the balance law is [7] ∂ Diffusive flux density j k Dα is expressed through the center of mass velocity where Diffusive flux j k Dα can be written using Fick's law [8] as where D α,m is diffusion coefficient.Density of flux of quantity ϕ we express as Density of production of ϕ is Local form of balance law for quantity ϕ can be written for the mixture as

Moist air
The moist air is mixture of dry air and superheated or saturated steam.For a description of the cooling tower flow, 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 flow.The density of the moist air is where the densities can be expressed using equation of state as Mass fraction of dry air and mass fraction of water vapour we can express as where Pressure can be expressed using Dalton's law Specific gas constant of mixture can be expressed by means of specific gas constants of both components Specific heat capacity of the mixture is expressed as Specific moisture indicates the weight of steam per kilogram of dry air Mass fractions can be expressed using specific moisture Equation of state of a mixture of dry air and water vapour is expressed using Dalton's law, the definition of specific humidity and the definition of the overall density as where M a is molecular weight of dry air, M v 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 The total energy of moist air is expressed as the sum of internal energy and kinetic energy 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 where t a = 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 This definition 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 define ad hoc internal energy of unsaturated humid air as The total energy of unsaturated moist air is expressed in the form where κ is Poisson coefficient

Governing equations
Continuity equation for dry air can be written as where D a,v is diffusion coefficient of dry air in the water vapour.We can express the continuity equation of water vapour in the form where σ v (x k , t) represents the density of source of water vapour and diffusion coefficient of water vapour in the air is D v,a = D a,v .The system of momentum equations can be written in the form where −ρgδ i3 represents the gravitational force acting on the flowing moist air, σ ik is the stress tensor in a flowing fluid and ζ represents a loss coefficient per meter.The energy equation is written in the form where σ q (x k , t) is density of heat source where only sensible part of heat source is considered because of the definition of total energy which is based on equation (25).

Quasi one-dimensional flow equations
The modification of governing equations for the case of quasi one-dimensional flow is performed for control volume A(z)dz.Dry air continuity equation has form 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 flow of the form where and vector of sources is The system of governing equations can be modified by summing of individual continuity equations to get mixture continuity equation.The water vapour continuity equation is modified to get water vapour mass fraction equation.
and modified vector of sources is

Computation of sources in fill zone
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 (x k , t) and density of heat source σ q (x k , t).The sources in the fill zone can be calculated using one dimensional model [1].

Balance laws
Due to the complexity of two phase flow occurring in the wet-cooling tower fill the one dimensional models of heat and mass transfer are used.These models are based on few assumptions which allow to create simplified one dimensional models.The first assumption talks about neglecting of the effects of horizontal temperature gradient in the liquid film, horizontal temperature gradient in air temperature 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 film 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 fill is based on balance laws.We have four variables in this problem: t a temperature of air, t w temperature of water, x specific humidity and ṁw water mass flow rate.Mass balance of the incremental step of the fill is given by d ṁw where ṁa is mass flow rate of dry air.The change in water mass flow rate can be expressed using mass transfer coefficient α m as where x (t w ) is saturated specific humidity at t w and dA is infinitesimal contact area.The energy balance can be written in the form where h 1+x is enthalpy of air water vapour mixture and h w is enthalpy of water.The change in total enthalpy can be evaluated using interface parameters similarly like in the case of mass balance where α is heat transfer coefficient and h v (t w ) is enthalpy of water vapour.

Model of Klimanek & Białecky [1]
To derive the system of ordinary differential 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 t w .The interface area can be expressed using variable z as dA = a q A f r dz, where a q is the transfer area per unit volume and A f r is the cross sectional area of the fill.We can derive equation for the change of the water mass flow rate from equation (38) To obtain the equation for the change of specific humidity in the fill we can substitute equation (42) into equation ( 37) we can get By using equation (39) after substitution of equation ( 42) and equation ( 40) we derive equation for the change of water temperature The application of the Lewis factor in the previous equations simplified the problem to find experimentally only mass transfer coefficient α m and calculate heat transfer coefficient α 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 (t w ) − x (t a ))) in the case of supersaturation.The system of ordinary differential equations has to be changed [1].The enthalpy of supersaturated air is 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 differential equations mentioned in this section represents the boundary value problem.We know air temperature t ai and specific humidity x i at air inlet and water temperature t wi and water mass flow rate ṁi on the opposite site of the fill because air and water flows in the opposite direction.There is additional unknown set of parameters in the system of equations, i.e. α m a q A f r .These parameters have to be solved by using experimentally obtained characteristics of the fill [1,2].

Fill zone simulation example
The calculation was performed for a fill of height H = 1 m.The air inlet mass flow rate is equal to ṁa = 14333 kg/s and water inlet mass flow rate is ṁw = 17200 kg/s.Inlet water temperature is t wi = 34.9• C. Air inlet temperature is t ai = 15.7 • C and specific humidity at air inlet is x i = 7.622 g/kg.The atmospheric pressure is p = 98100 Pa.Outlet water temperature is t wo = 25.7 • C.
The numerical solution of the fill zone model was performed using Runge-Kutta method combined with shooting method [2].The calculated distributions of temperature and humidity are depicted in figures 1 and 2. From the  mentioned figures is possible to observe supersaturation in the upper part of of the fill, it means the air temperature is equal to adiabatic saturation temperature and specific humidity is slightly higher then saturation humidity.The distribution of the density of mass source is shown in the figure 3 and the density of sensible part of heat source is in the figure 4, where density of mass source of water vapour can be expressed as EPJ Web of Conferences 02044-p.4

Cooling tower flow computation
Numerical solution is based on the (32), where because of the definition of total energy (25) the density of heat source should be calculated as

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 and corrector step is where Jameson's artificial diffusion AD(W n i ) [10] is where Vector of conservative variables is computed in new time level as Total pressure p 0 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.

Calculation example
Natural draft cooling tower of rotational hyperboloid shape with is selected [12]. .
Total pressure at the tower inlet was recalculated using previous equation to cooling tower inlet height 10m.In figures 5, 6, 7, 8 and 9 the distribution of chosen parameters are shown over the tower height.The pressure decay in the figure 5 is affected not only by gravity but also by losses in the rain and fill zones.There is a drop in the density distribution in the figure 6 which is connected with heat addition in the fill zone.The heat addition is affecting also the distribution of temperature in the figure 7 where the sharp increase of temperature in the fill zone is observed.The mass source in the fill zone leads to jump in the distribution of specific humidity as shown in the figure 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 flowing moist air.The effect of heat addition in the fill zone is perceptible also in the distribution of velocity.There is slight increase of velocity in the fill zone, see figure 9.The increase of velocity in fill zone have to compensate the mentioned drop in the density distribution.

Conclusions
The one dimensional model of natural draft cooling tower flow 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 flow 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 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.

Fig. 2 .Fig. 3 .
Fig.2.Distribution of specific humidity x, saturation specific humidity at air temperature x (t a ) and saturation specific humidity at water temperature x (t w )

Fig. 4 .
Fig. 4. Density of sensible part of heat source distribution

3 Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 6.Density of mixture distribution [13]r height is H t = 150m and fill zone is placed at the height of 11.5m.Calculation was performed for fill height H = 1m.Water inlet mass flow rate is ṁw = 17200 kg/s.Inlet water temperature is t wi = 34.9•C. Air inlet temperature is t ai = 15.7 • C and specific humidity at air inlet is x i = 7.622 g/kg.The atmospheric pressure is p = 98100 Pa.Outlet water temperature is t wo = 25.7 • C. Prescribed parameters allows to quantify fill requirement to cool water from t wi = 34.9•C to t wo = 25.7 • C. Outlet static pressure in the height H t is calculated using[13]