Modelling atmospheric effects on performance and plume dispersal from natural draft wet cooling towers

J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Contents lists available at ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Modelling atmospheric effects on performance and plume dispersal from natural draft wet cooling towers A. Chahine n, P. Matharan, D. Wendum, L. Musson-Genon, R. Bresson, B. Carissimo CEREA (Atmospheric Environment Teaching and Research Center, Joint Laboratory École des Ponts ParisTech-EDF R&D, Université Paris Est), EDF-R&D, 6 quai Watier, 78401 Chatou Cedex, France

art ic l e i nf o

a b s t r a c t

Article history: Received 27 June 2014 Received in revised form 16 October 2014 Accepted 5 November 2014

The study of cooling towers in realistic atmospheric conditions is important to assess their performance and plume dispersal as a function of wind speed and hygrometric conditions for their economical and public health consequences. In the present work the plume formation and its dispersal are modeled using a finite volume method on grid including details of the towers. The air mass flow and heat transfer are simulated inside and outside the tower. The equations for the transport of momentum, mass and liquid potential temperature are solved using the CFD software Code Saturne. In contrast to other studies dealing with the performance of cooling towers, the exchange processes inside the tower are represented with a source terms approach. The adiabatic expansion of the plume at the exit of the tower in the atmosphere is accounted for using the thermodynamical laws. The results of the simulation are compared to the measurement at Bugey nuclear power plant. The results of the model are shown to be in good agreement with the field measurements from point of view of the air flow structures, plume patterns and thermodynamical variables. Based on this reference simulation we study the variation of the cooling tower performance against the wind speed and quantify the effect of the ambient wind speed increase on the reduction of the natural draft. We also discuss the sheltering effect that the upstream tower can have on the other. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Modelling Cooling tower Plume Dispersal CFD

1. Introduction A natural draft wet cooling tower (NDWCT) is an important part of the cooling system in thermal power plants. Its use is based on circulating the warm water to be cooled inside the tower where mass and heat are exchanged between air and water and then they are ejected into the atmosphere via a circulation induced by a buoyancy mechanism. However, the NDWCTs is still a subject of many questions about their efficiency and socio-economical impacts. In fact, the rejected saturated air could be charged with pollutants that can affect the health of the neighbouring populations. The extent of the contaminated zone depends actually on the towers operating conditions, its design and the atmospheric conditions. For the consideration of the environmental impact, many efforts were deployed to investigate the plumes behavior and their effects in the proximity of towers. For example, a model for seasonal and annual prediction of the plume drift deposition and shadowing effect (SACTI model) was developed by the Argonne

n

Corresponding author. E-mail address: [email protected] (A. Chahine).

http://dx.doi.org/10.1016/j.jweia.2014.11.007 0167-6105/& 2014 Elsevier Ltd. All rights reserved.

National Laboratory of University of Illinois (Policastro et al., 1994). The model is composed of several submodels allowing analytical computation of plume length, plume rise and drift deposition from single and several NDCTs. The different submodels were validated on the basis of the field experiments, like Chalk Point dye tracer experiment where dye is added to the water inside NDWCT to be able to measure the drift deposition at different distance from the tower (Meyer, 1975). In addition, mathematical models for plume trajectory prediction were developed and they are based on the lagrangian plume puff trajectory (Hanna, 1975; Winiarski and Frick, 1976). Carhart et al. (1982) examined the performance of sixteen analytical models by comparison of model prediction with field data. His main conclusion is that the performance of the models is sensitive to some parameters, such as the dilution rate of the plume in the ambient air and the presence of the thermodynamic of the moisture. Many theoretical and experimental studies were conducted in order to better understand the cooling tower performance with regard to internal and external parameters. Williamson et al. (2008) examined the effect of all parameters related to the warm water and ambient air conditions on the NDWCT performance. His main conclusion is that the NDWCT performance is strongly influenced by the thermodynamical state of the ambient air. Other studies

152

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

were focused on the effect of the internal operating conditions, considering the NDWCT as an open thermodynamic system interacting with a windless atmosphere (Kaiser et al., 2005; Fisenko and Brin, 2007; Heidarinjad et al., 2009; Rahard, 2010). As the windless atmosphere is just an idealized state, numerical and experimental studies were carried out to understand the effect of wind on the flow structure inside and around the NDWCTs and on their performance (Su et al., 1999; Gao et al., 2007; Hyhlik, 2011). Via CFD and experiments, it was found that the performance and flow structure vary as the external wind speed varies (Al-Waked and Behnia, 2006; Al-Waked, 2010). Al-Waked and Behnia (2007) pointed out that there exists a critical wind velocity below and beyond which the performance decreases and increases respectively. For the NDWCTs resembling the Bugey1 NDWCT, this critical velocity is 7.5 ms
1. All the previous studies deal with the mass–heat transfer and air flow in the proximity of the tower. In the studies of the plumes emitted from cooling towers and its dispersal in the atmosphere, the NDWCTs are usually considered as filled obstacles (Meroney, 2006; Lucas et al., 2010). The mass rate and plume characteristics are fixed as a Dirichlet boundary conditions at the top of the NDWCT (Bouzereau et al., 2008). Following this approach, the external atmospheric conditions affect the plume and flow structure only outside the NDWCT. Lucas et al. (2010) investigated the effect of psychrometric ambient conditions on drift deposition in the proximity of the cooling tower, by considering the water from the cooling tower as discrete droplets. At this scale, to our knowledge all the available works on drift were undertaken with Lagrangian approach under a predefined wind and thermodynamic temperature conditions. In fact, under the realistic atmospheric stability state of the atmosphere, the hot plume at the exit of the cooling tower is subject to the adiabatic expansion and diffusion limiting the rise of the plume and influencing the droplet deposition. Moreover, at the exit of the cooling tower, the plume of saturated air undergoes condensation when it meets the atmosphere. Hence, the concentration of the droplets of water is very sensitive to the drift eliminator efficiency and the atmospheric condition at the top of the tower and the Lagrangian study of droplet cloud dispersal becomes difficult to perform. One advantageous alternative is the use of the droplets of liquid water as dispersed phase within the Eulerian approach. The fraction of liquid water in proximity of the tower is then determined from the thermodynamic equilibrium between the liquid and gaseous phases of water in the atmosphere (Bouzereau et al., 2008). This approach is promising especially when the plume rises to higher altitudes where the water phase distribution is governed by the microphysics laws. So, in the present work we deal with modelling air flow inside and outside the NDWCT and plume dispersal under realistic atmospheric conditions with an Eulerian approach. The atmosphere state is almost neutral over the first kilometer, followed by a stably stratified layer. The liquid water phase evolution is solved with the Eulerian approach and its phase change is represented by specific laws of microphysics. As the plume dispersals at the exit of the towers are adiabatic in the absence of the abrupt phase change, we use the potential temperature which is a conservative variable in the energy equation. The concept of potential temperature will be discussed at length in Section 4. In addition, the heat and mass transfer between the warm water and air inside the cooling tower will be represented by the global source term approach, without accounting for the heat and mass exchange mechanisms required when the liquid phase is treated by lagrangian approach. The momentum loss in the towers will be also

