Research Article  Open Access
Wind Load Characteristics and Action Mechanism on Internal and External Surfaces of SuperLarge Cooling Towers under WindRain Combined Effects
Abstract
By focusing on windrain twoway coupling algorithm, simulation iterations of wind field and raindrops in the world highest cooling tower (210m) in northwest China were carried out using continuous phase and discrete phase models based on CFD numerical simulation. Firstly, influence laws of 9 wind velocityrainfall intensity combinations on windinduced rainfall, raindrop additional force, and equivalent pressure coefficient on internal and external surface of the tower body were discussed. On this basis, speed flow line, turbulence energy strength, raindrop running speed, and track on the tower body in the windrain coupling field were disclosed. Finally, qualitative and quantitative contrastive analyses on wind pressure, rain pressure, and equivalent pressure coefficient on internal and external surfaces of the tower body were conducted under different working conditions. Thus, the most unfavorable windrain combination was identified. Calculation formulas of equivalent internal and external pressure coefficients of superlarge cooling towers were fitted from nonlinear least square method. Research results demonstrate that the 3D effect of equivalent internal and external pressure coefficients with considerations to windrain twoway coupling is more prominent. Particularly, there is strong transition on the windward region of the external surface and leeside region at bottom of internal surface. The quantity of caught raindrops on the structural surface is negatively related to wind velocity but is positively related to rainfall intensity. Rain load and rainfall coefficients on the external surface are significantly higher than those on the internal surface. Equivalent internal pressure coefficient has a sharp reduction on the leeside region under different working conditions. Besides, equivalent internal pressure coefficient of different meridians decreases with the increase of height. The maximum and minimum are 0.574 and 0.282, respectively. The proposed equivalent internal and external pressure coefficients of superlarge cooling tower can predict wind load under extreme climate conditions accurately.
1. Introduction
With universal uses of highcapacity and highparameter power units, cooling tower which is one of the important structures of fire/nuclear power plants is developing toward (super) large types [1, 2]. Compared to conventional cooling tower, superlarge cooling tower is a typical wind sensitive structure. Wind load is the control load of internal force design of tower body. Existing structural specifications of cooling towers all neglect impacts of rainfall on tower body. Nevertheless, structures have to bear combined effect of strong wind and rainstorm under extreme climatic conditions. Due to collaborative driving by wind force and gravity, the movement track of raindrops will deviate. Some will impact on external surface of the cooling tower strongly at the leading edge, while some will impact on inner wall strongly through the open tower top, resulting in significant changes of aerodynamic force distribution on surface. Meanwhile, rainstorm will deteriorate turbulence effect of turbulence wind where there is no rain. Under this circumstance, airflows close to tower body and in the tower are complicated, which changes movement track, additional force, and internal and external pressure of raindrops. For this reason, it has important theoretical significance and engineering value to study wind load distribution characteristics on internal and external surfaces of superlarge cooling tower under different wind velocities and rainfall intensities as well as the corresponding mechanism of action.
Researches concerning wind load of superlarge cooling towers mainly focus on stochastic behaviors of single wind load [3], extreme wind pressure distribution [4], static interference effects [5], and dynamic interference effect [6]. Chinese code [7] offers the mean wind pressure coefficient distribution curve on external surface of hyperbolic cooling tower and internal pressure coefficient of single tower (0.5), but it neglects variations of internal attractive force along height and meridians completely. Germany code [8] gives the mean pressure coefficient distribution pattern on external surface according to external surface roughness but still overlooks effects of 3D characteristics and rainfall on internal pressure coefficient. However, Shen et al. [9] simulated wind load on internal surface of cooling towers under coldhot air selfcirculating system and external wind field by using the realizable kε turbulence model and multiphase flow model. They found that internal pressure which was produced in external wind field changes dramatically with height and attitude. Zou et al. [10] analyzed 3D effects of wind pressure distribution in internal surface of tower body under different ventilation rates (100% and 30%) by rigid pressure test, finding nonuniform distribution of wind pressure on internal surface along height and circumferential directions. Associated results could solve values of wind pressures on internal and external surfaces of conventional cooling towers. Nevertheless, there are hardly researches about aerodynamic performances on internal and external surfaces of superlarge cooling towers under windrain combined effect, not to mention qualitative and quantitative comparisons of internal and external pressures under different wind velocityrainfall intensity combinations.
In addition, some researches [11–19] concerning windrain combined effect mainly focus on lowrise buildings, bridge cable, power transmission tower, and wind turbine. They only consider effects of windinduced rainfall but overlook reaction of rain to wind. Currently, there are rare researches that have discussed whether windrain coupling affects aerodynamic performances of internal and external surface of superlarge cooling towers under poor climate conditions.
In this study, by focusing on windrain twoway coupling algorithm, simulation iterations of wind field and raindrops in the world highest cooling tower (210m) in northwest China were carried out using continuous phase and discrete phase models based on CFD numerical simulation. On this basis, mechanism of action of internal and external pressures of the tower body under windrain combined effect was studied. Influence laws of different wind velocityrainfall intensity combinations on windinduced rainfall, raindrop additional force, and equivalent pressure coefficient on internal and external surface of the tower body were discussed. Finally, the most unfavorable windrain combination was identified. The calculation formulas of equivalent internal and external pressure coefficients of superlarge cooling towers were fitted based on nonlinear least square method.
2. WindRain TwoWay Coupling Algorithm
2.1. Rainfall Intensity
Classifications of rainfall intensity are shown in Table 1 [20–22]. Two rainfall intensity levels have different sampling times. Measurement results of the same rainfall differ significantly. Structure checking based on hourly rainfall intensity is safer than that based on daily rain intensity. Hourly rainfall intensity can intuitively reflect influences of instantaneous rainfall intensity on structural performances under extreme climate conditions which are mostly concerned in engineering.

