Reliability analysis of a large curved-roof structure considering wind and snow coupled effects

Ensuring the safety of large curved-roof structures subjected to blizzards in cold regions remains challenging. This paper assesses the reliability of a large curved-roof structure considering the coupled effects of wind and snow. To this end, a two-way coupled simulation considering the interaction between snow particles and turbulent wind is first carried out using a computational fluid dynamics scheme and a spectral representation method based on the wavenumber-frequency joint power spectrum and stochastic harmonic function; thereby, the uneven snow distribution and stochastic wind field on a large curved roof are revealed. The probability density evolution method is then adopted to conduct a stochastic response and reliability analysis of the structure. The simulation results reveal that a credible uneven distribution of snow on the large curved roof can be obtained using the suggested snowdrift simulation scheme


INTRODUCTION
Urbanization has received much attention from researchers worldwide in recent decades. With an increasing urban population, there has been considerable development of non-residential buildings. Many public buildings have been constructed with a large spatial configuration; e.g., industrial plants, exhibition halls, airport terminals, indoor gymnasiums and concert halls.
Compared with traditional buildings having normal spans, modern large-span spatial structures have more spacious and more open interior spaces and provide more possibilities for indoor activities. However, the increasing structural span means that the design of large-span spatial structures has to consider two issues. (i) The sensitivity of large structures to the gravity load is high, which is of particular concern for large-span spatial structures located in cold regions, where the snow load is a control load of the structural design [1] ; and (ii) Large structures have high structural flexibility and low natural vibration frequency and are thus prone to wind-induced damage under strong winds or winds with dominating low-frequency fluctuating components [2] . In fact, many instances of snow-induced collapse and wind-induced damage have been recorded for large-span spatial structures. For instance, in August 2005, the hurricane Katrina struck the city of New Orleans in the southern United States (with an average wind speed of 40.2 m/s), tearing and blowing over the ethylene propylene diene monomer rubber roof of the Louisiana Superdome. In April 2016, winds gusting at magnitude 13 (with an average wind speed of 39.9 m/s) hit the Changzhou area in Guangzhou, China, overturning the iron roof of the Guangwai Gymnasium in University Town. In January 2008, rarely heavy snowfall fell in Shanghai, China, and damaged and collapsed many houses, parking garages and warehouses. In February 2006, a blizzard hit Moscow, Russia, and snow accumulated continually on the roof of the Basmany Market on the outskirts of the city over a period of days and finally collapsed the roof, killing more than 50 people. In December 2010, strong winds and blizzards struck Minneapolis in the United States, resulting in the accumulation of more than 61 cm of snow on the roof of the home stadium of the Vikings football team, eventually destroying the roof.
The above-mentioned structural failure events suggest that there has been insufficient research on wind-induced damage, snow-induced damage and their coupled effects on large-span spatial structures. It is thus important to structural engineering practice to determine the global reliability of large-span spatial structures considering the coupled effects of wind and snow. Relevant studies on the distribution of snow and wind loads on a roof [3] , statistical analysis of structural responses [4] and reliability assessment [5] have been conducted. However, although the associated work has made great progress, there remain critical issues to address. (i) Most numerical simulations of the snow distribution on building roofs have been carried out on small models and only occasionally calibrated with the results of wind tunnel tests or field measurements [6,7] , and the numerical simulation of prototype models of building roofs has not received enough attention; (ii) Previous wind-induced snowdrift models are mostly based on the assumption of one-way coupled effects [3,8] , ignoring the reaction of moving snow particles to turbulent wind and its further effect on snow particle motion, and the simulation of snow-and wind-load distributions considering the two-way coupled effects of snow particles and turbulent wind needs to be addressed; and (iii) The vibration response analysis of large curved-roof structures subjected to snow and wind loads in previous studies has been largely based on the equivalent static wind load or frequency-domain analysis [9,10] , the argument for assessing the structural performance has been a deterministic index or second-order statistical moment, and the accurate reliability of structures has not been secured owing to a lack of sufficient probabilistic information of the responses.
The present study thus conducts the reliability assessment of a large curved-roof structure through time-history analysis considering the coupled effects of wind and snow. As a practical example of structural design, a large curved-roof structure located on the outskirts of a city in a cold region of China considering the design return period of 100 years is used. The contributions of the study are as follows.
• Wind and snow pressure coefficients of a large curved-roof structure are determined through computational fluid dynamics simulations considering two-way coupled effects of the snow particles and turbulent wind.
• A two-dimensional fluctuating wind field for time history analysis of a large curved-roof structure is simulated by adopting a refined spectral representation method combining the wavenumber-frequency joint power spectrum and stochastic harmonic function.
• A stochastic response analysis and performance evaluation of a large curved-roof structure subjected to an uneven snow distribution and fluctuating wind field are carried out using the probability density evolution method.
• The reliability of the large curved-roof structure under stochastic wind, snow and self-weight loads with a 100-year return period is assessed.
The rest of this paper is organized as follows. Section 2 simulates the wind pressure and snow pressure coefficients by introducing the two-way coupling of snow particles and turbulent wind and also simulates the fluctuating wind field using a spectral representation method based on the wavenumber-frequency joint power spectrum and stochastic harmonic function. Section 3 presents the results of stochastic response analysis of the structure based on the probability density evolution method (PDEM). Section 4 conducts a dynamic reliability assessment of the large curved-roof structure. Concluding remarks are given in Section 5.

