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

Authors: Hongjie Ling; Zhidong Wang; Shuai An; Guohuai Sun; Yangyue Yan
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

With the frequent occurrence of extreme weather, the rain load of structures under high wind speeds accompanied by heavy rainfall conditions has become a hot issue for research. In this paper, based on the discrete particle model, the raindrop impact load correction formula is proposed. A rapid calculation method for the rain load of structures is formed through experimental verification. And the forecasting program is developed independently to complete the calculation of rain load per unit flat plate under different wind speed and rainfall intensity combination states. The rain load correction coefficientΔCw is defined as the dimensionless coefficient of the rain load per unit area of the flat plate in a fully developed rain field driven by wind speed. The rain load is factored into the wind load calculation formula by way of the correction coefficient ΔCw to form the calculation formula for the rain load of structures. The results show that ΔCw has little correlation with wind speed but is closely related to rainfall intensity R. The calculation formula for ΔCw about rainfall intensity is fitted. When rainfall intensity R = 709.2 mm/h, ΔCw is approximately 22.3%. The research results in this paper provide theoretical and technical support for guiding the rapid calculation of rain load on structures.


---

## Body

## 1. 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. The 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–3] and the separation of the soil [4].The evolution of raindrops in the impact process is the key to explaining the impact force of raindrops. Adler [5] proposed a finite element method to simulate the raindrop impact on a deformable target. The 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 fitted 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. This finding indicated that the impact mechanisms of the bead and the liquid drops were completely different 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. The 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 field 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 difficult to generate and control the wind-rain coupling field in the laboratory environment. There are a few research papers on the quantitative test methods and results of the rain load on structures. Choi [11–13] has made major breakthroughs in the use of numerical simulations employing computational fluid dynamics (CFD) in wind-driven rain (WDR) research, conducted based on a steady-state 3D wind flow 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 fluid dynamics. The 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 difficult to calculate, which made it inconvenient for engineering use. There is an urgent need to develop an efficient and fast rain load forecasting method and form a simplified calculation formula to improve the efficiency of rain load calculation.According to international conventions, the maximum wind force near the center is shown in Table1. 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.Table 1
Typhoon scale.
TypeTropical depressionTropical stormSevere tropical stormTyphoonSevere typhoonSuper typhoonWind level6∼78∼910∼1112∼1314∼15≥16Wind speed (m/s)10.8∼17.117.2∼24.424.5∼32.632.7∼41.441.5∼50.9≥51Table 2
Rainfall intensity scale.
Rainfall typeLight rainModerate rainHeavy rainRainstormHeavy rainstormWeakMediumStrongExtremeDaily rainfall intensity (mm/24 h)0.1–1010–2525–5050–100200Hourly rainfall intensity (mm/h)<2.52.5–88–1616–3264100200709.2The research methods, achievements, and shortcomings of wind-driven rain load by domestic and foreign scholars are summarized in Section1. The 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 verified 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 flat plate under different rainfall intensity and wind speed combination states is carried out. The rain load correction coefficient ΔCw is obtained, and the formula of ΔCw and rainfall intensity R is fitted. The 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.


## 2. 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. The MP DSD was proposed by Marshall and Palmer in 1948. Its form is shown as equation (1):(1)ND=N0exp−ΛD.A shape factorμ is introduced into MP DSD to form Gamma DSD. It is shown as(2)ND=N0Dμexp−ΛD,where N(D) is the number of raindrops in a unit scale interval and unit volume, and the unit is m−3·mm−1; N0 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.Whenμ > 0, the curve bends upward; when μ < 0, the curve bends downward; and when μ = 0, it degenerates to MP DSD.The parameter selection of MP and Gamma DSD is shown in Table3. 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. The 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.Table 3
DSD parameter [20].
DSDμN0m−3mm−1Λmm−1MP080004.1R−0.21Gamma09057R0.1774.37R−0.17631.19×105R−0.3526.78R−0.17661.44×106R−0.8809.16R−0.176Figure 1
Distribution of raindrops with different rainfall intensities. (a)R = 100 mm/h. (b) R = 709.2 mm/h.
(a)(b)