In addition, the main purpose of this paper is to compare the influence of different wind speedrain intensity combinations on wind load characteristics and action mechanism on internal and external surfaces of superlarge cooling towers. Considering the possibility that wind speed and rainfall both reach a certain intensity under windrain combined effects, this paper chooses hourly rainfall intensity, rather than the rainfall intensity for minor durations with a very low probability of occurrence, and adopts the data based on the downpour grades under extreme conditions which may occur in China.
2.2. Raindrop Spectral Distribution
Rainfall spectrum approximately obeys negative exponential distribution. Common models [23] have Best Spectrum, MarshallPalmer Spectrum, and Gamma Spectrum. In this study, MarshallPalmer Spectrum was used:where D_{p} is the diameter of raindrop (mm), n(D_{p}) is concentration spectrum of raindrops with different diameters, N_{0} is concentration (=8,000), and λ is a scale parameter. λ can be expressed aswhere R is rainfall intensity (mm/h).
2.3. Final Velocity of Raindrops
When the falling height is over 20m, raindrops of almost all diameters can reach the final velocity [24]. Raindrops of diameter smaller than 2mm fall approximately as balls, but larger raindrops will suffer large deformations due to air resistance. The vertical final velocity of raindrops can be calculated by the following equation [25, 26]:where v(D_{p}) is the vertical final velocity of raindrops with a diameter of D_{p} (m/s).
2.4. WindRain TwoWay Coupling Algorithm
During downpours, the volume fraction of raindrops in air is far smaller than 10%. Therefore, raindrop can be simulated by DPM model. Simulated raindrops are added to the continuous phase as the second phase after stabilization of wind field computing for windrain twoway coupling [27, 28]. The motion equilibrium equation of raindrops in wind field iswhere is discrete phase particle velocity, is continuousphase fluid velocity, () is the drag force of unit particle mass, ρ_{p} and ρ are particle and fluid densities, and is interaction force between discrete phase and continuous phase.where μ is viscosity coefficient of fluid, D_{p} is particle diameter, and Re is relative Reynolds number.
After discrete phase of raindrops is considered, the basic control equation of wind continuous phase iswhere S_{m} is the second discrete phase mass which is added to the continuous phase, p is pressure, is tensor of stress, and is gravity. Among them, tensor of stress can be expressed aswhere I is unit tensor and the second term in the right is volume expansion.
2.5. Wall Collision Equation
The impacting process of raindrops onto tower walls complies with the law of conservation of momentum. The key to solve impact force lies in impact time. Evaporation, splash, and rupture of raindrops during the impacting process are neglected. It is believed the raindropstructure interaction follows Newton’s second law. According to theorem of momentum,where f(t) is impact force vector of single raindrop (N) and v is velocity vector of raindrop.
The impact force of raindrop onto the structure in unit time (F(τ)) is
If the falling raindrops are viewed as balls approximately, it has
Since raindrop diameter is generally lower than 6mm and the horizontal final velocity before the impact is relatively high, for calculation simplification, the impact time τ is determined:
Then, the impact force of raindrops to structure can be simplified as
3. Brief Introduction to the Project and Working Condition Setting
3.1. Brief Introduction to the Project
The superlarge reinforced concrete hyperbola indirect aircooling tower with natural ventilation built is 210m high. Elevation of the throat is 157.5m and elevation of the inlet mouth is 32.5m. Inner diameters of the throat and inlet mouth are 110m and 159m, respectively. The zero diameter is 180m. The tower body employs exponential thickening. The minimum thickness (0.37m) lies in the throat section and the maximum thickness (2.0m) is at the lower ring beam. The tower body is connected with circumferential plate foundation by 52 pairs of Xshaped columns. Xshaped columns have rectangular section (1.8m×1.2m). The circumferential foundation is a castinplace reinforced concrete structure: 12.0m (width)×2.5m (height). Main structural size and picture of the cooling tower are shown in Table 2.

