Solution of heat and mass transfer in counterflow wet

International Communications in Heat and Mass Transfer 36 (2009) 547–553

Contents lists available at ScienceDirect

International Communications in Heat and Mass Transfer j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / i c h m t

Solution of heat and mass transfer in counterflow wet-cooling tower fills☆ A. Klimanek ⁎, R.A. Białecki Institute of Thermal Technology, Silesian University of Technology, Konarskiego 22, 44-100 Gliwice, Poland

a r t i c l e

i n f o

Available online 17 April 2009 Keywords: Evaporative heat mass transfer Counterflow heat mass exchangers Wet-cooling towers

a b s t r a c t Model of heat and mass transfer in wet cooling tower fills is presented. The model consist of a set of four 1D ODEs describing the mass and energy conservation and kinetics with boundary conditions prescribed on opposite sides of the computational domain. Shooting technique with self adaptive Runge–Kutta step control is applied to solve the resulting model equations. The developed model is designed to be included in a large scale CFD calculations of a natural draft cooling tower where the fill is treated as a porous medium with prescribed distributions of mass and heat sources. Thus, the technique yields the spatial distributions of all flow parameters, specifically the heat and mass sources. Such distributions are not directly available in standard techniques such as Merkel, Poppe and e-NTU models of the fill where the temperature of the water is used as an independent variable. The method is validated against benchmark data available in the literature. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Wet cooling towers are huge structures whose aim is to cool industrial water. The cooling effect comes from heat convection and partial evaporation of the water into the air. Increase of the humidity of the air, and hence its density, induces buoyancy driven movement of the air. The evaporation takes place mainly in the fill, being a extended surface heat and mass exchanger, where thin film of the water is in contact with the air. CFD simulation is a well established technique of simulating coupled mass, heat and momentum transfer processes. The idea of using this technique in the context of cooling tower is not new [1–4]. The difficulties associated with the application of CFD to simulate cooling towers come mainly from the different geometrical scales arising in these problems. The height of cooling towers often exceeds hundred meters. To account for the interaction with the surrounding air, the dimension of the CFD model should be of the order of a kilometer. The fill, where around 90% of the heat transfer takes place occupies a very small fraction of the volume of the tower. Working with an exact geometry of the fill, whose channels are of the order of centimeters, would lead to prohibitive storage and execution time. A natural tendency is thus to treat the fill as a porous medium. This semi-continuous model frees the modeler from reproducing the detailed geometry of the elements of the fill. This in turn opens the possibility of using cells much larger than the thickness of the water film and channels of the fill [4,5].

☆ Communicated by W.J. Minkowycz. ⁎ Corresponding author. E-mail addresses: [email protected] (A. Klimanek), [email protected] (R.A. Białecki). 0735-1933/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.icheatmasstransfer.2009.03.007

The presence of evaporation and convective heat transfer in the fill should then be accounted for by defining a distribution of mass (water vapor), heat and momentum sources in the porous medium. The momentum source is readily determined by measuring the pressure drop on the fill. To get the distribution of the heat and mass sources a model of the mass and heat transfer in the fill should be used. As the vertical flow in the fill strongly dominates the transverse one, transport of mass and heat between the water film and the air can be described by a 1D model. As a result, the widely used models of the fill based on 1D countercurrent vertical plug flow approximation of both the water and the air. The Merkel [6], e-NTU [7] and the most accurate Poppe [8] methods are the most popular models of fill. Their comparison and rigorous derivation is presented in detail in [9] and [10] and will not be repeated here. The models are formulated as sets of differential equations and consist of mass and energy conservation equations accompanied with kinetics of mass and energy transfer. The oldest of these techniques, the Merkel method, simplifies the model by neglecting some less important features of the phenomena. The result is just one ordinary differential equation that can be integrated by hand. The youngest and most accurate Poppe's method solves the complete set of equations. The e-NTU method is based on the same critical assumptions as the Merkel method. Due to its simplicity, the eNTU is especially useful in the solution of the crossflow cooling tower fills [11] where the heat and mass transfer problem is described by partial differential equations. All these techniques have been devised to be implemented in 1D models of the cooling towers. The solution is obtained by following an iterative procedure to find the flow parameters that satisfy both the energy and the draft equations. Such models are, by definition, not able to reproduce the radial non-uniformities in the flow parameters.

