Analytical solution of the governing equations for heat and mass transfer in evaporative cooling process

Journal Pre-proof

Analytical solution of the governing equations for heat and mass transfer in evaporative cooling process ´ ´ , Jesus on ´ R. Ortiz-del-Castillo , Oscar M. Hernandez-Calder ´ Erika Y. Rios-Iribe , Marcos D. Gonzalez-Llanes , Eusiel Rubio-Castro , Maritza E. Cervantes-Gaxiola PII: DOI: Reference:

S0140-7007(19)30495-5 https://doi.org/10.1016/j.ijrefrig.2019.11.019 JIJR 4586

To appear in:

International Journal of Refrigeration

Received date: Revised date: Accepted date:

22 August 2019 21 October 2019 19 November 2019

´ ´ , Please cite this article as: Jesus Oscar M. Hernandez-Calder on ´ R. Ortiz-del-Castillo , ´ Erika Y. Rios-Iribe , Marcos D. Gonzalez-Llanes , Eusiel Rubio-Castro , Maritza E. Cervantes-Gaxiola , Analytical solution of the governing equations for heat and mass transfer in evaporative cooling process, International Journal of Refrigeration (2019), doi: https://doi.org/10.1016/j.ijrefrig.2019.11.019

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights     

An analytical solution for the Poppe’s model is developed using power series. The analytical solution is restricted to cooling processes with unsaturated air. The solution obtained is computationally more efficient. The experimental validation of the analytical solution is presented. The solution could be easily implemented in complex optimization schemes.

1

Analytical solution of the governing equations for heat and mass transfer in evaporative cooling process Jesús R. Ortiz-del-Castilloa, Oscar M. Hernández-Calderóna, Erika Y. Rios-Iribeb, Marcos D. González-Llanesa, Eusiel Rubio-Castroa, Maritza E. Cervantes-Gaxiolaa* a

Chemical Engineering Department, Facultad de

b

ABSTRACT In this work, the governing equations of the heat and mass transfer in evaporative cooling process are solved by using the power series method. The physical properties, including the Lewis factor, are considered constant along the cooling process. The water loss from water stream vaporization is taken into account. An iterative procedure is developed for calculating the expansion coefficients of the humidity ratio, the air enthalpy, and the number of transfer units. In all study cases, the power series solution is convergent for the heat and mass transfer equations, except for the number of transfer units equation. Thus, Gauss quadrature technique is implemented as an alternative method for determining the number of transfer units profile. As a comparison, the study cases are also solved numerically by the Dormand-Prince Runge-Kutta method. The numerical and analytical results are found to be in excellent agreement when the mass flow-rate ratio between water and dry air is low. The computational execution time of the analytical solution is 50 times faster than the numerical solution. Furthermore, the proposed technique was applied to a study case previously reported and the results were properly represented with an average error of 3%.

KEYWORDS: evaporative cooling process, power series solution, mass transfer, heat transfer *

Corresponding author: [email protected] Blvd. de las Américas y Josefa Ortiz de Domínguez s/n, Ciudad Universitaria, Culiacán, Sinaloa, México.

2

1. Introduction The main function of a cooling tower, which is a heat releasing equipment, is to remove waste heat from warm water to the air in the atmosphere. The process of heat release in a cooling tower implicates heat and mass transfer because it consists of convection caused by the contact of water droplets with the surrounding air and it also includes the evaporation of a small quantity water into moving air (Kotb, 2013). Cooling towers can be classified by means of air flow or by the movement of water, as mechanical draft or natural draft types, or as counter-flow or cross-flow types, respectively. In addition, they are regularly used in refrigeration, air conditioning industries, and power generation units (ASHRAE-Handbook, 2008). Besides, the cooling towers are designed by different methods; for example, the Merkel (1925), e-NTU (Jaber and Webb, 1989) and the Poppe methods (Poppe and Rögener, 1991). For any method, it is necessary to have the appropriate knowledge of the governing equations for heat and mass transfer to describe the thermal performance of the cooling tower. Particularly, the cooling tower theory was first given by Merkel (1925), who developed the first practical model for a cooling tower operation, which is obtained by neglecting the water loss of evaporation and assuming the Lewis factor to be one (Khamis Mansour and Hassab, 2014); therefore, Merkel’s model is not accurate enough and not suitable for real applications (Liu and Cai., 2002). Whereas Sutherland (1983) presented a more rigorous analysis of a cooling tower model than Merkel method and showed that designs obtained with the last one method are undersized between 5 and 15%. In addition, e-NTU method (effectiveness- Number of Transfer Units) is one simplified model used commonly for heat exchangers. It was adapted for cooling towers design by Jaber and Webb 3

(1989), and can be applied to both crossflow and counterflow wet cooling towers. This method requires that the heat and mass transfer equations are described by one-dimensional PDEs. Braun et al. (1989) developed an effectiveness model for cooling towers. In this case, they utilized the assumption of a linearized air saturation enthalpy. This model was useful for both process design and simulation. However, the correlation proposed by Braun et al. (1989) needs iterative computation to obtain results, therefore this method does not an option for complex optimization problems (Liu and Cai., 2002). Khamis Mansour and Hassab (2014) reported a simplified correlation for effectiveness-NTU of counterflow cooling tower; which is based on the Merkel theory and it was developed to be used by the process engineer or cooling tower designer, since the calculation of the cooling tower efficiency are done without complex numerical or analytical calculus. Khamis Mansour (2016) developed a new correlation for effectiveness-NTU of counterflow wet cooling coil to predict the thermal performance under both a unit and a non-unit Lewis factors. On the other hand, Poppe and Rögener (1991) developed a rigorous method to evaluate cooling tower performance; they determined the Lewis factor using Bosnjakovic’s formula (Bosnjakovic, 1965), taking into account the water lost by evaporation.