3.2. Working Conditions
This cooling tower locates in the Btype landform. With considerations to conventional working state of louver of cooling tower, the open effect of louver was considered according to 30% of ventilation rate [1, 2]. Influences of combinations of three wind velocities and three rain intensities on aerodynamic performance of the cooling tower were discussed. Weak wind, moderate wind, and strong wind were divided according to the maximum wind velocities in the return period of 10 years, 50 years, and 100 years. Rainfall intensity was classified into weak downpour, moderate downpour, and heavy downpour. Therefore, a total of 9 working conditions could be set (Figure 1).
4. Numerical Simulation of WindRain TwoWay Coupling
4.1. Establish the WindRain Field Model
To ensure the cooling tower in the rainfall region and full development of wake flow, dimension of calculation domain was set 3000m (wind length)1500m (wind width)600m (height). The bottom tower center was chosen as the origin of the coordinate system and the following wind direction is consistent with Xaxis. To consider both computing efficiency and accuracy, the cooling tower was divided into local and periphery windrain fields during gridding. The local windrain field contained the cooling tower model and was meshed by unstructured grid. The periphery windrain field had a regular shape and was meshed by highquality structured grid. The minimum grid size in the core region was 0.2m. The total number of grids in the overall model was over 18 millions and grid mass was higher than 0.40 (it is required to be higher than 0.1 in order to prevent negative volume). Both number and quality of grids can meet the calculation requirements. Overall gridding and local gridding are shown in Figure 2.
(a) Overall gridding
(b) Local gridding
Inlet of the calculation domain is the boundary of velocity inlet and the outlet is boundary of pressure outlet. Two side walls and top surface use symmetry boundaries and the cooling tower and ground adopt wall boundary. The overlapping face of local and periphery calculation domains are set the interface. The calculation domain of the windrain field and its boundary conditions are shown in Figure 3.
4.2. WindRain Field Coupling Calculation
Windrain twoway coupling simulation theory and discrete phase tracking method are applied to this superlarge Reynolds number structure, which proposes extremely high requirements on memory capacity of computer. Therefore, numerical calculation in this study is accomplished based on large compute server in the wind turbine aerodynamics highperformance computing center of the project group (Figure 4). The server processer used Intel(R) Xeon(R) CPU E52650 v3 @ 2.30GHz 2.30GHz (2 processor) and the installed memory was 256GB. The 64bit operating system was employed.
(a) Calculation center
(b) Calculation server
Numerical calculation used 3D single precision and separated solver. Flow velocity of fluid field was used as the absolute velocity. Air model was equivalent to the ideal incompressible fluid and the turbulence model chose the kω shearstress (SST) control equation. Inlet of the calculation domain used the wind profile model with an exponential index of 0.15. Wind velocities at 10m over the ground were set as three standard wind velocities in Section 3.2. The flow field was solved through coupling between seed and pressure by using the SIMPLEC algorithm. The solving formulation of convective term was the second order. During the calculation process, grid slope correction was set in order to increase hybrid grid computing effect. Calculation residue of the control equation was set as 1×1. Finally, the wind field was initialized for iterative computations. Simulation curves and theoretical curves of mean wind velocity and turbulence profile are shown in Figure 5. Results indicate that simulated results of mean wind velocity and turbulence profile conform to theoretical values, indicating that simulation standards of wind field meet engineering requirements.
(a) v0=20m/s
(b) v0=23.7m/s
(c) v0=25.3m/s
After stabilization of wind field solving, the discrete phase was added for iterative operation of windrain field coupling. Here, raindrops with six different diameters in the range of 1.0~6.0mm were used to simulate rainfall of continuous diameter distribution (Table 3). Number and volume occupancy of raindrops with different diameters were determined by MarshallPalmer spectrum in Section 2.2. Raindrops were “plane” released at the horizontal velocity of 0. Driven by wind, they kept consistent movement velocity with horizontal wind velocity in the location. The vertical release velocity was 5m/s. The collaborative effects of gravity and resistance make raindrops reach the final velocity in (3) to the maximum height range.