1

Near Lyon, France.

represented by a specific source term. The use of the global source term approach would allow easy understanding of the flow structure and plume dispersal under several atmospheric conditions. The approach with source term was already used at little extent by Hyhlik (2011) for studying cooling tower flow. The model will be used to simulate the realistic plume dispersal test case for validation purpose. The data and observations of plume dispersal from NDWCT of Bugey collected by EDF will be used. Some aspects on NDWCTs are given in Section 2, while in Section 3 we describe the cooling tower of Bugey site and the measurement protocols. In Section 4 the numerical description of the model and different theoretical schemes are presented. Next, the used methodology will be described, and the qualitative and quantitative analysis of the results will be made in Section 5. Sensitivity of the NDWCT behaviour to the external wind speed will be done in Section 6. Finally, we conclude with discussions and drawing conclusions concerning the results and methodologies in Section 7.

2. Background on NDWCTs A NDWCT is a thermal system that extracts heat from the warm water by evaporative and, to a lesser extent, convective processes. The heat extraction is ensured naturally by the buoyancy effect resulting from the temperature gradient between the indoor and outdoor air of the tower. In Fig. 1 are illustrated the main parts of NDWCT used in thermal and nuclear power plants. The warm water at the temperature Twi is sprayed from the water system in the spray zone. Just below, we find a fill of height ranging between 1 and 2 m. The fill allows the heat and mass exchange between air and water. It is the zone where the main mass and heat exchange occurs; about 82% of total heat is exchanged (Al-Waked and Behnia, 2006). Next, once the water has left the bottom of the fill it falls down over a distance of about 8–9 m in the rain zone before reaching the water bassin at the temperature Two. Thus, the efficiency of the cooling tower is characterized by the temperature decrease of the water δT w ¼ T wi
T wo . Usually, the performance of the cooling tower is defined by the following expression:

ξ¼

δT w T wi
T awb

ð1Þ

where Tawb represents the wet bulb temperature of the ambient air, it is the minimal temperature that water could reach at the end of its course in the rain zone. From the point of view of energy balance, the factor ξ gives an incomplete picture, as δT w represents only the sensible heat lost by the water, without including the latent heat due to the _ the rate of water exchanged evaporation of water. If we note δm and rejected into the atmosphere, the total heat gained by air and lost by water is given by this expression: _ v δQ ¼ m_w Cpw δT w þ δmL

ð2Þ

where m_w is the warm water flow rate (kg/s), Cpw (J/kg K) is the heat capacity of water and Lv (J/kg) is the latent heat of water vaporization. This extracted energy is also equal to the enthalpy gained by the air at the exit of the tower.

δQ ¼ m_a Cpm ðT ao
T ai Þ þ m_a ðC o
C i ÞLv

ð3Þ

where m_a , Tao and Tai are the air flow rate, the outlet and inlet air temperature respectively. The parameter Cpm represents the heat capacity of the mixture air and vapor. The symbols Co and Ci stand for the fraction of water vapour content at the outlet and inlet of the tower (kg/kg).

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

3. Description of the NDWCTs of Bugey and experimental protocols The nuclear power plant of Bugey is located in the Rhone valley in East central France and is equipped with four towers, equidistant from each other, as shown in Fig. 2. The height H of each tower is 127 m. Their bottom and top diameters are 104 m and 68 m. The diameter of their throat is 61 m. The distance between the center of each couple of towers is about 180 m. The air inlet has a height of 12 m, and it corresponds to the rain zone height. The rain zone is followed by a fill zone of height of 1.5 m. Experiments on Bugey NDWCTs were performed during March 1980 by EDF (Electricité de France). Both the shape of the plume and its thermophysical properties were investigated (Hodin, 1984). The measurements were taken at the ground level and at higher altitude, using a meteorological plane equipped with sensors in the wing. The plane crosses the plume horizontally over a certain

153

distance before it crosses it again at higher altitude. The relevant variables relative to the plume measured during the experiments were the temperature, vertical velocity and its fluctuations. Instead of measuring water droplets number through the plume, the concentration of liquid water was measured using two methods: Johnson–Williams (JW) and FSSP 100 (Forward Scattering Spectrometer) (Hodin, 1984). The variables related to the atmospheric state, such as wind velocity and hygrometry, are obtained from nearby radiosoundings.

4. Numerical modelling In this section we present the details concerning the modeled domain, the phenomenological equations governing the flow and plume dispersal. We also emphasise the microphysics and the behaviour laws used in the model.

4.1. Mesh description

Fig. 1. Main parts of NDWCT.

The flow is solved in a cubic domain of size Lx  Ly  Lz ¼ 3 6:8  8:0  2:5 km . The domain is discretized with hybrid hexaedral and prismatic cells in the directions x, y and z. Within each tower the mesh was refined in order to better solve the flow structure, heat and momentum transfer. The four cooling towers were created and meshed using the hyperboloid profile of 2 NDWCTs ðx2 =a2 Þ þ ððz
zt Þ2 =b Þ ¼ 1; where a2 and b2 are respectively equal to 930.25 and 3614.44 for the Bugey towers, and the throat altitude zt is 95 m. The mesh size in the horizontal direction is 5 m. The same mesh size is kept on the outside of the tower but the cells are gradually stretched. In the vertical direction, the mesh size along the height of the tower is 4 m, with refinement in the fill zone. Above the tower, the mesh is gradually stretched vertically. The overall number of cells used is 9  106 . Fig. 3 shows the mesh structure on surfaces near the towers.

