Numerical simulation and wind tunnel test for redistribution of snow on a flat roof

J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

Contents lists available at ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Numerical simulation and wind tunnel test for redistribution of snow on a flat roof Xuanyi Zhou, Luyang Kang, Ming Gu n, Liwei Qiu, Jinhai Hu State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China

art ic l e i nf o

a b s t r a c t

Article history: Received 16 July 2015 Received in revised form 16 March 2016 Accepted 17 March 2016

A real snowstorm often lasts for hours to several days. Simulating snow drifting on a roof through traditional transient method of computational fluid dynamics would be time consuming. Thus, a quasisteady method is proposed in this paper to simulate redistribution of snow on roofs. The process of snow drifting is divided into several phases according to the characteristics of meteorological data. The erosion/deposition of snow in each phase is regarded as a steady state. In the numerical simulation, the angle of repose of snow is taken into account. The boundaries of snowpack are updated based on the calculated results of change rate of snow depth in the previous phase. The proposed numerical method was applied to predict snow redistribution on a flat roof. A scaled test was conducted in a wind tunnel to validate the accuracy of the numerical simulation, in which a kind of high-density particle, silica sand, was used to model snow particles. The snow redistributions obtained from the numerical simulation agree well with those from the test results. As the wind tunnel test, the numerical simulation also predicted deposition on a small windward area of the flat roof. Then the transport rate of snow or silica sand on roof surface was investigated. The influence of the threshold velocity of snow and silica sand on prototype duration was also studied. Finally the snow redistribution on a flat roof during a real snowstorm was simulated in this study to show the practicability of the proposed method. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Redistribution of snow Numerical simulation Quasi-steady method Wind tunnel test Flat roof Meteorological data

1. Introduction Snow disaster challenges the environment and engineering for a long time. Snow drifting is associated with a significant number of snow disasters, which cause heavy losses to human beings. In the field of civil engineering, snow drifting may lead to an unbalanced snow distribution on roofs and ultimately collapse roof structures. Thus, accurate prediction of snow distribution on roof surface is vital to structural design. However, current load standards or codes are generally based on engineering judgement or limited field observations rather than theoretical foundation, which need the attention of researchers on snow engineering. With the development of computers, researchers began to employ computational fluid dynamics (CFD) to study the natural phenomenon of snow drifting. Most of these studies focused on the phenomenon around a model (Uematsu et al., 1991; Beyers et al., 2004; Tominaga et al., 2011a, 2011b). But only a small fraction of researches related with the snow drifting on roofs. Thiis and Ramberg (2008) employed a two-phase CFD simulation to simulate snow deposition on a curved roof. Zhou and Li (2010) simulated the redistributions of snow loads on the roof structures of n

Corresponding author.

http://dx.doi.org/10.1016/j.jweia.2016.03.008 0167-6105/& 2016 Elsevier Ltd. All rights reserved.

Beijing airport terminal using CFD technique. Although transient computation was performed in some of previous studies, it has to be pointed out that, the transient simulation on snow drifting is a very time-consuming work when the duration of snowstorm often lasts for many hours to several days. Another problem is that, the effects of snow surface changes caused by snow drifting were generally not considered in the previous studies. Wind tunnel test is also an important method to reproduce snow drifting phenomenon. Aside from the studies on the mechanism of snow drifting, researchers have also studied snow distribution in relation to engineering objectives and scaled models had to been employed in this kind of wind tunnel test. Isyumov and Mikitiuk (1990) employed bran as model particle to simulate the drift formation on the lower level of a two-level roof. Tsuchiya et al. (2002) studied the relationship between the snow drifting patterns on model roofs observed from field measurements and wind acceleration in the vicinity of roof surface. O'Rourke et al. (2004, 2005) analytically simulated snow drift loads on roofs and proposed a convenient method to calculate snow loads for roof structures. Flaga et al. (2009) predicted snow loads on two large-span roofs through wind tunnel tests, taking the theory of dispersion into account. Zhou et al. (2014) simulated the redistribution of snow load on a stepped flat roof in a wind tunnel using particles with different densities. However, no studies combined the method of

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

wind tunnel test with the CFD technique to simulate snow drifting on roofs. Regardless of their use of CFD technology and the wind tunnel test, most studies have been only concerned about the mechanism of snow drifting or snow drifting around simple models. Only a few of them were related with snow redistributions on roofs. CFD technique, through which a full scale model can be established, is expected to be an appropriate tool for simulating snow loads on roofs. Validation is always an essential part of CFD modeling of snow drifting. Field observation studies on snow drifting are very valuable. However, the results of full-scale measurements can only provide a qualitative validation to snowdrift pattern. In a real environment, wind velocities, wind directions or thermal conditions of roof cannot be controllable, which adds the uncertainties to the results of snow drifting on roof surface from field observations. Compared with field observation, the wind tunnel test can be performed under controllable conditions. Despite the use of a scaled model, the wind tunnel test is definitely a powerful method to reproduce snow drifting phenomenon with the development of similarity criteria theory. Thus, this study attempts to combine the CFD technique with the wind tunnel test to simulate snow drifting on roofs. Wind is sometimes accompanied by snowfall. Snow drifting usually occurs because of the wind, which usually changes the profile of snow on roofs. The current transient simulation on snow drifting using CFD would become a very time-consuming work since a snow drifting process often lasts for many hours to several days. To overcome these difficulties, a quasi-steady method was proposed in the CFD simulation, which considers the effects of temporal change on snow drifting on a flat roof. The angle of repose was taken into account in the CFD simulation to exactly replicate real snowpack on roofs. The wind tunnel test was performed to simulate the phenomenon of snow drifting on the flat roof, in which silica sand, a type of high-density particle, was employed as modeling particle herein. Then, a comparison was made between snow distributions on the roof obtained from CFD and wind tunnel test. The influence of threshold velocity of silica sand and snow on the duration of snow drifting was analyzed. Finally, the snow redistribution on a flat roof was simulated through the proposed numerical method in a real snowstorm.

2. Numerical model The turbulent air field is modeled by employing the Reynoldsaveraged Navier–Stokes (RANS) equations of motion, in which the Realizable k–ε model is adopted. These classical equations, which can be easily found in the other references, are not given in the study. The Eulerian description is also adopted for the simulation of snow drifting. The snow particles are treated as a continuous phase in the flow field. Characteristics of the snow phase are obtained by solving the convection–diffusion conservation equation of the snow phase. The commercial CFD software FLUENT is employed in this study with an additional user C code added to simulate the transportation of snow particles. The erosion and accumulation of snow particles on snow surface is treated as the boundary condition of the computational domain, in which some empirical formulas are employed. A snowstorm often lasts for many hours to several days. To improve the calculation efficiency, a quasi-steady method is proposed in this paper. The method attempts to compute the snow drifting on roofs during a real snowstorm. 2.1. Equations governing the snow phase The simulation of wind-driven snow is assumed to have a highly dilute air–snow system. Moreover, only one-way coupling

93