548

A. Klimanek, R.A. Białecki / International Communications in Heat and Mass Transfer 36 (2009) 547–553

Thus, the influence of wind onto the operation of the tower, switching off some zones of the fill cannot be accounted for. The later mode of operation is often implemented in the wintertime to prevent the water from freezing. A common feature of all standard fill models is that the independent variable of the ODEs is the temperature of the cooled water. This makes the models suitable for scaling calculations, where the temperature drop of the water is known. Generally, implementation of the standard formulations in multidimensional CFD models of the tower are difficult. First reason for this is that keeping the temperature drop of the cooling water constant, would require an adjustment of the fill height. This in turn would change the geometry of the computational domain. A numerical implementation of such option is very cumbersome. As already explained, to treat the fill as a porous medium, the spatial distributions of the mass and heat sources along the vertical coordinate of the fill should be known. Extracting these data from the existing fill models is difficult. The developed technique differs from the known approaches which do not use the natural, vertical coordinate of the fill as an independent variable. Thus, the distribution of the heat and mass sources can readily be evaluated. The set of ODEs is solved using a self adaptive integration technique, where the spatial step is determined automatically in such a way that the error of integration is kept in the prescribed interval. The boundary conditions needed to solve the system are readily imported from the CFD code. The complete set of boundary conditions necessary to solve the fill equations consist of the inlet flowrate and temperature of the water and air and the inlet air humidity. The problem belongs to the class of two point boundary problems, as some boundary conditions (inlet temperature and flowrate of the water) are known at the top of the fill while the inlet temperature and humidity of the air is known at the bottom of the fill. Such problems require iterative integration of the set of ODEs. The solver implemented in this paper is the shooting technique. 2. Heat and mass transfer in the counterflow wet-cooling tower fill Equations governing the heat and mass transfer in the fill link four dependent variables: temperature of the water Tw, temperature of the air Ta, humidity ratio of the air X, mass flowrate of water mw. The independent variable is the vertical position in the fill z. Vertical countercurrent plug flow of water film and air is assumed. The air at the water air interface has the temperature of the water. The bulk temperature of the air is Ta. The formulation of the governing equations depends on saturation level reached by the air. The unsaturated case is considered first. 2.1. Governing equations for unsaturated air The mass and energy balances on a differential slice dz of the fill perpendicular to the air and water flow take the form dmw dX = ma dz dz w

cw T w

dmw dh w dT + mw cw w = ma a dz dz dz

ð1Þ

ð4Þ

dA = aAz dz

where a is the transfer area per unit volume and Az is the cross sectional area of the fill, perpendicular to the air and water flow. The equation of kinetics of mass evaporation reads w dA w dmw = β Xs − X aAz = β Xs − X dz dz

ð5Þ

where β is the average mass transfer coefficient, Xw s is the humidity ratio of air at the air water interface. The air at this location is in saturated state and has the temperature of the water Tw. X is the humidity ratio of the air at the bulk air temperature. It is assumed that the absolute pressure p in the fill is constant (pressure drop is negligible), therefore for a given pressure the saturation humidity Xas depends on air temperature only a

Xs = 0:622

ps ðTa Þ p − ps ðTa Þ

ð6Þ

where ps is the saturation pressure of steam while p stands for the total pressure. The energy transferred from the water to the air due to evaporation of mass and heat convection is balanced with the increase of the enthalpy of the air ma

w dha = he β Xs − X + α ðTw − Ta Þ aAz dz

ð7Þ

where ma is the flowrate of the dry air, α is the average heat transfer coefficient and he is the enthalpy of the water evaporating to the air. This quantity is evaluated at the bulk water temperature Tw and can be expressed as w

he = r0 + cpv Tw

ð8Þ

where cpv is the specific heat of water vapor. Eqs. (1), (2), (5) and (7) comprise a set of four ordinary differential equations. The calorific quantities arising in the equations, i.e. cpv, cpa and r0 are well known temperature dependent quantities. Obtaining the mass β and heat transfer α coefficients is cumbersome as it requires time consuming and expensive experiments. The dimensionless measure of the relative rates of heat and mass transfer in an evaporative process is the Lewis factor Lef. This number is defined as the ratio of the heat transfer Stanton number St to the mass transfer Stanton number Stm. For air–water vapor system it takes the form Lef =