Fig. 2. Plant design of Bugey NDWCTs.

154

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Fig. 3. Top and front views of mesh sample inside and around the cooling towers.

Fig. 4. The velocity profile used for boundary conditions and initialization, based on nearby radio-sounding.

4.2. Boundary conditions and initialization

4.3. Governing equations

The ground and tower shells are considered as rough walls with roughness length of 0.1 m and 0.005 m respectively. Over the four lateral boundaries, the measured meteorological profiles of the velocity, temperature and humidity are applied. At the top of the domain a symmetry condition is applied for all variables. At the bottom, Dirichlet conditions are applied for all variables. These conditions are extrapolated from the measurements near ground level. Concerning the boundary condition for the vertical velocity w, the large scale field is considered with w ¼0 (no subsidence). During the experiment, the meteorological conditions were not exactly locally measured. So, we have spatio-temporally interpolated the meteorological data measured at other times and neighbouring sites. In the velocity profile, the effect of the geostrophic wind is accounted for by the change in the wind direction with altitude, as this can be seen in Fig. 4. Concerning the flow initialization, the flow is initialized with meteorological profiles in whole domain, except inside the towers where the horizontal velocity components are set to zero, and the vertical velocity is initialized to very small value in order to initiate the vertical air motion inside the towers. In Fig. 4 is shown the vertical profiles of the velocity components used for flow initialization and boundary conditions. The parameter representing the total number of water droplets is initialized to zero everywhere, except in the fill zone where n is set to 14,600 (Bouzereau, 2004).

The flow inside and outside the NDWCT is governed by momentum equations for turbulent flows. As for the water content and heat transfer, the equations of species and heat transfer are solved. For turbulent closure scheme, the k
ϵ model is used. All mean variables are noted without overbar, while the fluctuations are presented as primed variables. 4.3.1. Momentum equations The Averaged Navier stokes equations solved together with the k
ϵ equations are as follows:
 
 ∂ ∂ ∂p ∂ ∂ui ∂uj u i ¼
þ ρg i þ þ Sui ρ þ ρuj μt þ ð4Þ ∂t ∂xj ∂xi ∂xj ∂xj ∂xi   ∂ðρkÞ ∂ðρuj kÞ ∂ μt ∂k þ þ P þ G
ρϵ ¼ ∂t ∂xj σ k ∂xj ∂xj

ð5Þ

   ∂ðρϵÞ ∂ðρuj ϵÞ ∂ μt ∂ϵ ϵ þ þ C ϵ1 ðP þC ϵ3 GÞ
C ϵ2 ρϵ ¼ ∂t ∂xj σ ϵ ∂xj k ∂xj

ð6Þ

In the above equations, ui stands for the i component of the ! velocity vector u , with i ¼ 1; 2; 3 for x; y and z directions; gi is the gravity vector (0, 0,
9.81). The symbol μt is the turbulent diffusivity and it is calculated using the following expression:

μt ¼ ρC μ

k

2

ϵ

ð7Þ

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Table 1 Turbulent parameters fot the k
ϵ turbulent model. C μ1

C ϵ1

C ϵ2

C ϵ3

σk

σϵ

0.09

1.44

1.92

jGj þ G 2jGj

1

1.3

The parameter Sui is the source term due to the momentum loss inside the tower. The source terms P in the k equation represents the momentum transfer from the mean flow to the fluctuating part, while G represents the production or dissipation of k due to buoyancy and they are given by the following relation:
 ∂ui ∂uj ∂ui 2 ! ! P ¼ μt þ
ρk þ μt divð u Þ divð u Þ ð8Þ ∂xj ∂xi ∂xj 3

ρ 0 0 θ ug ð9Þ θv v i i where θv is the virtual potential temperature that will be disG¼


cussed in the next section. The turbulent parameters appearing in the momentum and k
ϵ equations are given in Table 1. 4.3.2. Continuity equation In the atmospheric flows, the anelastic approximation is usually used. This approximation states that the temporal variation of the air density is negligible compared to the other terms (Pielke, 1984). As the result, the continuity equation reduces to ∂ ðρu i Þ ¼ 0 ∂xi

ð10Þ

In this equation the density ρ of the mixture of air, liquid and vapour water phases is computed using the appropriate law that involves the fraction and density of each component of the mixture. For the gaseous phases (air and water vapor) the perfect gas law is used. Eq. (10) is used to derive the pressure Poisson equation for the pressure correction based on the predicted mass fluxes over the boundaries of the control volume. The pressure is initialized to the atmospheric pressure. Once the pressure is corrected, the predicted velocity is recomputed using the algorithm SIMPLE. For more details, the reader can refer to the manual of Code_Saturne (Archambeau et al., 2004). 4.3.3. Water conservation equation The transport of water inside and outside of the NDWCT is ensured by the equation of the conservation of mass and it is given for the gaseous and liquid phases by the equations for qv and ql:
   ∂ ∂ ∂ μt ∂qv ∂q ρ þui ð11Þ qv ¼ þ ρ v jE=C ∂t ∂xi ∂xi σ t ∂xi ∂t


   ∂ ∂ ∂ μt ∂ql ∂q ql ¼
ρ v jE=C ρ þui ∂t ∂xi ∂xi σ t ∂xi ∂t

ð12Þ

where ∂qv =∂t is the rate water evaporation or condensation within and outside the cooling towers. To avoid the term of phase change of water due to the evaporation and condensation the equation for the fraction of the total water content qw (kg/kg) (qw ¼ qv þql ) is indeed solved in the model. Its equation of transport is as follows:

   ∂ ∂ ∂ μt ∂qw ∂ql qw ¼ þρ ρ þui þ Sw ð13Þ ∂t ∂xi ∂xi σ t ∂xi ∂t sed where Sw is the source term of water. The distribution of water between liquid and vapour phase is given by thermodynamic and microphysics laws that will be discussed in the next section. The conservative variable qw is often used in micro-meteorological models (Bouzereau, 2004).