## 3. Rain Load Calculation Method

### 3.1. Computational Domain Space Division and Algorithm Flow

The computational domain is defined by the Cartesian coordinate system, and the computational space is divided by the coordinates of the diagonal points. The detailed implementation steps of the rain load calculation method for structures based on discrete particle model are as follows:(1)
As shown in Figure2, 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 intensityR, an appropriate DSD is selected to generate discrete raindrops with different diameters in the rainfall area (AB) which is generated, and the raindrop field Pini,xi,yi,zi,di,ui,vi,wi is initialized.(3)
The initialized wind fieldUw=UwU,V,W is given in the computational domain boundary (OE), where U, V, and W are the functions of time and space. Since the raindrops are discrete, the spatial occupation is much less than 10%, the influence of the rain field on the wind field is small, so the wind-driven rain one-way coupling algorithm is adopted in this paper.(4)
The 3D point cloud data nodes of the structure are loaded into the program. The 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)
Whent=t+Δt, the discrete particle model of wind-driven rain is solved. The raindrops in the rainfall area move under the combined action of the wind field and gravitational field. 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 intox, y, and z directions, and the load is integrated to obtain the rain load in each direction.Figure 2
Spatial division of computational domain.The detailed calculation process is shown in Figure3.Figure 3
Calculation flowchart.


### 3.2. Selection of DSD Model and Grid Discretization Method of Raindrop Space

The Gamma (3) distribution was chosen for the raindrop distribution function. It can be seen from Figures4 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 field 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.Figure 4
Gamma (3) DSD.Figure 5
Number of raindrops in unit space.According to DSD and rainfall intensity, the rain field is discretized into raindrop particle groups consisting of raindrops with different diameters and different numbers. The total number of raindrop particle groups in unit space is counted asN. 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. The nodes after spatial isometric segmentation are shown in Figure 6.Figure 6
Spatial isometric nodes.Figure 7
Schematic diagram of random movement of spatial nodes.The 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 Figure7, 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−Dmax/2, and Dmax is the maximum raindrop diameter.(2)
We introduce a uniform random numberξ, −1≤ξ≤1. The original spatial node coordinates are randomly moved in the spherical domain to obtain new space node coordinates, as shown in Figure 8.(3)Nodexnew,ynew,znew=Nodex+3ξ⋅r3,y+3ξ⋅r3,z+3ξ⋅r3.(3)
Nonrepeating random integers are generated according to the total number of raindropsN. The raindrops of different diameters and their corresponding quantities are randomly filled in the generated space nodes, as shown in Figure 9.Figure 8
Nodes after random movement.Figure 9
Spatial distribution of raindrops.


### 3.3. Unstructured Mesh Generation of Surfaces Based on Curst Algorithm [22]

Definition 1.
TheVoronoi diagram of a point set H is a collection of polygonal regions that contain only one point h in the point set. The 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.Definition 2.
The graph formed by connecting the points in the adjacent regions (regions with common edges) of theVoronoi diagram is called the Delaunay triangulation.Definition 3.
SupposeI is a point set, and H is the vortex set of Voronoi of I.
The complex surfaces generated by the multimodal function are shown in Figure10. It can be seen that the algorithm realizes the unstructured mesh generation of the complex surfaces.Figure 10
Diagram of unstructured mesh.