Overview of a large curved-roof structure
The considered large curved-roof structure serves as a terminal on the outskirts of a city in a cold region of China. A Google map of the terminal and its adjacent buildings is shown in Figure 1A. The main section of the terminal has a longitudinal span of 204 m (with a longitudinal dimension of 216 m), a transverse span of 75 m and a maximum overhang width of 7.5 m (with a transverse dimension of 90 m), which represents a typical large-span spatial grid structure; see the plan of the terminal in Figure 1B. The building has three stories, with the first and second stories having a reinforced-concrete frame structure and a height of 7.8 and 8.2 m, respectively, whereas the third story is a steel frame structure with a bidirectional curved and doublelayer reticulated shell roof, with the floor height being 21 m; see the elevation of the terminal in Figure 1C. The roof comprises plane trusses, restraint trusses, steel pipe columns and steel pipe diagonal braces and varies in elevation by approximately 5 m. All connections of the steel pipe columns and steel pipe diagonal braces with the bottom concrete frame column and the top steel roof are cast steel nodes, which are approximately hinged. The region of the terminal has a temperate continental climate, including a long cold winter with severe wind and snow conditions. According to Chinese Load Code for the Design of Building Structures (GB50009-2012), the average wind pressure and average snow pressure for a 100-year return period are 0.7 and 1.0 kN/m 2 , respectively.