4.3. Validity Verification
The mean wind pressure coefficients in the throat area on external surface of the cooling tower under three wind velocities were compared with existing Chinese and foreign relevant codes [7, 8] and measurement curves [29] (Figure 6). According to analysis, the angle between negative pressure extreme point and separation point on the mean wind pressure distribution curve in the throat area conforms to that of measurement curve under all three wind velocities. The negative pressure value in the leeside region is slightly higher than fire engineering codes and measured values but is lower than Germany code. Other numerical values of angles basically agree with codes and measurement results. To sum up, numerical simulation results are valid to some extent.
5. Contrastive Analysis of Results
5.1. Analysis of Wind Field
Threedimensional wind velocity flow lines and vorticity distribution of the cooling tower under three standard wind velocities before raindrops are added are shown in Figures 7 and 8, respectively. Clearly, the incoming flow runs through the tower body and splits in the windward side. It is accelerated along two external side walls of the tower body and forms different sizes of vortexes on the leeside surface. Some airflow enters into the tower body through louver. Meanwhile, inflows at tower top skim over and change partial rising airflow directions. Airflow generates reflux and forms a largesized vortex at the throat area. With the increase of wind speed, airflow runs at an accelerated velocity, which intensifies vortex falloff, dense distribution of wind velocity flow lines, and interaction between rising airflows and incoming flows at the tower top. Turbulence energy intensity is positively related to wind velocity, especially in middle and lower parts of the leeside surface. Peak mainly lies in top of windward region and throat area of the leeside region.
(a) v0=20m/s
(b) v0=23.7m/s
(c) v0=25.3m/s
(a) v0=20m/s
(b) v0=23.7m/s
(c) v0=25.3m/s
Threedimensional distribution patterns of external and internal pressure coefficients under 9 working conditions are shown in Figures 9 and 10.
(a) Working condition 1
(b) Working condition 2
(c) Working condition 3
(d) Working condition 4
(e) Working condition 5
(f) Working condition 6
(g) Working condition 7
(h) Working condition 8
(i) Working condition 9
(a) Working condition 1
(b) Working condition 2
(c) Working condition 3
(d) Working condition 4
(e) Working condition 5
(f) Working condition 6
(g) Working condition 7
(h) Working condition 8
(i) Working condition 9
Distribution pattern of external and internal pressure coefficients under all working conditions is basically consistent with numerical values. They are all highly symmetric around the wind axis.
External pressure coefficients on different crosssections of tower body under all working conditions all present a Vshaped variation trend with the expansion of circumferential angle and become stable close to the leeside region. Numerical values of external pressure coefficient are relatively stable close to windward region and leeside region except for gentle fluctuation in the crosswind region.
Influenced by airflow which runs through louver, the internal pressure coefficient at tower bottom jumps strongly along circumferential direction. The numerical value at the leeside surface is significantly lower than those in other regions. With the gradual growth of internal airflow from tower bottom, the circumferential random impact of tower body from middle and upper airflows is relatively weak. Therefore, the circumferential numerical value of internal pressure coefficient at middle and upper parts of the tower body basically remains constant.
5.2. Analysis of Rain Field
Raindrop movement tracking based on resultant velocity of particles in windrain field under 9 working conditions is shown in Figure 11. Raindrop density was processed by equal proportional coarsening. The following can be seen from Figure 11:
(a) Working condition 1
(b) Working condition 2
(c) Working condition 3
(d) Working condition 4
(e) Working condition 5
(f) Working condition 6
(g) Working condition 7
(h) Working condition 8
(i) Working condition 9
Raindrops on near wall of the cooling tower show complicated movement patterns. Some raindrops impact on the windward surface of the cooling tower under wind actions, while the remaining raindrops separate at two structural sides along the airflow. Only few adhere onto side walls of the cooling tower and most enter into the wake flow with winds.
Raindrops at upper front end of the cooling tower run through the outlet of the cooling tower and enter into the inside region and impact on middle and upper leeside regions of internal surface at a great velocity. Moreover, the stronger rainfall intensity is, the more raindrops caught by the internal wall will be.
With the increase of wind velocity, horizontal acting force of raindrops increases dramatically, which makes raindrops move along the wind direction at accelerating speeds. Tremendous raindrops skim over the tower top and move to rear of the tower body. As a result, number of raindrops which enter into the cooling tower drop sharply.
Impacts of number of raindrops with different diameters, velocity, and occupancy of speed on external and internal surfaces of the cooling tower under 9 working conditions are shown in Figures 12 and 13.
(a) Number of raindrops
(b) Occupancy of speed
(c) Horizontal velocity
(a) Number of raindrops
(b) Occupancy of speed
(c) Horizontal velocity
Diameter of raindrops caught by internal and external surfaces of the cooling tower mainly ranges within 3~6mm. Raindrops with 5mm diameter occupy the highest proportion. This is because, driven by the same wind intensity, movement speed of smaller raindrops increases quickly. Before reaching the height range of the cooling tower, raindrops have skimmed over the tower body in horizontal direction and entered into the wake flow region.
Number of raindrops caught by internal and external surfaces decreases with the increase of wind velocity. With the increase of rainfall intensity, number of raindrops on external surface is far higher than (about 10 times) that on the internal surface.
Occupancy of horizontal speed of raindrops which are caught by external surface shows basically the same distribution laws. It increases firstly and then decreases with the increase of final horizontal velocity. The mean horizontal velocity of raindrops which are caught by internal and external surfaces is significantly lower than the standard wind velocity. Moreover, the horizontal impact velocity of raindrops is positively related to diameter of raindrops.
Rain loads on internal and external surfaces of the cooling tower under 9 working conditions are calculated from (14). Meanwhile, the rainwind load ratios at different heights of internal and external surfaces of the cooling tower are given. Results are shown in Figure 14 and Tables 4–6.