155

4.3.4. Energy equation Usually, the equation of energy solved in CFD models is the equation for thermodynamic temperature T. However, in modelling flow and plume dispersal in the atmosphere it is convenient to use the conservative energy variable introduced by Betts (1973), that is the liquid potential temperature. In order to justify our choice, one can firstly consider the fate of a dry air parcel at the exit of the tower. The air parcel moves up toward lower pressure levels where it undergoes an expansion. Because the local motion of air parcel is quick and the thermal conductivity of the plume is weak, no energy heat is exchanged between air parcel and its neighbouring. So, the air parcel has been the subject of the adiabatic expansion. The resulting temperature after adiabatic expansion is the potential temperature θ and it is given by the thermodynamic law of perfect gas:
ð
R=C p Þ p θ¼T ð14Þ p0 where R is the perfect gas constant for air and Cp is the heat capacity of air, the pressure p0 is equal by definition to 105 Pa and T is the thermodynamic temperature at the level where the pressure is p. The equation of energy in terms of the potential temperature is then as expressed as follows:
   ∂ ∂ ∂ C p μt ∂ θ θ ∂qv j þ Sθ ðC p θÞ ¼
ρLv ρ þ ui ð15Þ ∂t ∂xi ∂xi σ t ∂xi T ∂t E=C In the presence of liquid water, the liquid potential temperature

θl is used instead. It is conservative both in dry or wet atmosphere. It corresponds to the potential temperature that the air parcel would have when all its liquid water content is evaporated. Its is given as a function of θ in the absence of the liquid water by

 L q L q ð16Þ θl ¼ θ exp
v l C θ 1
v l CpT CpT Combining Eqs. (15) and (16) leads to the equation of the transport of liquid potential temperature θl that is solved in the model

   ∂ ∂ ∂ C p μt ∂θl θ ∂ql ρ þ ui þ Sθ l ð17Þ ðC p θl Þ ¼ þ ρLv ∂t ∂xi ∂xi σ t ∂xi T ∂t sed In this above equation no source term related to the evaporation or condensation appears because of the definition of the θl. In addition to the potential temperature definitions introduced above, we infer the virtual temperature θv, which is the temperature that would have the dry air with the same density and pressure as those of humid air. Its expression is given by θv ¼ θð1 þ 0:608qv
ql Þ, and it is used in the turbulent buoyancy production of turbulent kinetic energy and dissipation in the Eqs. (5) and (6). 4.4. Modelling source terms The local heat, water and momentum exchanges between warm water and air are modeled using the source term approach. This approach is based on the energy, mass and momentum balances. 4.4.1. Heat source term As reported previously, most of the energy exchange between warm water and air inside the NDWCT is localized in the fill zone. So, in Eq. (17), the source term representing the volumic air heat gain Sθl must be added at the level of the cells located in the fill zone. The value of heat gain P t is calculated using the energy balance given by Eqs. (2) and (3). In Table 2 are given the energy balance for one dispersal situation obtained by measurement and the numerical model for evaporative cooling tower computation Teferi tool (Bourillot, 1983; Grange, 1994) and obtained by Bouzereau et al. (2008). The value of Sθl is computed by dividing P t by the volume of the fill zone.

156

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Table 2 The characteristics of NDWCTs of Bugey reported in Bouzereau et al. (2008) and computed with Teferi tool for the 12 March trial. w is the vertical velocity at the exit of the cooling tower, Qair is the air flow rate, δT is the temperature excess at the exit of towers and Pt is the heat power exchanged between warm water and ambient air. Method

w ðms
1 Þ

Q air ðm3 s
1 Þ

δT ð1CÞ

P t ðMwattÞ

Bouzereau Teferi

3.70 4.26

13434 15456

17.20 17.10

277.16 317.04

where the last two terms represent the evaporation and the condensation of the water droplets, and the new droplet generation by nucleation. As the liquid droplets are assumed spherical, the variables ql and n verify the following relation:

4.4.2. Momentum losses As the air flows from the bottom toward the exit of the towers, it loses its momentum in different zones of the tower. In the rain and spray zones, the loses are due to the drag exerted by water droplets on the air. Whereas in the fill zone, the loses are imposed by the geometry of the fill and draining water. Instead of considering the effect of individual tower elements on the air flow, we represent their global head loss effect by the added volumic source term Sui at the cells where momentum loss occurs in Eq. (4). In our study, we neglect the momentum loss outside the fill as it is negligible when compared to that inside the fill. The general expression of the momentum loss in the fill is given by the following expression (Williamson et al., 2008; Hyhlik, 2011): 1 Sui ¼
ρK h;v V 2 2

ð18Þ

where V is the velocity magnitude of air in the fill and K h;v are the momentum loss coefficients associated to the horizontal and vertical directions. As the flow is essentially vertical in the fill zone, the coefficients of momentum loss in the horizontal direction Kh is set to an infinite value. In the vertical direction, the value of Kv is set to the constant value suggested by Caudron (1991), Kv ¼27.5 m
1. 4.4.3. Total water source term For the water equation, the volumic source term Sw is added only in the fill zone of the cooling towers. The value of the source term is deduced from the water balance between the air at the cooling tower exit and the neighbouring ambient air, as expressed by the following relation: _ a ðC o
C i Þ Sw ¼ m

Eqs. (19) and (20) constitute the basis of the thermodynamic model for the phase change. If we note n the number of droplets, its equation of evolution is then

  
 ∂ ∂ 1 ∂ μt ∂n ∂n ∂n þ ui þ ð22Þ n¼ þ ∂t ∂xi ∂t e=c ∂t nuc ρ ∂xi σ t ∂xi

ð19Þ

_ a , Co and Ci standing respectively for the air flow rate at the with m tower exit, and the water concentration of air at the tower outlet and inlet. 4.5. Microphysics and thermodynamic laws The solved variable qw gives the total water content in each cell. For this reason a diagnostic equation for liquid water content ql is necessary. Both inside and outside the cooling tower, the liquid phase water is considered as a dispersed set of droplets. At local temperature T and pressure p, the liquid phase ql is determined on the basis of total water content qw and the water vapour content at the saturation qs ðT; PÞ in the following manner:

qw
qs ðT; pÞ if qw Z qs ðT; pÞ ð20Þ ql ¼ 0 Else: where the water content at the saturation is described by the following formula: 8 0:62197 es ðTÞ > > > qs ðT; pÞ ¼ > < p þ ð0:62197
1Þ es ðTÞ ð21Þ
 > 17:2694ðT
273:15Þ > > ðTÞ ¼ 6:1078 exp e > : s T
35:86

4 3

ρql ¼ π r 3 3 ρw n

ð23Þ

where r 3 is the volumic mean average radius of the droplet. The spectral distribution of the droplets radius, which is important in the calculation of sedimentation and evaporation-condensation terms, follows the log-normal law with a logarithmic standard deviation σr equal to 0.53 (Bouzereau et al., 2008). In contrast to the evaporative model where the liquid and gaseous phase are determined by heat and mass transfer processes, in the thermodynamic model the two phases are determined by their thermodynamic equilibrium for given local thermodynamic conditions of pressure and temperature. Concerning the subgrid variation of ql and n, the all or nothing subgrid condensation scheme is used (Bouzereau et al., 2008). The nucleation of the droplets is included to account for the formation of the new droplets. For the supersaturated air, water vapour deposits on cloud condensation nuclei (CCN) to form new water droplets. The rate of formation of the new water droplets is related to the concentration of the CCN, NCCN, and the existing water droplets concentration in the atmosphere n. The process of nucleation in the atmosphere is described in detail by Pruppacher and Klett (2000). In the literature many models of the NCCN are suggested and we adopt in our study the formulation of Cohard et al. (1998) and Cohard and Pinty (2000), where NCCN is expressed in terms of the deviation of partial pressure of water vapour from saturation state.

5. Methodology, plume structure and model validation 5.1. Methodology The plume dispersal from NDWCT of Bugey (March 12, 1980) is used as test case to validate the model. For all cells in the domain, either inside or outside the cooling tower, the transported variables are solved with the finite volume method, with boundary conditions on the simulated domain, set to the constant radio sounding data, and source terms applied within the towers. The steady state solution of the equations is obtained as a transient solution, with fixed boundary conditions (assuming steady large scale conditions over the period of integration) and forcing within the tower. The observed conditions are also used for initializing all fields. The Coriolis force is neglected in the Navier Stokes equations because the Rossby number is large for such a small domain. Its effect is however implicitly taken into account through the observed wind boundary conditions applied to the simulations. In our calculation the flow is driven by the prescribed inflow conditions and the corresponding pressure field is computed by the code. This approach is very often used in CFD and seems to be a good approximation for such domains. The computation convergence is checked by the variation with iteration of the microphysical variables in some control volume. At the end of each iteration, the physical time within the domain is variable. Fig. 5 shows that the convergence is reached before the number of 2000 iterations in the considered probe, ensuring good accuracy of the results.

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

157

Fig. 5. Variation of the real temperature T and velocity W at the exit of the towers with iterations.

Fig. 6. The simulated plume shape given by the iso-surface ql ¼ 0:05 g/kg colored with the local vertical velocity w (a), particle trajectories colored with ql (b) and the observed shape by photograph (c).

5.2. Qualitative features of the plume and wind flow In Fig. 6 are given the simulated and observed shapes of the plume. The simulated plume is given by the iso-surface of ql, ql ¼ 0:05 g/kg colored with the local vertical velocity, (Fig. 6a) and particle trajectory tracking (Fig. 6b). The simulated plume rises to the same altitude as that observed in the experiment, about 1500 m, before it falls down slightly afterwards, as this is shown by the negative values of the vertical velocity of the plume. At the lower altitude the plume spreads vertically with little tilt toward the North direction because of the weak wind speed. At higher altitude, the accumulated cloud is slightly directed toward the South due to the wind direction change, and afterward it is advected horizontally by the strong horizontal wind. The model reproduces well the general behaviour of the plume. Fig. 7 shows the flow structure on the horizontal plane located at height z ¼ 100 m. As the wind blows from the South-West to

North-East, the flow is decelerated upwind of the four towers due to the resistance that towers exerts on the flow. Moreover, we observe clearly the skirting effect of the towers T1, T2, T3 on the flow, generating flow acceleration on the lateral side of the towers located on the transverse direction of the wind. The skirting effect and flow acceleration are less pronounced in the North-West side of the tower T2 because of the sheltering effect of upwind tower T1. The same observation is valid for the tower T3, where the flow is less accelerated on the lateral sides because of the sheltering effects of the towers T4 and T1. 5.3. Quantitative validation of the model At the exit of the towers, the mean flow properties are calculated and compared to those reported by Bourillot (1981), and they are summarized in Table 3. The simulated mean vertical velocities w are close to those of Bourillot (1981), with a maximum difference of

158

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

for the liquid water, ql is represented by the horizontal mean across the plume at each altitude. Even if at the exit of the towers a good agreement is obtained, see Table 3, at higher altitude there is a discrepancy between simulation and measurements. This could be due to the uncertainty related to the profile of humidity deduced from linear interpolation in space and time. In fact, the temporal and spatial variation of the humidity does not follow a linear law. In the appendix are shown the different profiles of humidity measured and reported in the literature for the experiment of March 12. Another probable source of uncertainty is the measurement quality: no negligible difference between measurements performed with JW and FSSP methods is noticeable.

6. Wind speed effect on NDWCT behaviour Fig. 7. The horizontal cross section of the horizontal velocity u at z ¼ 100 m. The small arrow indicates the wind direction.

Table 3 Mean plume characteristics at the cooling tower exit obtained by Code_Saturne and Bourillot (1981). Variables

Tower T1

T2

T3

T4


1

w (ms ) Model Bourillot

4.52 3.74

4.53 3.64

4.56 3.74

T (C) Model Bourillot

24.13 24.22

24.07 23.66

24.36 23.66

24.19 24.22

ql (g/kg) Model Bourillot

0.54 0.68

0.55 0.68

0.54 0.68

0.54 0.68

4.539 3.64

