Research on Calculation Method of Rain Load on Structures Based on Discrete Particle Model

,


Introduction
Rainfall is a common natural phenomenon. In severe weather, such as typhoons and thunderstorms, raindrops achieve high movement speeds driven by high-speed wind. Raindrops form a rain load when they collide with the structure's surface at high speeds. Te high-speed impact of raindrops on the surface of the building produces large amounts of rain pressure, which has a great impact on the erosion of the building surface [1][2][3] and the separation of the soil [4].
Te evolution of raindrops in the impact process is the key to explaining the impact force of raindrops. Adler [5] proposed a fnite element method to simulate the raindrop impact on a deformable target. Te impact time was so short that load sensors generally could not achieve the necessary sampling frequency to accurately measure the impinging force. Nearing et al. [6] conducted an experiment to measure the raindrop impinging force using a piezoelectric transducer and derived a formula with ftted parameters for the mean peak force. Soto et al. [7] measured the raindrop impinging force using a piezoelectric quartz transducer and simultaneously captured the impinging process using a high-speed camera. Grinspan and Gnanamoorthy [8] found that the output signal of a transducer during bead and liquid drop impact tests was essentially the same when the mass of the bead was approximately 8.5 times than oil droplets and the height of droplet was higher than the bead. Tis fnding indicated that the impact mechanisms of the bead and the liquid drops were completely diferent and that the impact evolution of a liquid drop was complex. Abuku et al. [9] investigated the impact, absorption, and evaporation of raindrops on building facades. Te measurements showed that large drops with high impact speeds splash, whereas drops with small impact speeds and small impact angles bounce. Huang et al. [10] simulated the raindrop impinging process based on the Navier-Stokes equation but did not address the interaction between the substrate and raindrop.
Structures placed in the wind and rain feld are subjected to the combined action of wind load and rain load, and the proportion of rain load is obviously weaker than that of wind load. It is difcult to generate and control the wind-rain coupling feld in the laboratory environment. Tere are a few research papers on the quantitative test methods and results of the rain load on structures. Choi [11][12][13] has made major breakthroughs in the use of numerical simulations employing computational fuid dynamics (CFD) in winddriven rain (WDR) research, conducted based on a steadystate 3D wind fow pattern. And this method is used in the present study. Ke et al. [14] and Fu et al. [15,16] carried out studies on the numerical prediction of rain loads on structures such as high-rise buildings, power transmission towers, and cooling towers using continuous-phase and discrete-phase models of computational fuid dynamics. Te results showed that the rain loads were closely related to the rainfall intensity, and the proportion of rain loads relative to wind loads increased continuously as the rainfall intensity increased. However, the CFD method was time-consuming and difcult to calculate, which made it inconvenient for engineering use. Tere is an urgent need to develop an efcient and fast rain load forecasting method and form a simplifed calculation formula to improve the efciency of rain load calculation.
According to international conventions, the maximum wind force near the center is shown in Table 1. Rainfall intensity is an important index to describe the intensity of rainfall, which is represented by I. It is usually expressed by the rainfall of statistical time lengths such as year, month, day, and hour, so the units are mm/24 h, mm/h, etc. Rainfall intensity grades are shown in Table 2, according to 24-hour rainfall and 1-hour rainfall, which can be divided into light rain, moderate rain, heavy rain, rainstorm, and heavy rainstorm.
Te research methods, achievements, and shortcomings of wind-driven rain load by domestic and foreign scholars are summarized in Section 1. Te raindrop spectrum distribution function and characteristics are given in Section 2. In Section 3, the algorithms involved in rain load calculation based on the discrete particle model are introduced in detail. A rain droplet-wall collision correction model is proposed, and its correctness is verifed through experiments. A fast rain load forecasting method is formed, and a rain load forecasting program is written using MATLAB. In Section 4, based on the research in Section 3, the fast forecasting of rain load per unit fat plate under diferent rainfall intensity and wind speed combination states is carried out. Te rain load correction coefcient ΔC w is obtained, and the formula of ΔC w and rainfall intensity R is ftted. Te calculation results are in good agreement with the published literature. In Section 5, typical conclusions of the present research results are given, and the shortcomings of the study and the subsequent research directions are discussed.

