Numerical simulation of deposition of drifts and salt from multiple super-large seawater cooling towers
With the rapid development of coastal industry and the increasing shortage of fresh water resources, the secondary circulation of large seawater cooling tower is often used as the main cooling method in coastal power plants . The hot and humid air discharged from the seawater cooling tower contains small water droplets with large soluble salt content. They are carried out of the tower by the updraft to form drifts, which may cause corrosion on the surface of equipment pipes , increasing the need of anticorrosion. It is also possible to cause the surrounding vegetation dryout and soil salinization, affecting the local ecology – . Such salt deposition issue prohibits the application and generalization of seawater cooling towers. In 1972, a series of studies were carried out in the United States, including cooling tower drifts, and numerical models such as ballistic model, Gaussian model and K theory model were proposed – . In 1990, the National Laboratory Center for Environmental Impact Research summarized previous models and field tests data to form the SACTI calculation program . The SACTI model was used to predict the plume dispersion discharging from cooling tower under normal operation and three different kinds of wind speeds by Guo RuiPing et. al. . The results have shown that the plume length, plume height and plume radius are function of wind speed. The SACTI program still has drawbacks in simulating complex terrain. Therefore, when the buildings and terrain near the cooling tower are too complex, researchers tend to perform wind tunnel tests or field tests for a more realistic observation, notably of plume disperstion . In the recent decades, Computational Fluid Dynamics (CFD) has improved dramatically and also become a widely used method in the study of flow field around cooling tower and its pollutant dispersion considering complex buildings and terrain within the calculation domain. Its basic principle is to use discrete numerical method to solve the fluid physics. It has strong adaptability and is widely used in many engineering fields. The CFD method was used to study the influence of ambient wind on the exhaust steam of cooling tower by Takata and Meroney , and the simulation results were compared with field measurement results, which proved good consistency between the two. Shi Xuefeng et. al. compared the CFD model and the SACTI model predictions of drifts and salt deposition in a cooling tower case, in which the CFD model considered evaporation, condensation, and local circulation and buildings, resulting in that the peak value of salt deposition predicted by the SACTI model was about 3 times as much as that of the CFD model. By comparing the simulation results of fog plume trajectory with the field observation results, CFD is deemed more reliable in flow field prediction, and suitable for the study of drifts and salt deposition problems. In this study, CFD was applied to a group of six super-large seawater cooling towers in a power site to study the drifts characteristics in the cooling towers, the salt depostion and the resulting influence to the environment, considering realistic layout of the surrounding buildings and measured meteo data. The results provide technical reference for the selection of potential new sites, design of layouts, and for the assessment of the impact of cooling towers to the environment. 

There are six cooling towers arranged in a north–south direction, with lower buildings on the east side of each cooling tower. The center spacing of the six towers from north to south is 250 m, 290 m, 250 m, 290 m and 250 m respectively, as shown in Fig. . The main dimensions of the cooling towers are shown in Table . "S-wave" fill of 2.0m height is used. The salinity of the circulating seawater is 28‰ and the concentration ratio is 1.5. 