Simulation of wind pressure and snow pressure coefficients
Wind-induced snow drift is a phenomenon of interactions between the near-ground wind and snow particles at the surface of ground snow and floating in the air. The classical scheme of wind-induced snow drift simulation thus includes the simulation of the near-ground wind field and the simulation of snow drift motion.
The Reynolds averaged Navier-Stokes (RANS) equations are often used to model the near-ground wind field by considering flow properties such as the wind velocity, wind pressure, flow density, wall shear stress and turbulent characteristics. To consider the effect of snow particle motion on the air turbulence characteristics and the air phase-averaged momentum, an extended two-way coupled scheme for wind-induced snow drift simulation is used. In this scheme, different damping source terms are introduced into the RANS equations and the transport equations of the turbulent kinetic energy and energy dissipation rate of the standardturbulence model [11] : where ( = 1, 2, 3) denotes the three-dimensional spatial coordinates;¯( = 1, 2, 3) denotes the average threedimensional components of the wind velocity; = 1.225 kg/m 3 is the flow density;¯is the average pressure, which depends on the design conditions of the local site; t = 1.619 × 10 −5 is the flow kinematic viscosity; k is the production of the turbulent kinetic energy , which is assumed to equal the energy dissipation rate in this study; 1 = 1.44 and 2 = 1.92 are model constants; and k = 1.0 and = 1.3 are Prandtl constants relating to the turbulent kinetic energy and energy dissipation rate , respectively [12] . i , k , and are respectively the damping source terms in the above three equations, which relate to the mass concentration of the snow drift, the relative velocity between moving the snow particles and air, the characteristics of snow particles and other factors [11] : Here, = 3Φ /(2 s s ) is the representative area of the obstacle related to the drag force for a mass concentration of the snow drift Φ and the diameter s = 0.5 mm and density s = = 200 kg/m 3 of the snow particles. is the fluid volume within a unit volume; f = 5.0 × 10 3 is the aerodynamic drag coefficient of the moving obstacle; p = 0.9 is a model coefficient related to the particle speed in terms of the wind speed; and p = 1.0 is a coefficient in the expression for .
The simulated snowdrift is often quantified as the flux of snow erosion/deposition, which is mainly contributed by snow particles in both the saltation and suspension layers [13] : The flux of snow erosion and the flux of snow deposition are expressed as [14,15] where 0 = 7 × 10 −4 is the empirical coefficient of snow erosion/deposition; * t is the threshold friction velocity, which is considered to be lower ( * t = 0.07-0.25 m/s) for fresh, loose and dry snow and during snowfall and to be higher ( * t = 0.25-1.0 m/s) for older, wind-hardened, dense, and/or wet snow for which interparticle bonds and cohesion forces are strong [16] ; * is the friction velocity at the snow surface, which relates to the wind velocity at a specified height obtained from the simulation of the near-ground wind field and the roughness length; and f = 0.2 m/s is the snowfall velocity. The mass concentration of the snow drift Φ is expressed by a transport equation of the snow drift based on convection-diffusion theory [13] : where t = 1.32 is the coefficient of turbulent diffusion, denoting the action of turbulent wind upon the suspended snow particles. The third term on the left side expresses convection along the height ( ) direction and reflects the effect of snowfall on the mass concentration of the snow drift in the suspension layer.
In representing the uneven distribution of wind and snow loads, the wind pressure and snow pressure coefficients are defined as [17,18] where w ( , * ) is the wind pressure coefficient; (in units of Pa) is the simulated wind pressure at point on the roof; the snow-drift action time * is taken as 48 h herein because the lowest root-mean-square error (RMSE) of simulated points between time steps that satisfies the convergence criterion is obtained for this time step, where the RMSE values relative to the previous two time steps are both less than 0.01; ∞ , ∞ and ∞ are respectively the wind pressure, flow density and average wind velocity at the reference position in the far field of the airshed and taken as 700.16 Pa, 1.225 kg/m 3 , and 33.81 m/s (with the reference position being a 10-m height and the exponential wind profile being used as the boundary condition of the wind velocity inlet); s ( , * ) is the snow pressure coefficient; ℎ( , * ) is the simulated snow depth at point on the roof; and ℎ 0 is the ground snow depth, which is taken as 51 cm according to the average snow pressure for a 100-year return period in the design code.
The finite volume method is applied to solve the partial differential equations in Equations 1-3 and 10 using the finite element software ANSYS 14.0 FLUENT. The airshed modeling and grid division of the flow domain of the large curved roof subject to wind blowing in the most unfavorable direction (i.e., along the transverse span of the terminal building because this span is just one-third of the longitudinal span) are shown in Figure 2. The distances from the centroid of the horizontal projection plane of the terminal building to the wind inlet, outlet and top surface are set at 10 , 20 and 15 , where and are the transverse dimension ( = 90 m) and height ( = 37 m) of the terminal building, respectively, and the airshed blocking rate is 1.06%. Tetrahedral elements are used to generate an unstructured mesh. A suitable grid resolution for the airshed is ensured by controlling the growth ratio of the cells from the buildings to the surrounding airshed, and the total number of cells is 1,254,266.   conditions of wind and snow loads for a 100-year return period, the snow pressure on the roof has an obvious uneven distribution. The snow pressure coefficient is within the range [0.64, 1.02], where a value exceeding 1.0 indicates snow deposition and mainly occurs in the upper-right corner of the roof near the eave. However, most of the roof area, especially the area near the ridge line in which there is much snow erosion, is subject to snow erosion. The two-way coupled simulation and one-way coupled simulation were compared in the authors' previous study [19] . That study revealed that the two-way coupled simulation was more accurate relative to field observations (having 7.9% lower simulation error) and more efficient (having 37.1% lower computational cost) than the one-way coupled simulation.