St α = a a St m + Xcpv β cpa

ð9Þ

ð2Þ

where ha represents the enthalpy of the humid air per unit mass of dry air, which for unsaturated air can be expressed as a a ha = cpa Ta + X r0 + cpv Ta

quantity is evaluated at water w or bulk air a temperature. Analogous notation is used to distinguish the saturation humidity ratios determined at water and bulk air temperature. The conservation equations are accompanied with two relationships describing the kinetics of evaporation through the air–water interface of the water film. The area of the interface is

ð3Þ

with r0 standing for the specific enthalpy of evaporation at Tw = 0 °C. The superscript attached to the specific heats indicates whether the

If the value of the Lewis factor is known, only one of the pair of coefficients: the mass or heat transfer coefficient needs to be found experimentally. Bosnjakovic [12] published a formula defining the Lewis factor for air–water vapor systems as a function of the air humidity and saturation humidity ratio.

2=3

Lef = 0:866

w

Xs + 0:622 X w + 0:622 − 1 − 1 ln s X + 0:622 X + 0:622

ð10Þ


In the case of supersaturated air the Lewis factor is evaluated as

Substituting Eq. (5) into mass balance Eq. (1) produces βaAz Xsw − X dX = dz ma

2=3

ð11Þ

Differentiation of Eq. (3) yields the spatial derivative of the enthalpy of air dT dX dha a a a a = cpa + Xcpv + r0 + cpv Ta dz dz dz

549

ð12Þ

Lef = 0:866

w

Xs + 0:622 X w + 0:622 − 1 − 1 ln sa a Xs + 0:622 Xs + 0:622

In order to take into account the mist in the supersaturated air, the specific heat of humid air in Eq. (9) should be modified so that the Lewis factor is defined as Lef =

α β capa + Xsa capv + caw ð X − Xsa Þ

ð19Þ

Simple algebra analogous to those used in the unsaturated air brings Eqs. (15)–(19) to a set of four ordinary differential equations

Combining Eqs. (12), (7), (4) and (11), yields h i dTa a a w a = βaAz Lef ðTw − Ta Þ cpa + cpv X + cpv Tw − cpv Ta ðXs − X Þ dz h i a a ma cpa + cpv X ð13Þ

Introducing Eqs. (7) and (5) into the energy balance (2) the spatial derivative of the water temperature can be calculated from h i dTw a a w w = βaAz Lef ðTw − Ta Þ cpa + cpv X + r0 + cpv Tw − cw Tw ðXs − X Þ dz w mw cw ð14Þ

w dmw a = βaAz Xs − Xs dz

ð20Þ

w a βaAz Xs − Xs dX = dz ma

ð21Þ

dTa βaAz a w w =− cpa Lef ðTa − Tw Þ − Xs r0 + cpv Tw dz ma w a a a + cw Lef ðTa − Tw Þ X − Xs + Ta Xs − Xs a a w + Xs r0 + cpv Lef ðTa − Tw Þ + cpv Tw

dXsa a a a a a a a cpa + cw X + r0 + cpv Ta − cw Ta + Xs cpv − cw dTa

½

ð22Þ

2.2. Governing equations for supersaturated air The water content of the air gradually rises so that finally supersaturated air containing droplets of mist is produced. The water in contact with the supersaturated air is further evaporated and cooled. The reason for this is that the evaporation process is driven by the difference in the concentration of the water vapor on the air–water interface and in the bulk air. When the air is supersaturated the bulk air humidity ratio is the humidity of the saturated air evaluated at the bulk air temperature Xas and the force driving the a evaporation process is (Xw s − Xs ). The water will be cooled due to evaporation until the driving force is positive, i.e. until the local water temperature is higher than the bulk air temperature. The description of the process has to be modified in order to capture the physics of the supersaturated air. Eq. (5) is modified to w dmw a = β Xs − Xs aAz dz

ð15Þ

The enthalpy of the air in Eq. (3) is, for supersaturated air evaluated from a a a a a hast = cpa Ta + Xs r0 + cpv Ta + X − Xs cw Ta

ð16Þ

and its spatial derivative is