Raindrop Size Distribution (DSD)
Raindrop size can be described by MP DSD [17], lognormal distribution [18], Gamma distribution [19] and other functions, among which MP DSD and Gamma DSD are the most commonly used. Te MP DSD was proposed by Marshall and Palmer in 1948. Its form is shown as equation (1): (1) A shape factor μ is introduced into MP DSD to form Gamma DSD. It is shown as where N(D) is the number of raindrops in a unit scale interval and unit volume, and the unit is m − 3 ·mm − 1 ; N 0 is concentration, and the unit is m − 3 ·mm − 1 ; Λ is the scale parameter, related to rainfall intensity R (mmh − 1 ), and the unit is mm − 1 ; D is the diameter of raindrop, and its unit is mm.
Te parameter selection of MP and Gamma DSD is shown in Table 3. As can be seen from Figure 1, compared with MP DSD, due to the existence of the shape factor μ, Gamma DSD can more accurately express the characteristics of small-diameter raindrops, which are more common than large-diameter raindrops. Research [21] shows that when μ � 3, Gamma DSD is in good agreement with actual rainfall. Te calculation results of R � 100 mm/h and R � 709.2 mm/h rainfall are shown in Figure 1. Other rainfalls are not given, but the results are similar to Figure 1

Computational Domain Space Division and Algorithm
Flow. Te computational domain is defned by the Cartesian coordinate system, and the computational space is divided by the coordinates of the diagonal points. Te detailed implementation steps of the rain load calculation method for structures based on discrete particle model are as follows: (1) As shown in Figure 2, the computational domain is divided into four areas: (i) computational domain boundary (OE); (ii) rainfall area (AB); (iii) structure area; (iv) near-object area (CD). (2) According to rainfall intensity R, an appropriate DSD is selected to generate discrete raindrops with diferent diameters in the rainfall area (AB) which is generated, and the raindrop feld Since the raindrops are discrete, the spatial occupation is much less than 10%, the infuence of the rain feld on the wind feld is small, so the winddriven rain one-way coupling algorithm is adopted in this paper. (4) Te 3D point cloud data nodes of the structure are loaded into the program. Te surface unstructured meshes are generated based on the Curst algorithm; then, the coordinates of the mesh nodes on the surface of the structure and the normal vectors at the center point of the mesh are calculated.
(5) When t � t + Δt, the discrete particle model of winddriven rain is solved. Te raindrops in the rainfall area move under the combined action of the wind feld and gravitational feld. If raindrops move beyond the rainfall area, the excess raindrops are supplemented in the rainfall area and move randomly without overlapping. (6) It is determined whether the movement of raindrops exceeds the computational domain boundary, and raindrops beyond boundary disappear. It is determined whether raindrops in the computational domain reach the near object region (CD), and raindrops entering the near object region participate in the near-wall particle search. It is determined whether raindrops collide with structure; if a collision occurs, solve the collision model, then raindrops disappear after the collision. (7) When the raindrops collide with the structure, the spatial position and the impact load of the raindrops at the current moment are statistically counted. According to the normal vectors at different positions of structure, the raindrop impact load is decomposed into x, y, and z directions, and the load is integrated to obtain the rain load in each direction.

Selection of DSD Model and Grid Discretization Method of
Raindrop Space. Te Gamma (3) distribution was chosen for the raindrop distribution function. It can be seen from Figures 4 and 5 that the number of raindrops per unit space increases with the increase of rainfall intensity R, and the size of raindrops also transitions from small diameter to large diameter. However, the diameter of raindrops in the rain feld is generally less than 5.0 mm. Even if under extreme rainfall conditions with a rainfall intensity R � 709.2 mm/h, the number of raindrops larger than 5.0 mm is extremely small, and the distribution of raindrop diameter is mainly concentrated in the interval of 0.5 mm-3.0 mm.
According to DSD and rainfall intensity, the rain feld is discretized into raindrop particle groups consisting of raindrops with diferent diameters and diferent numbers. Te total number of raindrop particle groups in unit space is counted as N. In order to ensure the number of nodes generated in the space N′ ≥ N, the space is divided into equal distances in the x, y, and z directions. If it is not an integer, it is rounded backwards. Te nodes after spatial isometric segmentation are shown in Figure 6.
Te spatial distribution of raindrops has the characteristics of randomness, disorder, and nonoverlapping. In order to realize the spatial distribution of raindrops to meet the requirements of "three properties," the random movement criteria of the spatial nodes are introduced as follows: (1) As shown in Figure 7, on the basis of the generated ordered spatial nodes, the coordinates of each node are taken as the center of the sphere, the radius of the sphere is r � Δd − D max /2, and D max is the maximum raindrop diameter. (2) We introduce a uniform random number ξ, − 1 ≤ ξ ≤ 1. Te original spatial node coordinates are randomly moved in the spherical domain to obtain new space node coordinates, as shown in Figure 8. (3) Nonrepeating random integers are generated according to the total number of raindrops N. Te raindrops of diferent diameters and their corresponding quantities are randomly flled in the generated space nodes, as shown in Figure 9.