loppers and

r ger (2005a) expressed

the equation for the tower performance calculation using the Poppe method, which was developed for unsaturated and supersaturated air into cooling towers. Heidarinejad et al. (2009) solved numerically the governing equations on wet-cooling towers divided into three zones (spray, packing and rain) and investigated the effect of the Lewis factor on the performance prediction of wet-cooling towers using the Bosnjakovic equation. The results were validated against the experimental measurements and concluded that the combination of three zones in the model increased the consistency between model and prototype resulting in a more accurate solution. Kloppers and Kröger (2005b) analyzed the derivation of heat and mass transfer 4

equations in counterflow wet cooling towers in detail using the Merkel and Poppe methods. Besides. Zheng et al. (2012) used the Poppe method to analyze the cooling tower performance under conditions of saturated and unsaturated air and the governing equations were solved with the finite difference method. In particular, Hernández-Calderón et al. (2014) used the orthogonal collocation technique to solve the Poppe method equations for heat and mass transfer in counterflowing wet-cooling towers. Their methodology was applied to eight study cases and the results were compared with those obtained from the describing equations integrated using the Dormand–Prince method. They concluded that the accuracy was similar between both techniques, however the orthogonal collocation requires less computational time. In addition, Liu et al. (2018) proposed a onedimensional sectional model for a counter-flow cooling towers and they constructed an experimental prototype in order to compare the results, the error were of 3.59% to predict the outlet water temperature and 9.54% to predict the wet-bulb temperature of the outlet air. Alternatively, Llano-Restrepo and Monsalve-Reyes (2017) provided a simpler mathematical derivation of a steady-state one-dimensional continuous differential air-water contactor (CDAWC) model that describes the material and energy balances with the objective to simulate counterflow wet-cooling towers and obtain values of the volumetric mass transfer coefficient from experimental measurements of the cooling tower performance. They demonstrated that the equations of CDAWC model are approximately equivalent to some models previously reported, for instance Poppe Model. On the other hand, in several works (Serna-González et al., 2010; Kintner-Meyer and Emery, 1995; Söylemez, 2001; Söylemez, 2004) the design of cooling towers were optimized economically based on Merkel and e-NTU methods. The results obtained were sub-optimal because the methods considered the Lewis factor as a unit. It is because of in both methods the 5

water lost by evaporation is neglected. Therefore, Rubio-Castro et al. (2011) reported a mixed integer non-linear programming (MINLP) model to obtain the optimal design for counterflow cooling towers based on the Poppe method and considered a detailed geometric design. They point out that a large number of variables and the non-linearities related to the numerical method used to solve the differential equations are some of the disadvantages of using the Poppe method for optimization problems. In addition, Rao and More (2017) proposed a self-adaptive Jaya algorithm for the optimal design of a mechanical draft cooling in six different examples to find the best possible design from an economic point of view. They concluded that the proposed algorithm was superior to the Merkel and Poppe methods; as well as ABD and basic Jaya algorithms in terms of optimal results, convergence and computational time. It should be noted that the rigorous representation of transport phenomena in cooling towers has great practical importance. It is due to this process equipment are too used in the industrial activity. However, this rigorous representation must be efficient and demands low computational efforts in order to design optimal cooling towers. Hence, in this work, the governing equations of the heat and mass transfer in evaporative cooling process are solved by using the power series method with the purpose of using analytical solutions to get faster solutions than in previous works, with a simple computational codification and an excellent representation of the transport phenomena. 2. Mathematical model The equations that describe the evaporative cooling process of the Poppe method are taken and adapted from Poppe and Rögener (1991) and Kröger (2004); the equations are obtained from the mass balance for the control volume (Fig. 1 and Fig. 2). A control volume in the packing of a counterflow wet-cooling tower is shown in Fig. 1, and an air-side control volume of the fill illustrated is shown in Fig. 2. The aforementioned equations are given below: 6

mw *  w  w ma dw  * * dTw ima  ima   Lef  1 ima  ima   w*  w  iv*    w*  w  c p , w Tw  T0 

(1)

 c p , w Tw  T0   w*  w  dima mw    c p,w 1 * * dTw ma  ima  ima   Lef  1 ima  ima   w*  w  iv*    w*  w  c p ,w Tw  T0    

(2)

c p,w

c p,w dNTU  * * dTw ima  ima   Lef  1 ima  ima   w*  w  iv*    w*  w  c p ,w Tw  T0 

(3)

This quotient indicates the relative rates of heat and mass transfer in an evaporative process. When the air is unsaturated, it can be determined by (taken from Kloppers and Kröger, 2005b):

 w*  0.622    w*  0.622   Lef  0.8652/3   1 / ln   w  0.662     w  0.662  

(4)

7

mw  dmw

ma 1  w  dw 

iw  diw

ima  dima

dz

mw

ma 1  w 

iw

ima

Fig. 1. Control volume of the counterflow fill.

mw  dmw

ma 1  w  dw 

iw  diw

ima  dima

dmw

mw

ma 1  w 

iw

ima

dz

dQC  h Tw  Tdb  dA

Fig. 2. Air-side control volume of the fill.

8

3. Analytical solution In this section the analytical solution of evaporative cooling process from series expansion of powers is addressed. This is obtained through the next subsections: the process of series expansion of powers in order to get the solutions of transfer mass and energy equations, the methodology to determine the undetermined coefficients W0 and I 0 and the procedure to calculate the number of transfer units (NTU). Furthermore, the descriptions of considered assumptions are given in Section I of Supplementary Material. 3.1 Power Series Solution of the Equation of Mass Transfer * * * Since the saturated air enthalpy is given by ima  c p,a Tw  T0   w  iv and the water vapour * enthalpy at water temperature is iv  v  c p ,v Tw  T0  , both equations can be combined to give * ima   c p ,a  w*  c p ,v  Tw  T0   w*  v . By applying the variable change Tw   x  