Under all working conditions, rain load on external surface of the cooling tower decreases firstly and then increases with the increase of height, reaching the valley at 0.23~0.31H and the peak at tower bottom or tower top.
Under all working conditions, rain load on internal surface of the cooling tower increases with the increase of height. There is no raindrop at places below 0.69H and most rain loads are concentrated in the height range of 0.90~1.0H. The ratio is about 95%.
Rainwind load ratio at different heights is very small. The highest rain load is only 3.655 of wind load, which occurs in the height range of 0.90~1.0H under working condition 3. Besides, rain load on external surface is significantly higher than that on internal surface under all working conditions.
Rain load in internal and external surface under different wind velocities climbs up significantly with the increase of rainfall intensity. Given the fixed rainfall intensity, rain loads on internal surface as well as middle and lower parts of external surface decline as wind velocity increases. In middle and lower parts of the cooling tower, number of raindrops influences rain load more significantly than wind velocity. However, accelerating movement of raindrops on middle and upper parts of the cooling tower increases rain load dramatically.
To display location and number of raindrops at typical heights and angles as well as corresponding pressure coefficients more clearly, 3D distributions of raindrops and raininduced pressure coefficients on internal and external surfaces under all working conditions are shown in Figures 15 and 16.
(a) Working condition 1
(b) Working condition 2
(c) Working condition 3
(d) Working condition 4
(e) Working condition 5
(f) Working condition 6
(g) Working condition 7
(h) Working condition 8
(i) Working condition 9
(a) Working condition 1
(b) Working condition 2
(c) Working condition 3
(d) Working condition 4
(e) Working condition 5
(f) Working condition 6
(g) Working condition 7
(h) Working condition 8
(i) Working condition 9
Under all working conditions, raindrop impact positions mainly concentrate on windward region of external surface and leeside region of internal surface. Driven by airflow vortexes, there are few raindrops adhered on leeside region of external surface and windward region of internal surface.
There are more raindrops caught on internal and external surface of the cooling tower under working condition 3 than the remaining working conditions. Number of caught raindrops drops sharply with the increase of wind velocity and increases gradually with the increase of rainfall intensity. Meanwhile, number of caught raindrops on external surface is far higher than that on internal surface.
Under all working conditions, raininduced external pressure coefficient mainly concentrates on 60° range of two sides of the windward surface, and it is basically 0 in the remaining ranges. The maximum value is 0.184, which is observed in the range of 0.15~0.23H under working condition 3.
Under all working conditions, raininduced internal pressure coefficient mainly concentrates in the meridian range of 0.9H~1.0H and 90° range of two sides of the leeside region. It is basically 0 in the remaining ranges.
5.3. Analysis of Equivalent Pressure Coefficient
For quantitative comparison of wind and raininduced pressure distribution on the cooling tower under different working conditions, the equivalent internal and external pressure coefficients are defined. Specific calculation steps are as follows: ① rain loads at different monitoring points of internal and external surfaces are converted into raininduced pressure. ② Calculate the ratios between rain pressure at monitoring points and wind pressure at the corresponding reference heights, that is, raininduced internal and external pressure coefficients. ③ Calculate the sum of vectors of raininduced coefficient and wind pressure coefficient, that is, equivalent internal and external pressure coefficients under windrain collaborative effects.
Comparison curves of equivalent internal and external pressure coefficients of 6 internal and external typical sections under different working conditions are exhibited in Figures 17 and 18.
(a) 0.27H (56m)
(b) 0.51H (104m)
(c) 0.73H (152m)
(d) 0.80H (168m)
(e) 0.87H(184m)
(f) 0.95H (200m)
(a) 0.27H (56m)
(b) 0.51H (104m)
(c) 0.73H (152m)
(d) 0.80H (168m)
(e) 0.87H (184m)
(f) 0.95H (200m)
Distribution laws and numerical values of equivalent external pressure coefficients of the section at the same height are consistent under different working conditions except for small differences between windward region and leeside region. The external pressure coefficient decreases firstly, increases, decreases, and finally stabilizes from the windward region to leeside region.
Under the same working condition, there is a small difference in numerical values of equivalent external pressure coefficients on 6 typical sections. However, they are all highly symmetric around the wind axis. The maximum negative pressure presents an inverted Vshaped variation with the increase of height. The negative pressure close to the throat area is about 1.5.
Threedimensional distribution effect of internal pressure coefficient is remarkable. Under the same working condition, equivalent internal pressure coefficients on 6 typical sections are different, but they are all symmetric around the wind axis. The circumferential distribution laws of equivalent internal pressure coefficients are basically consistent under different working conditions.
Equivalent internal pressure coefficients of all sections decrease in the leeside region of the cooling tower (which is in the blue dotted framework). This is mainly because airflow runs through the louver and impacts on internal surface of the leeside region. Influenced by the incoming wind, raindrops will impact on upper internal surface of the cooling tower. Moreover, raindrops influence internal pressure coefficient more than the incoming wind.
Equivalent internal and external pressure coefficient curves of four typical meridians (0°, 75°, 120°, and 180°) of the tower body are shown in Figures 19 and 20.
(a) 0°
(b) 75°
(c) 120°
(d) 180°
(a) 0°
(b) 75°
(c) 120°
(d) 180°
The equivalent external pressure coefficient of the 0° meridian is basically stable, while those of the 75° and 180° meridians increase firstly and then decrease with the increase of height. The equivalent external pressure coefficient of the 120° meridian remains basically constant with the increase of height, but the numerical value decreases above the throat area.
Distribution trends of internal pressure coefficients of all four meridians are basically consistent. They all decrease gradually with the increase of height. The maximum numerical value is 0.574, which is at inlet height of 75° meridian. The minimum numerical value is 0.282, which is at the tower top of 180° meridian.
Compared to other meridians, the equivalent internal pressure coefficient of internal surface bottom of the 180° meridian decreases sharply, while the equivalent internal pressure coefficient of the tower top declines slightly. This is mainly because bottom equivalent internal pressure coefficient is influenced by impacts of airflow that runs through the louver on internal wall of the cooling tower, whereas the top equivalent internal pressure coefficient is caused by raindrop impact load on internal surface of the leeside region.
5.4. TwoDimensional Fitting Surface of Equivalent Pressure Coefficient
For the convenience of engineering designers to get equivalent internal and external pressure coefficients of superlarge cooling towers, the calculation formula of equivalent external pressure coefficient is fitted based on the principle of nonlinear least square method by using meridian height and circumferential angle as the objective function. The cooling tower is divided into N_{1} sections uniformly along the circumferential direction and divided into N_{2} sections uniformly along the meridian direction. Let N=N_{1}N_{2}. The specific definition of the formula iswhere I is the N1 matrix which covers only elements valuing 1, is the N1 matrix that circulates by N_{2} times that using N_{1} angles as the circulation unit, Z is the N1 matrix that circulates by N_{1} times that using N_{2} same heights as the circulation unit, is multiplying corresponding elements of the matrix, is dividing corresponding elements of the matrix, is the n power of corresponding matrix elements, exp () is the matrix of exponent that uses each element of the matrix in the return brackets e as the base number, and is the N1 matrix that uses wind pressure coefficient of N_{1} circumferential angles as the unit and changes by N_{2} times along the meridian height. a_{i}, b_{i}……o_{i} (i=1, 2,…9) are fitting coefficients. The fitting coefficients under the most unfavorable working condition (condition 3) are shown in Table 7.