Unstructured Mesh Generation of Surfaces Based on Curst
Algorithm [22] Defnition 1. Te Voronoi diagram of a point set H is a collection of polygonal regions that contain only one point h in the point set. Te distance from any position to h in the region is shorter than the distance from this position to all other points in the point set.   Figure 10. It can be seen that the algorithm realizes the unstructured mesh generation of the complex surfaces.

Mathematical Model of Raindrop Movement in Wind
Field. Choi [12] proposed the motion equation of raindrop particles in wind-driven rain, as shown in the following equation: where R e is the Reynolds number; ρ a is the air density, 1.235 kg/m 3 ; ρw is the raindrop density, 1000 kg/m 3 ; μ is the dynamic viscosity coefcient; g is the gravitational acceleration; U, V, and W are the components of wind speed in x, y, and z directions; m is the mass of raindrop; C D is the drag coefcient of raindrops, related to its diameter D. Te research of Van Mook [23] shows that when the diameter of a raindrop is larger than 6 mm, the surface tension of the     Advances in Civil Engineering raindrop cannot maintain its shape, and the phenomenon of fragmentation occurs during the movement.
Substituting equations (5) and (6) into (3), we get Te acceleration, velocity, and displacement of raindrops can be obtained by solving (8) using self-programming with the fnite element method. According to the experimental results of the drag coefcient of equivalent raindrop diameter 0.1∼5.8 mm given by Gunn and Kinzer [24], Equation (8) is obtained by ftting the experimental results. It can be seen from Figure 11 that the ftting results of equation (8) are in good agreement with the experimental values.
Te distribution of the rain felds driven by the uniform wind with a horizontal wind speed of 30 m/s is shown in Figures 12 and 13. It can be seen from Figure 12 that the horizontal terminal velocity of particles with diferent diameters all reach the wind speed under the drive of the uniform wind feld. It can be seen from Figure 13 that the raindrop particles move under the combined action of the wind feld and gravitational feld. Te small-diameter raindrops accelerate rapidly in the horizontal direction and slowly in the vertical direction, located at the top. While the large-diameter raindrops accelerate slowly in the horizontal direction and rapidly in the vertical direction, located at the bottom. Terefore, particles with diferent diameters in the rain feld exhibit delamination and difusion in the vertical direction. Advances in Civil Engineering near the structure, and the corresponding velocity is V � {V 1 , V 2 ,. . ., V n }. Rain drops acting on the surface of the structure produce an instantaneous impact load. A projection is made to the three normal vectors at this impact location. Te rain load on the structure is obtained by integrating all individual impact loads at the same moment. Te specifc calculation steps are as follows:

Near-Wall Collision
(1) At time t, the maximum raindrop velocity V max � {V 1 ,V 2 , . . ., V n } is obtained among all raindrops in the near object area. (2) As shown in Figure 14, the solid wall point set S is taken as the center of the circle, a circle is drawn with the radius R � V max •∆t, and the raindrop point set P k � {p i , . . ., p j }, P k ∈ P is calculated in the corresponding space. If P k � Φ (empty set), it means that there is no raindrop impact at the current position, k � 1, 2, . . ., m. (3) Te set of raindrops near the s k solid wall point is solved as P k . Te minimum Euclidean distance is calculated between the raindrops in the set and s k . Te raindrop point with the minimum Euclidean distance is determined to collide with the solid wall point.

Collision Model Modifcation of Raindrop and Wall.
Assume that after the raindrop hits the structure vertically, the velocity in the vertical direction is zero. Te impact force of a single spherical raindrop on a solid wall can be calculated by the impulse theorem using size, mass, fnal velocity, and impact time.
where τ is the impact time.