1.22 ms
1. Concerning the mean temperature at the exit of the towers, we find a good agreement between simulation and results of Bourillot (1981). This demonstrates that the heat source term approach used is correct and the temperature field within the towers is well solved. Moreover, the simulated liquid water is sufficiently close to those of Bourillot (1981). So, the processes of the condensation, evaporation and the microphysics of the water are valid. Fig. 8 shows the comparison with aircraft measurements of vertical profiles of the thermodynamic temperature T, vertical velocity w and the liquid water ql. As the measured temperature at each altitude is the mean across the plume, the measurements are presented together with the minimum and maximum simulated thermodynamic temperature deduced from the solved liquid potential temperature. During the measurement, the plane flew at different levels passing through the plume and clear air. However at the time (1980s) there were no onboard recording of the horizontal position of the measurement relative to the towers. To overcome this lack of information, in the postprocessing of the simulation results, we have taken for each horizontal layer of cells the minimum and maximum value of w, T and ql. This procedure has been applied to all levels in the simulation (and not only for the altitude of each flight level), as we do not have an accurate height of the plume and limit of visibility. Within the plume, itself, there is a maximum and minimum of w, ql, and T that are characteristic of the plume. Above the plume the selection of max and min is continued until the top of the domain. At all altitudes, the measured mean temperature is between the simulated minimum and maximum. This shows that the choice of the liquid potential temperature as thermal variable is consistent. Concerning the vertical velocity w, the simulated minimum and maximum vertical velocity are between the measurement points. As

To evaluate the effect of the wind speed on the tower exit characteristics, the average variables at the exit of the tower are calculated. The original measured wind profile, given in the appendix, for which the model is validated, is weighted by the factor α ranging between 0.16 and 4.5 to obtain the modified wind profile. To highlight the effect of wind speed on the tower behaviour, the velocity values of the profiles at the height of 150 m are used as reference velocities in the graphics that follows. In the measured profile case, the reference velocity at z ¼ 150 m is equal to 1.65 ms
1. Fig. 9 shows the variation of the vertical velocity w and the temperature T at the tower exit. The velocity w at the exit of all towers decreases as the wind speed increases for sufficiently high wind speed. For the tower 3, located downstream the three others, we see that there exist two critical wind speeds Vc1 and Vc2, where the velocity w variation changes, with the values of Vc1 and Vc2 equal to 1.65 and 3.3 ms
1 respectively. For wind speed less than Vc1, the velocity W decreases with wind speed. When the wind speed reaches Vc1 the velocity W at the exit of the tower increases with the wind speed going up to Vc2 before it begins to increase. Concerning the temperature profiles, we observe that for high wind speed, the mean temperature at the exit of the tower increases as the wind speed increases. We observe also the presence of the same two critical wind speeds. Fig. 10 shows the variation of the liquid water ql and total water qw rates at the exit of the towers. It follows the same tendency as that of the vertical velocity. The air flow rate at the exit and tower performance are presented in Fig. 11. The performance ξ, defined by Eq. (1), is deduced by equating the enthalpy loss of water, given by Eq. (24) and the enthalpy gain of air, given by (25) and dividing by T wi
T ai : _ v Lv _ w Cpw ðT wi
T wo Þ
δm Qw ¼ m

ð24Þ

_ a Cpa ðT ao
T ai Þ Qa ¼ m

ð25Þ

The final expression of performance is as follows: _ a Cpa ðT ao
T ai Þ δm_ v Lv m þ ξ¼ _ _ w Cpw ðT wi
T ai Þ m w Cpw ðT wi
T ai Þ m

ð26Þ

_ v is the simulated water rate at the exit of the tower, where δm Lv is the latent heat of the water vaporization (Lv ¼ 2501 kj/kg), Cpw _ a is the simulated is the heat capacity of water (Cpw ¼4200 J/kg K), m air flow rate, Cpa is the heat capacity of air (Cpa ¼1017.24 J/kg K). The parameters related to the warm water to be cooled are set _ w ¼ 17265 kg=s, temperature inlet to: warm water flow rate m Twi ¼44.1 1C. The performance decreases with the wind speed for sufficiently high wind speed, except for the tower 3 presenting a complex variation. The air flow decreases as the wind speed increases in all towers, except for the tower 3 presenting the same two critical wind speeds previously discussed. The decrease of the air flow rate

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

159

Fig. 8. The simulated and measured properties of the plume, (a) T (1C), (b) w (m/s) and (c) ql (g/kg).

Fig. 9. The mean plume velocity Wðm=sÞ and temperature Tð1CÞ at the exit of the towers for different reference wind speeds.

at the exit of the towers with the wind speed is due to the ambient flow structure at the exit. When the wind speed increases, at the top of the tower, vortices occur and reduce the exit area of the

tower. The two critical velocities that appear in pronounced manner in the tower 3 could be explained by the fact that the tower is located in the sheltering zone, so it would behave to some

160

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Fig. 10. Same as Fig. 9 but for liquid phase ql and total water qw flow rates.

Fig. 11. Air flow rate ma (kg/s) at the exit of the towers and their performance, computed with Eq. (26), as a function of the wind speed.

extent like a singular tower for which the performance decreases with the wind speed, before it increases as the wind speed increases more (Al-Waked and Behnia, 2007). In addition, the non-unicity of the critical velocity could be attributed to the complexity of the flow and multi-towers interaction. The air flow rate has the same evolution as the vertical velocity w, ql and the total water qw. This is consistent with the expected

results, as for all wind speeds the amount of water released inside the towers as source term is the same and the rejected water into the atmosphere is directly related to the air flow rate. Concerning the correlation between the air flow rate and the mean temperature at the exit of the towers, we see that the air flow rate varies with wind speed in opposite manner compared with the temperature. This is due to the fact that for all wind

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

speeds, the heat released in the tower is the same, and the air become hotter as the air flow rate is less. For all towers, the performance is computed using the model results and standard data of the cooling tower of Bugey. To better understand the effect of the wind speed on the cooling tower behaviour, the flow structure is examined. Fig. 12

161