Control function The operating conditions including loads, meteo and the mass flowrate of the circulate water is assumed stable, the airflow field is assumed steady around the towers. The function of flow movement is shown in reference . The porous media model of Fluent was adopted in the fill area, where the resistance loss coefficient of the fill is shown in Eq. ( ). Discrete Phase Model (DPM) was used in rain zone, where the air resistance to raindrops was analyzed according to Eqs. ( ) to ( ). MERKEL model was used to analyze the heat exchange process between fill and rain zone. 1 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$_{f} =  V^{M} _{a} }}{{ V^{2} H_{f} }}$$ ξ f = 2 A p V M γ a ρ V 2 H f 2 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$_{p} }}{dt} = F_{D} ( { - _{p} } ) + _{x} (_{p} -  )}}{{_{p} }}$$ d u → p dt = F D V → - u → p + g → x ( ρ p - ρ ) ρ p 3 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$F_{D} =  d_{p}^{2} }} {} }}{24}$$ F D = 18 μ ρ p d p 2 C D Re 24 4 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$C_{D} = }}( {1 + b_{1} {}^{{b_{2} }} } ) +  {} }}{{b_{4} + {} }}$$ C D = 24 Re 1 + b 1 Re b 2 + b 3 Re b 4 + Re where u p is particle velocity vector, m/s, V is air velocity vector, m/s, g x is acceleration vector of particle, m/s 2 , ρ p is density of raindrops, kg/m 3 , d p is raindrop diameter, m, C D is raindrop resistance coefficient, and bi is experience coefficient. The commercial computational fluid software FLUENT was used to solve the three-dimensional numerical model of cooling tower. For natural-draft cooling tower, buoyancy is one of the main forces in momentum equation, so buoyancy effect was considered. Air flow, heat transfer, raindrop evaporation and water vapor diffusion need to be solved coupled. The discrete equation is solved by the separation of variables method, and the decoupling of velocity and pressure is performed by the SIMPLEC algorithm. QUICK format was adopted in discretization of velocity and temperature field. The second boundary conditions of outers, the resistance characteristics of the fill, and the evaporation of the raindrop are calculated by the secondary interface of FLUENT software, with user-defined function (UDF) enabled. Grid establishment GAMBIT grid division software was used to establish a three-dimensional model, as shown in Fig. , with calculation range of 12500 m × 4500 m × 2500 m, divided into near field and far field, where the near field is divided into inside-tower and outside-tower, the outside-tower is divided into the perimeter and top, the inside-tower is divided into above-throat, below-throat, fill area, bottom of the tower and air inlet area. The grid of tower outside is sparse, and the fill area is the densest, with a total grid of about 8 million. In order to calculate the trajectory of particles in different wind bearings, a set of grids were established with interval of 15° starting from 0° where the wind blows parrallely to the towers row, such that the distribution of salt deposition in all bearings and distances around the tower group could be conducted. Boundary condition The computational domain is divided into two parts: inside and outside the cooling tower. For the boundary of the outer region, the bottom is an adiabatic boundary. When there is no wind, the other surfaces are pressure outlet boundaries. In case of ambient wind, the inlet surface is a velocity inlet boundary, and the other surfaces are pressure outlet boundaries. The cooling tower shell is arranged as an adiabatic boundary, and the inlet and outlet of the tower are arranged as the internal boundary. The fill is set to porous medium area, where the resistance coefficient is set according to the experiment experience. Meteorological parameters and operating parameters are shown in Table . The inject surface of the raindrop is the filltop. The initial velocity, temperature and amount of raindrop are set in the DPM. which can be calculated by operating parameters. The motion characteristics of floating droplets under different particle size distribution conditions were studied. The drifts spectrum is shown in Table . In this study, the floating water rate of the water collector is 5 parts per million (5 × 10 –6 ). The motion trajectory of a drift is calculated according to Eqs. ( ) to ( ). The evaporation of drifts is neglected when calculating the salt deposition. The annual or monthly distribution of salt deposition is closely related to the distribution of deposited quantity of drifts at each wind direction and wind speed Fig. shows the annual wind rose diagram, in which 1 m/s represents the wind speed interval of 0 ~ 1 m/s, 2 m/s represents the wind speed interval of 1 ~ 2 m/s, etc. The ambient wind speed is correlated with the height, as shown in Eq. ( ), the wind shear index is related to site and atmospheric stability category, here for example in neutral category is 0.19 . The UDF function is used to setup the wind speed vertical profile. To characterize the influence of ambient humidity, air moisture content and temperature were set in the pressure outlet boundary or velocity inlet boundary. 5 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_{w,y} = v_{w} ( {^{{}} }}} )^{0.19}$$ v w , y = v w y h ref 0.19 where v w,y is wind speed measured y meters above the ground, m/s, v w is wind speed measured at reference height, m/s, h ref is the reference height, 10 m. Data statistical method For each grid, regular lattices of 50 m × 50 m  on ground-level are set-up to count the quantity of deposited drifts, such that for each wind direction, each wind speed and each released drifts size, the total deposited mass of drifts are calculated, and afterward the deposition of salt in each bearing and distance can be conducted, as Eq. ( ). 6 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}_{{}} { = }_{j = 1}^{j = N} {}_{{{}_{j} }} }}{{50  50  _{{{}_{j} }} }}  3600  1000  {}_{{{}_{j} }}  {}_{{{}}} }$$ q i = ∑ j = 1 j = N m d j 50 × 50 × ρ d j · 3600 · 1000 · f d j · f wind where q i is the amount of droplets in the i region, mm/h; m dj is the mass of droplets with particle diameter of d j falling in this region, kg/s; ρ dj is the density of droplets with particle diameter of d j , kg/m 3 ; f dj is the frequency of droplets with particle diameter of d j , kg/kg; f wind is the wind frequency of the corresponding direction in the region. 