Simulation of the fluctuating wind field
The wind pressure coefficient for the roof shown in Figure 4A is attained in steady computational fluid dynamics (CFD) simulation, and it is essentially an average value of the real wind pressure field. In reality, carrying out a transient CFD simulation is time-consuming. A practical approach is to decompose the wind field into a deterministic average wind component and a stochastic fluctuating wind component through Reynolds decomposition. The deterministic average wind component can be simulated using the RANS equations (see Equations 1-3), whereas the simulation of the stochastic fluctuating wind component often resorts to a spectral representation.
The elevation variation of the roof is far less (by two orders of magnitude) than the plan size, and the threedimensional spatial fluctuating wind field can thus be approximately presented as a two-dimensional fluctuating wind field. Moreover, to reduce the computational cost, the same sub-regions of the roof used for the wind pressure coefficient distribution are adopted, where the center point is taken as the representative point and other points in the same sub-region are assumed to have the same fluctuating wind as the representative point. Considering the low efficiency of the conventional scheme of spectral representation methods in wind field simulation, the spectral representation method based on the wavenumber-frequency joint power spectrum and stochastic harmonic function is used in this study [20,21] : where ( , , ) is the fluctuating wind velocity at point ( , ) on the roof plane; (W−F) is the wavenumberfrequency joint power spectrum; , and are the numbers of discretized points in the wavenumber domains comprising and dimensions and the frequency domain, respectively, all of which are set at 50 in this study; ( , , ) denotes the coordinates of the -th point in the wavenumber-frequency domain; is the representative volume of the -th point in the wavenumber-frequency domain; and , ( = 1, 2, 3, 4) denotes mutually independent random variables uniformly distributed on [0, 2 ]. The discretized scheme of the wavenumber-frequency domain is expressed as [21] =  Equation 14 from the spectral representation method based on the classical wavenumber-frequency joint power spectrum, which is referred to as the wavenumber-frequency joint power spectrum and stochastic harmonic function based spectral representation method. It is also noted from Equations 19 and 20 that an exponential formula is used in the discretized scheme of the wavenumber-frequency domain so as to retain dense point sets close to the base point and sparse point sets far from the base point.
This study uses the Davenport spectrum and coherence function model to derive the wavenumber-frequency joint power spectrum as [21] (W−F) where is the mean wind velocity at the roof height corresponding to the 100-year return period, which is calculated as 40.35 m/s; * is the friction velocity, which is calculated as 2.680 m/s for the geomorphic conditions of the local site; is the circular frequency; and and are the decay parameters along and directions, respectively, both of which are taken as 10 [22] .
To validate the accuracy of the wavenumber-frequency joint power spectrum and stochastic harmonic function based spectral representation method, the simulated fluctuating wind velocities at three representative points on the roof, namely the center points of the sub-regions Z9 (labeled P1), Z10 (labeled P2) and Z4 (labeled P3), are considered. Representative time histories of the fluctuating wind velocities at P1, P2 and P3, and their auto-power spectra and coherence functions are shown in Figures 5 and 6. It is seen that the auto-power spectra of the simulated fluctuating wind velocities at the different spatial points and their coherence functions are in good agreement with the targets; see the Davenport spectrum and the related coherence function, which indicates that the wavenumber-frequency joint power spectrum and stochastic harmonic function based spectral representation method can be used to accurately simulate the two-dimensional fluctuating wind field for a large curved-roof structure.