and

*

substituting the power series for w , the saturated air enthalpy can be expressed by the equation



* ima  c p ,a  c p ,v  m0 Wm* x m 

  x    T     0

v

 m0

Wm* x m

(5)

By regrouping terms in Equation (5), it can be written * ima    T0  c p ,a    T0  c p ,vW0*  vW0*

(6)

 c p ,a    T0  c p ,vW1*  c p ,v W0*  vW1*  x 

    T0  c p ,vWm*  c p ,v Wm*1  vWm*  x m m2

* ma

Thus, if i





* m m 0 m

I x then

I 0*    T0  c p ,a    T0  c p ,vW0*  vW0*

(7.a)

I1*  c p ,a    T0  c p ,vW1*  c p ,vW0*  vW1*

(7.b)

I m*    T0  c p ,vWm*  c p ,vWm*1  vWm*

(7.c)

m2

9





*  ima   m0 Am x m then it is necessary that Also, if ima   m0 I m x m and Lef ima 



Am  Lef  I m*  I m 



Regarding the term 1  Lef

  w  w i *

* v

m0

(8)

it is possible to express that

   w iv*    Wm*  Wm  x m  v  c p ,v  x    T0    m0 

w

*

(9)

By expanding the product in Equation (9) and then regrouping in powers of x , it is possible to obtain the following equation:

w

 w  iv*  v    T0  c p ,v  W0*  W0 

*







  v    T0  c p ,v  Wm*  Wm   c p ,v Wm*1  Wm 1  x m m 1



Thus, if 1  Lef



  w  w i   *

* v

 m0

Bm x m , then



B0  1  Lef v    T0  c p ,v  W0*  W0 



(10)

(11.a)





Bm  1  Lef v    T0  c p ,v  Wm*  Wm   c p ,v Wm*1  Wm1 

m 1

(11.b)

On the other hand, by rewriting the term  w*  w Tw  T0  as

w

*

   w  Tw  T0     Wm*  Wm  x m   x    T0   m0 

(12)

By grouping in power of x , Equation (12) can be write by: 

 w  w  T *

* * *   m w  T0     T0  W0  W0      T0  Wm  Wm    Wm 1  Wm 1   x

(13)

m 1

Therefore, if  w*  w c p , w Tw  T0   Ck x m , then it is necessary that the following conditions n

m0

are met

10

C0    T0  c p,w W0*  W0  Cm    T0  c p , w Wm*  Wm    c p ,w Wm*1  Wm1 

(14.a) m 1

(14.b)

Also it is possible to propose that m  mw * w  w   w,in   wout  w    w*  w   ma  ma 

(15)

And since w   m0 Wm x m and W0  wout , is obtained 

 m   mw * w  w    w,in  Wm x m   Wm*  Wm  x m  ma m 1  ma  m0

(16)

By substituting W0  mw,in ma and Wm  Wm into Equation (16), it is simplified to the form:

mw *      w  w    Wm x m    Wm*  Wm  x m   ma  m0   m0 

(17)

where for m  1 . Expanding the product in Equation (17), it is obtained that   m  mw * w  w   W j Wm* j  Wm j   x m  ma m0  j 0 





(18)

Thus, if c p ,w  mw ma  w*  w   m0 Dm x m then

Dm 



m

c j 0

W j Wm* j  Wm j 

p,w

m0

(19)

By making the variable change Tw   x   in the derivative dw / dTw and substituting

w   m0 Wm x m , the derivative is expressed as: n

 dw  Pm x m dTw m0

(20)

11

where Pm   m  1Wm1 /  with m  0 . Therefore, Equation (1) is transformed to the following equation 

P x

m

m

m0

  



Dm x m

m0 

(21)

E xm m0 m

where

Em  Am  Bm  Cm

m0

(22)

By rearranging Equation (21), the following is obtained 



m

Pj Em j xm  Dm x m

m 0 j 0

As



m j 0

(23)

m0

Pj Em j  Dm , it is possible to obtain a recursive solution for the coefficients Wk through

this equation. Notice that the superscript m  1 is the highest superscript used for the coefficient

Wk in Equation (23). Specifically, the coefficient Wm1 appears in the coefficient Pm , which can be evaluated by the following equation

Pm 

Dm 



m 1

P Em j

j 0 j



m 1

(24)

E0

and subsequently using the expression Wm1   Pm  m  1 for obtaining the coefficient Wm1 . 3.2 Power Series Solution of the Equation of Heat Transfer Although in Equation (2) it could have been used an expansion in power series of normalized temperature for air humidity ratio and the enthalpy, is preferred using a macroscopic energy balance applied on the system enclosed between an arbitrary level and the top of the cooling tower, this is

12

ima 

mw,in ma

iw,in  ima ,out 

mw iw  0 ma

(25)

Considering that the specific heat capacity is constant along the cooling tower and that the water loss is by evaporation from water-stream to air-stream, then ima 

m  c p , w Tw,in  T0   ima ,out   w,in   wout  w   c p , w Tw  T0   0 ma  ma 

mw,in

(26)

Introducing the expansions w   m0 Wm x m and ima   m0 I m x m into Equation (26) and the 



variable change Tw   x   , then 

I m 1

m mx 

 m  c p , w   T0    w,in  Wm x m  c p , w  x    T0   0 ma m 1  ma 

mw,in

(27)

where wout  W0 and ima ,out  I 0 . Expanding the product in Equation (27) and rearranging in power of x thereby 