The operating conditions including loads, meteo and the mass flowrate of the circulate water is assumed stable, the airflow field is assumed steady around the towers. The function of flow movement is shown in reference . The porous media model of Fluent was adopted in the fill area, where the resistance loss coefficient of the fill is shown in Eq. ( ). Discrete Phase Model (DPM) was used in rain zone, where the air resistance to raindrops was analyzed according to Eqs. ( ) to ( ). MERKEL model was used to analyze the heat exchange process between fill and rain zone. 1 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$_{f} =  V^{M} _{a} }}{{ V^{2} H_{f} }}$$ ξ f = 2 A p V M γ a ρ V 2 H f 2 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$_{p} }}{dt} = F_{D} ( { - _{p} } ) + _{x} (_{p} -  )}}{{_{p} }}$$ d u → p dt = F D V → - u → p + g → x ( ρ p - ρ ) ρ p 3 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$F_{D} =  d_{p}^{2} }} {} }}{24}$$ F D = 18 μ ρ p d p 2 C D Re 24 4 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$C_{D} = }}( {1 + b_{1} {}^{{b_{2} }} } ) +  {} }}{{b_{4} + {} }}$$ C D = 24 Re 1 + b 1 Re b 2 + b 3 Re b 4 + Re where u p is particle velocity vector, m/s, V is air velocity vector, m/s, g x is acceleration vector of particle, m/s 2 , ρ p is density of raindrops, kg/m 3 , d p is raindrop diameter, m, C D is raindrop resistance coefficient, and bi is experience coefficient. The commercial computational fluid software FLUENT was used to solve the three-dimensional numerical model of cooling tower. For natural-draft cooling tower, buoyancy is one of the main forces in momentum equation, so buoyancy effect was considered. Air flow, heat transfer, raindrop evaporation and water vapor diffusion need to be solved coupled. The discrete equation is solved by the separation of variables method, and the decoupling of velocity and pressure is performed by the SIMPLEC algorithm. QUICK format was adopted in discretization of velocity and temperature field. The second boundary conditions of outers, the resistance characteristics of the fill, and the evaporation of the raindrop are calculated by the secondary interface of FLUENT software, with user-defined function (UDF) enabled. 

GAMBIT grid division software was used to establish a three-dimensional model, as shown in Fig. , with calculation range of 12500 m × 4500 m × 2500 m, divided into near field and far field, where the near field is divided into inside-tower and outside-tower, the outside-tower is divided into the perimeter and top, the inside-tower is divided into above-throat, below-throat, fill area, bottom of the tower and air inlet area. The grid of tower outside is sparse, and the fill area is the densest, with a total grid of about 8 million. In order to calculate the trajectory of particles in different wind bearings, a set of grids were established with interval of 15° starting from 0° where the wind blows parrallely to the towers row, such that the distribution of salt deposition in all bearings and distances around the tower group could be conducted. 

The computational domain is divided into two parts: inside and outside the cooling tower. For the boundary of the outer region, the bottom is an adiabatic boundary. When there is no wind, the other surfaces are pressure outlet boundaries. In case of ambient wind, the inlet surface is a velocity inlet boundary, and the other surfaces are pressure outlet boundaries. The cooling tower shell is arranged as an adiabatic boundary, and the inlet and outlet of the tower are arranged as the internal boundary. The fill is set to porous medium area, where the resistance coefficient is set according to the experiment experience. Meteorological parameters and operating parameters are shown in Table . The inject surface of the raindrop is the filltop. The initial velocity, temperature and amount of raindrop are set in the DPM. which can be calculated by operating parameters. The motion characteristics of floating droplets under different particle size distribution conditions were studied. The drifts spectrum is shown in Table . In this study, the floating water rate of the water collector is 5 parts per million (5 × 10 –6 ). The motion trajectory of a drift is calculated according to Eqs. ( ) to ( ). The evaporation of drifts is neglected when calculating the salt deposition. The annual or monthly distribution of salt deposition is closely related to the distribution of deposited quantity of drifts at each wind direction and wind speed Fig. shows the annual wind rose diagram, in which 1 m/s represents the wind speed interval of 0 ~ 1 m/s, 2 m/s represents the wind speed interval of 1 ~ 2 m/s, etc. The ambient wind speed is correlated with the height, as shown in Eq. ( ), the wind shear index is related to site and atmospheric stability category, here for example in neutral category is 0.19 . The UDF function is used to setup the wind speed vertical profile. To characterize the influence of ambient humidity, air moisture content and temperature were set in the pressure outlet boundary or velocity inlet boundary. 5 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_{w,y} = v_{w} ( {^{{}} }}} )^{0.19}$$ v w , y = v w y h ref 0.19 where v w,y is wind speed measured y meters above the ground, m/s, v w is wind speed measured at reference height, m/s, h ref is the reference height, 10 m. 