### 3.4. 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:(4)md2xdt2=18cdReπμDU−dxdtmd2ydt2=18cDReπμDV−dydtmd2zdt2=−mg1−ρaρw+18cDReπμDW−dzdtwhere Re is the Reynolds number; ρa is the air density, 1.235 kg/m3; ρw is the raindrop density, 1000 kg/m3; μ is the dynamic viscosity coefficient; 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; CD is the drag coefficient of raindrops, related to its diameter D. The research of Van Mook [23] shows that when the diameter of a raindrop is larger than 6 mm, the surface tension of the raindrop cannot maintain its shape, and the phenomenon of fragmentation occurs during the movement.(5)Re=ρaDμU−dxdt2+V−dydt2+W−dzdt2,(6)m=πD36ρw.Substituting equations (5) and (6) into (3), we get(7)d2xdt2=3Cdρa4DρwU−dxdtU−dxdt2+V−dydt2+W−dzdt2,d2ydt2=3Cdρa4DρwU−dydtU−dxdt2+V−dydt2+W−dzdt2,d2zdt2=3Cdρa4DρwU−dzdtU−dxdt2+V−dydt2+W−dzdt2−g1−ρaρw.The acceleration, velocity, and displacement of raindrops can be obtained by solving (8) using self-programming with the finite element method. According to the experimental results of the drag coefficient of equivalent raindrop diameter 0.1∼5.8 mm given by Gunn and Kinzer [24], Equation (8) is obtained by fitting the experimental results. It can be seen from Figure 11 that the fitting results of equation (8) are in good agreement with the experimental values.(8)CD=0.7565−0.0926D+0.0196D2−0.411lnDD.Figure 11
Fitting diagram of raindrop drag coefficient.The distribution of the rain fields driven by the uniform wind with a horizontal wind speed of 30 m/s is shown in Figures12 and 13. It can be seen from Figure 12 that the horizontal terminal velocity of particles with different diameters all reach the wind speed under the drive of the uniform wind field. It can be seen from Figure 13 that the raindrop particles move under the combined action of the wind field and gravitational field. The 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. Therefore, particles with different diameters in the rain field exhibit delamination and diffusion in the vertical direction.Figure 12
Velocity distribution of particles.Figure 13
Diameter distribution of particles.


### 3.5. Near-Wall Collision Particle Search Algorithm

After the unstructured grid of the structural surface is divided, the solid wall point cloud datasetS = {s1, s2, s3, …, sm} is formed by the center positions of m grid nodes, and the corresponding normal vector ψ = {ψ1, ψ2,…, ψm}. The point cloud dataset P = = {p1, p2, p3, …, pn} is composed of n raindrops near the structure, and the corresponding velocity is V = {V1, V2,…, Vn}. 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. The rain load on the structure is obtained by integrating all individual impact loads at the same moment. The specific calculation steps are as follows:(1)
At timet, the maximum raindrop velocity Vmax = {V1,V2, …, Vn} is obtained among all raindrops in the near object area.(2)
As shown in Figure14, the solid wall point set S is taken as the center of the circle, a circle is drawn with the radius R = Vmax∙∆t, and the raindrop point set Pk = {pi, …, pj}, Pk ∈ P is calculated in the corresponding space. If Pk = Φ (empty set), it means that there is no raindrop impact at the current position, k = 1, 2, …, m.(3)
The set of raindrops near thesk solid wall point is solved as Pk. The minimum Euclidean distance is calculated between the raindrops in the set and sk. The raindrop point with the minimum Euclidean distance is determined to collide with the solid wall point.(4)
The rain load is calculated throughFd=∑i=1mFi⋅ψi, where the normal vector is ψi=ψix,ψiy,ψiz.(5)
t = t + ∆t, after updating the computational domain rain field, repeat step (1).Figure 14
Schematic diagram of near-wall collision particle.


### 3.6. Collision Model Modification of Raindrop and Wall