The comparison curves between twodimensional actual values and fitting values of equivalent external pressure coefficients under the most unfavorable working condition (condition 3) are shown in Figure 21. Numerical values of scattered points are actual equivalent external pressure coefficients of the cooling tower. The corresponding numerical values of the hook surface are equivalent external pressure coefficients which are gained by twodimensional fitting formula. It can be seen from Figure 21 that equivalent external pressure coefficient of the cooling tower has obvious twodimensional features along the meridian and circumferential directions. Viewed from overall distribution of fitting hook surface, the variation laws along meridian and circumferential directions basically conform to distribution of actual equivalent external pressure coefficient except for slight difference (<10%) in the region of maximum negative pressure. The comparison results demonstrate that the proposed twodimensional fitting formula can be used as the valuing reference of equivalent external pressure coefficient of superlarge cooling towers.
Similarly, the calculation formula of equivalent internal pressure coefficient of superlarge cooling towers is fitted:where Z is meridian height of the tower body, is circumferential angle of tower body, is fitting value of equivalent internal pressure coefficient that uses Z and as the target function, and a_{i}, b_{i}……o_{i} (i=1, 2,…9) are fitting coefficients. The fitting coefficients under the most unfavorable working condition (condition 3) are shown in Table 8.

Comparison between twodimensional actual values and fitting values of equivalent internal pressure coefficients of cooling tower under the most unfavorable working condition (condition 3) is shown in Figure 22. Numerical values of scattering point are actual equivalent internal coefficients of the cooling tower. Corresponding numerical values of the hook surface are equivalent internal pressure coefficients which are fitted by the twodimensional fitting formula. It can be seen from Figure 22 that equivalent internal pressure coefficient of the cooling tower decreases gradually with the increase of height, but the numerical value in the bottom leeside region decreases significantly. According to overall distribution of fitting hook surface, the variation laws of meridian direction and circumferential direction comply basically with actual equivalent internal pressure coefficients. The numerical values envelop real numerical values and the error rate is smaller than 10%. The comparison reveals that the proposed twodimensional fitting formula can be used as the valuing reference for equivalent internal pressure coefficient of superlarge cooling towers.
6. Conclusions
(i)Wind field and wind pressure coefficients of cooling towers have distinct threedimensional effect before the rain field is added. Horizontal force of raindrops which are added subsequently increases with the growth of wind velocity. Tremendous raindrops skim over the tower top and separate at two sides of the tower body and finally enter into the wake flow region. Only few raindrops adhere on external surface and enter into the tower and impact on the internal wall with airflow.(ii)Number of raindrops which are caught by internal and external surfaces is negatively correlated with wind velocity but is positively related with rainfall intensity. The number of raindrops caught by internal and external surfaces under working condition 3 (wind velocity=20m/s and rainfall intensity=200mm/h) is higher than those under the remaining working conditions. Moreover, raindrops with 5mm diameter occupy the highest proportion on internal and external surfaces. Horizontal velocity of raindrops on external surface is mainly in the range of 3~10m/s and the horizontal velocity of raindrops on internal surface is mainly in the range of 1~3m/s. Number of caught raindrops on external surface is far higher than that on internal surface.(iii)Under different working conditions, raindrop impact positions mainly concentrate in the range of 60° at two sides of the external windward surface and the range of 90° at two sides of internal leeside region near the top of cooling tower. Rain load and raininduced pressure coefficient of external surface are significantly higher than those of internal surface. The maximum raininduced external pressure coefficient is 0.184, which is observed near the bottom of cooling tower under working condition 3. The maximum raininduced internal pressure coefficient is 0.0038, which is observed near the top of cooling tower under working condition 3.(iv)Under different working conditions, equivalent external pressure coefficient at section of the same height differs slightly in the windward region and leeside region. The maximum negative pressure values and negative pressure values of the leeside region increase firstly and then decrease with the increase of height. The maximum negative pressure in the throat area is about 1.5. Equivalent internal pressure coefficient under different working conditions drops dramatically in leeside region. Equivalent internal pressure coefficient of different meridians decreases gradually with the increase of height. The maximum and minimum are 0.574 and 0.282, respectively.(v)The twodimensional fitting formulas of equivalent internal and external pressure coefficients which are proposed based on nonlinear least square method can predict load value under similar extreme climate conditions. The maximum error rate is controlled lower than 10%.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project is jointly supported by the National Natural Science Foundation (U1733129; 51761165022; 51208254), Jiangsu Province Outstanding Natural Science Foundation (BK20160083), China Postdoctoral Science Foundation (2013M530255; 1202006B), and Jiangsu Qinglan Project.
References
 S. T. Ke, Y. J. Ge, L. Zhao, and Y. Tamura, “Stability and Reinforcement Analysis of Superlarge Exhaust Cooling Towers Based on a Wind Tunnel Test,” Journal of Structural Engineering, vol. 141, no. 12, p. 04015066, 2015. View at: Publisher Site  Google Scholar
 S. T. Ke, J. Liang, L. Zhao, and Y. J. Ge, “Influence of ventilation rate on the aerodynamic interference for two IDCTs by CFD,” Wind and Structures, An International Journal, vol. 20, no. 3, pp. 449–468, 2015. View at: Google Scholar
 J.F. Zhang, Y.J. Ge, and L. Zhao, “Influence of latitude wind pressure distribution on the responses of hyperbolodial cooling tower shell,” Wind and Structures, An International Journal, vol. 16, no. 6, pp. 579–601, 2013. View at: Publisher Site  Google Scholar
 X. Chen, L. Zhao, S. Cao et al., “Extreme Wind Loads on SuperLarge Cooling Towers,” Journal of the International Association for Shell & Spatial Structures, vol. 57, pp. 49–58, 2016. View at: Google Scholar
 H.J. Niemann and H.D. Köpper, “Influence of adjacent buildings on wind effects on cooling towers,” Engineering Structures, vol. 20, no. 10, pp. 874–880, 1998. View at: Publisher Site  Google Scholar
 S. T. Ke and Y. J. Ge, “The influence of selfexcited forces on wind loads and wind effects for superlarge cooling towers,” Journal of Wind Engineering & Industrial Aerodynamics, vol. 132, pp. 125–135, 2014. View at: Publisher Site  Google Scholar
 GB/T 501022014, Code for Design of Cooling for Industrial Recirculating Water , The Ministry of Construction of China, Beijing, China, 2014 (Chinese).
 VGBR610Ue, VGBGuideline: Structural Design of Cooling TowerTechnical Guideline for The Structural Design, Computation And Execution of Cooling Towers, Essen: BTR Bautechnik Bei Kuhlturmen, 2005.
 G. Shen H, C. Zhang S, and N. Sun B, “Numerical simulation of wind load on inner surface of large hyperbolic cooling tower,” Journal of HaErBin Institute of Technology, vol. 43, no. 4, pp. 104–108, 2011. View at: Google Scholar
 Y. F. Zou, H. W. Niu, and Z. Q. Chen, “Threedimensional effect of wind load on the single tower of the special large cooling tower, and its design value,” Journal of HuNanUniversity, vol. 32, no. 1, pp. 76–82, 2015. View at: Google Scholar
 X. Fu, H.N. Li, and G. Li, “Fragility analysis and estimation of collapse status for transmission tower subjected to wind and rain loads,” Structural Safety, vol. 58, pp. 1–10, 2016. View at: Publisher Site  Google Scholar
 Z. Wang, Y. Zhao, F. Li, and J. Jiang, “Extreme dynamic responses of mwlevel wind turbine tower in the strong typhoon considering windrain loads,” Mathematical Problems in Engineering, vol. 2013, Article ID 512530, 13 pages, 2013. View at: Publisher Site  Google Scholar
 D. Xin, H. Li, L. Wang, and J. Ou, “Experimental study on static characteristics of the bridge deck section under simultaneous actions of wind and rain,” Journal of Wind Engineering & Industrial Aerodynamics, vol. 107108, pp. 17–27, 2012. View at: Publisher Site  Google Scholar
 M. A. Goudarzi and S. R. SabbaghYazdi, “Effects of modeling strategy on computational wind pressure distribution around the cooling tower's,” Wind and Structures, An International Journal, vol. 14, no. 1, pp. 81–84, 2011. View at: Publisher Site  Google Scholar
 M. Bennett, V. Kodakalla, and V. Gupta, “Vibration mitigation measures in cable stayed bridges,” Journal of Molecular Structure, vol. 996, no. 13, pp. 64–68, 2011. View at: Google Scholar
 L. Y. Wang and Y. L. Xu, “Active stiffness control of windraininduced vibration of prototype stay cable,” International Journal for Numerical Methods in Engineering, vol. 74, no. 1, pp. 80–100, 2010. View at: Google Scholar
 B. Blocken, G. Dezsö, J. van Beeck, and J. Carmeliet, “Comparison of calculation models for winddriven rain deposition on building facades,” Atmospheric Environment, vol. 44, no. 14, pp. 1714–1725, 2010. View at: Publisher Site  Google Scholar
 Q.C. Zhang, W.Y. Li, and W. Wang, “Static bifurcation of rainwindinduced vibration of stay cable,” Wuli Xuebao/Acta Physica Sinica, vol. 59, no. 2, pp. 729–734, 2010. View at: Google Scholar
 B. Blocken and J. Carmeliet, “A review of winddriven rain research in building science,” Journal of Wind Engineering & Industrial Aerodynamics, vol. 92, no. 13, pp. 1079–1130, 2004. View at: Publisher Site  Google Scholar
 GB/T 285922012, Grade of Precipitation , Standards Press of China, Beijing, China, 2012 (Chinese).
 F. Tauro, J. Selker, N. van de Giesen et al., “Measurements and Observations in the XXI century (MOXXI): innovation and multidisciplinarity to sense the hydrological cycle,” Hydrological Sciences Journal, vol. 63, no. 2, pp. 169–196, 2018. View at: Publisher Site  Google Scholar
 X. P. Wu, Numerical Study on The Effects of Wind And Rain on LowRise Buildings, Zhejiang University, 2008.
 G. M. Mcfarquhar and R. List, “The Raindrop Mean Free Path and Collision Rate Dependence on Rainrate for ThreePeak Equilibrium and MarshallPalmer Distributions,” Journal of the Atmospheric Sciences, vol. 48, no. 48, pp. 1999–2004, 2010. View at: Google Scholar
 R. Gunn and G. D. Kinzer, “The terminal fall velocity for water droplets in stagnant air,” Journal of the Atmospheric Sciences, vol. 6, no. 4, pp. 243–248, 1949. View at: Google Scholar
 J. S. Marshall and W. M. K. Palmer, “The distribution of raindrops with size,” Journal of the Atmospheric Sciences, vol. 5, no. 4, pp. 165166, 1948. View at: Publisher Site  Google Scholar
 E. C. Rigby, J. S. Marshall, and W. Hitschfeld, “The development of the size distribution of raindrops during their fall,” Journal of the Atmospheric Sciences, vol. 11, no. 5, pp. 362–372, 2010. View at: Google Scholar
 S. Liu, S. H. Huang, and Q. S. Li, “3D numerical simulation of winddriven rain on bridge deck sections based on eulerian multiphase model,” Engineering Machanics, vol. 34, no. 4, pp. 63–71, 2017. View at: Google Scholar
 E. C. Douvi and D. P. Margaris, “Aerodynamic performance investigation under the influence of heavy rain of a NACA 0012 airfoil for wind turbine applications,” International Review of Mechanical Engineering, vol. 6, no. 6, pp. 1228–1235, 2012. View at: Google Scholar
 T. Sun F and M. Zhou L, “Without ribs the elliptic wind pressure distribution of the cooling tower full size measurement and wind tunnel study,” Journal of Air Dynamics, vol. 12, no. 4, pp. 12–17, 1983. View at: Google Scholar
Copyright
Copyright © 2018 Shitang Ke et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.