I m 1

m

 m   x m   c p , w w,in  c p , w   T0 W1  x   c p , w   T0 Wm  c p , wWm 1  x m  0 ma m2  

(28)

In order, to satisfy the Equation (28), the following relations have to be met

I1  c p , w

mw,in ma

 c p ,w   T0 W1

I m  c p ,w Wm1  c p ,w   T0 Wm

(29.a) m2

(29.b)

Notice that the coefficient I m 1 can be evaluated immediately after that the coefficient Wm1 are determined. 3.3 Determination of the coefficients W0 and I 0 The boundary conditions established in the mathematical formulation of this work does not allow to evaluate analytically the coefficients W0 and I 0 , but that it is necessary to satisfy the following

13

conditions



n m0

Wm  win and



n

I  ima ,in , because generally the inlet air-stream conditions

m0 m

are known. Therefore, the implementation of an iterative procedure in order to satisfy all constrains is unavoidable. For this purpose, the Newton-Raphson method was chosen for its quick convergence; and also, a procedure for obtaining the seed values of the coefficients W0 and

I 0 was developed. 3.3.1 Relation between the coefficients W0 and I 0 Applying an energy balance on the system enclosed between the base and top of the cooling tower, it can be expressed by

ima ,in 

m  c p ,w Tw,in  T0   ima ,out   w,in   wout  win  c p ,w Tw,out  T0   0 ma  ma 

mw,in

(30)

As W0  wout and I 0  ima ,out then ima ,in 

m  c p , w Tw,in  T0    w,in  win  c p ,w Tw,out  T0   I 0  W0c p , w Tw,out  T0  ma  ma 

mw,in

(31)

Introducing some simplifying notation as following

W0  a  bI 0

(32)

where  mw,in   win  c p , w Tw,out  T0   ima ,in  c p , w Tw,in  T0   m  a a c p , w Tw,out  T0 

b

1

c p , w Tw,out  T0 

(33.a)

(33.b)

That means that the coefficients W0 and I 0 are not independent from each other. 3.3.2 Procedure for the calculation of the seed values of coefficients W0 and I 0 14

Rearranging Equation (2) and subsequently evaluating at Tw  Tw,in (top of the cooling tower), is obtained the next equation:

dima dTw

c p,w  Tw Tw ,in

Making

Lef  i

the

* ma , out



mw,in *  Lef  ima ,out  ima ,out   1  Lef ma 



 ima ,out   1  Lef

variable

w

* out

 wout  i

* v ,out

Tw   x  

change

and

w

* out

 w

* out

 wout  iv*,out    wout  c p , w Tw,in  T0 

employing

the

series

(34)

expansions

di ma dTw  1  m1 mI m x m1 , w   m0 Wm x m and ima   m0 I m x m , Equation (34) is rewritten as 

I1



n

c p,w 

n



mw,in *  Lef  ima ,out  I 0   1  Lef  ma



* Lef  ima , out  I 0   1  Lef

w

* out

w

* out

 W0  iv*,out  

*  W0  iv*,out   wout  W0  c p , w Tw,in  T0 

(35)

where di ma dTw  1 I1 , w  W0 and ima  I 0 at x  0 (top of the cooling tower). Considering that the variation of enthalpy of air with respect to water temperature is constant along the cooling tower

dima ima ,out  ima ,in ima ,in  I 0   dTw Tw,in  Tw,out 

(36)

And because di ma dTw  1 I1 then I1  ima ,in  I 0 . Thus, Equation (35) is only function of the coefficients W0 and I 0 , and it is satisfied approximately. If the coefficient I 00 satisfies exactly the expression I1  ima ,in  I 0 then

ima ,in  I



0 0

c p,w 



mw,in * 0  Lef  ima ,out  I 0   1  Lef ma 



* 0 Lef  ima , out  I 0   1  Lef



w

* out

 W00  iv*,out  

 wout*  W00  iv*,out   wo*ut  W00  c p,w Tw,in  T0 

(37)

Where the coefficients W00 and I 00 must satisfy the energy balance given by Equation (32), this is W00  a  bI 00

(38)

then, the Equation (37) can be simplified to the following quadratic equation 15

 ima ,in  a  G ' G    ima ,in  a  H ' bG ' H  I 00  bH '  I 00   0     2

(39)

where

G   c p,w

mw,in

H   c p,w

mw,in



*  Lef ima ,out  1  Lef ma 



w

* out

 a  iv*,out  

(40.a)



 Lef  1  Lef biv*,out   ma 



* G '  Lef ima ,out  1  Lef



w

* out

(40.b)

*  a  iv*,out   wout  a  c p ,w   T0 

(40.c)



H '  Lef  1  Lef biv*,out  bc p ,w   T0 

(40.d)

Here the Equation (39) is a quadratic expression and hence have two roots. However, only the minor root has a physical meaning. Therefore, the coefficients W00 and I 00 can be considered proper seed values of the coefficients W0 and I 0 . 3.3.3 Evaluation of the Coefficients W0 and I 0 . Because that the seed value W00 satisfies approximately the boundary conditions and







 m0

Wm  win

I  ima ,in , an iterative procedure for refining an initial solution was performed. Let us

m0 m

n

consider the function g  Wm  win , then the root of g is the coefficient W0 . m0

n

Applying the Newton-Raphson method to g  Wm  win , it is possible to obtain the following m0

recurrence formula

W0k 1  W0k 

g W0k  dg dW0

k 0

(41)

W0k

16

where W0k 1 with k  0 is the value of the coefficient W0 obtained in the iteration k and

 dg

dW0 W k is the derivative of the function g with respect to the coefficient W0 evaluated in 0

W0  W0k . The derivative of g is computed numerically by the Euler forward difference method.