PDEM
There is randomness associated with the fluctuating wind field, and the dynamical system of the large curvedroof structure subjected to the coupled effects of the wind and snow is thus a typical random vibration problem, which needs to be solved using probabilistic methods. In this study, the PDEM, which has been widely adopted in recent years, is used to efficiently and accurately estimate the statistical information of the structural response [23] .
Without loss of generality, we write the equation of motion of an -degree-of-freedom system as where M is the mass matrix; f (·) is the nonlinear restoring force vector; and X, X and X are the acceleration, velocity and displacement of the structure relative to the ground. F(·) is the stochastic external load, where is a stochastic vector representing the randomness associated with the external loads.
As regards the physical quantities Z = =1 T of interest (e.g., the displacement, velocity, acceleration, stress and internal force), the augmented system (Z, ) is probability preserved because all the random factors are involved. The generalized probability density evolution equation (GDEE) is thus [23] Θ ( , , ) + where is a sample of the stochastic vector and is the velocity of the physical quantity of interest. It is noted that the GDEE reveals the intrinsic connections between a stochastic dynamical system and its deterministic counterpart. To solve the GDEE, an initial condition is given by where (×) is the Dirac delta function; 0 is the initial time; and 0 is a deterministic initial value.

Solving Equation 25 under the initial condition given as Equation 26
, the instantaneous probability density function (PDF) of ( ) is obtained as where Ω Θ is the sample space of the stochastic vector .
The numerical procedure for solving the GDEE is as follows. Step 1: Select the representative point set P @ = 1, , 2, , L, s, | = 1, 2, L, pt in the domain W Q , where pt is the total number of selected points. The representative sub-domains, as a partition of W Q , are determined using Voronoi cells [24] .
Step 2: For the prescribed Q = , solve Equation 23 using a deterministic time integration method to evaluate ( , ), where = D , = 0, 1, 2, L, and D is the time step in the numerical integration.

Finite element modeling and analysis
Using the finite element analysis platform SAP2000 v19, the deterministic structural dynamic analysis involved in Step 2 in Section 3.1 is carried out. In the concrete frame-steel roof hybrid structure, the lower concrete structure has higher rigidity than the top large-span steel roof structure, and the structural internal force and displacement are generally small under wind and snow loads. To simplify the calculation, therefore, only the steel roof at the top of the terminal building is modeled using finite elements, and the bottom of the steel pipe column and steel diagonal bracing adopts hinge supports, of which the number is 58. For the steel pipe columns, steel diagonal bracings, upper chords, lower chords and diagonal rods of the steel roof, frame elements are used, and the section size of each element is selected according to the actual design size. The steel roof support shell adopts massless shell elements for roof load transfer. By doing so, the finite element model of the large curvedroof structure includes a total of 4241 nodes, 16,073 frame elements and 3083 massless shell elements, with a total structural mass of 2924.6 t. The three-dimensional finite element model is shown in Figure 7A.
To clarify the basic dynamic characteristics of the large curved-roof structure, structural modal analysis is carried out using the Rayleigh-Ritz vector method. It is found that the first four orders of the modal frequencies are 0.854, 0.767, 0.553 and 0.463 s, and they sequentially correspond to the modal characteristics of horizontal vibration along the axis, horizontal vibration along the axis and torsional and vertical vibrations along the axis. The first mode is shown in Figure 7B. It is seen that the basic natural vibration period of the structure is higher than 0.25 s for the basic natural vibration period of flexible structures according to Chinese Technical Regulation for Spatial Grid Structures (JGJ7-2010); therefore, the wind-induced vibration problem should be considered in the structural analysis. Additionally, the natural vibration periods of the structure are distributed in a narrow range; i.e., the structure has densely distributed frequencies. Thereby, the contribution of each mode to the dynamic response of the structure is equally important, and ignoring the contribution of high-order modes will inevitably lead to large computational errors. Moreover, the overall vertical stiffness of the structure is high, and it is predicted that the vertical deflection of nodes of the roof under the loads will be smaller than that for a flexible roof of the same size.
As a practical situation, the large curved-roof structure is simultaneously subjected to the wind load comprising average and fluctuating wind components, the snow load unevenly distributed on the roof and the self-weight load. Structural analysis should consider a certain combination of the effects of these loads. In this study, they are in weighted linear combination according to Chinese Load Code for the Design of Building Structures (GB50009-2012). Figure 8A and B show the steady member axial stress and node vertical deflection of the large curved-roof structure under representative wind, snow and self-weight loads (at a time instant of 300 s).
It is seen that under the representative wind, snow and self-weight loads, the axial stress acting on members is in the range of -198 to 154 MPa, and the lower roof chord and roof inclined rod experience higher stresses as same as the range in the structural level, whereas the upper roof chord and roof support column experience lower stresses in the range of -141 to 106 MPa. Under the representative wind, snow and self-weight loads, the node vertical deflection is in the range of -71.5 to 16.5 mm and the deflection is a maximum at the nodes in the mid-span area of the roof, whereas the overhanging part of the roof has a small upward deflection. Therefore, the member axial stress of the large curved-roof structure reaches a maximum of 198 MPa, and the node vertical deflection reaches a maximum of 71.5 mm. Considering that the strength of the structural member is 345 MPa, the maximum member axial stress is approximately 0.574 times the maximum allowable stress. The threshold of the node vertical deflection is 1/250 of the horizontal span of the structure according to the design provision (i.e., 300 mm), and the maximum node vertical deflection is approximately 0.238 times the threshold. The member axial stress is closer to its limit value than the node vertical deflection, and the risk of the member yielding is thus higher than that of node deflection, which is in agreement with the results of structural modal analysis.