Assume that after the raindrop hits the structure vertically, the velocity in the vertical direction is zero. The impact force of a single spherical raindrop on a solid wall can be calculated by the impulse theorem using size, mass, final velocity, and impact time.(9)∫0τftdt+∫v0mdv=0Fτ=1τ∫0τftdt=mvτ=16τρπD3v,where τ is the impact time.Taking the numerical simulation of impact load with the raindropD = 3.0 mm, v = 50 m/s as an example, the high-speed impact duration curve and pulse form are shown in Figure 15. It can be seen as follows:(1)
The entire impact process is completed in a very short time of 1e − 4 s. The impact is initially characterized as a rigid impact load. After a rapid climb to the maximum value, the raindrop produces a flexible 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 classification of (2), four pulse forms can be given: Rigid body impact, pulse time τ = 0.5D/v, impact load Frig=1/3ρπD2v2; Continuous impact, pulse timeτ = 0.5D/v, impact load Fcon=1/6ρπD2v2; the peak impact load of raindrops are calculated as the impulse load Fp=Cr⋅Frig; The average load of raindrops are used as the impulse load Fm=Cp⋅Fp.Figure 15
Raindrop impact duration curves withD = 3.0 mm, v = 50 m/s.The maximum impact loads are shown in Figure16, generated by raindrops with different diameters impacting the solid wall at the final falling velocity [24]. The correction coefficient Cr = 0.7481, and the corrected results are in good agreement with the numerical calculation results.Figure 16
Statistical diagram of raindrop impact load.


### 3.7. 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 Figure17) and parameters are as follows: (1) high-speed camera (Photron): sampling frequency of 10 kHz, shooting photo resolution of 1024∗1000. (2) piezoelectric sensor with a maximum range of 5 N, a measurement accuracy of 0.5%FS, and a sampling frequency of 128 kHz; (3) pipette with a precision of 0.5 μl/1 μl; (4) dropper; (5) graduated scale; (6) fill light.Figure 17
Test instrumentation. (a) High-speed camera. (b) Piezoelectric sensor. (c) Pipette. (d) Dropper.
(a)(b)(c)(d)The size of the dropper aperture directly governs the size of the raindrop volume. Twelve droppers of different 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 different dropper diameters; (4) slowly rotating the top knob of the pipette, each rotation of 1 frame, discharge a certain volume of liquid V1 (that is, the accuracy of the pipette). The 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 V1; (6) each raindrop is repeated at least 3 times to ensure the repeatability of the generated raindrops.Table 4
Volume calibration table for raindrops.
VolumeV0 (μL)4910121516202225313942Equivalent diameterD (mm)1.972.582.672.843.063.133.373.483.633.904.214.31A schematic diagram of the raindrop drop motion test is given in Figure18, 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 final velocity. The 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. The photo of the test site is shown in Figure 19.Figure 18
Schematic diagram of the test.Figure 19
Testing site.An overlay of the program automatically recognizing the shape and spatial location of raindrops at different moments is shown in Figure20, in which the red dots are filmed raindrops, and the black dots are the raindrops identified by the image. In order to improve the accuracy of the calculation and eliminate the influence of local anomalies, the first 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. The time interval of adjacent photos is 0.0001 s, and the falling time is 0.002 s. The linear regression method of least squares is used to calculate the end velocity of raindrops. The end velocities of raindrops with different diameters falling from different heights are shown in Figure 21. The impact load duration curve recorded by the piezoelectric sensor from the raindrop is shown in Figure 22. The maximum impact load is extracted from the raindrop impact load duration curve. The maximum impact loads Frig and Fp 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 Fp have a more consistent law, and the error is less than 10%. This verifies the correctness of the maximum load correction formula for raindrops Fp=Cr⋅Frig. Therefore, the raindrop pulse load in this paper is selected as Fp.Figure 20
Overlay of raindrop photos taken by high-speed camera.Figure 21
Comparison graph of raindrop end velocity.Figure 22
Impact load duration curve ofH = 93 cm, V0 = 39 μL.Figure 23
Comparison chart of maximum impact load.


## 3.1. Computational Domain Space Division and Algorithm Flow

