Numerical Simulation of Wind-Driven Rain Based on the Eulerian Model

The wind-driven rain model based on the Eulerian multiphase model is the most widely used in engineering applications of green buildings. In this study, a real-world case of wind-driven rain was simulated using a wind-driven rain solver based on the Eulerian multiphase model developed by OpenFOAM. Separation of flow was observed at the front edges of the top of the building in the numerical simulation results, and the separated flow lines reconnected in the middle region of the roof. This type of corner flow with separation produces an accelerated horizontal airflow and increases the catch ratio; hence, there is often a larger catch ratio at the windward facades near the roof. The simulation results show that the Eulerian model used in this study possesses adequate capability and accuracy for wind-driven rain simulation.


INTRODUCTION
Sustainable utilization of water resource has been a urgent topic around world (Zhang et al., 2021a;Zhang et al., 2021b). Wind-driven rain is a natural phenomenon in which raindrops fall with a nonzero horizontal velocity due to the influence of the wind field to form oblique rain. This oblique rain affects the water and heat exchange properties and material durability of building facades when acting on buildings and is a significant source of moisture for building facades. Its destructive consequences can take many forms, including moisture accumulation in porous materials leading to water seepage, frost damage, moisture-induced salt migration, weathering and discoloration, and cracking of structures due to thermal and moisture gradients. Runoff following wind-driven rain is also a contributing factor to the appearance of surface contamination patterns on facades. Moreover, it has been noted in the literature that the combination of wind and rain has a non-negligible effect on the load effect on the building structure as additional rain loads are generated by raindrops striking building facades . Therefore, the study of the wind-driven rain phenomenon is of great significance to the performance design of buildings. Raindrops are constantly in a state of non-equilibrium motion due to gravity and viscous drag during their descent as well as the variation of wind speed, resulting in uncertainty in the velocity of raindrops; this is one of the complexities of the wind-driven rain problem.
Early wind-driven rain studies were based on a semi-empirical approach combining a theoretically derived form with field measurements (Lacy, 1971;Rydock, et al., 2005). The standard meteorological data measured at weather stations were wind speed, wind direction, and horizontal rainfall, and wind-driven rainfall was not usually measured. Therefore, researchers have established semi-empirical relationships between wind-driven rain loads and meteorological parameters (wind speed, wind direction, and horizontal rainfall). The development of the relationships has been guided by experimental observations in which both free wind-driven rain and wind-driven rain on buildings increase in approximate proportion to wind speed and horizontal rainfall. These formulae are semi-empirical in that there is a theoretical basis for their form, but free parameters are chosen to fit the experimental data. The many limitations of this approach, including the difficulties in fully revealing the intrinsic complexity of the wind-driven rain problem associated with using experimental and semiempirical methods, have led to the increased use of numerical methods in recent research efforts. Choi (1993) introduced the first Euler-Lagrange method to solve the raindrop trajectory and conducted a wind-driven rain study on two-dimensional and three-dimensional buildings. He proposed a numerical method to simulate wind-driven rain around buildings. This model is made up of two components. The flow around the building is calculated by solving k − ε of the two-equation turbulence model. The motion of the raindrops in the wind field is calculated by considering the forces acting on the droplets, then deriving the raindrop trajectories by solving the equations of motion. Wind-driven rain around rectangular buildings with different aspect ratios was investigated. The raindrop trajectories for various sizes of droplets were determined, and their velocities and directions when they struck the buildings at different locations were calculated. The amount of rainfall that fell on different parts of building facades and the corresponding intensities were obtained. The effects of wind speed, rainfall intensity, and building geometry on wind-driven rain were also investigated. In a subsequent study, Choi (1999) proposed another method to calculate the impact of wind-driven rain on building facades. This method enables the quantitative numerical simulation of wind-driven rain. Blocken and Carmeliet (2000a) and Blocken and Carmeliet (2000b) further developed the work of Choi by examining the transient variation of wind-driven rain over time and space during rainfall while considering meteorological conditions that could vary at any time. They incorporated numerical techniques for wind-driven rain simulation into a method for estimating wind-driven rain loads on a practical new building envelope based on buildings, which yields accurate estimates of spatial and temporal distribution, as well as meteorological data for the building site. This numerical method was also applied to several arithmetic sequences from a validated low-level test. The calculated results showed that the local distribution of driven rain loads on the building envelope can be accurately predicted with this numerical method. Ooi et al. (2020) used Lagrangian particle tracking methods to assess the effectiveness of different mitigation methods under the wind-driven rain. The research showed that the effectiveness of horizontal and vertical barriers was dependent on the size of the raindrops, with vertical barriers being more effective for smaller sized raindrops and horizontal barriers being more effective for larger sized raindrops.
The discrete model in the Euler-Lagrangian system proposed by previous scholars was once the main method for wind-driven rain studies. However, the disadvantages of this method are obvious: only the relative values at the physical level can be obtained, various field maps such as velocity and pressure fields of low buildings in continuous phase wind fields cannot be obtained similarly, and the pre-processing calculation and post-processing work of the simulation are time-consuming. Huang and Li (2010) conducted the first numerical simulation study of buildings under the coupling effect of wind and rain using a multiphase model in the Eulerian-Eulerian system. They proposed a new numerical simulation method for wind-driven rain based on the Eulerian multiphase model and verified its accuracy for wind-driven rain coupled with the measured data of windward facades of isolated low-rise buildings reported by Blocken et al. (2013). This model divides the raindrops into several different phases according to their particle size, and, when combined with the wind phase, flow control equations can be established for each phase in the model. The coupling effect of the interacting forces between wind and rain and the gravitational effect of raindrops are incorporated into the model in the form of a term at the right end of the equation, thereby simplifying the setting of the model's boundary conditions. Kubilay et al. (2013) simulated the wind-driven rain field of a high-rise tower based on the Eulerian multiphase model, then compared and analyzed the calculated results with the field measured data to validate the model, and found that the calculation time required by the Eulerian multiphase model was reduced about 10 times compared with that of the earlier Lagrangian model. Compared with the earlier Eulerian-Lagrangian method, the biggest advantage of wind-driven rain simulation using the Eulerian multiphase model is that the setting of rain phase boundary conditions is simplified; only the volume fraction and the corresponding velocity components of the rain phase are required to be set at the boundary, and the wind-driven rain parameters are relatively easier to obtain (Kubilay et al., 2015a;Kubilay et al., 2015b;Pettersson, 2016). This model also has better method robustness, and the computational time can be reduced significantly. Wang et al. (2020) analyzed the wind-driven rain distribution characteristics of low-rise buildings under typhoons using the Eulerian multiphase model and field measurement, verifying the reliability of the model and the important effects of wind speed, wind direction, and horizontal rainfall intensity. He also found that the wind-driven rain effect received by the windward facade of the building under extreme meteorological conditions is higher than that of the building roof. In recent years, the numerical simulation method based on the Eulerian multiphase model has become extensively used in the field of wind-driven rain research.
In this study, the wind-driven rain solver based on the Eulerian multiphase model developed by OpenFOAM was used to analyze the model solution process and applied to the simulation of a wind-driven rain instance. The numerical results were then compared with experimentally measured results to verify the accuracy and effectiveness of the model and to lay the foundation for future research and development.