dX a dhast a a a a a a a dTa s = r0 + cpv Ta − cw Ta + Xs cpv − cw + cpa + Xcw dz dTa dz a

+ cw Ta

dX dz

ð18Þ

ð17Þ

The derivative dXas /dTa is a known function of saturation temperature. The explicit expression can be obtained by expressing in (6) the saturation pressure as a function of temperature and differentiating the result.

dTw w w w a = βaAz r0 + cpv Tw − cw Tw Xs − Xs dz a a a a a w + Lef ðTw − Ta Þ cpa + cw X − Xs + cpv Xs = cw mw

½

ð23Þ

3. Solution of the governing equations Under the assumptions made when deriving the equations the mass and heat transfer problem in the fill is described by a set of four ordinary differential equations explicit in the derivatives of the unknown functions. The problem is well defined, as the number of boundary conditions is equal to the number of unknowns. Should all the boundary condition be defined at one point (bottom or top of the fill), distribution of the unknown functions: Ta, Tw, X and mw could be obtained by any standard ODE solver. However, as already mentioned, the problem at hand belongs to the class of two point boundary value problems, for which some of the boundary conditions are known at one boundary point while the remaining at the other. The boundary conditions known at the bottom of the fill (z = 0) are the air temperature Ta and humidity ratio X. At the top of the fill (z = H) the water temperature Tw and mass flowrate mw are known. In order to integrate the equations, the missing boundary conditions at the bottom of the fill need to be guessed and iteratively adjusted so that the known boundary conditions at the top of the fill are satisfied. The shooting method [13] is used to solve two point boundary value problem applying multidimensional globally convergent Newton–Raphson method. The differential equations are integrated using the adaptive step size control Runge–Kutta technique. The solver automatically adjusts the integration step in such a way that the error is kept within a prescribed interval. Before the governing equations are solved, the value of the heat α or mass β transfer coefficient has to be determined. The only feasible method of getting this quantities is an experiment. An elegant and robust technique of retrieving the coefficient from the measurements

550


Fig. 3. Water mass flowrate distribution in the fill.

Fig. 1. Air and water temperature distributions in the fill.

is the inverse analysis. This approach is, however, seldom used in the context of the cooling towers. Traditionally, the mass transfer coefficient is included in Merkel's number defined as Z Me =

βaAz dz mw

ð24Þ

which is determined by solving the governing equations with the experimentally determined water outlet temperature. The Merkel number for various flow conditions is usually expressed as a function of air and water mass flowrates. Sometimes more complex empirical expressions are employed to take into account the influence of inlet water temperature and fill height. The mass transfer coefficient β or the product βaA z can be calculated from the experimentally determined Merkel number. This, however, needs to be done iteratively since the integral on the right hand side of Eq. (24) is not known a priori. The initial guess for the iterative procedure can be ðβaAz Þj = Me

mwi H

ð25Þ

that the left hand side of Eq. (24) can be determined by appending to the integrated set a relationship ðβaAz Þj dMej = dz mw

ð26Þ

The product (βaAz)j is adjusted until the obtained Mej number approaches the experimental value Me. Since the water mass flowrate is not changing much, the initial guess (25) is a good approximation and the convergence is obtained after few iterations. 3.1. Numerical implementation Integration of the set from z = 0 to z = H yields the values of the four dependent variables at each intermediate integration step. The initial values of the unknown boundary conditions at z = 0 are selected as follows. The water outlet temperature is taken from experimental approximation given in [14] Two =

Twi + 2Twb + Tai 4

ð27Þ

where Me is the Merkel number known from experiments. The governing equations are solved using the initial guess of (βaAz)j, so

where Twb is the wetbulb temperature of the inlet air and the indices i and o indicate the inlet and outlet parameters, respectively. This

Fig. 2. Distributions of humidity and saturation humidity ratios in the fill.

Fig. 4. Merkel number distribution in the fill.


Fig. 5. Distributions of heat and mass sources in the fill.

Fig. 6. Air and water temperature distributions.

approximation is used for high inlet air temperatures only. For low air temperatures, the outlet water temperature is calculated using the energy balance (2) assuming that dmw T + Twi ao ≈0; Tao ≈ ai ; Xo ≈Xs dz 2

ð28Þ

where Xao is the saturation humidity ratio evaluated at the s approximated air outlet temperature Tao. The outlet water temperature can be then calculated as Two = Twi −