The computational domain is defined by the Cartesian coordinate system, and the computational space is divided by the coordinates of the diagonal points. The detailed implementation steps of the rain load calculation method for structures based on discrete particle model are as follows:(1)
As shown in Figure2, 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 intensityR, an appropriate DSD is selected to generate discrete raindrops with different diameters in the rainfall area (AB) which is generated, and the raindrop field Pini,xi,yi,zi,di,ui,vi,wi is initialized.(3)
The initialized wind fieldUw=UwU,V,W is given in the computational domain boundary (OE), where U, V, and W are the functions of time and space. Since the raindrops are discrete, the spatial occupation is much less than 10%, the influence of the rain field on the wind field is small, so the wind-driven rain one-way coupling algorithm is adopted in this paper.(4)
The 3D point cloud data nodes of the structure are loaded into the program. The 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)
Whent=t+Δt, the discrete particle model of wind-driven rain is solved. The raindrops in the rainfall area move under the combined action of the wind field and gravitational field. 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 intox, y, and z directions, and the load is integrated to obtain the rain load in each direction.Figure 2
Spatial division of computational domain.The detailed calculation process is shown in Figure3.Figure 3
Calculation flowchart.


## 3.2. Selection of DSD Model and Grid Discretization Method of Raindrop Space

The Gamma (3) distribution was chosen for the raindrop distribution function. It can be seen from Figures4 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 field 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.Figure 4
Gamma (3) DSD.Figure 5
Number of raindrops in unit space.According to DSD and rainfall intensity, the rain field is discretized into raindrop particle groups consisting of raindrops with different diameters and different numbers. The total number of raindrop particle groups in unit space is counted asN. 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. The nodes after spatial isometric segmentation are shown in Figure 6.Figure 6
Spatial isometric nodes.Figure 7
Schematic diagram of random movement of spatial nodes.The 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 Figure7, 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−Dmax/2, and Dmax is the maximum raindrop diameter.(2)
We introduce a uniform random numberξ, −1≤ξ≤1. The original spatial node coordinates are randomly moved in the spherical domain to obtain new space node coordinates, as shown in Figure 8.(3)Nodexnew,ynew,znew=Nodex+3ξ⋅r3,y+3ξ⋅r3,z+3ξ⋅r3.(3)
Nonrepeating random integers are generated according to the total number of raindropsN. The raindrops of different diameters and their corresponding quantities are randomly filled in the generated space nodes, as shown in Figure 9.Figure 8
Nodes after random movement.Figure 9
Spatial distribution of raindrops.


## 3.3. Unstructured Mesh Generation of Surfaces Based on Curst Algorithm [22]

Definition 1.
TheVoronoi diagram of a point set H is a collection of polygonal regions that contain only one point h in the point set. The 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.Definition 2.
The graph formed by connecting the points in the adjacent regions (regions with common edges) of theVoronoi diagram is called the Delaunay triangulation.Definition 3.
SupposeI is a point set, and H is the vortex set of Voronoi of I.
The complex surfaces generated by the multimodal function are shown in Figure10. It can be seen that the algorithm realizes the unstructured mesh generation of the complex surfaces.Figure 10
Diagram of unstructured mesh.