The general procedure for the calculation of the coefficients W0 and I 0 is presented in Fig. S2 in Section III of Supplementary Material. 3.4 Determination of number of transfer units The function to calculate the number of transfer units ( NTU ) is obtained analytically by the power series method, and is determined numerically by the Gauss quadrature method. 3.4.1 Expansion in Power Series of the Equation of Number of Transfer Units Let us consider NTU   m0 U m x m , which is an expansion for NTU in power of the 

normalized temperature of water. The first derivative of NTU with respect to Tw is given by

dNTU    Qm x m dTw m 0

(42)

where Qm  (m  1)U m1  for m  0 . As c p , w is constant, then c pw   m0 Fm x m where F0  c pw 

and Fm  0 for m  1 . Thus, Equation (3) can be expressed by the following equation 

 Qm x

m

m0

  



F xm

m0 m 

Em x m

(43)

m 0

By rearranging the equation above, we get 

m

 Q E

m0 j 0

j



m m m  j x   Fm x

(44)

m0

17

m

To satisfy the Equation (44), the following equation has to be met:

Q E j 0

j

m j

 Fm Thus, the

evaluation of Qm can be obtained by the equation below

Qm 

Fm 



m 1

Q j Em j

j 0



m1

(45)

E0

This calculation provides the coefficient U m1 through U m1   Qm / (m  1) . Certainly, this procedure allows to determine the unknown coefficients of the expansion NTU   m0 U m x m , n

but it does not guarantee the converge of the series solution. 3.4.2 Determination of the Number of Transfer Unit Using Gauss Quadrature Method The various techniques of quadrature try to determine with the minimal error the integral: b

A   f ( x)dx

(46)

a

The Gauss method (1814) (Lanczos, 1988; Davis and Rabinowitz, 1975) is one of the most efficient intergration method because it is based in Legendre polynomials roots  k , leading us to the following data points in [a,b]: xk 

ba a b k  2 2

k  [1,1],

k  1, 2,..., N

(47)

where the  k are the zeros of Pn ( k ) . Then, an interpolating polynomial G( x) is constructred, which approximates to f ( x) in [a,b] and allows to give a value close to Equation (46): b

 a

ba G ( )d 2 1 1

f ( x)dx 

(48)

So the corresponding interpolating polynomial G( x) of N  1 degree, adopts the form:

18

G( x)  y1 p1 ( x)  y2 p2 ( x) 

 y N p N ( x)

(49)

with the basic property G( xk )  yk  f ( xk ) . Therefore, Equation (48) expresses the Gaussian quadrature formula for simple points: b



f ( x)dx 

a

ba N k f (k ) 2 k 1

(50)

where k are the weight factors, and they are calculated by the following integral Cheney and Kincaid (2012) 1 N

k    1 j 1 j k

x  j

k   j

dx

(51)

Combining Equations (1) and (3), is obtained that

dw dTw

dNTU  mw * dTw  w  w ma





(52)



As dTw   dx , dw dTw  Pm x m and c p ,w  mw ma  w*  w   m0 Dm x m then m 0



m dNTU  c pw  m0Pm x   dx  Dm xm 

(53)

m0

Integrating from x  1 (base of the cooling tower) to x  z with 0  z  1 (any arbitrary level of the cooling tower), it is possible to obtain the function NTU ( z ) , which is the number of transfer units profile between the base and any arbitrary level of the cooling tower. Thus z

NTU ( z )   1

 c p , w  m0Pm x m 





D xm m0 m

dx

(54)

where NTU (1)  0 . Applying the Gauss quadrature rule to Equation (54), it can be rewritten as:

19

m z  1 N k  c p , w  m0Pm xk NTU ( z )    2 k 1  Dm xkm 

(55)

m0

where xk 

z 1 z 1 k  2 2

(56)

The maximum value of NTU occurs when z  0 , and such value is equal to NTU max  NTU (0) . 4. Results and discussion To demonstrate the accuracy degree of the obtained solution by power series a theoretical case and an experimental case are addressed. 4.1 Theoretical case For this scenario, five study cases reported by Rubio-Castro et al. (2011) were considered. To know the inherent error of the analytical solution, as well as to compare the computer processing time, the Equations (1)-(3) were numerically solved using the Dormand-Prince Runge-Kutta method (ODE45 MATLABTM solver, Matlab, 2015). This methodology required the implementation of the bisection method for obtaining wout , whose value is necessary for the calculation of mw along the cooling tower. The bisection method was initialized in the range of  win ,W00  to reduce the number of iterations required to meet the chosen tolerance (in this work,

104 ). Thus, Table 1 shows the operation configurations for the evaporative cooling process analyzed. In all cases, the calculation of the following variables: Ta ,out , ima ,out , wout and NTU , obtained from analytical solution, have a percentage error less than 0.5% with respect to the numerical methodology. Note that the runtime of analytical solution was 50 times faster than the runtime of ODE45 MATLAB solver. Fig. 3a-d show the absolute humidity profile, air enthalpy, and NTU along the tower; and a high degree of accuracy is observed. The NTU profile was 20

obtained by developing the power series and by numerical integration with Gaussian quadrature rule, and in most cases, quadrature rule exhibited minor deviation from the behavior obtained by the numerical methodology. Also, three sets of simulations of the evaporative cooling process were performed. In the first set, cooling water ranged from 50 °C to 20 °C using air with a dry bulb temperature of 22 °C and a wet bulb temperature of 12 °C under different inlet water mass and air mass flow rate ratios, mw,in / ma  0.40, 0.65, 0.90, 1.15, and 1.40 were considered. The high degree of accuracy was