exists between the air and snow phase. The air phase is treated as unaffected by the airborne snow, and the snow transport is mainly controlled by the air flow and roof surface conditions. The equations that govern the snow phase and adopted in the paper were referenced from Uematsu et al. (1991) and Tominaga et al. (2011a)
  ∂ðρs f Þ ∂ðρs f uj Þ ∂ ∂ρ f ∂  þ ¼ νt s þ  w f ρs f ð1Þ ∂t ∂xj ∂xj ∂x3 ∂xj where uj is the wind velocity which is calculated through the equations governing the air phase, νt is the turbulent kinematic viscosity, ρs is the density of snow, f is the snow volume fraction, ρsf is the drifting snow concentration, and wf is the settling velocity of snow particles. Snow fraction f on the ground and inlet is set as zero. For the upper boundary and outlet of the flow domain, snow fraction f fits the relation ∂f/∂n ¼0. On the snow surface, the boundary conditions of snow fraction f are determined by Eq. (2). ρm is the mixture density, given by ρm ¼ f ρs þ ð1  f Þρ. The terms in the equation represent various transport processes of snow particles: the unsteady term and convection term on the left-hand side, and the diffusive term and additional convection term on the right-hand side. For the transportation of snow is supposed to be steady in each phase, the unsteady term is set as zero in the current study. An assumption is made that the turbulent diffusion flux of drift density is similar to the momentum diffusion of flow. Thus the turbulent diffusion coefficient is νt. This assumption is applied in most previous snowdrift studies (Tominaga et al., 2011a). The second term on the right-hand expresses the convection effect by gravitational sedimentation on snow particles in terms of snowfall settling velocity wf. 2.2. Erosion and accumulation on snow surface Whether erosion or accumulation occurs on the snow surface is mainly determined by friction velocity and snow surface property, which is represented as threshold friction velocity ut . Once the friction velocity u* exceeds the threshold friction velocity ut , snow drifting would occur. Snow particles can enter the computational domain, which results in erosion. If the friction velocity u* is lower than the threshold friction velocity ut , then the snow would leave the computational domain, depositing on the snow surface. The erosive or accumulative snow fluxes on the snow surface as the boundary conditions are calculated using the formula (Naaim et al., 1998; Beyers et al., 2004) as follows: qdep ¼ ϕwf qero ¼

u2t  u2 ; u2t

Aero ðu2  u2t Þ;

u o ut u 4 ut

ð2Þ

where Φ is the snow concentration, which can be obtained from the numerical results of snow field. Aero is a coefficient that represents the bonding strength of snowpack. A value of Aero ¼7  10  4 is employed for the present analysis (Naaim et al., 1998; Beyers et al., 2004). The total snow mass flux qtotal on the snow surface is expressed as qtotal ¼ qero þqdep

ð3Þ

2.3. Quasi-steady simulation method of snow drifting on roof surface Snow erosion/deposition caused by a long period of snowstorm could result in the change of snow cover boundaries on roof surface. This change affects the flow field around roofs. The algorithm described in Sections 2.1 and 2.2 presents only one-way coupling between air and snow phases. Another problem is that a snowstorm often lasts for many hours to several days. Thus, tremendous computing time would be consumed if traditional transient

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

simulation method is applied. The methods of altering boundaries of computational domain were adopted in previous research considering the effects of the change of snow depth during unsteady calculation (Sundsbø, 1998; Beyers et al., 2004). However, these methods are appropriate only for simple conditions such as snowdrift around a cube, especially for snowdrift with short duration. The methods are unsuitable for simulating snow drifting on building roofs in a snowstorm that could last several days. Meanwhile, when snow accumulates on the roof surface, snow particles reside at the edge of roof with a certain angle of repose, which can affect snow drifting on the roof surface. This feature distinguishes snow drifting on the roof surface from that on the ground. However, no study has considered this feature when simulating snow drifting on roof surface. A practical quasi-steady method is proposed in this paper to simulate snow drifting on roofs in a long period of snowstorm. This method attempts to balance computational efficiency and calculation accuracy. Based on the meteorological data provided by meteorological department, the process of snow movement on roof surface is divided into several phases. In each phase, the steady and one-way coupling method described in Sections 2.1 and 2.2 is used to calculate snow drifting on roof surface. And the angle of repose of snowpack is simulated in the initial profile of snow cover on roof surface while establishing computational domain. The boundaries of snow cover is updated according to the calculated results of change rate of snow depth in the previous phase, taking into account the effects of change of snow depth during a snowstorm. Specific procedures are as follows: (1) Based on the meteorological data provided by a meteorological department, the duration T of certain snow drifting is divided into n phases. The duration of phase i is represented by Δti. The snow particles are supposed to be paved uniformly on the roof surface during the precipitation. The erosion/ deposition of snow in each phase caused by snow drifting is regarded as a steady state. An example using the real meteorological data is given in Section 4.5. (2) The steady velocity field of phase i is obtained by solving the governing equations of air. During the simulating process for velocity field, meteorological data are used. According to the velocity field in phase i, the friction velocity u* on roof surface can be calculated. (3) Based on the calculation results of Step 2, steady and one-way coupling method (Sections 2.1 and 2.2) is employed to solve the governing equation of snow phase. The change rate of snow depth per unit time during phase i can be calculated by the equation Δhi ¼ qi; total =ρs . (4) The initial snow depths of phase iþ1 are calculated according to the following equation: hi þ 1 ¼ hi þ Δhi U Δt i þ di

ð4Þ

where hi and hi þ 1 represent initial snow depth of phase i and iþ1, respectively, and di is the increased snow depth as a result of snowfall during phase i. The angle of repose of snowpack is simulated in the initial profile of snow cover on roof surface while establishing computational domain. Snow boundaries are updated and the domain is regrided for the computation of the next phase after the change rate of snow cover profile is obtained (Wang et al., 2013). To be specific, in each phase the snow depths are calculated by Eq. (4), where change rate of snow depth per unit time obtained from the steady computation is used. The boundary of the roof with snowpack is the sum of roof height and snow depth, which determines the coordinates of the nodes of the boundary of roof in the next phase. The coordinates of the nodes are then

used to generate the curve of the snow surface. In the software for grid generation, the approximate smooth curves can be generated through a series of separated nodes automatically. Finally, the flow field is regrided according to the new boundaries using quadrilateral cells. If the number of nodes is large, the boundary of snow surface could be thought to be smooth sufficiently. (5) Steps 2–4 are repeated until the final phase n, when the final distribution of snow on roof surface is obtained. The purpose of this quasi-steady method is to simulate a specific process of snow drifting on roof surface, through which the snow distribution on roofs after a snowstorm can be obtained based on meteorological data. Although one-way coupling method is employed in each phase, the influence of change of snow depth on flow field is considered by altering the profile of snow cover in the next phase. That takes into account the effects of motion of snow particles on the flow field to some extent. 2.4. Simulation parameters Snow loads on a flat roof are important for most current load codes and standards (ISO4355, 2013; EN19 91-1-3, 2003; NBC2005, 2005; ASCE7-10, 2010; GB50009-2012, 2012). Thus, this paper chooses a flat roof as its research object. The span of the studied flat roof is 12.0 m, and the height is 3.0 m. Fig. 1 shows the dimensions of the flat roof. In commercial software FLUENT, the prototype flat roof is modeled in two dimensions using the numerical method proposed in the study. Thus, subscript p is used to denote physical quantities related to the prototype roof. Moreover, the scaled model of the roof is employed in the wind tunnel test. Thus, subscript m is employed to denote physical quantities related to the scaled model below. Physical properties of snow particles adopted in the numerical simulation are shown in Table 1. This study does not measure the physical properties of natural snow particles. Their values are determined by some related literature (Kuroiwa et al., 1967; Tominaga et al., 2011a). The angle of repose of the pulverized snow varies from 45° to 50° in the range of  35 °C to  3.5 °C. It approaches 90° near the melting point (Kuroiwa et al., 1967). In this study it is supposed that the temperature is below 0 °C during snowfall and that the angle of repose of snow is 50°. The profile of neutral equilibrium atmospheric boundary layer proposed by Yang et al. (2009) is employed at the inlet for approaching wind profile
 u z ln þ1 ð5Þ U ðz Þ ¼ z0 κ k¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C 1 lnðz þ z0 Þ þ C 2

ð6Þ

repose angle of snow θ=50°

H=3.0m

94

L=12.0m Fig. 1. Dimensions of flat roof.

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

∂u ∂z