## 3.4. 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:(4)md2xdt2=18cdReπμDU−dxdtmd2ydt2=18cDReπμDV−dydtmd2zdt2=−mg1−ρaρw+18cDReπμDW−dzdtwhere Re is the Reynolds number; ρa is the air density, 1.235 kg/m3; ρw is the raindrop density, 1000 kg/m3; μ is the dynamic viscosity coefficient; 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; CD is the drag coefficient of raindrops, related to its diameter D. The research of Van Mook [23] shows that when the diameter of a raindrop is larger than 6 mm, the surface tension of the raindrop cannot maintain its shape, and the phenomenon of fragmentation occurs during the movement.(5)Re=ρaDμU−dxdt2+V−dydt2+W−dzdt2,(6)m=πD36ρw.Substituting equations (5) and (6) into (3), we get(7)d2xdt2=3Cdρa4DρwU−dxdtU−dxdt2+V−dydt2+W−dzdt2,d2ydt2=3Cdρa4DρwU−dydtU−dxdt2+V−dydt2+W−dzdt2,d2zdt2=3Cdρa4DρwU−dzdtU−dxdt2+V−dydt2+W−dzdt2−g1−ρaρw.The acceleration, velocity, and displacement of raindrops can be obtained by solving (8) using self-programming with the finite element method. According to the experimental results of the drag coefficient of equivalent raindrop diameter 0.1∼5.8 mm given by Gunn and Kinzer [24], Equation (8) is obtained by fitting the experimental results. It can be seen from Figure 11 that the fitting results of equation (8) are in good agreement with the experimental values.(8)CD=0.7565−0.0926D+0.0196D2−0.411lnDD.Figure 11
Fitting diagram of raindrop drag coefficient.The distribution of the rain fields driven by the uniform wind with a horizontal wind speed of 30 m/s is shown in Figures12 and 13. It can be seen from Figure 12 that the horizontal terminal velocity of particles with different diameters all reach the wind speed under the drive of the uniform wind field. It can be seen from Figure 13 that the raindrop particles move under the combined action of the wind field and gravitational field. The 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. Therefore, particles with different diameters in the rain field exhibit delamination and diffusion in the vertical direction.Figure 12
Velocity distribution of particles.Figure 13
Diameter distribution of particles.


## 3.5. Near-Wall Collision Particle Search Algorithm

After the unstructured grid of the structural surface is divided, the solid wall point cloud datasetS = {s1, s2, s3, …, sm} is formed by the center positions of m grid nodes, and the corresponding normal vector ψ = {ψ1, ψ2,…, ψm}. The point cloud dataset P = = {p1, p2, p3, …, pn} is composed of n raindrops near the structure, and the corresponding velocity is V = {V1, V2,…, Vn}. 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. The rain load on the structure is obtained by integrating all individual impact loads at the same moment. The specific calculation steps are as follows:(1)
At timet, the maximum raindrop velocity Vmax = {V1,V2, …, Vn} is obtained among all raindrops in the near object area.(2)
As shown in Figure14, the solid wall point set S is taken as the center of the circle, a circle is drawn with the radius R = Vmax∙∆t, and the raindrop point set Pk = {pi, …, pj}, Pk ∈ P is calculated in the corresponding space. If Pk = Φ (empty set), it means that there is no raindrop impact at the current position, k = 1, 2, …, m.(3)
The set of raindrops near thesk solid wall point is solved as Pk. The minimum Euclidean distance is calculated between the raindrops in the set and sk. The raindrop point with the minimum Euclidean distance is determined to collide with the solid wall point.(4)
The rain load is calculated throughFd=∑i=1mFi⋅ψi, where the normal vector is ψi=ψix,ψiy,ψiz.(5)
t = t + ∆t, after updating the computational domain rain field, repeat step (1).Figure 14
Schematic diagram of near-wall collision particle.


## 3.6. Collision Model Modification of Raindrop and Wall