For each grid, regular lattices of 50 m × 50 m  on ground-level are set-up to count the quantity of deposited drifts, such that for each wind direction, each wind speed and each released drifts size, the total deposited mass of drifts are calculated, and afterward the deposition of salt in each bearing and distance can be conducted, as Eq. ( ). 6 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}_{{}} { = }_{j = 1}^{j = N} {}_{{{}_{j} }} }}{{50  50  _{{{}_{j} }} }}  3600  1000  {}_{{{}_{j} }}  {}_{{{}}} }$$ q i = ∑ j = 1 j = N m d j 50 × 50 × ρ d j · 3600 · 1000 · f d j · f wind where q i is the amount of droplets in the i region, mm/h; m dj is the mass of droplets with particle diameter of d j falling in this region, kg/s; ρ dj is the density of droplets with particle diameter of d j , kg/m 3 ; f dj is the frequency of droplets with particle diameter of d j , kg/kg; f wind is the wind frequency of the corresponding direction in the region. 

Flow field Figure shows the flow characteristics of inside and outside of cooling towers when the ambient wind direction is parallel to the tower row. When the wind speed is 1 m/s, the hot air from the first cooling tower outlet on the windward side is obviously deflected, while the hot air from other cooling towers outlets is not obviously deflected. When the wind speed is 3 m/s, the outlet air deflection of the first and second cooling towers on the windward side is large, while the outlet air deflection of other cooling towers is relatively small. When the wind speed is greater than 5 m/s, the air from the cooling tower outlets has a large deflection, while the leeward side is slightly smaller. When the wind direction is parallel to the tower row, the cooling tower on the windward side can play a protection role and reduce the deflection range of hot air from the leeward side cooling tower outlet. Figure shows the pathline of hot and humid air from cooling tower outlet under the wind speed of 5 m/s with different wind directions. As can be seen from the velocity distribution of outlet air flow in Fig. , when the wind direction is parallel to the tower row, the overlapping of hot air at the outlet leads to the larger velocity of drifts. From the trajectories, it can be seen that when the wind direction is parallel to the tower row, are supported by the air flow to move in a higher region. When the wind direction is 45 degrees away from the tower row, the air flow is not superimposed, at this time, the drifts is obviously slower, and are more likely to fall, resulting in salt deposition. When the wind direction is at some certain angle, the hot air from outlet of the six cooling towers does not overlap, and strongly deflects under the act of wind. Drifts trajectory with different particle sizes Figure shows the escape tracks of drifts with different particle sizes. When the partical diameter is small, the drifts barely drop. When the particle diameter is about 350 μm or larger, the escaped drifts drop to the ground within an obviously shorter distance. When the particle diameter is 550 μm or larger, drifts cannot escape from the cooling tower outlet. The interaction between air and raindrop is related to raindrop diameter and relative velocity. If the air velocity inside the tower is low or the diameter of the raindrops is large, the raindrops will not be able to float out of the tower. Theoretically, only when Eq. ( ) is greater than 0, raindrops are able to float out of the tower. Therefore, the relationship between the maximum particle diameter of floating droplets and the wind speed of the fill section is shown in Fig. . When the wind speed of the fill section is 2 m/s, the maximum particle diameter allowing drifts to float out is about 500 μm, which is in good agreement with the above numerical analysis. It is proved that the numerical model can predict the motion trajectory of floating droplets well. Distribution of salt deposition According to a NUREG report, salt deposition of 10 kg/(ha month) or more in any single month of plant-growing seasons may lead to leaf damage of many species. Figure shows the predicted distribution of salt deposition from the group of six simultaneously-operating super-large seawater cooling towers, the deposits mass is accumulated to each month. It can be seen that the highest peak of salt deposition happens in July and is about 8.5 kg/(ha month), located at NNW and about 2 km away from the site. The peak of salt deposition in January and February is the smallest, about 1.2 kg/ (ha month), located around NNW and ENE, 2 km away from the site. It can be seen that the peak of salt deposition is generally larger in summer (June, July and August), followed by spring (March, April and May), and smaller in winter (December, January and February). When salt deposition exceeds 200 kg/(ha·year), it is considered that it will cause serious damage to plants, and preventive measures should be taken. Figure shows the predicted distribution of salt deposits during the year. The annual peak salt deposition is about 38.8 kg/(ha·year), which occurs at NNW, about 2–3 km away from the site. 