ε ¼ C 1=2 μ k ðz Þ

ð7Þ

where the von Karman constant is 0.42. Constants in Eq. (6) are fitted by the profile of turbulent kinetic energy obtained from the wind tunnel test below. k(z)¼ 1.2[U(z)I(z)]2, where U(z) is the wind velocity at the height of z, I(z) is average turbulent intensity at the heights of z, C1 ¼  0.064, C2 ¼0.588 and Cμ ¼ 0.09. The descriptions of approaching wind profile in the wind tunnel test can be found in detail below. The characteristics of wind velocity during snowfall period in Shenyang and Nanjing, which are regarded as the representative cities of northern and southern China respectively, are analyzed in this paper. Based on the meteorological data provided by the China Meteorological Data Sharing Service System, Fig. 2 shows the probability distribution of wind velocity at 10 m height during snowfall period in Shenyang and Nanjing from 1980 to 2014. The average wind velocities during snowfall at 10 m height of Shenyang and Nanjing are 3.7 and 3.4 m/s, respectively. The probability of wind velocity being lower than 6.4 m/s is 86% and 95%. As an example of engineering application, the initial snow depth on prototype roof in the numerical simulation was set as 50 cm and the wind velocity at a height of 10 m was set as 6.4 m/s. Friction velocity on the ground snow surface u* was 0.35 m/s, and roughness height z0 was 0.00426 m. Friction velocity and roughness height were obtained by fitting the measured data obtained from the wind tunnel test according to Eq. (5). Thus, the wind velocity at prototype roof height is U(H)¼5.4 m/s. Given that a snowstorm could last from many hours to several days, the duration of wind Tp was set as 10.3 days. Hence, the process of snow drifting was divided into three phases, and the time duration of each phase was 2.4, 5.4, and 2.5 days. The wind duration of numerical simulation for each phase was calculated by the results from wind tunnel test. For the wind duration of the model in the wind tunnel test is fixed, after the numerical simulated mass transport rate Qp and experimental results Qm was obtained, the wind duration of prototype could be determined Table 1 Physical properties of natural snow particle and silica sand.

according to Eq. (11). From Table 6, we could see the proportion of Qp from the numerical simulation is different from those of wind tunnel test. Hence, though in wind tunnel test, the wind duration was the same for each phase, the values of wind duration of prototype for each phase were different according to Eq. (11). In this paper, the transport rates from numerical simulation are found to be close to those from empirical formula used by O'Rourke et al. (2005), which can be found in Table 6. Thus the prototype time calculated from this similarity relation could be reasonable. For the precipitation could occur several times during these three phases, the snow was supposed to be driftable for the entire simulated duration. Three phases are represented by I, II, and III. Subscript S (Start) represents the initial time of a certain phase. Subscript E (End) represents the end of a certain phase. Table 2 presents the major parameters in the numerical simulation. Table 3 summarizes the computational conditions in the CFD simulation. Structured progressive mesh was employed as grid scheme. The growth factor of grid was 1.1, and the total number of quadrilateral grid cell was 18.4 thousand. Fig. 3 shows the influence of grid size on the calculation results of initial friction velocity. The minimum size of the grid near the snow surface is H/15, H/20 and H/30, where H is the roof height. It can be seen from the figure that friction velocity is almost the same for different minimum grid size near the whole snow surface, which confirmed that the prediction results did not change significantly with finer grids. In this paper, the minimum grid size is set as H/20.

3. Description of wind tunnel test The wind tunnel test is an effective way to study snow redistribution on roof surface. Given that wind tunnel provides an artificially controlled platform to conduct snow drifting test, this offers a big advantage compared with field measurement, which is Table 2 Parameters of numerical simulation and wind tunnel test. Description

Name of particle

Snow particle

Silica sand

Diameter, dp (mm) Particle density, ρp (kg m  3) Bulk density, ρb (kg m  3) Threshold friction velocity, ut (m s  1) Angle of repose, θ (deg) Settling velocity, wf (m s  1)

0.15 250 150 0.20 50 0.2

0.20 2784 1436 0.26 34 0.6

95

Numerical simulation Wind tunnel test

Geometric scale ratio 1:1 Initial depth of snow, S0 50 cm Wind velocity at roof height, U H (m/s) 5.4 0.35 Friction velocity, u (m/s) 0.00426 Roughness length, z0 (m) 13 Turbulent intensity at roof height, I H (%) 2.4dþ 5.4d þ2.5d Wind duration

Fig. 2. Probability distribution of wind velocity during snowfall period.

1:25 20 mm 7.0 0.46 0.00017 13 1 min þ1 minþ1 min

96

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

Table 3 Computational conditions for flow field. 16L(x)  30H(z), L: building length; H: building height (here H ¼3.0 m, L ¼12.0 m) The minimum grid width is 0.05H Profile of velocity and turbulence intensity is fitted by wind tunnel test data according to the form proposed by Yang Yi et al. (2009) Downstream boundary Zero gradient condition (outflow) Upper surface of computational domain Zero gradient condition is used for all variables, except that the normal components of velocity with respect to the boundaries are set to zero (symmetry) Building and ground surface boundary Standard wall functions (wall) Scheme for convection terms Second order Computational domain Grid discretization Inflow boundary

Q

0.6 Wind

H/15 H/20 H/30

Friction velocity (m/s)

0.5

0.4

0.3

0.2

Fig. 4. Transport rate of snow particle on flat roof.

Front

Middle

Rear

0.1

Relative position X/H Fig. 3. Comparison of initial friction velocity for different grid size.

influenced by many complicated factors such as meteorological conditions. A scaled test was conducted in the wind tunnel for this flat roof to validate the accuracy of numerical simulation. The test is to simulate the redistribution of snowpack on the flat roof when it is subjected to wind with a constant velocity. An application example of a real snowstorm, according to meteorological record, would be given later. In this section, similarity parameters are first introduced. Test parameters and the process of wind tunnel test are then described in detail. 3.1. Similarity parameters Silica sand, which is a high-density particle, was used to model snow particles on the surface of roof in the wind tunnel test. Some researchers considered that high-density particles can reproduce the phenomenon of wind-induced redistribution of snow better than low-density particles (Kind, 1986; Zhou et al., 2014). A number of similarity requirements should be satisfied in the wind tunnel test on snow drifting on roof surface, such as the similarity of terrain, flow field, ejection process, particle trajectory, and deposition pattern (Zhou et al., 2014). In the wind tunnel test, most similarity parameters can be determined if the type of modeled particle and geometric scale ratio are fixed. The uncertain parameters are test velocity and wind duration. Thus, the following content focuses only on the similarity parameters of wind velocity and wind duration. In the snowdrift test, the similarity parameter of wind velocity proposed by Anno (1984) was employed     U ðH Þ U ðH Þ ¼ ð8Þ ut m ut p where U(H) is the wind velocity at roof height, and u*t is the threshold velocity of particles. Eq. (8) guarantees that the shear stress on snow surface is similar with respect to the threshold

velocity of particles. Thus, similar erosion/deposition characteristics could be obtained by the wind tunnel test and numerical simulation. Based on the similarity of dimensionless volume of a leeward snowdrift between a model and its prototype, Anno (1984) proposed a similar criterion related to the measured results. Some other researchers also used Eq. (9) to determine the wind duration of prototype, such as Smedley et al. (1993) and Beyers and Harms (2003) for outdoors modeling of snowdrift. This relationship was validated by Anno (1984) through using two different materials for the same experiment. Smedley et al. (1993) found that this similarity parameter were close between the prototype and scaled modeling. Hence we adopted this relationship in the study to calculate the wind duration of prototype " # " # tQ η tQ η ¼ ð9Þ ρb L 2 m ρb L2 p where Q is mass transport rate, η is the obstruction’s collection coefficient of snow particles, ρb is the bulk density of particles, and L is characteristic length. t on the right-hand side is the duration of prototype storm, whereas the t on the left-hand side is the duration of model test. The snow transport rate on roof surface for the flat roof is defined as the product of bulk particle density and the integration of particle depth change per unit width per unit time along the entire roof span. The calculation formula is as follows: Z ð10Þ Q ¼ ρb ΔhðxÞdx L