ma ðh − hai Þ mw cw ao

ð29Þ

The outlet water mass flowrate can be estimated as ao mwo = mwi − ma Xs − Xi

ð30Þ

At every integration step the integration procedure tests whether the air is unsaturated or supersaturated. Depending on the condition (X N Xas ) either equations for unsaturated or supersaturated air are invoked. As the integration is adaptive, the values of the unknown functions are known at some not evenly distributed points. In order to obtain the values of the dependent variables at any desired position z, the original results are interpolated using cubic splines.

H = 1.2 m. The cross-sectional area of the fill is taken Az = 1 m2. The air and water mass flowrates are chosen to be equal to ma =mw = 3.0 kg/s. The Merkel number for this flowrates is Me = 1.8613. In the computations performed in this study, the product βaAz is determined from the Merkel number using the iterative procedure described in section 3. The air inlet temperature is Tai = 30.0 °C, inlet humidity ratio Xi = 0.00262 kg/kg dry air and the inlet water temperature is Twi = 37.0 °C. The results of the analysis are presented in Figs. 1–5 where the nodes indicate the values produced by the adaptive step size control integration technique (calc.) as well as the spline interpolation (int.). As can be seen in Fig. 1 the temperature profiles cross. The temperature of the air flowing upwards is decreasing in the bottom part of the fill and starts to increase in the upper part. In some situations the air leaving the fill can even be cooler than the air entering the fill. The temperature of the water steadily decreases so that the cooling process continues even over the crossing point of the air and water temperatures. This is possible because the inlet air is very dry so that the evaporative heat transfer dominates over the convective one. Fig. 2 shows the humidity and saturation humidity distributions. The curves do not cross which indicates that the air does not become supersaturated and only the first set of equations is used. The

3.2. Numerical example The numerical example of fill analysis presented below shows a particular situation occurring in wet cooling towers under hot and very dry air conditions. The analysis was performed for a fill of height

Table 1 Comparison of the results obtained with Poppe method by Kloppers [14] and presented approach.

Two/°C Me Tao/°C Xo/kg kg− 1 Xas /kg kg− 1 mwo/kg s− 1 hao/kJ kg− 1

551

Kloppers [14]

Proposed approach

Relative error %

21.41 1.5548 26.71 0.02789 0.02718 12170.5 96.303

21.41 1.5548 26.70 0.02767 0.02718 12174.2 96.301

0.00 0.00 0.02 0.80 0.01 0.03 0.00

Fig. 7. Distributions of humidity and saturation humidity ratios.

552


distribution of the mass flowrate of the water is plotted in Fig. 3. Merkel's number distribution is depicted in Fig. 4. The sources of mass and heat required for the CFD simulation and calculated from the spatial distribution of flow variables are presented in Fig. 5. Since the air is cooled in the bottom part of the fill, the sources are negative in that part. The mass sources decrease at the bottom part of the fill and increase while approaching the top.

4. Validation of the method The presented method is validated against benchmark results acquired from [14]. The results are obtained by application of the Poppe method to the analysis of a full scale cooling tower. The Merkel number in this case is the sum of the Merkel numbers of all three heat and mass transfer regions, namely the spray, fill and rain zones. The number amounted Me = 1.5548. The inlet air temperature is Tai = 15.45 °C, inlet humidity ratio Xi = 0.008127 kg/kg dry air, air mass flowrate ma = 16672.19 kg/s, water inlet temperature Twi = 40.0 °C and water inlet mass flowrate mwi = 12,500.0 kg/s. The ambient pressure was p = 84100 Pa. In the calculations performed by Kloppers [14] the outlet water temperature is known Two = 21.41 °C and the Merkel number according to Poppe theory is calculated. In the computations performed in this study, the product βaAz is determined from this Merkel number using the iterative procedure described in Section 3. As the result, the calculated water temperature is determined and compared with the known value. The results of the comparison are presented in Table 1. As can be seen results obtained from both methods are in good agreement with each other. The slight discrepancy can be due to the application different integration techniques. Another source of discrepancy may be different functions representing temperature dependence of specific heats and saturation pressure. The spatial distributions of the dependent variables for this example are shown in Figs. 6–9. As can be seen in Fig. 7 the humidity and saturation humidity ratios cross somewhere close to z = 2 m. The crossing point indicates that the air became supersaturated and the equations for supersaturated air are applied. It is visible that the adaptive step size control integration technique requires more steps in the vicinity of the saturation region increasing the accuracy in that region, whereas the steps can be much longer in other parts of the integration space. Presented validation of heat and mass transfer model conducted on fill computations and the entire cooling tower analysis show good agreement between data obtained with the Poppe