Figure shows the flow characteristics of inside and outside of cooling towers when the ambient wind direction is parallel to the tower row. When the wind speed is 1 m/s, the hot air from the first cooling tower outlet on the windward side is obviously deflected, while the hot air from other cooling towers outlets is not obviously deflected. When the wind speed is 3 m/s, the outlet air deflection of the first and second cooling towers on the windward side is large, while the outlet air deflection of other cooling towers is relatively small. When the wind speed is greater than 5 m/s, the air from the cooling tower outlets has a large deflection, while the leeward side is slightly smaller. When the wind direction is parallel to the tower row, the cooling tower on the windward side can play a protection role and reduce the deflection range of hot air from the leeward side cooling tower outlet. Figure shows the pathline of hot and humid air from cooling tower outlet under the wind speed of 5 m/s with different wind directions. As can be seen from the velocity distribution of outlet air flow in Fig. , when the wind direction is parallel to the tower row, the overlapping of hot air at the outlet leads to the larger velocity of drifts. From the trajectories, it can be seen that when the wind direction is parallel to the tower row, are supported by the air flow to move in a higher region. When the wind direction is 45 degrees away from the tower row, the air flow is not superimposed, at this time, the drifts is obviously slower, and are more likely to fall, resulting in salt deposition. When the wind direction is at some certain angle, the hot air from outlet of the six cooling towers does not overlap, and strongly deflects under the act of wind. 

Figure shows the escape tracks of drifts with different particle sizes. When the partical diameter is small, the drifts barely drop. When the particle diameter is about 350 μm or larger, the escaped drifts drop to the ground within an obviously shorter distance. When the particle diameter is 550 μm or larger, drifts cannot escape from the cooling tower outlet. The interaction between air and raindrop is related to raindrop diameter and relative velocity. If the air velocity inside the tower is low or the diameter of the raindrops is large, the raindrops will not be able to float out of the tower. Theoretically, only when Eq. ( ) is greater than 0, raindrops are able to float out of the tower. Therefore, the relationship between the maximum particle diameter of floating droplets and the wind speed of the fill section is shown in Fig. . When the wind speed of the fill section is 2 m/s, the maximum particle diameter allowing drifts to float out is about 500 μm, which is in good agreement with the above numerical analysis. It is proved that the numerical model can predict the motion trajectory of floating droplets well. 

According to a NUREG report, salt deposition of 10 kg/(ha month) or more in any single month of plant-growing seasons may lead to leaf damage of many species. Figure shows the predicted distribution of salt deposition from the group of six simultaneously-operating super-large seawater cooling towers, the deposits mass is accumulated to each month. It can be seen that the highest peak of salt deposition happens in July and is about 8.5 kg/(ha month), located at NNW and about 2 km away from the site. The peak of salt deposition in January and February is the smallest, about 1.2 kg/ (ha month), located around NNW and ENE, 2 km away from the site. It can be seen that the peak of salt deposition is generally larger in summer (June, July and August), followed by spring (March, April and May), and smaller in winter (December, January and February). When salt deposition exceeds 200 kg/(ha·year), it is considered that it will cause serious damage to plants, and preventive measures should be taken. Figure shows the predicted distribution of salt deposits during the year. The annual peak salt deposition is about 38.8 kg/(ha·year), which occurs at NNW, about 2–3 km away from the site. 

A three-dimensional numerical simulation with CFD method was adopted to analyze the air flow and drifts trajectories from a group of six super-large seawater cooling towers. Distribution of drifts and salt depostion was couducted, and the impact of salt depostion to the environment was assessed. Main conclusions are as follows: When the wind direction is parallel to the tower row, the cooling tower on the windward side can play a shielding role, reducing the deflection amplitude of the hot air at the outlet of the cooling tower on the leeward side. The hot air at the outlet overlap to form as one airflow, which enhances dispersion, therefore flattening the deposition peak of drifts and salt. When the particle diameter of drifts is greater than 550 μm, they cannot float out of the cooling tower. Under normal operation conditions of the 6 super-large seawater cooling towers, the salt deposition reaches its maximum in summer, while minimum in winter. The monthly peak of salt deposition is predicted less than 8.5kg/(ha month), which would not cause damages to leaves. The annual peak of salt deposition is predicted 38.8 kg/(ha·year) less than 200 kg/(ha·year), which would not cause serious damage to plants. 