According to Eq. (10), the transport rate is actually the total particle mass blown off the roof per unit time and unit width (Fig. 4). Based on the transport rate Q, a dimensionless parameter for duration ratio, which resembles that which was proposed by Anno (1984), is defined in this paper. " # " # tQ tQ ¼ ð11Þ ρb L 2 m ρb L2 p where the characteristic length L herein is the span of flat roof (Fig. 1). This relation was based on the similarity of dimensionless

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

volume of a snowdrift. In this paper, the dimensionless volume is the amount of snow blown off the roof normalized by standard volume ρbL2. When the value of left side for model is fixed, to obtain the same pattern of snow redistribution, the wind duration of prototype would change for different values of the simulated transport rate of prototype, Qp. 3.2. Description of wind tunnel test The scaled model of snow redistribution on the flat roof was tested in the TJ-1 Boundary Layer Wind Tunnel at Tongji University. The TJ-1 wind tunnel is an open straight-circuit type, and its working section is 1.8 m wide and 1.8 m high. During the test for threshold velocity, the baseboard of the working section was paved with a thin layer of particles. The wind velocity was increased until a large number of particles began to move. The wind velocity then gradually decreased. The wind velocity at 1.0 m from the baseboard was recorded when the particles ceased to move. This velocity is considered the threshold velocity. Thus, the threshold friction velocity can be calculated using Eq. (5). This threshold velocity for silica sand is 0.26 m/s. Physical properties of silica sand are given in Table 1. The test region layout of the wind tunnel is presented in Fig. 5. The testing baseboard was paved with a piece of 16 mm thick wooden board to fix the model in place. The front edge of the board had a sharp transition, which prevented the flow separation at the front edge. The windshield was installed at a distance of 20 cm ahead of the model. At the beginning of the test, the

97

windshield was lifted, and the wind velocity was gradually increased. The windshield was maintained at the erect position until the wind velocity at roof height reached test velocity. The string, which was tied to the windshield, was then pulled to make the windshield fall into the board groove. The depth of the board groove is 16 mm, which is equal to the thickness of the windshield. That makes the test baseboard smooth when the wind shield was put down. In this way, the flow field would not be affected by the windshield. The purpose of wind shield is to keep the particles on the roof surface at rest during the increase of wind velocity and to ensure the stability of the wind field throughout the following test. As this test aims to simulate snow drifting on 2D roofs, two side plates with a sharp front edge were installed on both sides of the model to remove the 3D effect of the model (Figs. 5 and 6). A sand fence made of airing gauze mesh was installed in the downwind area located at a distance of 2.0 m from the model to collect the testing particles blown away from the model surface. A tailormade snow scale was applied to record the depth of the particles on the roof surface. This scale was employed to measure the snow depth at the central axis of the roof. The snow scale had both horizontal and vertical scales at a precision of 0.5 mm. The bottom of the snow scale was sharpened to reduce the disturbance of the particles during insertion into the particle layer for measurement. The snow scale is shown in Fig. 7. During the measurement, the scale was slowly inserted into the particle layer. Photographs were taken for subsequent data reading and processing. The average wind velocity profile, turbulent intensity, and turbulent kinetic energy are shown in Fig. 8(a)–(c), respectively. In

plan view

side wall of wind tunnel baseboard side plate sand fence

roof model

wind shield groove

wind side plate

side wall of wind tunnel camera

vertical view

side plate

wind shield wind

baseboard

roof model

groove

During the process of increase of wind velocity side plate

baseboard

roof model

wind groove

wind shield After the wind velocity is stable Fig. 5. Layout of test region: (a) plan view and (b) vertical view.

98

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

these figures, the scattered points are the measured values, and the curves were fitted which are also used for the numerical simulation. The turbulent kinetic energy is calculated from the test wind velocity and test turbulent intensity, k(z) ¼1.2[U(z)I(z)]2. The influence of sand fence on flow field is considered in Fig. 8. According to the similarity parameter of wind velocity in Eq. (8), the wind velocity of test at roof height should be 7.0 m/s. The turbulent intensity at roof height is about 13%. The wind tunnel test is divided into three phases with symbols I, II, and III to denote each phase. Similar to the numerical simulation, subscript S (Start) is employed to denote the initial time of one phase, and subscript E (End) is used for the end. The duration of wind in each phase tm is one minute.

Geometric scale ratio is 1:25 in the wind tunnel test, and the initial snow depth on roof surface of the test is 20 mm. The major parameters in the wind tunnel test are also given in Table 2. According to the physical properties of modeled and prototype particles shown in Table 1 and the information given in Table 2, the similarity parameters for model and prototype could be obtained (Table 4). The Reynolds number of aerodynamic roughness height of test u3t =2g ν maintains the fully rough flow of the snow particles. With respect to dynamic similarity, it satisfied the requirement of more than 30 (Kind, 1986). The similarity paraρu2 was in the meter of test of aerodynamic roughness height ρ Hg p

same order of prototype. However, it was larger than the prototype. Different from snowdrift in open terrain, strong separation of flow appears at the windward edge of roof. Thus, the aerodynamic roughness height caused by particle saltation is not the major influence factor of flow field. Therefore, the requirement of simiρu2 could be relaxed. With respect to the larity parameter of ρ Hg p

similarity in particle ejection process, Kind and Murray (1982) and Kind (1986) indicated that density ratio ρp/ρ should be larger than 600. The density of silica sand employed by the test was large enough to meet this requirement. The densimetric Froude number of the test was different from the prototype. However, Anno (1984) Table 4 Similarity parameters of CFD prototype and wind tunnel test scale model.

Fig. 6. Side plates on both sides of model.

Description

Similarity parameter

Snow (CFD)

Silica sand (Wind tunnel test)

Reynolds number of aerodynamic roughness Aerodynamic roughness height

u3t =2gv 4 30

—

61.8

ρu2 ρp Hg

2.04E  05

7.92E  05

0.134

0.015

ρ

Ejection condition

u2t

ðρp  ρÞ gdp U ðH Þ ut

Ejection condition

ρp =ρ Z 600

Ejection condition

Particle trajectory

ρP U 2 ðH Þ ρp  ρ gH wf U ðH Þ tQ ρb L 2

Particle trajectory Time

27

27

— 1.0

2272.7 41.7

0.037

0.086

Phase I: 0.050 Phase II: 0.057 Phase III: 0.012

0.050 0.057 0.012

Fig. 7. Snow scale.

10

8

6

4

2

0 0.6

10 Exp Fitted

8 Height (h/H)

Exp Fitted Height (h/H)

Height (h/H)

8

10

6 4

1.6

0 4%

6 4 2

2

0.8 1.0 1.2 1.4 Wind velocity (U/UH)

Exp Fitted

8% 12% 16% Turbulent intensity

20%

0 0.3 0.6 0.9 1.2 1.5 Turbulent kinetic energy (k/kH)

Fig. 8. (a) Experimental and fitted mean wind velocity profiles. (b) Experimental and fitted turbulent intensity profiles. (c) Experimental and fitted turbulent kinetic energy profiles.

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

conducted a contrast test of different particles and pointed out that Froude number had a minimal influence on the final pattern of redistribution of snow. Thus, the similarity parameter of ρ ρ ρ ð p Þ u2t

ρ

2