CONTROL EQUATIONS FOR WIND-DRIVEN RAIN
In this study, the three-dimensional steady-state Reynoldsaveraged Navier-Stokes (RANS) equations were used for the wind field control equations, and the renormalization group (RNG) k − ε model was used for the turbulence model (Yakhot et al., 1992). In the Eulerian multiphase model, the rain phase is simulated based on the continuum assumption, and the wind and rain phases are considered to be unidirectionally coupled, i.e., only the effect of the wind phase on the rain phase is accounted for; the effect of the rain phase on the wind phase is ignored. Every rain phase corresponds to a certain range of raindrop sizes within which the raindrops have similar wind-driven rain behavior. For each rain phase, its continuity equation and corresponding momentum equation are as follows: where d is the diameter of the raindrop, α d is the volume fraction of the rain phase d, u d,j is the velocity component of the rain phase d in the j direction, u i is the wind speed component in the i direction, ρ w is the density of water, μ a is the dynamic viscosity of air, g is the acceleration due to gravity, C d is the drag coefficient of the raindrop, and Re R is the relative Reynolds value calculated based on the relative velocities of the wind and rain phases. A horizontal line on a variable indicates that Reynolds averaging is performed for this variable.

CHARACTERIZATION PARAMETERS OF WIND-DRIVEN RAIN
The intensity of wind-driven rain is usually expressed in terms of the specific catch ratio η d (d) for horizontal rainfall intensity and the catch ratio η. Specific catch ratio η d (d) is related to a specific rain phase d, whereas catch ratio η is related to the sum of the rain phases corresponding to all droplet sizes in the entire raindrop spectrum. The specific catch ratio η d (d) and the catch ratio η can be determined using the following equations: where R wdr is the wind-driven rain intensity; R h is the horizontal rainfall intensity (R h 10mm/h in the calculations of this study); f h (R h , d) is the size distribution of raindrops that pass through the horizontal plane, which is simulated in this study based on the raindrop distribution spectrum given by Best (1950) and |V n (d)| is the size of the rain phase velocity perpendicular to the building facade.