provided by the power series expansion for w and ima (see Fig. 4a-d). However, the degree of accuracy for the profile of the NTU was affected by the mw,in / ma ratio; in fact, in the methodology of power series for NTU function at mw,in / ma  1 the degree of accuracy was good, in contrast for high values of that ratio the convergence was not achieved. On the other hand, Gauss quadrature method to evaluate NTU was always consistent, showing the following values of percentage error with respect to the numerical solution: 0.4124%, 0.8910%, 1.8873%, 3.0846%, and 1.5780%. The runtime of the analytical solution was approximately 50 times faster than the numerical methodology. In the second set of simulations, cooling water ranged from 50 °C to 20 °C using air with a dry bulb temperature of 22 °C and a wet bulb temperature of 7 °C under different inlet water and air mass-flow rate ratios, mw,in / ma  0.40, 0.65, 0.90, 1.15, and 1.40 were considered. The results (see Fig. S3a-d in Section III of Supplementary Material) were similar to those obtained in the first set of simulations, maintaining approximately the same ratio of computational runtime between analytical and numerical solution of 50. Also, the Gauss quadrature method to evaluate NTU was always consistent again, presenting the following percentage error values with respect to the numerical solution: 0.4394%, 0.6247%, 1.1336%, 3.1607%, and 5.0411%. 21

In the third set of simulations, the temperature range of the cooling water was changed, now ranging from 50 °C to 30 °C, and using air with a dry bulb temperature of 22 °C and a wet bulb temperature of 20 °C under different inlet water and air mass flow rate ratios, mw,in / ma  0.55, 0.95, 1.35, 1.75, and 2.15 . As a result, the value of 50 for the ratio between

runtime of analytical solution and the numerical methodology remained the same, with a high degree of accuracy of the solution in power series for w and ima (see Fig. S4a-d in Section III of Supplementary Material). In the case of the power series solution for NTU, this was not always convergent, but using the Gaussian quadrature was acceptable, with the following percentage errors: -0.1310%, -0.0835% -0.0340%, 0.4088%, and 3.4102%.

22

Results

Data

Table 1. Comparison of the analytic and numerical results obtained when solving the governing equations of heat and mass transfer in evaporative cooling process for five study cases reported by Rubio-Castro et al. (2011).

Process Variables

Case 1 Numerical Power Series Method Method

Case 2 Numerical Power Series Method Method

Case 3 Numerical Power Series Method Method

Case 4 Numerical Power Series Method Method

Case 5 Numerical Power Series Method Method

Tw,in

38.8866

29.5566

45.4517

24.1476

42.9877

Tw,out

20

20

20

15

25

Ta ,in

22

17

22

22

22

Twb,in

12

12

7

12

12

mw,in

29.9843

60.0479

22.1726

59.2602

31.0874

ma

43.2373

71.2273

31.4714

85.9841

35.8909

win

0.004682

0.006722

0.000182

0.004682

0.004682

ima ,in

34049.00

34138.77

22605.53

34049.00

34049.00

Ta ,out

28.5043

28.4974

23.4086

23.4078

30.4359

30.4267

21.2761

21.2722

30.7737

30.7805

wout

0.024163

0.024162

0.017756

0.017755

0.027029

0.027024

0.015627

0.015628

0.027573

0.027566

NTU

90428.43 2.360174

68756.88 1.680199

99773.86 2.050207 2.060613 0.036052

61129.69 4.406334

61128.54 4.407704 4.402474 0.030773

101545.08 1.404098

2.001074

68755.34 1.678701 1.678900 0.027927

99795.96 2.060704

CPU time (s)

90419.24 2.354766 2.359112 0.041459

101535.10 1.404784 1.402892 0.029488

% Error

ima , out

1.800249

1.88661

1.899575

1.821984

Ta ,out

-0.0023%

-0.0003%

-0.0030%

-0.0007%

0.0023%

wout

-0.0032%

-0.0015%

-0.0179%

0.0024%

-0.0244%

-0.0102% -0.2292% -0.0450%

-0.0022% -0.0892% -0.0773%

-0.0222% -0.5094% -0.0044%

-0.0020% 0.0134% -0.1053%

-0.0098% 0.0488% -0.0859%

ima ,out NTU

23

Fig. 3. Results for the governing equations of Poppe method obtained by the power series method and the Dormand-Price method, for five case studies reported by Rubio-Castro et al. (2011) (a) Profile of the air humidity ratio. (b) Profile of the moist air enthalpy. (c) NTU profile using the power series method; (d) NTU profile using the Gauss quadrature method.

24

Fig. 4. Results for the governing equations of Poppe method obtained by the power series method and the Dormand-Price method, using the following process parameters: Tw,in  50 C , Tw,out  20 °C , Tdb,in  22 °C , Twb,in  12 °C and mw,in / ma  0.40,0.65,0.90,1.15, and 1.40 kg water/kg dry-air. (a) Profile of the air humidity ratio. (b) Profile of the moist air enthalpy. (c) NTU profile using the power series method; (d) NTU profile using the Gauss quadrature method.

25