ðH Þ and ρ p ρ UgH could be relaxed in practical application. With p reference to the similarity of particle trajectory, the ratio of settling velocity of silica sand and reference velocity was larger than that of prototype snow particles. However, the settling velocity of snow was between 0.2 m/s and 0.5 m/s considering the inhomogeneity w of snow particles. Thus, U ðHf Þ of silica sand was within the range of gdp

w

the prototype. The value of U ðHf Þ ranged from 0.037 to 0.093. The wind velocity and duration of the wind tunnel test were determined according to Eqs. (8) and (11). Thus, these two similarity parameters of wind velocity and time were the same between the model and the prototype.

4. Results and discussion

99

4.1. Characteristics of flow field and friction velocity of snow surface Figs. 9 and 10 show the vectors of wind velocity around the flat roof at the start time of Phase I. The angle of repose of snow in Fig. 9 was set as 90°, and the flow separated strongly at point A at the front edge of the windward side, forming a large-scale vortex attached to the roof surface. At this time, the direction of wind velocity close to the snow surface was opposite to the approaching flow, except the leeward region of roof. The results with repose angle of 50° are shown in Fig. 10. The separation points at the front edge of snow surface under approaching flow could be points B and C. However, due to a smaller repose angle, the strength of flow separation significantly weakened compared with the condition of the right angle seen in Fig. 10b. The direction of wind velocity near snow surface is approximately parallel to the snow surface, which is different from that shown in Fig. 9. Meanwhile, the vortex cannot be observed evidently from the velocity vector diagram (Fig. 10a) at the front edge of the windward side. Thus, the angle of repose of snow must be taken into account appropriately in the

CFD method is employed to analyze the characteristics of flow field around the flat roof and the distribution of friction velocity on roof surface in this section. Then the results from the numerical simulation and the wind tunnel test are compared. Given that the threshold friction velocity of silica sand and snow directly affects the amount of particles blown off, the influence of threshold friction velocity on duration of prototype is analyzed using the proposed numerical method. Finally the proposed numerical method is applied at the end of this section to predict the snow redistribution on a flat roof during a real snowstorm, which shows the practicability of this method.

Fig. 10. Velocity vectors around snow surface (repose angle as 50°, phase Is, prototype): (a) whole roof and (b) front edge.

Fig. 9. Velocity vectors around snow surface (repose angle as 90°, phase Is, prototype): (a) whole roof and (b) front edge.

Fig. 11. Velocity vectors around snow surface (phase IIs, prototype).

100

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

phase Ιs

0.5

Fig. 12. Velocity vectors around snow surface (phase IIIs, prototype).

Friction velocity (m/s)

phase ΙΙs phase ΙΙΙs

0.4

0.3

0.2 Front

Middle

Rear

0.1 0

1

2 Relative position X/H

3

4

Fig. 14. Comparison of friction velocity among front, middle and rear regions (prototype).

Table 5 Average friction velocity of snow surface.

Fig. 13. Snow concentration around flat roof (unit: g/m3, prototype): (a) phase Is, (b) phase IIs and (c) phase IIIs.

numerical simulation of snow redistribution on the roof surface. Otherwise, a great difference of the flow field near the roof surface would be generated. This great difference could lead to a serious deviation between the simulated results and actual situation. Therefore only results with repose angle of 50° simulated are analyzed below. Velocity vector diagrams of Phase II and Phase III at start time are presented in Figs. 11 and 12, respectively. The characteristics are similar to those of Phase I. Thus, an analysis of these characteristics is not given. Fig. 13 presents the snow concentration around the roof at the start time of every phase. The snow concentration near snow surface in the initial profile of snow cover was between 20 g/m3 and 60 g/m3. Snow concentration was slightly larger than 60 g/m3

Region

Phase Is (m/s)

Phase IIs (m/s)

Phase IIIs (m/s)

Front Middle Rear

0.345 0.253 0.243

0.241 0.237 0.245

0.177 0.218 0.238

only in a small area at the end edge of the snow surface. Fig. 13 (a) shows that snow concentration near the snow surface increases gradually from the windward edge to the leeward edge. The wind velocity close to the snow surface shown in Fig. 10 indicates that the snow particles moved downwind to the end edge of the roof. The entrainment process by wind resulted in the largest concentration at the end edge of roof. The snow concentration at the end edge of roof reached its maximum. The numerical simulation reflects the snow drifting on the roof under approaching flow. A comparison of friction velocity on the snow surface of the three phases that was simulated numerically for the prototype is presented in Fig. 14. The friction velocity of each phase was obtained from the snow cover pattern at the start time of each phase. At Phase Is, the friction velocity at the front edge of roof was large, which led to the significant erosion of snow in that region. The change in the profile of snow cover caused by erosion also altered the flow field around front edge of roof. This condition explains why the friction velocity at the windward edge changed drastically. The flow field at the middle and rear roof was relatively stable, and friction velocity changed smoothly. To further analyze the characteristics of friction velocity in different regions, the snow surface was divided into three regions: front, middle, and rear. The division was according to the change features of friction velocity. The proportions of three regions to the roof span are 2/8, 5/8, and 1/8 (Fig. 14). Table 5 shows the average friction velocity of each region. The average friction velocity at the front region dropped rapidly, and at Phase IIs and Phase IIIs, it had decreased by 30% and 27%, respectively, compared with the previous phase. At Phase IIIs, the average friction velocity of front region was 0.177 m/s, which was already smaller than the threshold friction velocity 0.20 m/s. The average friction velocity of the middle region decreased slightly. In each phase, the decrease amplitude was no more than 0.02 m/s, which was 6% and 8% of the average friction velocity of the previous phase. At Phase IIIs, the average friction velocity was 0.218 m/s, which is close to the

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

101

Table 6 Transport rate of snow or silica sand on flat roof. Phase I

Phase II

Phase III

3.31 (100%)

1.64 (50%)

0.76 (23%)

Wind tunnel test 10  2 kg m  1 s  1

1.74 (100%)

1.98 (113%) 0.42 (24%)

Eq. (14) 10  4 kg m  1 s  1

4.94 (100%) 2.92 (59%)

1.78 (36%)

Eq. (15) 10  4 kg m  1 s  1

2.20 (100%)

0.79 (36%)

Numerical simulation 10  4 kg m  1 s  1

1.30 (59%)

threshold friction velocity of 0.20 m/s. The average friction velocity at the rear region was almost unchanged throughout the entire wind duration. The velocity remained about 0.24 m/s. Due to the formation of a mound in the front region of snow surface (Fig. 15), the friction velocity decreased gradually as snow drifting continued.

The results of the wind tunnel test are first analyzed. Fig. 15 (a) reflects the erosion/deposition status on the entire flat roof surface. The particles at the front region of roof were severely eroded during Phases I and II. The degree of erosion of modeled particles in Phase I was relatively uniform in the middle and rear region, whereas the erosion in the middle roof in Phase II was no longer uniform. The particles in the middle and rear region of roof in Phase II were severely eroded. The total change rate of particle distribution pattern in Phase III decreased significantly compared with that in Phases I and II. A few particles accumulated in the middle region in Phase III compared with that in Phase II, except most modeled particles were eroded on the flat roof. Fig. 15 (b) shows the enlarged view of the snow distribution. The numerical simulation results are then analyzed. The snow distribution on the roof surface obtained from the numerical simulation was similar to that of the wind tunnel test in Phase I. Snow on the entire roof surface was eroded, and the erosion in the middle region was relatively uniform. However, in Phase I, the amount of eroded particles obtained from the numerical simulation was slightly larger than that from the wind tunnel test in the middle roof. Moreover, the amount of eroded particles obtained from the numerical simulation was slightly smaller than that from the wind tunnel test in the front and rear regions. In Phase II, the numerical simulation also predicted the non-uniform erosion in the middle region. Moreover, the amount of eroded snow obtained from the numerical simulation was larger than that from the wind tunnel test, whereas the opposite effect was found at the front region. In Phase III, snow was only eroded slightly in the middlerear region. The total change of snow distribution pattern on roof surface was already very small. Consistent with the results from the wind tunnel test, the numerical simulation also predicted the deposition of particles on the roof surface. However, deposition occurred in the front region during the numerical simulation, and deposition occurred in the middle region close to the front region during the wind tunnel test, as shown in Fig. 15(b). The deposition of the numerical simulation occurred earlier than in the wind tunnel test. Deposition occurred in Phase III of the wind tunnel test, whereas deposition occurred in Phase II in the numerical simulation.