BOUNDARY CONDITIONS
The velocity profile of the wind speed at the entrance boundary is determined according to the logarithmic law using the following equation: where U(z) is the average wind speed at height z; u p ABL is the frictional velocity of the atmospheric boundary layer; κ is the von Karman constant, which usually takes the value of 0.42; and z 0 is the aerodynamic roughness length (z 0 0.49m in the calculations of this study). In this study, u p ABL 0.24U ref , where U ref is the reference wind speed at a height of z 0 . The profiles of turbulent kinetic energy k and turbulent dissipation rate ε at the entrance boundary are determined using the following equations: where I u is the turbulence intensity along the flow lines. As for the solid wall boundary, the standard wall function (Launder and Spalding, 1983) is used for the simulation, and the wall roughness is corrected by the gravel roughness (Cebeci and Bradshaw, 1977). The relationship between the gravel roughness height k s and the roughness constant C s is determined using the following equation (Blocken et al., 2007): where E is the empirical constant and takes the value of 9.793. The rain phase boundary conditions of the Eulerian model are simpler to set than those of the Lagrangian method, requiring only the volume fraction and initial velocity of the raindrops at the entrance boundary and top boundary. The horizontal velocity of the raindrops is assumed to be equal to the wind speed at the entrance boundary and top boundary; therefore, the relative horizontal velocity of the wind phase and the rain phase at the entrance boundary are zero. The vertical velocity of the rain phase is the vertical landing terminal velocity V t (d) of the raindrops, which can be determined using Stokes' law, and the phase volume fraction of the raindrops is determined using the following equation: The two sides of the computational domain are set as symmetric boundaries, and the exit boundary is set as the pressure exit boundary.

MATERIALS AND METHODS
A real-world case of wind-driven rain was simulated using a wind-driven rain solver based on the Eulerian multiphase model developed by OpenFOAM (Kubilay et al., 2015a;Kubilay et al., 2015b). The geometric shape of the calculations was as shown in Figure 1. There were two rectangular buildings of 6 m in length and 2 m in width within the computational domain. The height of the upstream building was 2.17 m, the height of the downstream building was 4.29 m, and the distance between the two buildings was 2 m. To obtain calculation results with adequate accuracy, a locally scaled grid was used. The grid was placed on both the top and sides of the windward facades of the building, with a minimum grid size of 0.004 m and a total grid size of about three million. The constant wind phase flow field around the building was the first to be simulated in the numerical simulation. A coupled pressure-velocity approach was used for the wind field simulation with the SIMPLE algorithm as the solution algorithm, and the second-order discrete format was used for both the convective and viscous terms in the wind field control equations. In this study, three different inflow velocity conditions were simulated, corresponding to wind speeds of U ref 1, 3, and 5 m/s, respectively. After obtaining the constant solution for the wind phase flow field, it was used as an input condition for solving the rain phase field. The criterion for assessing convergence when solving for the wind and rainphase fields was that the residual values were less than 10 −5 . In the rain phase field, a total of 17 rain phases with different diameters (diameters of 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 3.0, 4.0, 5.0, and 6.0 mm) were simulated. The specific catch ratio of each rain phase was determined individually using Eq. 3, and then the catch ratio was determined using Eq. 4.