Fig. 9. Merkel number distribution in the fill.

method and the proposed approach. This confirms that the proposed approach is properly implemented and equivalent to the Poppe analysis. 5. Conclusions The proposed method of solution of heat and mass transfer process in wet cooling tower fills was shown to be equivalent to the Poppe method of analysis, known as an accurate computational technique. The comparison of both methods showed very good agreement between results in the validation computations. The applied adaptive step size control integration technique adjusts the step size in the regions of steep gradients and therefore it should be a method of choice as it increases the accuracy of the integration with a small additional computational effort. The proposed approach gives insight to the physics of the process revealing the spatial distributions of the dependent variables. The distributions may be required when multidimensional CFD codes are used to analyze the cooling towers. In such situation the spatial distributions of heat and mass sources are required to be introduced into the continuity, species convection– diffusion and energy equations by projection on the numerical grid of the CFD code [4,5]. Acknowledgement The work of Adam Klimanek was financed by the Polish Ministry of Science and Higher Education as a research grant in years 2007–2010. The contribution is fully acknowledged. References

Fig. 8. Water mass flowrate distribution in the fill.

[1] N. Williamson, S. Armfield, M. Behnia, Numerical simulation of flow in a natural draft wet cooling tower — the effect of radial thermofluid fields, Applied Thermal Engineering 28 (2008) 178–189. [2] D. Radosavljevic, The numerical simulation of direct-contact natural-draught cooling tower performance under the influence of cross-wind, PhD thesis, University of London, London, UK, 1990. [3] M.N.A. Hawlader, B.M. Liu, Numerical study of the thermal-hydraulic performance of evaporative natural draft cooling towers, Applied Thermal Engineering 22 (1) (2002) 41–59. [4] A. Klimanek, R.A. Białecki, On a numerical model of a natural draught wet-cooling tower, Archives of Thermodynamics 29 (4) (2008) 63–72. [5] A. Klimanek, R.A. Białecki, Modelling of flow, heat and mass transfer in natural draught wet-cooling towers, Proceedings of XVI Mdzinarodna Vdecka Konferencja — Aplikacia Experimentalnych a Numerickych Metod v Mechanike Tekucin, vol. 1, Zilina-Terchova, Slovakia, 2008, pp. 73–78. [6] F. Merkel, Verdunstungskühlung, 70, VDI-Zeitschrift, 1925, pp. 123–128. [7] H. Jaber, R.L. Webb, Design of cooling towers by the effectiveness-NTU method, Journal of Heat Transfer 111 (1989) 837–843.

A. Klimanek, R.A. Białecki / International Communications in Heat and Mass Transfer 36 (2009) 547–553 [8] M. Poppe, H. Rögener, Berechnung von Rückkühlwerken, VDI-Wärmeatlas, 1991, pp. Mi 1–Mi 15. [9] J.C. Kloppers, D.G. Kröger, A critical investigation into the heat and mass transfer analysis of counterflow wet-cooling towers, International Journal of Heat and Mass Transfer 48 (2005) 765–777. [10] D.G. Kröger, Air-cooled Heat Exchangers and Cooling Towers, PennWell Corporation, Tulsa, Oklahoma, 2004. [11] J.C. Kloppers, D.G. Kröger, A critical investigation into the heat and mass transfer analysis of crossflow wet-cooling towers, Numerical Heat Transfer, Part A 46 (2004) 785–806.

553

[12] F. Bosnjakovic, in: Perry L. Blackshear Jr. (Ed.), Technical Thermodynamics, Holt, Rinehart and Winston, New York, 1965. [13] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, Cambridge University Press, New York, 1992. [14] J.C. Kloppers, A critical evaluation and refinement of the performance prediction of wet-cooling towers, Doctoral thesis, Department of Mechanical Engineering, University of Stellenbosch, South Africa, 2003.