4.2. Redistribution of snow on roof surface

4.3. Transport rate of snow on roof surface

Snow redistribution on the roof surface at the end of each phase from both the numerical simulation and the wind tunnel test is shown in Fig. 15. Dimensionless snow depth is the ratio of snow depth at the end of each phase on the roof surface to the initial snow depth S0, which is 50 cm for the prototype and 20 mm for the model roof. The amount of erosion predicted by the numerical simulation for snow redistribution at the front region of roof surface was smaller than that of the wind tunnel test. The results from the numerical simulation and wind tunnel test were relatively consistent in the middle and rear regions of the roof. For the middle region, a longer distance from the front region corresponds to increased erosion.

Snow transport rate on roof surface is the quantitative expression of snow drifting. It is an important indicator to determine snow loads on roofs. In this section, the snow transport rate, which is obtained by the numerical simulation and wind tunnel test, is analyzed in detail. Table 6 shows the calculated results of snow transport rate obtained from both the numerical simulation and wind tunnel tests using Eq. (10). The values of snow transport rate obtained by two ways were normalized using the transport rate in Phase I as reference. The snow transport rate from the numerical simulation was the result at the start time in each phase, whereas it was the average value in each phase for the wind tunnel test.

phase IIE numerical results

Normalized snow depth S/S0

1.0

0.9 phase IIIE wind tunnel results Phase IIE Wind tunnel results Phase IIIE Wind tunnel results Phase IE Numerical results Phase IIE Numerical results

0.8

0.7

Wind Front

Middle

0.6 0.5

1.0

1.5 Relative position X/H

2.0

2.5

Detail with enlarged scale Fig. 15. Snow distribution on flat roof: (a) whole roof and (b) detail with enlarged scale.

102

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

As seen from Table 6, as snow drifting proceeded, snow transport rate obtained from the numerical simulation decreased gradually, and the rate became nearly half compared with that of the previous phase. Unlike in the numerical simulation, the transport rate of roof obtained by the wind tunnel test did not decrease in Phase II. Particle transport rate in Phase II was close to that in Phase I and rose slightly. The transport rate on the roof from the wind tunnel test after entering Phase III began to drop dramatically. Interestingly, the sharp decrease of transport rate coincided with the particle deposition in both the wind tunnel test and the numerical simulation. According to Section 4.2, snow deposition occurred in Phase II during the numerical simulation, and particle deposition occurred in Phase III during the wind tunnel test. Transport rate declined after deposition. However, future studies should be performed to study their correlation. Meantime, almost the same decrease in the amplitude of transport rate was eventually achieved in the wind tunnel test and numerical simulation. The transport rate in Phase III was about 1/4 of Phase I both in the wind tunnel test and numerical simulation. Snow transport rate obtained by the numerical simulation were further compared with the calculated results obtained through an empirical formula to verify the proposed quasi-steady numerical method to a certain extent. O'Rourke et al. (2005) proposed an analytical method for calculating drift snow loads on a roof; this method is applied by the current American Society of Civil Engineers ASCE standard (ASCE7-10, 2010). The method employed by O'Rourke et al. (2005) to calculate snow transport rate on a roof was adopted to compare numerical simulation results. However, the method of O’Rourke was modified slightly for the research object is different. In the design of snow fences for highways, Tabler (1994) presented the following relation for the transport rate for snow along an open field: Q ðU 10 Þ ¼ 4:27  10  6 ðU 10 Þ3:8

ð12Þ

where the snow transport rate Q is fully developed for snow drifting, and U10 is the wind velocity at 10 m for approaching flow. Full-scale measurements obtained by Takeuchi (1980) suggest that the transport rate becomes fully developed for an upwind fetch of about 210 m. At smaller fetches, the reduction factor for snow transport rate on the ground appears toffi be proportional to the pffiffiffiffiffiffiffiffiffiffiffiffiffi square root of the fetch distance F=210, where F is the upwind fetch length. Thus, the formula used by O'Rourke et al. (2005) to calculate snow transport rate on a roof surface is as follows:
 L 0:5 Q ðU 10 Þ ¼ 4:27  10  6 ðU 10 Þ3:8 ð13Þ 210 where L is the span of upwind roof. O'Rourke et al. (2005) studied stepped and gable roofs. L is the span of an upper-level roof for stepped roofs and upwind roof for gable roofs. The upwind roof is the physical part of the roof where erosion occurs. It acts as a snow source area for leeward roof. Almost the entire flat roof was eroded in this study. Moreover, the entire roof surface can be regarded as a snow source area. Hence, L could be the entire span of the flat roof when applying Eq. (13) to calculate snow transport rate. Due to snow conditions and location, the minimum fetch distance for fully developed snow transport could vary over a wide range. Tabler (1994, 2003) gave this distance as 3 km in Wyoming, which is considerably different from the 210 m measured by Takeuchi (1980). In this paper, the snow is supposed to be fresh, which could have a lower minimum fetch before the snow profile is fully developed. Hence minimum fetch length of 210 m is used herein for the calculation of snow transport rates on roofs (O'Rourke et al., 2005). The influence of roof profile on transport rate is not considered when Eq. (13) is employed for snow drifting on a roof surface.

Unlike the situation of an open field, the wind velocities of approaching flow are different at different heights even if no interference is seen from other buildings. Wind velocity has to be modified to apply Eq. (13) to the study. The modification of wind velocity uses the friction velocity on the snow surface of a roof. The modified formula is as follows: 

u;r 3:8 L 0:5 Q ðU 10 Þ ¼ 4:27  10  6 U 10 ð14Þ 210 u;g where the wind velocity of approaching flow at 10 m height above the ground surface is U10 ¼6.4 m/s, u;r is the average friction velocity on the entire snow surface of roof, and u;g ¼ 0:35 m=s is the friction velocity on the ground (namely u in Table 2). For the three phases of numerical simulation, the revised wind velocity was 5.1, 4.4, and 3.9 m/s, respectively. U 10 uu;r ;g The snow transport rate of the prototype calculated by the modified Eq. (14) is shown in the third line of Table 6. The results of snow transport rate were also normalized by using transport rate at the start time of Phase I. The numerical simulation results were in the same order as the results of the empirical formula. The numerical simulation results were relatively smaller. The snow transport rate on the roof surface calculated by the empirical formula also decreased gradually with time, and the decline amplitude was about 40%. The decline of empirical calculation was slightly smaller compared with 50% decline from the numerical simulation. Combined with the results in Section 4.2, a mound of snow formed near the front edge of the roof, which blocked flow downwind (Fig. 15). Hence, the friction velocity on snow surface decreased and the transport rate also declined. Snow redistribution was more significant in the early stage of snow drifting. According to the square-root relationship used by O'Rourke et al. (2005), When roof span is larger than 210 m, the reduction factor would be larger than 1.0. Whereas, Tabler (1994, 2003) used a better model to reflect the increase of transport rate with fetch distance, where the upper limit for the reduction factor is 1.0. Thus, the following equation is also used to compare the simulated results. The snow transport rate calculated by the modified Eq. (15) is shown in the last line of Table 6.
 u;r 3:8 Q ðU 10 Þ ¼ 4:27  10  6 U 10 ð1  0:14L=210 Þ ð15Þ u;g For the fetch distance of flat roof is fixed, the reduction factor is the same in each phase. Hence, the decline amplitude is identical to those from Eq. (14). Transport rates from Eq. (15) are less than half as those from Eq. (14), and they are closer to the numerical simulation compared with those from Eq. (14). However, it should be pointed out that the minimum distance for fully developed snow transport could influence the results significantly. In this paper for fresh snow, the minimum distance is set as 210 m, which was also used by O'Rourke et al. (2005). Just as described previously, this distance could vary over a wide range. 4.4. Parameter analysis The threshold friction velocity of modeled particles is difficult to accurately measure in the wind tunnel test. Moreover, threshold friction velocity, which is an important property of snow particles, is affected by many factors. Thus, its value generally varies within a certain range. However, the value of the threshold friction velocity directly affects the results of both the numerical simulation and wind tunnel test. Hence, the effects of threshold velocity of snow and silica sand on the duration of the prototype are analyzed in this section. According to the similarity parameter of wind velocity in Eq. (8), under the condition of the fixed wind velocity of the test, the

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