4.2 Experimental case The proposed analytical solution was validated with the experimental data reported by Braun et al. (1989) and compared with respect to the results obtained by the correlation of Khamis Mansour and Hassab (2014). In Table 2 the comparison of the experimental data with the values reported by Khamis Mansour and Hassab (2014) and the values obtained in this work are given. Here, it is possible to see that the proposed analytical solution presents a minor average error (approximately 3%) than the obtained by Khamis Mansour and Hassab (2014), whose average error is equal to 11%. Besides, it is important to mention that the proposed analytical solution allows to get the profiles along the packed cooling tower and not only the outlet conditions. In this sense, Fig. 5a-c presents the humidity, enthalpy of humid air and NTU profiles. When comparing the experimental results with respect to the correlation of Khamis Mansour and Hassab (2014) and the analytical solution of this work (see Fig. 5c), it is observed that the correlation underestimates the required NTUs (except in Case 2); in contrast to the analytical solution, whose estimated NTUs are very similar to the experimental ones. Failures to take into account this underestimation of the cooling tower design, operation problems could occur, specifically, a lower cooling capacity. On the other hand, for all cases, it is observed that the humidity and enthalpy profiles of the analytical solution are very similar to those obtained by the numerical integration of the system of differential equations of the Poppe model (see Fig. 5a-b). Thus, in the case of cooling tower under real operating conditions, Poppe model adequately describes its behavior, and also guaranteed that the analytical solution obtained (despite considering the constant thermophysical and transport properties) is very close to the exact solution of Poppe model. Finally, it is found that the humidity and enthalpy profiles strongly exhibit a linear behavior with respect to temperature, which establishes that the first and second

26

terms of the humidity and enthalpy solutions in power series are the ones that primarily determine the phenomenological behavior of the cooling process. It is important to indicate that for the above comparison the determination of coefficients W0 and

I 0 was not required; it is because all operating conditions are known. Therefore, the coefficients of power series were directly and recursively determined using the information of each experiment.

Table 2. Comparison between the experimental results by Braun et al. (1989) and the results calculated with the analytical solution in this work and the correlation reported Khamis Mansour and Hassab (2014). Case

mw,in / ma

NTU Braun

NTU Khamis

21.11

0.651

1.992

1.782

26.00

21.11

0.651

1.992

24.22

21.11

1.061

1.645

34.50

26.22

21.11

1.061

5

38.78

29.33

26.67

6

38.78

29.33

26.67

Tw,in

Tw,out

Twb,in

(°C)

(°C)

(°C)

1

31.22

23.88

2

41.44

3

28.72

4

Error of NTU Water NTU by (this work) lost Khamis (%) (%)

Error of NTU by this work (%)

-10.54

2.0379

1.12

2.30

2.016

1.20

1.9137

2.25

-3.93

1.485

-9.73

1.7208

0.70

4.61

1.645

1.453

-11.67

1.6398

1.21

-0.32

0.797

1.841

1.628

-11.57

1.8985

1.40

3.12

0.806

1.831

1.634

-10.76

1.9100

1.49

4.31

27

Fig. 5. Results by analytical solution using the parameters of Table 2: (a) Profile of the air humidity ratio. (b) Profile of the moist air enthalpy. (c) NTU profile versus experimental points.

28

5. Conclusions An analytical power-series solution for the problem of heat and mass transfer in evaporative cooling process was carried out. The variation of water flow-rate along the cooling process was taken into account and described as a series expansion in power of the normalized water temperature. The iterative procedure developed for the calculation of the coefficients of series expansions for the humidity ratio, air enthalpy, and NTU was simple to implement computationally. In all cases, the series solutions for the heat and mass equations were convergent and quite close to the numerical solution via Dormand-Prince method. The series solution for the equation of NTU was convergent for mw,in / ma  1 . Gauss quadrature integration method was easily applied as an alternative method for calculating the NTU profile. In general, the assumption of considering the physical properties and the Lewis factor as constant was an acceptable approximation. A comparison between the analytic solution and the numerical solution for the calculation of NTU gave  1 % percentage error at mw,in / ma  1 , and in the other cases no greater than 6 % . The validation of experimental data allows to affirm that the analytical solution developed in this work is a tool trustworthy that can be used to design cooling towers with a rigorous representation of transport phenomena. Finally, the computational implementation of this analytical expression is simple therefore it can be employed in complex optimization schemes. Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments The authors acknowledge the financial support by SEP-CONACYT (CB-2013/221249).

29

Nomenclature

Am

Coefficient of Eq. (8) from expansion in power series of the equation of mass transfer

Bm

Coefficient of Eq. (11) from expansion in power series of the equation of mass transfer

Cm

Coefficient of Eq. (14) from expansion in power series of the equation of mass transfer

cp

Specific heat capacities, J/kg K

cp

Arithmetic mean of specific heat capacities, J/kg K

Dm

Coefficient of Eq. (19) from expansion in power series of the equation of mass transfer

Em

Coefficient of Eq. (43) from expansion in power series of the equation of heat transfer

Fm

Coefficient of Eq. (44) from expansion in power series of the equation of number of transfer units

g

Auxiliary function of Eq. (41)

i

Enthalpy, J/kg

I0

Coefficient for power series

k

Index

Lef

Lewis factor

m

Mass flow rate, kg/s

NTU

Number of transfer units 30

Pm

Coefficient of Eq. (24) from expansion in power series of the equation of mass transfer

Qm

Coefficient of Eq. (45) from expansion in power series of the equation of number of transfer units

S

* m

Fit coefficients of series expansion of saturated humidity ratio

t

Temperature, °C

T

Temperature, K

w

Humidity ratio through the cooling tower, kg water vapor/kg dry air

W0

Coefficient for power series

x

Normalized temperature, K

Subscript 0

Reference

a

Air

in

Inlet

j

Index

k

Index

m

Mass transfer or index

n

Degree of polynomial solution

ma

Mixture of air and water vapor

out

Outlet

v

Vapor

w

Water

wb

Wet-bulb 31

db

Dry-bulb

top

Outlet air-stream

bottom

Inlet air-stream

Superscript *

Saturation at water temperature

Greek



Temperature range, Tw,out  Tw,in , K

k

Legendre polynomials roots

v

Latent heat of vaporization, J/kg

  Tw,in

Inlet water temperature, K