Assume that after the raindrop hits the structure vertically, the velocity in the vertical direction is zero. The impact force of a single spherical raindrop on a solid wall can be calculated by the impulse theorem using size, mass, final velocity, and impact time.(9)∫0τftdt+∫v0mdv=0Fτ=1τ∫0τftdt=mvτ=16τρπD3v,where τ is the impact time.Taking the numerical simulation of impact load with the raindropD = 3.0 mm, v = 50 m/s as an example, the high-speed impact duration curve and pulse form are shown in Figure 15. It can be seen as follows:(1)
The entire impact process is completed in a very short time of 1e − 4 s. The impact is initially characterized as a rigid impact load. After a rapid climb to the maximum value, the raindrop produces a flexible 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 classification of (2), four pulse forms can be given: Rigid body impact, pulse time τ = 0.5D/v, impact load Frig=1/3ρπD2v2; Continuous impact, pulse timeτ = 0.5D/v, impact load Fcon=1/6ρπD2v2; the peak impact load of raindrops are calculated as the impulse load Fp=Cr⋅Frig; The average load of raindrops are used as the impulse load Fm=Cp⋅Fp.Figure 15
Raindrop impact duration curves withD = 3.0 mm, v = 50 m/s.The maximum impact loads are shown in Figure16, generated by raindrops with different diameters impacting the solid wall at the final falling velocity [24]. The correction coefficient Cr = 0.7481, and the corrected results are in good agreement with the numerical calculation results.Figure 16
Statistical diagram of raindrop impact load.


## 3.7. 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 Figure17) and parameters are as follows: (1) high-speed camera (Photron): sampling frequency of 10 kHz, shooting photo resolution of 1024∗1000. (2) piezoelectric sensor with a maximum range of 5 N, a measurement accuracy of 0.5%FS, and a sampling frequency of 128 kHz; (3) pipette with a precision of 0.5 μl/1 μl; (4) dropper; (5) graduated scale; (6) fill light.Figure 17
Test instrumentation. (a) High-speed camera. (b) Piezoelectric sensor. (c) Pipette. (d) Dropper.
(a)(b)(c)(d)The size of the dropper aperture directly governs the size of the raindrop volume. Twelve droppers of different 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 different dropper diameters; (4) slowly rotating the top knob of the pipette, each rotation of 1 frame, discharge a certain volume of liquid V1 (that is, the accuracy of the pipette). The 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 V1; (6) each raindrop is repeated at least 3 times to ensure the repeatability of the generated raindrops.Table 4
Volume calibration table for raindrops.
VolumeV0 (μL)4910121516202225313942Equivalent diameterD (mm)1.972.582.672.843.063.133.373.483.633.904.214.31A schematic diagram of the raindrop drop motion test is given in Figure18, 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 final velocity. The 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. The photo of the test site is shown in Figure 19.Figure 18
Schematic diagram of the test.Figure 19
Testing site.An overlay of the program automatically recognizing the shape and spatial location of raindrops at different moments is shown in Figure20, in which the red dots are filmed raindrops, and the black dots are the raindrops identified by the image. In order to improve the accuracy of the calculation and eliminate the influence of local anomalies, the first 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. The time interval of adjacent photos is 0.0001 s, and the falling time is 0.002 s. The linear regression method of least squares is used to calculate the end velocity of raindrops. The end velocities of raindrops with different diameters falling from different heights are shown in Figure 21. The impact load duration curve recorded by the piezoelectric sensor from the raindrop is shown in Figure 22. The maximum impact load is extracted from the raindrop impact load duration curve. The maximum impact loads Frig and Fp 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 Fp have a more consistent law, and the error is less than 10%. This verifies the correctness of the maximum load correction formula for raindrops Fp=Cr⋅Frig. Therefore, the raindrop pulse load in this paper is selected as Fp.Figure 20
Overlay of raindrop photos taken by high-speed camera.Figure 21
Comparison graph of raindrop end velocity.Figure 22
Impact load duration curve ofH = 93 cm, V0 = 39 μL.Figure 23
Comparison chart of maximum impact load.


## 4. Analysis of Calculation Results of Unit Plate Rain Load