shows the flow structure on the horizontal plane at the exit of the towers for values of α equal to 0.5, 1.5 and 4. We observe clearly that for the small wind speed, the flow at the exit is symmetrically distributed; this structure is close to the case of the windless atmosphere (Al-Waked and Behnia, 2006). For α equal to 1.5, the symmetry disappears completely in three towers, but it remains in

Fig. 12. Horizontal cross sections of (top to bottom) vertical velocity w (ms
1), turbulent kinetic energy k (m2 s
2 ), temperature difference DT (1C) between the simulated field and the reference temperature far from the towers exit level, and liquid water ql (kg/kg) fields at the exit of the towers for three reference wind speeds (α equal to 0.5, 1.5 and 4 from left to right). The wind direction at the top of the towers is from South-West.

162

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

tower 3 as it is subject to the sheltering effect and the wind speed is not sufficiently high to disturb completely the flow at the exit. For higher velocity, the symmetry disappears even in tower 3.

As the wind speed increases, the vertical velocity in the upwind part of the exit area is lowered and the high vertical velocity at the exit of the tower is pushed toward the downwind part. This is due

Fig. 13. Fields of the vertical velocity w, temperature difference DT, turbulent kinetic energy k and liquid water content ql on the vertical plane crossing tower 3 for three reference wind speeds (α equal to 0.5, 1.5 and 4 from left to right).

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

163

Fig. 14. The potential temperature and humidity profiles used for boundary conditions and initialization.

to the increase in the turbulence intensity with the wind speed around the exit, as this can be seen in Figs. 12 and 13. The shape of the velocity field at the exit of the tower suggests that the air circulating along the tower is composed of two counter-rotating vortices, and this can also be confirmed by the shape of the liquid water concentration just above the tower exits. The vertical structure of the flow is presented in Fig. 13 for different reference wind speeds (α equal to 0.5, 1.5 and 4). For α equal to 0.5, the plume rises almost vertically with high vertical velocity w, equivalently high mass flow rate, and low turbulence at the exit of the tower. Concerning the area with large negative vertical velocity downwind the plume top visible for α ¼0.5, we also observed it, to a lesser extent, for α ¼1.5 and also in Fig. 6 with the observations (α ¼ 1, by definition) and the model trajectories. We think that this is an overshoot of the plume, which goes higher than its equilibrium position and then returns. Some oscillations can be noticed on Fig. 6, characterizing the gravity wave under stable zone at high altitude (Fig. 14). As the wind speed is low (α ¼ 0.5), the plume charged with high humidity rises up vertically with higher speed and the effect is stronger. Also at higher altitude, where the plume reaches its maximum height, it becomes heavier because of the condensation of the water vapour increasing its negative vertical velocity. This is highlighted with the ql contours at the bottom of Fig. 13 which shows this effect for α ¼0.5 and 1.5, with the maximum ql region at the top of the plume. For α ¼4, the wind is too strong to see the effect. As the wind speed increases, the plume is tilted and the vertical velocity is reduced, inducing consequently the reduction in the air flow rate and tower performance. At high wind speed, the air flow at the tower exit is surrounded by a zone of high turbulence, limiting air and water extraction from the towers.

7. Discussion and conclusions By using realistic atmospheric conditions observed on the site of the Bugey nuclear power plant (4 NDWCT) we have successfully simulated the cooling towers thermodynamical circulations and the plume dispersal, using the open source CFD software Code Saturne. For these simulations, the heat and mass transfer between warm water and cooling air within the towers are represented with a global source term approach. The water phase change inside and outside the towers is represented by the thermodynamic and microphysics laws,

using the liquid potential temperature, total water content and the number of droplets as prognostic variables. Comparisons of the model results and observations at the cooling tower exit and obtained with aircraft flights showed that the plume dispersal is well reproduced. For example we have good comparison of photograph of the plume with iso-surface of simulated liquid water content. In addition vertical profiles of temperature, vertical velocity, and liquid water content have been shown to agree well with the measurements. Some conclusions can also be drawn concerning the methodology used and the physical processes taken into account:

 The use of the liquid potential temperature as the thermal



variable for plume dispersal at the scale of the atmospheric boundary layer is suitable to take into account the different relations between the microphysics of water in the atmosphere, the adiabatic expansion of the plume cloud, the evaporation and condensation of droplets and the atmosphere stratification in terms of the temperature and pressure. The use of the global source term approach appears as an alternative for study of plume dispersal from cooling tower when the spatial scale of plume modelling is large enough.

There can be a strong influence of the tower on each other, for example in our case the tower 3 is strongly influenced by the wake from the other upstream towers. A sensitivity study of the cooling tower for different wind speeds, and the analysis of the flow field near the tower exits has revealed how the plume characteristics at the exit and tower performance are influenced by the wind speed. As the wind speed increases, the air flow rate within the towers and their performance decreases for the towers configuration dealt within this study. By this study, we have quantified the effect of strong wind on reducing the natural draft from the towers. For example the performance goes from 0.455 for a 1.24 ms
1 wind speed to 0.438 for a 7.40 ms
1 wind speed. Although our study is focused only on one configuration of four cooling towers and one wind direction the fact that the model results compare well with the observations, together with the fact that our previous studies have shown good agreement for the same site and different meteorological situations, will allow us to apply our tool to other sites and configurations. In addition, as our study deals with the realistic NDWCT of Bugey, our results about cooling tower performance could not be further quantitatively

164

A. Chahine et al. / J. Wind Eng. Ind. Aerodyn. 136 (2015) 151–164

Fig. 15. Vertical profiles of humidity at different times and sites.

compared to those available in the literature considering in most cases an individual tower. In future we intend that the model will be further validated to better understand its limit and scope. Moreover, this model could be used to study the dispersal of the possible contaminants present in the plume in order to assess the level of the underlying pollution and risk exposure.

Acknowledgments Many thanks are for Frank David for his fruitfull advices, his help and for the Teferi tool that he has provided us. We would also like to thank Eric Gilbert for his help and informations provided about the NDWCTs of Bugey. Our thanks are also for referes for their helpful comments.