References ASHRAE-Handbook, 2008. HVAC Systems and Equipment, Cooling equipment and components, Cooling towers, Atlanta, GA, USA. Bosnjakovic, F., 1965. Technische thermodynamik, fourth ed. Theodor Steinkopff, Dresden. Braun, J., Klein, S., Beckman, W., Mitchell, J., 1989. Methodologies for optimal control of chilled water systems without storage. ASHRAE Trans. 95, 652. Cheney, E.W., Kincaid, D.R., 2012. Numerical Mathematics and Computing, seventh ed. Brooks Cole, Boston. Davis, P.J., Rabinowitz, P., 1975. Methods of numerical integration. Academic Press, New York. Heidarinejad, G., Karami, M., Delfani, S., 2009. Numerical simulation of counter-flow wet cooling towers. Int. J. Refrig. 32, 996. DOI: 10.1016/j.ijrefrig.2008.10.008 32

Hernández-Calderón, O.M., Rubio-Castro, E., Rios-Iribe, E.Y., 2014. Solving the heat and mass transfer equations for an evaporative cooling tower through an orthogonal collocation method. Comput. Chem. Eng. 71, 24. DOI: 10.1016/j.compchemeng.2014.06.008 Jaber,

H.,

Webb,

R.L.,

1989.

Design

of

Cooling

Towers

by

the

Effectiveness-NTU

Method.

J.

Heat

Transfer.

111,

837.

https://doi.org/10.1115/1.3250794 Khamis Mansour, M., Hassab, M.A., 2014. Innovative correlation for calculating thermal performance of counterflow wet-cooling tower. Energy. 74, 855. https://doi.org/10.1016/j.energy.2014.07.059 Kintner-Meyer, M., Emery, A., 1995. Cost-optimal design for cooling towers. ASHRAE J. 37, 46. Kloppers, J.C., Kröger, D.G., 2005b. A critical investigation into the heat and mass transfer analysis of counterflow wet-cooling towers. Int. J. Heat Mass Transf. 48, 765. https://doi.org/10.1016/j.ijheatmasstransfer.2004.09.004 loppers

.C.

r ger D. .

a. Cooling ower Performance Evaluation Merkel Poppe and e-NTU Methods of Analysis. J. Eng. Gas Turb.

Power. 127, 1. https://doi.org/10.1115/1.1787504 Kotb, A., 2013. Determination of optimum height for counter flow cooling tower. Asian Journal of Applied Science and Engineering. 2, 33. Kröger, D.G., 2004. Air-Cooled Heat Exchangers and Cooling Towers. PennWell Corp., Tulsa, Oklahoma, USA. Lanczos, C., 1988. Applied analysis. Dover Publications, New York. Liu, L., Cai., W., 2002. A Universal Engineering Model For Cooling Towers, International Refrigeration and Air Conditioning Conference. Paper 625. Liu, Y., Yu, J., Chen, L., Jin, S., 2018. Theoretical investigation and experimental verification of a mathematical model for counter-flow spray separation tower. Int. J. Heat Mass Transf. 120, 316. https://doi.org/10.1016/j.ijheatmasstransfer.2017.12.034 Llano-Restrepo, M., Monsalve-Reyes, R., 2017. Modeling and simulation of counterflow wet-cooling towers and the accurate calculation and correlation of mass transfer coefficients for thermal performance prediction. Int. J. Refrig. 74, 47. https://doi.org/10.1016/j.ijrefrig.2016.10.018 Mansour, M.K., 2016. Practical effectiveness-NTU model for cooling and dehumidifying coil with non-unit Lewis Factor. Appl. Therm. Eng. 100, 1111. DOI:10.1016/j.applthermaleng.2016.02.096 Matlab, 2015. Release 2015b, Release 2015b ed. The MathWorks, Inc., Natick, Massachusetts, United States. Merkel, F., 1925. Verdunstungskühlung. Z. Ver. Dtsch. Ing. 70, 123. Poppe, M., Rögener, H., 1991. Berechnung von Rückkühlwerken. VDIWärmeatlas, Mi 1- Mi 15.

33

Rao, R.V., More, K.C., 2017. Optimal design and analysis of mechanical draft cooling tower using improved Jaya algorithm. Int. J. Refrig. 82, 312. https://doi.org/10.1016/j.ijrefrig.2017.06.024 u io-Castro E.

erna- on le M. Ponce-Ortega, J.M. im ne - uti rre A.

. Optimal Design of Cooling Towers, in: Hossain, P.M.M. (Ed.),

Heat and Mass Transfer - Modeling and Simulation. InTech. Serna-González, M., Ponce-Ortega, J.M., Jiménez-Gutiérrez, A., 2010. MINLP optimization of mechanical draft counter flow wet-cooling towers. Chem. Eng. Res. Des. 88, 614. https://doi.org/10.1016/j.cherd.2009.09.016 Söylemez, M.S., 2001. On the optimum sizing of cooling towers. Energy Convers. Manag. 42, 783. https://doi.org/10.1016/S0196-8904(00)00148-5 Söylemez, M.S., 2004. On the optimum performance of forced draft counter flow cooling towers. Energy Convers. Manag. 45, 2335. https://doi.org/10.1016/j.enconman.2003.11.023 Sutherland,

J.W.,

1983.

Analysis

of

Mechanical-Draught

Counterflow

Air/Water

Cooling

Towers.

J.

Heat

Transfer.

105,

576.

https://doi.org/10.1115/1.3245624 Zheng, W.-Y., Zhu, D.-S., Zhou, G.-Y., Wu, J.-F., Shi, Y.-Y., 2012. Thermal performance analysis of closed wet cooling towers under both unsaturated and supersaturated conditions. Int. J. Heat Mass Transf. 55, 7803. https://doi.org/10.1016/j.ijheatmasstransfer.2012.08.006

34