The rain load forecast program is independently written using MATLAB. The rain load calculation per unit area of the flat plate under different combinations of rainfall intensity and wind speed is carried out.As shown in Figure24, the horizontal distance between the flat plate and the rain field satisfies the conditions for the sufficient development of the rain field driven by the wind field. When the raindrops accelerate to the position of the flat plate, they reach the same speed as the wind speed, as shown in Figure 25. In order to eliminate the phenomenon of stratification and diffusion during the movement of the rain field caused by gravity, the gravity field is ignored in the calculation. It can be seen from Figure 26 that the raindrop particles of different diameters in the rain field are fully mixed without stratification and diffusion. The calculation conditions are shown in Table 5, with a time step of 0.01 s and a duration of 50 s.Figure 24
Distribution of rain load calculation domain per unit area of flat plate.Figure 25
Distribution of rain field velocity.Figure 26
Distribution of particle diameters.Table 5
Calculation conditions.
Research objectRainfall intensityR (mm/h)Wind speedUw (m/s)Condition numberFlat plate (S = 1 m2)10, 50, 100, 200, 300, 400, 500, 600, 709.210∼55,∆v = 590In order to study the statistical characteristics of the flat 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 fitting, as shown in Figure27.(10)fx=12πσe−x−μ2/2σ2,where f(x) is the probability density function; μ is the expected number; σ is the standard deviation.Figure 27
Statistical probability distribution of raindrop impact loads. (a)R = 709.2 mm/h-v = 30 m/s. (b) R = 709.2 mm/h-v = 55 m/s.
(a)(b)All the rain loads for the calculated working conditions are averaged to plot Figure28. It can be seen that the rain load satisfies a quadratic relationship with the wind speed for the same rainfall intensity. The 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 coefficient of variation Cv=σ/u has been introduced. It can be seen from the distribution diagram of the variation coefficient 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 flat plate decreases with the increase of rainfall intensity.Figure 28
Statistical chart of rain loads.Figure 29
Variation coefficient of rain loads.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 unified consideration, and a correction term is introduced, which is defined as follows:(11)ΔCw=2FdρaUw2S,where Fd is the rain load on the structure; S is the windward area.It can be seen from Figure30(a) that ΔCw 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 ΔCw, is used as a variable to fit the rain load correction coefficient formula as follows:(12)ΔCw=a⋅Rb,where a = 0.01206, b = 0.4488, RMSE = 0.003771.Figure 30
Proportion diagram of rain load. (a) Transient pulse load superposition. (b) Momentum average rain load.
(a)(b)The high-rise trip,let tower studied by Wang et al. [25] belongs to a typical large facade configuration. The number of raindrop particles in the large scale range is huge, and the calculation efficiency of commercial software is low. Three rainfall intensities of R = 64 mm/h, R = 100 mm/h and R = 200 mm/h are given by Wang. The flat plate studied in this paper can be considered as the windward side of a rectangular body, which has a similar structural form to the high-rise triplet tower. The calculation results are shown in Figure 30(a). The calculation results in this study are in good agreement with those in the open literature, indicating that the rain load calculation method established in this paper is correct.It is noted that rain load in this paper adopts the transient impact pulse load superposition algorithm, which is different 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 Figure30(b) that the rain load obtained by the quantitative average algorithm is significantly less than the superposition of transient impact load.Therefore, the calculation formula for the combined wind and rain load on the structure can be expressed as(13)F=12Cw+ΔCwρaUw2S,where Cw is the wind resistance coefficient, taken as 1.0.


## 5. Conclusions and Discussions

A rapid calculation method for rain load on structures is proposed and a forecast program is independently developed in this study. Then, the calculation of the rain load per unit area of flat plate and the statistical analysis of the calculation results are completed. The main conclusions are as follows:(1)
The self-developed rain load forecasting program for structures can realize the prediction of wind-driven rain fields and the calculation of rain loads on structures.(2)
The rain load on the structure obeys the Gaussian distribution, and the degree of dispersion decreases with the increase in rainfall intensity.(3)
The 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)
The rain load correction coefficientΔCw is closely related to the rainfall intensity R, but has little to do with the wind speed.The 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. The fluctuating wind will inevitably cause the pulsation of the rain field, and the structural rain load characteristics in the pulsating rain field are also one of the issues worthy of study [26, 27].



---

*Source: 2107987-2022-11-02.xml*