Appendix A In Fig. 14 are shown the profiles of the potential temperature and humidity used for flow initialization and boundary conditions. In Fig. 15 are illustrated the vertical variation of the measured humidity profiles in Satolas and Loyettes meteorological stations together with the resulting interpolated profile and that used for validation of the model Gedeon (Louis et al., 1984). We see clearly the high variability of the humidity with space and time which could explain the discrepancy between the simulated and measured vertical profile. References Al-Waked, R., 2010. Crosswinds effect on the performance of natural draft wet cooling towers. Int. J. Therm. Sci. 49, 218–224. Al-Waked, R., Behnia, M., 2006. CFD simulation of wet cooling towers. Appl. Therm. Eng. 26, 382–395. Al-Waked, R., Behnia, M., 2007. Enhancing performance of wet cooling tower. Energy Convers. Manag. 48, 2638–2648. Archambeau, F., Méchitoua, N., Sakiz, M., 2004. Code_saturne: a finite volume code for the computation of turbulent incompressible flows—industrial applications. Int. J. Finite 1, 1–62. Betts, A., 1973. Non-precipitating cumulus convection and its parametrization. Q. J. R. Meteorol. Soc. 99, 178–196.

Bourillot, C., 1981. Teferi et soutef: Modéles numériques pour le calcul des performances des réfrigérants atmosphériques humides. Technical Report, EDF-DER HT/31-81.06. Bourillot, C., August 1983. Teferi, numerical model for calculating the performance of an evaporative cooling tower. Technical Report EPRI Report CS-3212-SR, Electric Power Reserach Institute, California, USA. Bouzereau, E., 2004. Représentation des nuages chauds dans le modéle météorologique mercure: applications aux panaches d'aéroréfrigérants et aux précipitations orographiques (Ph.D. thesis), Ecole doctorale des sciences de l'environnement d'ile-de-france. Bouzereau, M., Musson Genon, L., Carissimo, B., 2008. Application of a semi-spectral cloud water parameterization to cooling tower plumes simulations. Atmos. Res. 43, 120–133. Carhart, R.A., Policastro, A., Ziemer, S., 1982. Evaluation of mathematical models for natural draft cooling tower plume dispersion. Atmos. Environ. 16, 67–83. Caudron, L., 1991. Les réfrigérants atmosphériques industriels. Collection de la Direction des Études et Recherches d'Électricité de France. Cohard, J., Pinty, J., 2000. A comprehensive two-moment warm microphysical bulk scheme. Q. J. R. Meteorol. Soc. 126, 1815–1842. Cohard, J., Pinty, J., Bedos, C., 1998. Extending twomey's analytical estimate of nucleated cloud droplet concentration from ccn spectra. J. Atmos. Sci. 55, 3348–3357. Fisenko, F., Brin, A., 2007. Simulation of a cross-flow cooling tower performance. Int. J. Heat Mass Transf. 50, 3216–3223. Gao, M., Sun, F., Wang, K., Shi, Y., Zhao, Y., 2007. Experimental research of heat transfer performance of natural draft counter flow wet cooling tower under cross-wind conditions. Int. J. Therm. Sci. 47, 935–941. Grange, J.L., September 1994. Calculating the evaporated water flow in a wet cooling tower. In: 9th IAHR Cooling tower and Spraying Pound Symposium. Brussels, pp. 20–23. Hanna, S.R., 1975. Predicted and observed cooling tower plume rise and visible plume length at the John E. Amos power plant. Atmos. Environ. 10, 1043–1052. Heidarinjad, G., Karami, M., Delphani, S., 2009. Numerical simulation of counterflow wet-cooling towers. Int. J. Refrig. 32, 996–1002. Hodin, A., 1984. Etude experimentale des panaches de réfrigérants atmosphériques du bugey, premier bilan des mesures physiques effectuées entre le 25 fevrier et le 13 mars 1980. Technical Report, EDF. Hyhlik, T., 2011. CFD computational of natural draft cooling tower flow. In: Colloquium Fluid Dynamics. Kaiser, A., Lucas, M., Viedma, A., Zamora, B., 2005. Numerical model of evaporative cooling processes in a new type of cooling tower. Int. J. Heat Mass Transf. 48, 986–999. Louis, F., Biscay, P., Saab, A., 1984. Modélisation tridimensionnelle de l'impact atmosphérique des aéroréfrigérants de grande puissance, validation du code gedeon. Technical Report, EDF. Lucas, M., Martinez, P., Ruiz, J., Kaiser, A., Viedma, A., 2010. On the influence of psychrometric ambient conditions on cooling tower drift deposition. Int. J. Heat Mass Transf. 53, 594–604. Meroney, R., 2006. CFD prediction of cooling tower drift. J. Wind Eng. Ind. Aerodyn. 94, 463–490. Meyer, J.H., 1975. Chalk point surface weather and ambient atmospheric profile data, first intensive test period 15–19 December. Applied Physics Lab., Johns Hopkins University, PPSP-CPCTP-4REV (and Environmental Systems Corporation document PPSP6-CPCTP-8). Pielke, R.A., 1984. Mesoscale Meteorological Modelling. Academic Press, New York, NY, p. 612. Policastro, A., Dunn, W., Carhart, R., 1994. A model for seasonal and annual cooling tower impacts. Atmos. Environ. 28, 379–395. Pruppacher, H., Klett, J., 2000. Microphysics of Clouds and Precipitation. Kluwer Academic Publishers, Boston, MA. Second Revised and Enlarged Edition with an Introduction to Cloud Chemistry and Cloud Electricity. Rahard, T., 2010. Modélisation de l'aérodynamique d'une tour d'aéroréfrigérant humide par chainage d'outil de CFD (Master's thesis). EDF-ENSMA. Su, M., Tang, G., Fu, S., 1999. Numerical simulation of fluid flow and thermal performance of a dry-cooling tower under cross wind condition. J. Wind Eng. Ind. Aerodyn. 79, 289–306. Williamson, N., Armfield, S., Behnia, M., 2008. Numerical simulation of flow in a natural draft wet cooling tower- the effect of radial thermofluid fields. Appl. Therm. Eng. 28, 178–189. Winiarski, L., Frick, W., 1976. Cooling Tower Plume Model. Technical Report, USEPA Ecological Research Series, EPA-600/3-76-100, USEPA, Corvalis, Oregon.