Advances in Civil Engineering
Taking the numerical simulation of impact load with the raindrop D � 3.0 mm, v � 50 m/s as an example, the highspeed impact duration curve and pulse form are shown in Figure 15. It can be seen as follows: (1) Te entire impact process is completed in a very short time of 1e − 4 s. Te impact is initially characterized as a rigid impact load. After a rapid climb to the maximum value, the raindrop produces a fexible deformation, and the load appears to fall slowly. (2) It is divided into two types by way of impulse load: given the pulse time, the impact load is calculated; given the impact load, the pulse time is calculated. (3) According to the classifcation of (2), four pulse forms can be given: Rigid body impact, pulse time τ � 0.5D/v, impact load F rig � 1/3ρπD 2 v 2 ; Continuous impact, pulse timeτ � 0.5D/v, impact load F con � 1/6ρπD 2 v 2 ; the peak impact load of raindrops are calculated as the impulse load F p � C r · F rig ; Te average load of raindrops are used as the impulse load F m � C p · F p .
Te maximum impact loads are shown in Figure 16, generated by raindrops with diferent diameters impacting the solid wall at the fnal falling velocity [24]. Te correction coefcient C r � 0.7481, and the corrected results are in good agreement with the numerical calculation results.

Single Raindrop Impact Test.
In order to realize the test of the spatial falling motion of raindrops and the whole process of the impact load acting on the surface of the structure, the instruments and equipment used in the test (see Figure 17) and parameters are as follows: (1) high-speed camera (Photron): sampling frequency of 10 kHz, shooting   Te size of the dropper aperture directly governs the size of the raindrop volume. Twelve droppers of diferent diameters were selected for the test. After calibration, the volume of raindrops was obtained from 4 μl to 42 μl, with equivalent diameters of 1.97 mm to 4.31 mm, as detailed in Table 4. Droplet volume calibration method: (1) choose an appropriate pipette; (2) pipette aspirate the volume of test liquid; (3) pipette with diferent dropper diameters; (4) slowly rotating the top knob of the pipette, each rotation of 1 frame, discharge a certain volume of liquid V 1 (that is, the accuracy of the pipette). Te liquid gradually converges into raindrops at the mouth of the dropper by surface tension; (5) when the droplets happen to drop, record the number of rotating frames K, that is, the volume of raindrops V � K V 1 ; (6) each raindrop is repeated at least 3 times to ensure the repeatability of the generated raindrops.
A schematic diagram of the raindrop drop motion test is given in Figure 18, where the dropper of the pipette drops calibrated raindrops from two heights of H � 0.93 m and H � 1.86 m, respectively. Under the joint action of gravity, air buoyancy, and surface tension raindrops fall and impact on the piezoelectric sensor with a certain fnal velocity. Te piezoelectric sensor is used to record the entire duration of the impact load. High-speed photography technology is used to record the process before the impact and the impact crushing process. Te photo of the test site is shown in Figure 19.
An overlay of the program automatically recognizing the shape and spatial location of raindrops at diferent moments is shown in Figure 20, in which the red dots are flmed raindrops, and the black dots are the raindrops identifed by the image. In order to improve the accuracy of the calculation and eliminate the infuence of local anomalies, the frst 20 photos of the raindrop impact structure are selected, and the distance l from its center point to the top of the photo is calculated. Te time interval of adjacent photos is 0.0001 s, and the falling time is 0.002 s. Te linear regression method of least squares is used to calculate the end velocity of raindrops. Te end velocities of raindrops with diferent diameters falling from diferent heights are shown in Figure 21. Te impact load duration curve recorded by the piezoelectric sensor from the raindrop is shown in Figure 22. Te maximum impact load is extracted from the raindrop impact load duration curve. Te maximum impact loads F rig and F p are calculated from the raindrop mass and the end velocity before impact using (9) and plotted to form Figure 23. It can be seen that the test results and the theoretical calculation F p have a more consistent law, and the error is less than 10%. Tis verifes the correctness of the maximum load correction formula for raindrops F p � C r · F rig . Terefore, the raindrop pulse load in this paper is selected as F p .