Statistical analysis of the stochastic responses of structures
In this study, there is a random vibration problem due to the randomness associated with the fluctuating wind field. Using the PDEM, a total of 300 samples of the fluctuating wind field are simulated and used for the stochastic response analysis and reliability assessment of the large curved-roof structure. For the purpose of illustration, the stochastic responses of sub-region Z10 are addressed. Time histories of the means and standard deviations of the axial stress of Member #14294 and the vertical deflection of Node #3019 in this sub-region are presented in Figure 9A-D, respectively. The simulated fluctuating wind field is a stationary state, and the time histories of statistical moments within just the first 300 s are thus presented to clearly show the initial non-stationary stage of the structural responses.
It is seen that the means and standard deviations of the member axial stress and node vertical deflection quickly reach a steady state; i.e., within 50 periods of the basic natural vibration of the structure. Additionally, the average wind component of the stochastic wind load has a suction effect on the roof, which partially offsets the structural response relating to the self-weight load and the unevenly distributed snow load, and the average wind component thus reduces the overall response of the structure. Meanwhile, the fluctuating wind component has an appreciable dynamic effect on the responses of the structure. Furthermore, the probability densities of the axial stress of Member #14294 at 100, 200, and 300 s and in the time interval [150, 160] s are presented in Figure 10. It is seen that the probability density of the structural response in the stage of the steady state only slightly changes and approximately follows a normal distribution. This can be explained in that the simulation of the stochastic wind field is a Gaussian process and the dynamic analysis using finite element software reveals the linear state of the structure.  Figure 10 also shows that although the mean of the member axial stress is far less than the maximum allowable stress, the member still faces the possibility of yielding because the axial stress varies appreciably. It is thus necessary to carry out reliability assessment in terms of the member axial stress in determining the safety of the structure.

PDEM-based reliability method
Without loss of generality, the associated reliability of the first-passage problem of structural dynamical systems is defined as where Pr{×} is the probability operator of a random event; ( , ) is the random event of physical quantities of concern such as the stress, internal force and deflection; Ω is the safe domain; and is the time duration corresponding to (·). When the structural system has multiple modalities of failure, the first-passage problem is expressed as where is the number of random events and { } =1 denotes the safe thresholds of the random events { ( , )} =1 . For the convenience of description, we refer to the single-passage level as a type-B barrier herein.
We construct an equivalent extreme-value event as a random variable: The reliability function in Equation 29 can then be defined alternatively as ( ) = Pr eq ( , ) ≤ 0 (32) Equation (32) can be rewritten as where eq ( , ) is the PDF of eq ( , ).
Therefore, the introduction of the equivalent extreme-value event criterion (see Equation 31) reduces the highdimension integral associated with the reliability of systems into the one-dimension integral [25] .
To derive the PDF of the equivalent extreme-value event, a pseudo-random process is first defined as which satisfies the conditions where denotes the generalized time and 0 , denote the initial and final time instants of the pseudo-random process, respectively.
Theoretically, arbitrary forms of the pseudo-random process ( ) can be assumed if they satisfy the conditions shown in Equation 35. This study uses the form ( ) = eq ( , )sin( ), where = 2.5 , = 1.
It is observed that the probabilistic information of the equivalent extreme-value event eq ( , ) is completely included in the final time instant of the pseudo-random process ( ), and they share the same randomness denoted by . Therefore, the augmented system ( ( ), ) is a probability conservation system that is governed by the equation for the generalized probability density evolution given as Equation 25; thereby, the PDF of the pseudo-random process (i.e., ( , )) can be derived. One then has