50 Precipitation of snow (mm)

Threshold velocity of silica sand (u*t)m (m/s) 0.24 0.26 0.28

7

6

5

10 Precipitation of snow Wind velocity at roof height

40

8

30

6

20

4

10

2

Simulated value 0

0 1

2

3

4 0.16

0.18 0.20 0.22 Threshold velocity of snow (u*t)p (m/s)

0.24

4 5 Date (day)

0.8

0.24 0.26 0.28

30

6th day 7th day

0.4

0.2

20

4th day 5th day

0.6 Snow depth (m)

Wind duration (prototype, unit: day)

1st day 2nd day 3rd day

40

7

6

Fig. 18. Meteorological data of a real snow storm.

Fig. 16. Influence of threshold velocity of silica sand and snow on wind velocity at roof height.

Threshold velocity of silica sand (u*t)m (m/s)

Wind velocity at roof height (m/s)

8 Prototype wind velocity at roof height (m/s)

103

Wind direction

Simulation value 10

0.0 0

0 0.16

0.18 0.20 0.22 Threshold velocity of snow (u*t)p (m/s)

wind velocity of the prototype at roof height is inversely proportional to the threshold friction velocity of silica sand and proportional to the threshold friction velocity of snow. The influence of threshold velocity of silica sand and snow on wind velocity at roof height is given in Fig. 16. The threshold friction velocity of silica sand and snow greatly influence the wind velocity of the prototype at roof height (Fig. 16), which would affect the results of snow transport rate on roofs. The wind duration of prototype is determined by Eq. (11). For the wind duration of the wind tunnel test is fixed, the wind duration of the prototype also changes according to the snow transport rate of prototype and that of model. If the prototype Qp is large, short duration of prototype then is obtained to reach the certain level of erosion/accumulation; whereas if the prototype Qp is small due to a lower wind velocity, longer duration of prototype is needed to reach the same level of erosion/accumulation. The change ensures that the similarity parameter of wind duration is the same in the numerical simulation and the wind tunnel test. The influence of threshold velocity of snow and silica sand on the wind duration of the prototype can be obtained through the proposed numerical method (Fig. 17). The prototype time duration of wind decreased as the threshold friction velocity of silica sand decreased or the threshold friction velocity of snow

2 Relative position X/H

3

4

Fig. 19. Snow distribution on the flat roof.

0.24

Fig. 17. Influence of threshold velocity of sand and snow on wind duration.

1

particles increased. Given that the threshold friction velocity of silica sand decreased or the threshold friction velocity of snow particles increased, the prototype wind velocity at roof height would increase. This increase caused the snow on the flat roof to become more erodible, and snow transport rate Q increased. According to Eq. (11), the wind duration in the wind tunnel test, which was required to reach the same pattern of snow cover, would be reduced. The threshold friction velocity of silica sand and snow were both 0.24 m/s when the wind velocity of prototype at roof height increased to 7.0 m/s. Moreover, the wind duration of the prototype was 4.6 days, which is only 45% of the result in the previous simulation. Therefore, the snow on roofs may accomplish the redistribution in less time when the wind velocity of prototype increases. 4.5. Case study Meteorological data can be easily applied by using the proposed numerical method to predict the redistribution of snow on a roof surface. A case study was performed to simulate the snow redistribution on a flat roof in a real snowstorm. The duration of the real snowstorm is 7 days. The daily snowfall and wind velocity are shown in Fig. 18. Total snowfall in the 7 days is 94.5 mm water equivalent. The total snowfall is equivalent to 63.0 cm of snow depth when the snow density is 150 kg/m3 for new snowfall.

104

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

test and numerical simulation validated the accuracy of the proposed numerical method to a certain degree. The conclusions of the study are as follows:

Table 7 Exposure coefficients for flat roofs (Meløysund et al., 2007). References

Otstavnov and Rosenberg Lutes Taylor O’Rourke et al. Høibø Løberg Com. Eur. Comm

Exposure coefficient Sheltered

Semi-sheltered

Windswept

0.98 0.90 0.80 0.76 0.82, 0.62 – 0.90

0.72 0.60 – 0.57 – 0.55 0.74

0.46 0.30 – 0.55 – 0.27 0.58

Maximum snowfall occurs in the first day, and the water equivalent of snowfall is 45 mm. The average wind velocity at the height of the roof is 5.6 m/s, which is 6.8 m/s when converted to the wind velocity at a height of 10 m. The maximum wind velocity of 7.46 m/s at roof height occurs in the third day. This velocity is 9.0 m/s when converted to the wind velocity at a height of 10 m. The physical properties of snow used herein are the same as those in the previous numerical simulation given in Table 1. The research object is the same flat roof as described above. The wind direction was assumed to be parallel to the span. According to the meteorological data, the snowstorm was divided into seven phases; one day corresponds to one phase. In the current simulation, the precipitating snow is supposed to pave uniformly with a repose angle on the snowpack surface of previous days. Fig. 19 shows the snow distribution obtained from the numerical simulation. Snowfall gradually accumulated, and the snow depth on the flat roof gradually increased. However, the accumulation of snow on the roof surface became unevenly distributed under wind action. Although the largest snowfall occurred in the first day, the distribution of snow on the roof was almost uniform and a certain angle of repose can be seen at the edge of the roof since the wind velocity was small. The maximum wind velocity appeared on the third day. The snow on the roof began to exhibit significant uneven distribution. Severe erosion occurred in the middle-rear part of the roof compared with that in the front region. The trend became more apparent in the following days. The average snow depth was 52.9 mm when the snowstorm ended. The depth was approximately 84.0% of the total snowfall. That means only 84.0% of snow accumulated on the roof surface, and 16.0% of fallen snow was blown off the roof under wind action. In the real environment, the wind direction and wind velocity are always changing with time. Hence it is difficult to obtain daily snow distribution on a flat roof under a stable wind environment. However, some papers measured the snow loads on roofs through field observations and presented them as the exposure coefficients of snow loading for flat roofs, which are the ratio of roof snow loads to ground snow loads. The measuring results could be seen in Table 7 (Meløysund et al., 2007). For the case study in this paper, the final exposure coefficient is 0.84, in which the wind velocity is relatively not very high and the simulated duration is about 7 days. The exposure coefficient of 0.84 seems somewhat reasonable when compared with the results from the field observations.