Analysis of Calculation Results of Unit Plate Rain Load
Te rain load forecast program is independently written using MATLAB. Te rain load calculation per unit area of the fat plate under diferent combinations of rainfall intensity and wind speed is carried out. As shown in Figure 24, the horizontal distance between the fat plate and the rain feld satisfes the conditions for the sufcient development of the rain feld driven by the wind feld. When the raindrops accelerate to the position of the fat plate, they reach the same speed as the wind speed, as shown in Figure 25. In order to eliminate the phenomenon  Advances in Civil Engineering of stratifcation and difusion during the movement of the rain feld caused by gravity, the gravity feld is ignored in the calculation. It can be seen from Figure 26 that the raindrop particles of diferent diameters in the rain feld are fully mixed without stratifcation and difusion. Te calculation conditions are shown in Table 5, with a time step of 0.01 s and a duration of 50 s.
In order to study the statistical characteristics of the fat plate subjected to the impact load of raindrops, 4800 sample points are selected, and the interval is divided into 100 parts. After probability statistics, normalization was carried out, and it was found that it obeys the Gaussian distribution. A Gaussian distribution function is used for ftting, as shown in Figure 27.
where f(x) is the probability density function; μ is the expected number; σ is the standard deviation. All the rain loads for the calculated working conditions are averaged to plot Figure 28. It can be seen that the rain load satisfes a quadratic relationship with the wind speed for the same rainfall intensity. Te rain load increases with the increase of rainfall intensity for the same wind speed. At the same wind speed, the rain load increases with the increase in rainfall intensity. In order to describe the dispersion degree of rain load, the coefcient of variation C v � σ/u has been introduced. It can be seen from the distribution diagram of the variation coefcient in Figure 29 that due to the increase in rainfall intensity, the concentration of raindrops per unit volume increases, and the load acting on the plate is more uniform. As a consequence, the degree of dispersion of the rain load on the fat plate decreases with the increase of rainfall intensity.
In order to simplify the calculation and facilitate the use of the engineering community, the raindrop impact load is integrated into the wind load for unifed consideration, and a correction term is introduced, which is defned as follows: where F d is the rain load on the structure; S is the windward area. It can be seen from Figure 30(a) that ΔC w has little correlation with wind speed and is strongly related to rainfall intensity R. In order to simplify the formula and improve its  practicability, the rainfall intensity R, which plays a leading role in ΔC w , is used as a variable to ft the rain load correction coefcient formula as follows:      It is noted that rain load in this paper adopts the transient impact pulse load superposition algorithm, which is diferent from the long-term momentum averaging algorithm. Since the rain load is a typical discrete problem, the momentum averaging ignores the pulse characteristics of the rain load. It can be seen from Figure 30(b) that the rain load obtained by the quantitative average algorithm is      14 Advances in Civil Engineering signifcantly less than the superposition of transient impact load. Terefore, the calculation formula for the combined wind and rain load on the structure can be expressed as where C w is the wind resistance coefcient, taken as 1.0.

Conclusions and Discussions
A rapid calculation method for rain load on structures is proposed and a forecast program is independently developed in this study. Ten, the calculation of the rain load per unit area of fat plate and the statistical analysis of the calculation results are completed. Te main conclusions are as follows: (1) Te self-developed rain load forecasting program for structures can realize the prediction of wind-driven rain felds and the calculation of rain loads on structures.
(2) Te rain load on the structure obeys the Gaussian distribution, and the degree of dispersion decreases with the increase in rainfall intensity. (3) Te calculation formula for including the rain load into the wind load in the form of a correction factor is proposed, and the calculation formula for the correction factor is given. (4) Te rain load correction coefcient ΔC w is closely related to the rainfall intensity R, but has little to do with the wind speed.
Te rain load calculation method established in this paper is the ultimate load superimposed by the transient rain load, which is suitable for the calculation of the local transient ultimate load and rain pressure of the structure. If the low-frequency response problem of rain load to the motion of marine structures is studied, then the long-term momentum averaging algorithm should be used for rain load. Since rain load is a typical discrete problem, the momentum averaging algorithm ignores the impulse characteristics of rain load. How to reasonably increase the impulse characteristics on the basis of the momentum averaging algorithm requires further in-depth research. Te fuctuating wind will inevitably cause the pulsation of the rain feld, and the structural rain load characteristics in the pulsating rain feld are also one of the issues worthy of study [26,27].

Data Availability
Te relevant data designed in this paper are calculated by self-programming according to the numerical algorithm in this paper. At the end of the paper, the corresponding reference 20 is applied in the reference, and the data are from Rain pressure distribution of three•power connected tall buildings with wind•rain coupled environment. J. Nanj. U. Aeron. Astron., 52(05), 825-834."

Conflicts of Interest
Te authors declare that they have no conficts of interest.