RESULTS AND DISCUSSION
The normalized velocity distribution and flow line distribution of the wind field at different reference velocities were as shown in Figure 2. It can be seen that the normalized velocity at the various locations of the wind field was affected by the reference velocity. The greater the reference velocity, the greater the normalized velocity at the corresponding location of the wind field. The reason for this is that the parameter used here for the normalized velocity is the reference velocity U ref with a corresponding height of z 0 0.49m. Therefore, it is obvious that the velocities at most of the locations, especially those higher than z 0 , were greater than the reference velocity U ref . As a result, the normalized velocity increased with the reference velocity U ref . On the other hand, flow separation was observed at the front edges of the top of both buildings in all three working conditions, as seen in the velocity and flow line diagrams (Figure 2), and the separated flow lines reconnected in the middle region of the roof. The main structures of the flow fields were very similar: there was a vortex with a scale comparable to the width of the building before the upstream building; there was also a vortex with a scale comparable to the width of the building between the two buildings; and there was another vortex of a larger scale after the downstream building due to flow separation. It is evident from the results that different reference velocities had little effect on the vortex structure of the flow field regardless of the size or location of the vortex. Figure 3 shows the catch ratio distribution of the building facades at different reference velocities. It can be seen for all wind velocity cases that the catch ratio of the upstream and downstream building facades at the same height was smaller than that of the upstream building at the same location, and the reason for this can be identified by analyzing the flow characteristics of the wind field and rain phase field. It can be seen from the velocity cloud diagrams and flow line diagrams in Figure 2 that the vortex generated in the area between the two buildings extended mainly in the vertical direction, whereas the vortex before the upstream building extended mainly in the horizontal direction. For this reason, the catch ratio of the windward facades of the downstream building was smaller than that of the upstream building at the same height. Figure 4 shows the catch ratio distributions of the windward building facades at different reference velocities. It can be seen that all of the catch ratio distributions were high at the top and low at the bottom. This is mainly because the airflow was accelerated near the roof, where corner flow is formed, resulting in greater horizontal velocity and a greater catch ratio, whereas further down the building, the horizontal wind speed became smaller due to obstruction by the building, giving a smaller catch ratio.
The wind-driven rain model developed based on OpenFOAM provided numerical simulation results that model the case study well. Separation of flow was observed at the front edges of the top of the building, and the separated flow lines reconnected in the middle region of the roof. This type of corner flow with separation produces an accelerated horizontal airflow and increases the catch ratio; hence, there is often a larger catch ratio at the windward facades near the roof. Compared with previous work, our work simplified the rain-phase boundary conditions, made it easier to obtain the wind-driven rain parameters, and significantly reduced the computational time. In addition, our simulation results can be a guide for the design of buildings. There was also some limitation in our numerical models. As with most numerical models, our work involves simplifying boundary conditions and material properties during the analysis process, which ignores some factors that may affect the accuracy of the results.

CONCLUSION
In this study, a numerical simulation of a specific case was performed using a wind-driven rain model developed based on OpenFOAM. The main conclusions are as follows: 1) Separation of flow was observed at the front edges of the top of the two buildings, and the separated flow lines reconnected in the middle region of the roof. Different reference velocities had little effect on the vortex structure of the flow field regardless of the size or location of the vortex. 2) The vortex generated in the area between the two buildings extended mainly in the vertical direction, whereas the vortex before the upstream building extended mainly in the horizontal direction. As a result, the catch ratio of the windward facade of the downstream building is smaller than that of the upstream building at the same height.
3) The airflow was accelerated near the roof, where corner flow was formed, resulting in greater horizontal velocity and a greater catch ratio, whereas further down the building, the horizontal wind speed became smaller due to obstruction by the building, giving a smaller catch ratio. 4) The simulation results suggest that the Eulerian model utilized in this study is capable and accurate enough to simulate winddriven rain.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.