5. Conclusion The study proposed a quasi-steady numerical method to simulate snow redistribution on a roof surface, and meteorological data were used in this method. The method was applied to a flat roof. A scale model test was conducted for this flat roof in a wind tunnel, in which silica sand was employed to model snow particles. The comparison between the results from the wind tunnel

(1) A quasi-steady numerical method is proposed to calculate snow redistribution on the roof surface. The process of snow drifting on the roof surface is divided into several phases according to the characteristics of meteorological data. In each phase, the erosion/deposition of snow caused by snow drifting is regarded as a steady state. The influence of the change in the snow profile on flow field is considered by updating the snow boundaries in the next phase. The proposed method has higher efficiency compared with the traditional transient method of CFD. (2) Snow on the flat roof mainly eroded under wind action. As snow drifting continued, a mound of snow formed near the front edge, thereby blocking flow downwind. Friction velocity gradually decreased, and the transport rate declined. Snow redistribution on flat roof was more significant in the early stage of snow drifting. (3) Both the numerical simulation and wind tunnel test predicted deposition in a small area on flat roof. The deposition in the numerical simulation occurred at the front region of the flat roof, whereas the wind tunnel test predicted that the deposition was in the middle region and close to the front region. The deposition during the numerical simulation occurred earlier than in the wind tunnel test. (4) The magnitude of threshold velocity of snow and silica sand could significantly affect the wind duration of the prototype. The increase of threshold friction velocity of snow and the decrease of threshold velocity of silica sand both reduced the wind duration of prototype. In return, the prototype time became longer. In the current wind tunnel test, only snow drifting without snowfall has been modeled, which means the proposed method is suitable for snow cover which is firstly accumulated and then redistributed. During snowfall, there is no limitation to the upwind fetch for snow transport on building roofs, the snow drifting could be different. In the future, further studies will be conducted to study the difference between snow drifting on building roofs with snowfall and that without snowfall.

Acknowledgments This project is jointly supported by the National Natural Science Foundation of China (51278368, 51478359), the Ministry of Science and Technology of China (SLDRCE14-B-10) and the Fundamental Research Funds for the Central Universities, which are gratefully acknowledged.

References American Society of Civil Engineers, 2010. Minimum Design Loads for Buildings and Other Structures, ASCE/SEI 7–10, USA. Anno, Y., 1984. Requirements for modeling of a snowdrift. Cold Reg. Sci. Technol. 8 (3), 241–252. Beyers, J.H.M., Harms, T.M., 2003. Outdoors modelling of snowdrift at SANAE IV Research Station, Antarctica. J. Wind Eng. Ind. Aerodyn. 91 (4), 551–569. Beyers, J.H.M., Sundsbø, P.A., Harms, T.M., 2004. Numerical simulation of threedimensional, transient snow drifting around a cube. J. Wind Eng. Ind. Aerodyn. 92 (9), 725–747. EN1991-1-3, 2003. Euro code 1—Actions on Structures—Parts 1–3: General Actions —Snow Loads. BSI, London, pp. 1–57. Flaga, A., Kimbar, G., Matys, P., 2009. A new approach to wind tunnel similarity criteria for snow load prediction with an exemplary application of football

X. Zhou et al. / J. Wind Eng. Ind. Aerodyn. 153 (2016) 92–105

stadium roof. In: Proceedings of the 5th European & African Conference on Wind Engineering, Florence, Italy. International Organization for Standardization, 2013. Bases for Design of Structures —Determination of Snow Loads on Roofs, third edition, ISO4355, Switzerland. Isyumov, N., Mikitiuk, M., 1990. Wind tunnel model tests of snow drifting on a twolevel flat roof. J. Wind Eng. Ind. Aerodyn. 36 (2), 893–904. Kind, R.J., Murray, S.B., 1982. Saltation flow measurements relating to modeling of snowdrifting. J. Wind Eng. Ind. Aerodyn. 10 (1), 89–102. Kind, R.J., 1986. Snowdrifting: A review of modeling methods. Cold Reg. Sci. Technol. 12 (3), 217–228. Kuroiwa, D., Mizuno, Y., Takeuchi, M., 1967. Micromeritical properties of snow. In: Proceedings of the International Conference on Low Temperature Science: Physics of Snow and Ice, vol. 1 (2), pp. 751–772. Meløysund, V., Lisø, K.R., Hygen, H.O., Høiseth, K.V., Leira, B., 2007. Effects of wind exposure on roof snow loads. Build. Environ. 42 (10), 3726–3736. Ministry of Construction of the People's Republic of China, 2006. Load Code for the Design of Building Structures (GB50009-2012). China Building Industry Press, Beijing. Naaim, M., Naaim-Bouvet, F., Martinez, H., 1998. Numerical simulation of drifting snow: erosion and deposition models. Ann. Glaciol. 26, 191–196. National Research Council of Canada, 2005. National Building Code of Canada, Canada. O'Rourke, M., DeGaetanob, A., Tokarczyk, J.D., 2004. Snow drifting transport rates from water flume simulation. J. Wind Eng. Ind. Aerodyn. 92 (14–15), 1245–1264. O'Rourke, M., DeGaetanob, A., Tokarczyk, J.D., 2005. Analytical simulation of snow drift loading. J. Struct. Eng. 131 (4), 660–667. Smedley, D.J., Kwok, K.C.S., Kim, D.H., 1993. Snowdrifting simulation around davis station workshop, antarctica. J. Wind Eng. Ind. Aerodyn. 50, 153–162. Sundsbø, P.A., 1998. Numerical simulations of wind deflection fins to control snow accumulation in building steps. J. Wind Eng. Ind. Aerodyn. 74 (98), 543–552. Tabler, R.D., 1994. Design Guidelines for the Control of Blowing and Drifting Snow. Strategic Highway Research Program, National Research Council.

105

Tabler, R.D., 2003. Controlling blowing and drifting snow with snow fences and road design. National Cooperative Highway Research Program Project. NCHRP Project 20-7(147). Takeuchi, M., 1980. Vertical profile and horizontal increase of drift snow transport. J. Glaciol. 26, 481–492. Thiis, T.K., Ramberg, J.F., 2008. Measurements and numerical simulations of development of snow drifts on curved roofs. In: Proceedings of the Snow Engineering VI, Whistler, Canada. Tominaga, Y., Okaze, T., Mochida, A., 2011a. CFD modeling of snowdrift around a building: an overview of models and evaluation of a new approach. Build. Environ. 46 (4), 899–910. Tominaga, Y., Mochida, A., Okaze, T., 2011b. Development of a system for predicting snow distribution in built-up environments: combining a mesoscale meteorological model and a CFD model. J. Wind Eng. Ind. Aerodyn. 99 (4), 460–468. Tsuchiya, M., Tomabechi, T., Hongo, T., Ueda, H., 2002. Wind effects on snowdrift on stepped flat roofs. J. Wind Eng. Ind. Aerodyn. 90 (12–15), 1881–1892. Uematsu, T., Nakata, T., Takeuchi, K., Arisawa, Y., Kaneda, Y., 1991. Threedimensional numerical simulation of snowdrift. Cold Reg. Sci. Technol. 20 (1), 65–73. Wang, W, Liao, Haili, Li, Mingshui, 2013. Numerical simulation of wind-induced roof snow distributions based on time variable boundary. J. Southwest Jiaotong Univ. 05, 851–856 þ 967 (in Chinese). Yang, Y., Gu, M., Chen, S., et al., 2009. New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering. J. Wind Eng. Ind. Aerodyn. 97 (2), 88–95. Zhou, X., Li, X., 2010. Simulation of snow drifting on roof surface of terminal building of an airport. Disaster Adv. 3 (1), 42–50. Zhou, X., Hu, J., Gu, M., 2014. Wind tunnel test of snow loads on a stepped flat roof using different granular materials. Nat. Hazards 74 (3), 1629–1648.