Reliability assessment of a structure
Using the PDEM and the associated reliability method, the PDFs and cumulative density functions (CDFs) of the equivalent extreme-value axial stress of members in the sub-regions Z1, Z8, Z10 and Z17 are derived, as shown in Figure 11. It is seen that all the curves related to the sub-regions under investigation are far below the maximum allowable stress (i.e., 345 MPa), indicating safety in these regions. It is also seen that among the four representative sub-regions under investigation, the sub-region Z10 has the largest member axial stress, followed by Z1, Z8 and Z17 in order. In view of Figure 4, the sub-region Z1 features lower wind suction and higher snow pressure and is thus expected to have appreciable member axial stress, whereas the other  three sub-regions feature almost similar wind suction and snow pressure and thus, compared with Z8 and Z17, the sub-region Z10 has a greater contribution from the fluctuating wind component to the member axial stress because this area is located in the rolling region of the curved roof. Furthermore, all the sub-regions are analyzed and their reliabilities are shown in Figure 12. Additionally, the PDF and CDF of the equivalent extreme-value axial stress of members of the roof structure are presented in Figure 13, Showing a safe structure in the global sense.
In summary, under stochastic wind, snow and self-weight loads with a 100-year return period, the equivalent extreme-value axial stresses of members in each sub-region and in the global structural system are both far less than the maximum allowable stress, and the large curved roof structure is deemed to be sufficiently safe.

CONCLUSIONS
Aiming at the reliability assessment of large curved-roof structures considering the coupled effects of wind and snow with a 100-year return period, this study conducted pertinent works, including a snowdrift simulation using computational fluid dynamics, a fluctuating wind field simulation using the wavenumber-frequency joint power spectrum and stochastic harmonic function based spectral representation method, and a stochastic response and dynamic reliability analysis of a structure using the probability density evolution method. The main remarks and conclusions of the study are listed as follows.
• Snowdrift simulation considering the two-way coupling of snow particles and turbulent wind provides a credible uneven distribution of snow on a large curved-roof structure.
• The wavenumber-frequency joint power spectrum and stochastic harmonic function based spectral representation method can be used to efficiently and accurately simulate the two-dimensional fluctuating wind field of a large curved-roof structure.
• The overall vertical stiffness of the structure is large, and the vertical deflection of nodes of the presented roof under loads is thus smaller than that of a flexible roof of the same size, and there is a higher risk of structural member yielding than structural node deflection.
• The average wind component of the stochastic wind load has a suction effect on the roof, which partially offsets the structural response for the self-weight load and the unevenly-distributed snow load and reduces the overall response of the structure. Meanwhile, the fluctuating wind component has an appreciable dynamic effect on the responses of the structure.
• The probability density of the structural response in the stage of a steady state changes only slightly and approximately follows a normal distribution because the simulation of the stochastic wind field is a Gaussian process and the dynamic analysis using finite element software reveals the linear state of the structure.
• Under stochastic wind, snow and self-weight loads having a 100-year return period, the member axial stress of the structural system is far less than the maximum allowable stress, and the large curved-roof structure is thus deemed sufficiently safe.
This study simulated snowdrift considering the two-way coupling of snow particles and turbulent wind. However, only steady simulations based on non-thermal and constant meteorological conditions were involved, and the simulation results may differ from reality. Therefore, unsteady simulations considering thermal effects such as roof heating and varying meteorological conditions including the wind velocity, wind direction, temperature and humidity need to be conducted in the near future. Additionally, a field measurement of the snow distribution on a large curved roof will be performed to validate the simulation results of the present study.