Probabilistic assessment of high-rise buildings with fluid viscous dampers under hurricane wind loads

Abstract

Current wind design codes incorporate turbulence through gust factors and rely on historical wind data, including tropical cyclones. While generally conservative, standard code wind profiles and spectra do not fully reproduce the vertical distribution and dynamic characteristics of hurricane winds, particularly in the supergradient region near the eyewall, and can sometimes underestimate tail risks, low-probability, high-impact events, as observed during Hurricane Otis in Acapulco (2023). This study probabilistically evaluates wind-induced vibrations in high-rise buildings with different lateral resisting systems equipped with fluid viscous dampers (FVDs), under non-tropical storm and tropical cyclone conditions. Along-wind loads were modeled in the time domain as stationary, multidimensional stochastic processes and analyzed using one million Monte Carlo simulations and Incremental Dynamic Analysis on the DelftBlue supercomputer. Statistical distributions of responses, bivariate dependence via copulas, and fragility curves were obtained. Results show that wind type, structural deformation mode, and damper properties significantly affect response distributions, correlation structures, and failure probabilities. FVDs effectively reduce structural dynamic response, improving serviceability, while increased shear stiffness further reduces fragility. Modeling hurricane winds as non-tropical storms can overestimate damper effectiveness. These findings provide insights for refining wind codes and designing high-rise buildings that remain safe and functional under extreme events.

Graphical abstract

Keywords

bivariate analysis / fluid viscous dampers / fragility curves / high-rise buildings / hurricane wind loads / Monte Carlo simulation

1 Introduction

High-rise buildings in coastal urban areas face increasing challenges due to extreme wind hazards associated with hurricanes, which threaten not only structural integrity but also urban functionality and safety. While hurricanes are expected to intensify as a consequence of climate change, current wind design codes often rely on simplified models that represent turbulent wind as equivalent static loads derived from non-tropical storms (NTS). This simplification may lead to underestimation of wind forces and overestimation of structural resilience in hurricane-prone regions. To address this, effective vibration control systems such as fluid viscous dampers (FVDs) offer a promising approach to mitigate dynamic responses of tall buildings under these complex wind loads. This study aims to probabilistically assess the performance of high-rise buildings equipped with FVDs under realistic hurricane wind conditions, providing insights for engineers and urban planners seeking to improve the resilience of urban infrastructure.

Recent studies indicate that, although the frequency of tropical cyclones (TC) may not increase, their intensity and severity are expected to rise due to global warming. Holland and Bruyère [1] reported a substantial regional and global increase of 25%–30% for each 1 °C of anthropogenic warming in the proportion of Category 4–5 hurricanes since 1975. Table 1 summarizes some of the most intense hurricanes over the past 15 years, highlighting their maximum intensity, landfall strength, locations, and associated economic losses (compiled from official reports and media sources). Some recent events, such as Hurricane Otis (2023), have demonstrated that exceptionally intense storms can retain considerable strength at landfall, causing multi-billion-dollar losses.

For practical engineering purposes, the mean wind velocity within the atmospheric boundary layer is typically modeled using empirical expressions for NTS, as included in wind design codes. These codes incorporate gust factors and are based on historical wind data, including TC, generally providing conservative estimates for structural design. However, standard formulations do not fully reproduce the supergradient wind region near the eyewall, where tangential wind speeds within the boundary layer can exceed those above it, as revealed by extensive field measurements using global positioning system (GPS) dropwindsondes and Doppler radar [2]. Observations from major hurricanes have shown complex boundary layer structures and strong radial variations of wind speed and turbulence that cannot be captured by traditional gradient wind or logarithmic models [3]. This behavior has been further explained through analytical models that incorporate multi-field parameter correlations and height-dependent eddy viscosity in the typhoon boundary layer, reproducing the non-exponential and supergradient characteristics of the vertical wind profile [4]. As a result, extremely rare, high-intensity events may exceed the reference wind speeds specified in design codes. For example, the maximum sustained wind speed recorded in Acapulco during Hurricane Otis was 270 km/h [5], whereas the Mexican design code [6] specifies 130 km/h for a 200 year return period. Probabilistic design maps should therefore continue evolving to better capture these tail risks, that is, low-probability, high-impact events, which is especially relevant in cities such as Acapulco and Cancun.

For the marine hurricane boundary layer, Vickery et al. [7] proposed a model for mean wind velocity based on six years of hurricane dropsonde data, while Snaiki and Wu [8] developed a semi-empirical model for landfalling hurricanes using radar and GPS dropsonde data from 1995 to 2012. Both models account for the supergradient wind region, unlike those derived for NTS. For non-tropical conditions, Petersen et al. [9] found that the Kaimal et al. [10] spectrum closely matches atmospheric turbulence, while the Von Kármán spectrum [11] better fits wind-tunnel data. However, none of these spectra reproduce the spectral energy characteristics of hurricane winds. According to Schroeder and Smith [12], Yu et al. [13], and Jung and Masters [14], hurricane winds exhibit greater energy at lower frequencies. To address this, Yu et al. [13] and Li et al. [15] proposed spectral models for TC based on field measurements from hurricanes in the US and China, respectively.

Traditional lateral-load-resisting systems, such as moment frames or braced frames, often provide insufficient stiffness for tall buildings subjected to strong winds or earthquakes. Alternatives include the use of shear wall cores or tube-in-tube systems, but these may significantly increase construction costs. FVDs provide an effective and economical alternative to enhance energy dissipation and reduce structural responses. These devices dissipate mechanical energy through the motion of a piston in a fluid-filled cylinder, converting kinetic energy into heat and reducing vibration amplitudes.

In the probabilistic analysis of wind-induced structural response, lognormal distributions are commonly assumed [16]. Chen et al. [17], however, analyzed wind-induced vibrations in structures with FVDs under non-tropical winds using a Gumbel distribution for wind field randomness, concluding that stochastic variability must be explicitly considered in FVD design. Zhou et al. [18] further demonstrated that higher vibration modes significantly influence wind-induced responses in high-rise buildings with FVDs. Sun et al. [19] and Li et al. [20] investigated hybrid and variable damping systems to improve FVD efficiency in tall buildings, though none of these studies addressed TC.

Regarding TC conditions, Li and Li [21] monitored a supertall building during 14 TC and reported satisfactory performance, although none reached hurricane intensity. Most FVD-related studies address seismic or non-tropical wind loads, while research on TC has focused on other slender structures [22,23]. Therefore, this study aims to bridge this gap by evaluating the probabilistic performance of high-rise buildings with FVDs under hurricane wind loads. The analysis involves fitting univariate distributions to simulated responses, assessing bivariate dependencies using copulas, and deriving fragility curves. The results provide insights for improving vibration control strategies and enhancing the reliability of tall buildings under extreme wind events.

2 Hurricane wind loads

A TC is a rapidly rotating low-pressure system that forms over warm tropical oceans, drawing energy from heat and moisture fluxes at the sea surface. Its structure comprises a calm central eye surrounded by the eyewall, where intense updrafts of moist air release latent heat, and outer spiral rainbands producing heavy rainfall (Fig. 1). The storm’s maximum wind speeds occur near the radius of maximum winds r_{\mathrm{m}}, typically located at the edge of the eyewall, while r denotes the radial distance from the hurricane center to the exposed building of interest. As the cyclone moves inland, the lack of oceanic heat and moisture causes a rapid decay in intensity.

TC originate as tropical disturbances driven by converging winds, humidity, and warm sea-surface temperatures. With increasing organization and a closed circulation, they evolve from tropical depressions to tropical storms, and finally to hurricanes when 1-min sustained winds exceed 119 km/h The Saffir–Simpson hurricane wind scale classifies these systems into five categories based on wind speed and central pressure, as summarized in Table 2 [24,25]. Typical wind field diameters range from 500 to 1000 km, although some events such as Typhoon Tip (1979; 2200 km), Hurricane Sandy (2012; 1380 km), and Hurricane Marge (1951; 1130 km) reached much larger scales [24].

2.1 Mean wind velocity

For a NTS, a good mathematical model to describe the wind speed profile in strong wind conditions is the logarithmic law, which is widely used in practical wind load codes. For greater accuracy at high-altitude above ground, Deaves and Harris [26] improved the logarithmic law to approximate the mean wind velocity throughout the atmospheric boundary layer as follows:

(1)\begin{array}{rl}\overline{U}\left(\mathit{\text{z}}\right)=& \frac{u^{\ast }}{\kappa }[\ln \left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_0}\right)+5.75\left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_{\text{g}}}\right)-1.875{\left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_{\text{g}}}\right)}^2\\ & -\frac{4}{3}{\left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_{\text{g}}}\right)}^3+\frac{1}{4}{\left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_{\text{g}}}\right)}^4],\end{array}

where \kappa \approx 0.4 is the Von Kármán constant, z is the height above ground, z₀ is the roughness length, u^{\ast } is the friction velocity and {\mathit{\text{z}}}_{\text{g}} is the gradient height; the latter two variables are given by

(2)u^{\ast }=\sqrt{\frac{{\tau }_{\text{s}}}{{\rho }_{\text{a}}}}\approx \overline{U}\left({\mathit{\text{z}}}_{\text{r}}\right)\left[\frac{\kappa }{\ln \left({\mathit{\text{z}}}_{\text{r}}/{\mathit{\text{z}}}_0\right)}\right],

(3){\mathit{\text{z}}}_{\text{g}}=\frac{u^{\ast }}{6f_{\text{c}}},

where {\tau }_{\text{s}} is the surface shear stress at the ground surface, {\rho }_{\text{a}} is the air density, {\mathit{\text{z}}}_{\text{r}} is the reference height, which is usually equal to 10 m, \overline{U}\left({\mathit{\text{z}}}_{\text{r}}\right)is the mean wind velocity at {\mathit{\text{z}}}_{\text{r}}, and f_{\text{c}} is the Coriolis parameter, which is given by

(4)f_{\text{c}}=2\mathrm{\Omega }\sin \varsigma ,

where \mathrm{\Omega }=7.27\times {10}^{-5}\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s} is the angular velocity of the Earth and \varsigma is the latitude of the site of interest in decimal degrees.

The commonly used logarithmic and power-law profiles in engineering practice, such as the one shown in Eq. (1), may be inappropriate for modeling the mean wind velocity within the atmospheric boundary layer of a hurricane, as they do not account for the super-gradient wind region, where the maximum wind occurs. The super-gradient wind region in hurricanes refers to a portion of the boundary layer where tangential winds (winds moving around the center of the storm) exceed the speeds predicted by the gradient wind model. This phenomenon typically occurs near the surface and results from complex interactions within the boundary layer, including friction and vertical mixing of air. It is commonly observed in hurricanes and plays a key role in understanding and modeling the wind dynamics and profiles in these storms.

The empirical wind profile proposed by Vickery et al. [7] is widely used and has been extensively validated for landfalling hurricanes. In this study, we employed the semi-empirical formulation of Snaiki and Wu [8], which is based on extensive field measurements over both land and ocean, accounts for the super-gradient wind region near the hurricane eye, includes the effect of wind direction through the inflow angle, and incorporates the height of maximum wind considering inertial stability and surface roughness. This profile was selected because it provides a realistic representation of the hurricane boundary layer, capturing variations in wind speed with height and radial distance that are not reflected in conventional code-based design wind speeds at 10 m. While the 10 m design wind speed for a given return period serves as a reference for structural design, it does not reflect the full vertical profile or the radial variation of the height of maximum wind, which tends to increase from the storm center toward the outer regions in real observations. Fang et al. [27] demonstrated that conventional code-suggested profiles may underestimate TC-induced wind loads along the building height and emphasized the benefits of height-dependent outer envelope and composite design wind profiles. By using the Snaiki and Wu model, we can evaluate wind-induced responses at different heights and radial distances, ensuring that high-rise buildings are accurately assessed under extreme hurricane winds.

According to Snaiki and Wu [8], the mean wind velocity at height z for marine and landfalling hurricane boundary layer can be computed by the following semi-empirical expression:

(5)\overline{U}\left(\mathit{\text{z}}\right)=\frac{u^{\ast }}{\kappa }\left[\ln \left(\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_0}\right)+{\eta }_0\sin \left(\frac{\mathit{\text{z}}}{\delta }\right)\exp \left(-\frac{\mathit{\text{z}}}{\delta }\right)\right],

where the constant {\eta }_0 is determined by the condition {{\displaystyle \frac{\mathrm{\partial }\overline{U}}{\mathrm{\partial }\mathit{\text{z}}}}|}_{\mathit{\text{z}}=\delta }, which results in {\eta }_0\approx 9.026, and \delta is the height of the maximum wind in the super-gradient wind region, which is given by

(6)\delta =\exp \left[a_1\text{}\ln \left(I\right)+a_2\text{}\ln \left(R_{\text{o}}\right)+a_3\right],

where I is the inertial stability that represents the ability of the TC to resist changes or disturbances in its rotational balance, R_o is the modified surface Rossby number that describes the relative importance of inertial forces to Coriolis forces for the rotating eyewall; and \left(a_1,a_2,a_3\right) are empirical constants. For marine hurricane boundary layer, a_1=-0.4807, a_2=-0.05, and a_3= 4.0221, whereas for landfalling hurricanes, a_1=-0.2452, a_2=-0.05, and a_3=5.6149. Regarding the inertial stability and modified surface Rossby number, both parameters can be computed by the following expressions:

(7)I=\sqrt{\left(f_{\text{c}}+\frac{2{\nu }_{\text{g}}}{r}\right)\left(f_{\text{c}}+\frac{{\nu }_{\text{g}}}{r}+\frac{\text{δ}{\nu }_{\text{g}}}{\text{δ}r}\right)},

(8)R_{\text{o}}=\frac{{\nu }_{\text{g}}}{{\mathit{\text{z}}}_{0}I},

where r is the radial position of the building with respect to the center of the storm (Fig. 1) and {\nu }_{\mathrm{g}} is the gradient wind speed, which represents the theoretical wind speed where the flow is influenced by the balance of pressure gradient force, Coriolis force and centrifugal force. Accordingly, the gradient wind speed and its derivative with respect to r are given by Ref. [28]:

(9)v_{\text{g}}=-\frac{\left|f_{\text{c}}\right|\text{}r}{2}+\sqrt{\frac{{\left(f_{\text{c}}\text{}r\right)}^2}{4}+\frac{r}{{\rho }_{\text{a}}}\left|\frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}\right|},

(10){\displaystyle \frac{\text{δ}{\nu }_g}{\text{δ}r}}=-{\displaystyle \frac{\left|f_{\text{c}}\right|}{2}}+{\displaystyle \frac{{\displaystyle \frac{f_{\text{c}}{}^{2}r}{4}}+{\displaystyle \frac{1}{2{\rho }_{\text{a}}}}\left|{\displaystyle \frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}}\right|+{\displaystyle \frac{r}{2{\rho }_{\text{a}}}}{\displaystyle \frac{\text{d}}{\text{d}r}}\left|{\displaystyle \frac{\text{δ}{p}_{\mathrm{s}}}{\text{δ}r}}\right|}{\sqrt{{\displaystyle \frac{{\left(f_{\text{c}}r\right)}^2}{4}}+{\displaystyle \frac{r}{{\rho }_{\text{a}}}}\left|{\displaystyle \frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}}\right|}}},

where p_{\mathrm{s}} is the surface pressure of the hurricane at a radial distance r, {\displaystyle \frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}}is the radial pressure gradient and {\displaystyle \frac{\text{d}}{\text{d}r}}\left|{\displaystyle \frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}}\right| is its derivative. These parameters can be computed by the following expressions [28,29]:

(11)p_{\text{s}}=p_{\text{cs}}+\mathrm{\Delta }{p}_{\text{s}}\exp \left[-{\left(\frac{r_{\text{m}}}{r}\right)}^b\right],

(12)\frac{\text{δ}{p}_{\mathrm{s}}}{\text{δ}r}=\frac{\mathrm{\Delta }{p}_{\text{s}}b}{r}{\left(\frac{r_{\text{m}}}{r}\right)}^b\exp \left[-{\left(\frac{r_{\text{m}}}{r}\right)}^b\right],

(13)\frac{\text{d}}{\text{d}r}\left|\frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}\right|=\frac{1}{\left|{\displaystyle \frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}}\right|}{\left(\frac{\text{δ}{p}_{\text{s}}}{\text{δ}r}\right)}^2\frac{1}{r}\left[b{\left(\frac{r_{\text{m}}}{r}\right)}^b-b-1\right],

where \mathrm{\Delta }{p}_{\text{s}} is the central pressure deficit, that is, the pressure difference between the center of the storm and the outermost closed isobar, b is the Holland’s scale parameter from the statistical models performed by Vickery and Wadhera [30], and r_{\text{m}} is the radius of maximum wind near the edge of the eyewall. Accordingly, these parameters are given by

(14)\mathrm{\Delta }{p}_{\text{s}}=p_{\text{es}}-p_{\text{cs}},

(15)b=1.833-0.326\sqrt{f_{\text{c}}{r}_{\text{m}}},

(16)r_{\text{m}}=\exp \left(3.015-6.291\times {10}^{-5}\mathrm{\Delta }{p}_{\text{s}}{}^2+0.0337\varsigma \right).

Based on Eq. (14), p_{\text{cs}} is the central pressure of the storm, which can be obtained from Table 2; and p_{\text{es}}=1013\text{hPa} is the external pressure of the storm, that is, the atmospheric pressure outside the hurricane over the ocean. It is important to mention that Eqs. (15) and (16) are relationships fitted to empirical data [30], which require a particular use of units: Eq. (15) requires that r_m be input in meters, while Eq. (16) is used to obtain r_m in kilometers but requires \mathrm{\Delta }{p}_{\text{s}} to be input in hPa and \zeta in decimal degrees.

As a building gets closer to the eye wall of a landfalling hurricane, the mean wind velocity increases. To illustrate this phenomenon, Fig. 2 compares the mean wind velocity for the city of Cancun, México, over terrain representative of coastal areas exposed to the open sea, such as beachfront hotels. For this comparison, Eq. (1) was used to represent the NTS profile, and Eq. (5) for hurricane winds, varying the radial distance from 10 to 200 km (Fig. 1). Although the Mexican wind design code [6] actually employs a power-law profile, for illustrative purposes we represent the code-based profile using a corrected logarithmic law (Eq. (1)), which is typically used in other international wind design codes. This facilitates comparison with the Snaiki and Wu [8] profile, which accounts for the supergradient region near the eyewall. The supergradient profile produces higher wind velocities at elevations typical of tall buildings, whereas the power-law profile of the Mexican code generally results in slightly lower velocities. Based on the Mexican wind design code, three return periods were considered: 10-year (category 1 hurricane, \overline{U}(10\mathrm{m})=107.82\mathrm{k}\mathrm{m}/\mathrm{h}, \mathrm{\Delta }{p}_{\text{s}}=28.44\text{hPa}); 50-year (category 3 hurricane, \overline{U}(10\mathrm{m})=164.16\mathrm{k}\mathrm{m}/\mathrm{h}, \mathrm{\Delta }{p}_{\text{s}}=59.93\text{hPa}); and 200-year (category 4 hurricane, \overline{U}(10\mathrm{m})=201.85\mathrm{k}\mathrm{m}/\mathrm{h}, \mathrm{\Delta }{p}_{\text{s}}=84.43\text{hPa}). Figure 2 demonstrates that while code-based profiles consider return periods associated with hurricanes, standard formulations do not fully capture the supergradient wind region, and in some cases can underestimate tail risks, that is, low-probability, high-impact events. Additionally, as the building gets closer to the eye wall, the supergradient wind region becomes more evident, with the maximum wind speed occurring at lower heights.

2.2 Power spectral density function

The power spectral density (PSD) function describes how wind energy is distributed across different frequencies. In other words, it is a mathematical representation of wind speed fluctuations (turbulence) in the frequency domain. For NTS, there are many mathematical models used to represent the energy across various wind frequencies. For example, the spectra from Solari [31], Kaimal et al. [10], and Davenport [32] do not adequately reproduce wind power at low frequencies, while the Von Kármán [11] spectrum underestimates wind power at high frequencies. On the other hand, the PSD function contained in the engineering science data unit (ESDU) [33] adequately reproduces the wind power for NTS at low frequencies, high frequencies, and in the inertial sub-range [34,35].

For TC, the PSD function for the along-wind velocity can be computed by the following expression proposed by Yu et al. [13]:

(17)S_{\text{u}}\left(\mathit{\text{z}},f\right)=\frac{{\sigma }_{\text{u}}{}^2}{f\mathrm{\Phi }}\left[\frac{p_1{\hat{f}}^2+p_2\hat{f}+p_3}{{\hat{f}}^3+q_1{\hat{f}}^2+q_2\hat{f}+q_3}\right],

where

(18)\hat{f}=\frac{f\mathit{\text{z}}}{\overline{U}\left(\mathit{\text{z}}\right)},

(19)\mathrm{\Phi }={\left[2.72-0.25{\log }_{10}\left({\mathit{\text{z}}}_0\right)\right]}^2,

(20){\sigma }_{\text{u}}=u^{\ast }\sqrt{\mathrm{\Phi }}=u^{\ast }\left[2.72-0.25{\log }_{10}\left({\mathit{\text{z}}}_0\right)\right].

In Eq. (17), f is the frequency, \hat{f} is the normalized frequency, \sqrt{\mathit{\Phi }} is the turbulence ratio and {\sigma }_{\text{u}} is the standard deviation of the along-wind velocity. According to Li et al. [36], Eqs. (19) and (20) are empirical expressions associated with hurricane winds. Concerning the empirical constants \left(p_1,p_2,p_3\right) and \left(q_1,q_2,q_3\right), these were proposed by Yu et al. [10] based on data of six hurricanes along the south-eastern coast of the United States: Bonnie (1998), Gordon (2000), Isidore (2002), Lili (2002), Ivan (2004) and Ike (2008). For marine hurricane boundary layer, p_1=-0.00598, p_2=0.1544, p_3=0.00001055, q_1=0.4458, q_2=0.06486, and q_3=0.00009754; whereas for landfalling hurricanes, p_1=-0.9999, p_2=3.112, p_3=0.0001159, q_1=18.64, q_2=1.188, and q_3=0.00335.

These empirical constants, originally calibrated using wind measurements at 10 m height by Yu et al. [13], are applied uniformly along the height of the structure in this study. Detailed PSD measurements at heights above 10 m for landfalling hurricanes are currently scarce, particularly over land. Therefore, this approach, which follows common practice in Refs. [37,38], combines surface-based PSD parameters with the mean wind velocity profile variation with height (Eq. (5)). Although this introduces some uncertainty, the main objective, to capture the greater energy content at low frequencies characteristic of hurricane wind spectra compared to non-hurricane spectra, is preserved, since the mean wind speed variation with height is explicitly accounted for in the simulations.

Figure 3 compares the PSD functions at different heights for NTS [7,8,27–29] and hurricanes (Eq. (17)) in Cancun for a 200-year return period (Fig. 2(c)). Turbulent energy is higher at low frequencies for hurricane winds, a phenomenon that becomes more pronounced as height above the ground increases. The radial distance of the hurricane, r, does not significantly affect the energy content across the frequency range. For the Von Kármán PSD, the integral length scale L_{\text{u}}\left(\mathit{\text{z}}\right) is evaluated at the height z of interest, following empirical relationships (e.g., ESDU [33]), allowing the spectrum to reflect local turbulence characteristics and be consistently compared with other PSD functions at different heights.

Recent advances in wind spectral modeling have highlighted the importance of accurately capturing both the variability and stochastic nature of turbulent energy in TC winds. Field-validated analytical models, such as that proposed by Li et al. [15], provide insights into the spectral characteristics of typhoon winds over the sea surface, including non-classical energy cascades in the eyewall region. Building upon these foundations, Liu et al. [39] developed stochastic PSD models for typhoon and non-typhoon winds based on long-term structural health monitoring data, revealing statistical correlations among spectral parameters and their dependence on mean wind speed, and enabling the random simulation of wind spectra under varying conditions. Additionally, Fang et al. [40] introduced a probabilistic gust factor model that accounts for the non-Gaussian attributes, time-varying mean wind speed, and dispersion of gusts, further allowing Monte Carlo-based estimation of site-specific wind hazards. Together, these studies reinforce the rationale for employing semi-empirical or data-driven PSD formulations in this work, ensuring that the probabilistic simulations performed capture the enhanced low-frequency energy, stochastic variability, and gust characteristics of hurricane winds, as reflected in the response spectra presented in Fig. 3.

2.3 Stochastic simulation

Based on Fig. 4, the hurricane will generate a drag force on the building as it impacts its façade. This force results from the combined effect of wind pressure on the windward side and suction on the leeward side. Accordingly, the fluctuating drag force in time domain is given by

(21)F_{\text{D}}\left(\mathit{\text{z}},t\right)=\frac{1}{2}{\rho }_{\text{a}}A{C}_{\text{D}}U{\left(\mathit{\text{z}},t\right)}^2,

where A is the projected area, that is, the effective area of the building facing the wind, {\rho }_{\text{a}} is the air density, t is the time, C_{\text{D}} is the drag coefficient and U\left(\mathit{\text{z}},t\right) is the along-wind velocity, which is given by the following expression for a stationary wind assumption:

(22)U\left(\mathit{\text{z}},t\right)=\overline{U}\left(\mathit{\text{z}}\right)+u\left(\mathit{\text{z}},t\right),

where \overline{U}\left(\mathit{\text{z}}\right) is the mean wind velocity within the hurricane boundary layer for a 10-min averaging time, which can be computed by using Eq. (5), and u\left(\mathit{\text{z}},t\right) is the turbulence or fluctuating component of the along-wind velocity, which is related to the PSD function of Eq. (17).

According to Dyrbye and Hansen [41], the PSD function of along-wind velocity shown in Eq. (17) can be transformed into wind drag force PSD function using the aerodynamic admittance function, which adjusts the magnitude of fluctuating wind loads applied to the building based on the influence of the structure’s surfaces on wind pressures. Therefore, the PSD function of the wind drag force is given by

(23)S_{\text{FD}}\left(\mathit{\text{z}},f\right)={\left[{\rho }_{\text{a}}{C}_{\text{D}}A\overline{U}\left(\mathit{\text{z}}\right)\right]}^2\chi {\left(\mathit{\text{z}},f\right)}^2{S}_{\text{u}}\left(\mathit{\text{z}},f\right),

where \chi (\mathit{\text{z}},f{)}^2 is the aerodynamic admittance function. There are several empirical formulas based on experimental evidence for the calculation of \chi {\left(\mathit{\text{z}},f\right)}^2, which are included in wind design codes. However, a widely used expression in the scientific literature is the one proposed by Vickery [42]:

(24)\chi {\left(\mathit{\text{z}},f\right)}^2={\left[1+{\left(\frac{2f\sqrt{A}}{\overline{U}\left(\mathit{\text{z}}\right)}\right)}^{4/3}\right]}^{-2}.

Based on the improved spectral representation method [43,44], the drag force can be simulated in time domain as a set of N homogeneous Gaussian multidimensional stochastic processes:

(25)F_{\text{D}}\left({\mathit{\text{z}}}_j,t\right)=\sum _{k=1}^N\sum _{n=1}^m\left|H_{jk}\left(f_n\right)\right|\sqrt{2\mathrm{\Delta }f}\cos \left[2\text{π}{f}_{n}t+{\varphi }_{kn}\right],

where j=1,2,3,\dots ,N, N is the total number of stories in the buildings, mis the total number of frequencies that composes the PSD function, \mathrm{\Delta }f is the sampling frequency of the PSD function, {\varphi }_{kn} represents independent random phase angles uniformly distributed over the interval \left[0,2\text{π}\right], and H_{jk}\left(f_n\right) represents the value in row j and column k of the lower triangular matrix resulting from applying the Cholesky decomposition to the cross-spectral density matrix.

For the purpose of applying Eq. (25), the cross-spectral density matrix represents the relationship between the spectral components of drag forces at different heights. In this context, a coherence function is needed to correlate the PSD function from Eq. (23) along the height of the building. The empirical coherence function proposed by Davenport [45] is commonly used to simulate NTS, and has also been applied for hurricanes [37,38]. The Krenk coherence function [46] provides a more accurate representation of along-wind turbulence for structures, correcting limitations of Davenport’s function such as overestimation of correlation at low frequencies and large separations. It is often used to evaluate wind-induced vibrations on buildings and other structures [41], and is given by

(26)Co{h}_{jk}\left(f\right)=\left(1-\frac{f^{\ast }{C}_{\text{z}}{r}_{\text{z}}}{\overline{U}\left({\mathit{\text{z}}}_j\right)+\overline{U}\left({\mathit{\text{z}}}_k\right)}\right)\exp \left(-\frac{2f^{\ast }{C}_{\text{z}}{r}_{\text{z}}}{\overline{U}\left({\mathit{\text{z}}}_j\right)+\overline{U}\left({\mathit{\text{z}}}_k\right)}\right),

where C_{\mathrm{z}}=5 is the vertical decay constant according to Hansen and Krenk [47], r_{\text{z}} is the vertical separation between point j and point k, and f^{\ast } is the modified frequency given by

(27)f^{\ast }=\sqrt{f^2+{\left(\frac{\overline{U}\left({\mathit{\text{z}}}_j\right)+\overline{U}\left({\mathit{\text{z}}}_k\right)}{2\text{π}\left[L_{\text{u}}\left({\mathit{\text{z}}}_j\right)+L_{\text{u}}\left({\mathit{\text{z}}}_k\right)\right]}\right)}^2},

where L_{\text{u}}\left(\mathit{\text{z}}\right) is the integral length scale characterizing the correlation of along-wind turbulence across the building height. It can be computed using empirical formulations contained in wind design codes, such as ESDU [33]. Huergo [34], Huergo et al. [35], and Huergo et al. [48] provide a detailed description of how to apply Eq. (25) using Krenk’s coherence function for NTS. In this study, L_{\text{u}}\left(\mathit{\text{z}}\right) was computed following the ESDU [33] formulation, which is based on strong wind measurements in a neutral atmosphere. Although these formulations were originally developed for non-tropical strong winds, they provide a reasonable approximation for hurricane winds in the absence of comprehensive field data. Accordingly, the integral length scale is given by

(28)L_{\text{u}}\left(\mathit{\text{z}}\right)=\frac{A_k^{3/2}{\left[{\sigma }_{\text{u}}\left(\mathit{\text{z}}\right)/u^{\ast }\right]}^3\mathit{\text{z}}}{2.5K_{\text{z}}^{3/2}{\left(1-\mathit{\text{z}}/{\mathit{\text{z}}}_{\text{g}}\right)}^2\left(1+5.75\mathit{\text{z}}/{\mathit{\text{z}}}_{\text{g}}\right)},

where

(29)A_k=0.115{\left[1+0.3151{\left(1-\frac{\mathit{\text{z}}}{{\mathit{\text{z}}}_{\text{g}}}\right)}^6\right]}^{2/3},

(30)K_{\text{z}}=0.19-\left(0.19-K_{\text{o}}\right)e^{-B_k{\left(\mathit{\text{z}}/{\mathit{\text{z}}}_{\text{g}}\right)}^{N_k}},

(31)K_{\text{o}}=\frac{0.39}{R_i^{0.11}},

(32)R_i=\frac{u^{\ast }}{f_{\text{c}}{\mathit{\text{z}}}_0},

(33)B_k=24R_i{}^{0.155},

(34)N_k=1.24R_i{}^{0.008}.

In the calculation of the integral length scale L_{\text{u}}\left(\mathit{\text{z}}\right), the parameters u^{\ast }, {\mathit{\text{z}}}_{\text{g}}, and f_{\text{c}} are the wind friction velocity, the gradient height, and the Coriolis parameter, defined in Eqs. (2)–(4), respectively. Additionally, R_i is the dimensionless Richardson number, which expresses the ratio between the potential and kinetic energy of the air.

For non-tropical winds, the standard deviation of the along-wind velocity, {\sigma }_{\mathrm{u}}, is given by the following expression proposed in ESDU [33]:

(35){\sigma }_{\text{u}}\left(\mathit{\text{z}}\right)=\frac{7.5\eta u^{\ast }{\left[0.538+0.09\ln \left(\mathit{\text{z}}/{\mathit{\text{z}}}_0\right)\right]}^p}{1+0.156\ln \left(R_i\right)},

where

(36)\eta =1-\frac{6f_{\text{c}}\mathit{\text{z}}}{u^{\ast }},

(37)p={\eta }^{16}.

It is important to note that across-wind effects, particularly vortex shedding, can dominate the response of tall buildings, especially those with aspect ratios, H/\sqrt{BD}, over 3. Previous studies [34,35,48,49] have shown that across-wind responses can exceed along-wind responses when the building’s natural frequencies approach the vortex shedding frequency. However, those analyses were based on PSDs derived from wind-tunnel experiments simulating a standard atmospheric boundary layer, not hurricane winds. In this study, we focus on along-wind loads under hurricane conditions, which generally governs the response for the high-rise buildings considered. This limitation is acknowledged here to clarify the scope of the analysis, while recognizing that across-wind effects remain important for future investigations.

3 Coupled shear-flexural model with nonlinear fluid viscous dampers

A N-story tall building, featuring any type of lateral resisting system and equipped with FVDs, can be idealized as shown in Fig. 5, where f_{\text{a},j} and c_{\mathrm{F}\mathrm{V}\mathrm{D}j} are the amplification factor and damping coefficient of the jth FVD, respectively. The model characterizes the overall flexural stiffness of the building through a flexural beam idealization (shear wall) and the overall effective shear stiffness via a multi-degree-of-freedom representation of a shear building. By coupling these stiffness components in parallel using axially rigid links, the model accurately captures the lateral deformation behavior associated with various lateral resisting systems, including moment-resisting frames, braced frames, dual systems (shear wall-frame buildings), tube-in-tube configurations, among others.

The equation of motion related to the mechanical system in Fig. 5 is given by

(38)\begin{array}{rl}& \frac{\overline{m}}{EI}\frac{{\mathrm{\partial }}^{2}u\left(\overline{\mathit{\text{z}}},t\right)}{\mathrm{\partial }{t}^2}+\frac{\overline{c}}{EI}\frac{\mathrm{\partial }u\left(\overline{\mathit{\text{z}}},t\right)}{\mathrm{\partial }t}+\frac{1}{H^4}\frac{{\mathrm{\partial }}^{4}u\left(\overline{\mathit{\text{z}}},t\right)}{\mathrm{\partial }{\overline{\mathit{\text{z}}}}^4}\\ -& \frac{{\alpha }^2}{H^4}\frac{{\mathrm{\partial }}^{2}u\left(\overline{\mathit{\text{z}}},t\right)}{\mathrm{\partial }{\overline{\mathit{\text{z}}}}^2}=\frac{\overline{f}\left(\overline{\mathit{\text{z}}},t\right)}{EI},\end{array}

where \overline{m} is the mass per unit length, EI is the global flexural stiffness, \overline{c} is the total damping coefficient per unit length (structural + aerodynamic + FVDs), H is the total height, u\left(\overline{\mathit{\text{z}}},t\right) is the dynamic displacement at the nondimensional height \overline{\mathit{\text{z}}}(ranging from zero at the base to one at the rooftop) and at time t, \overline{f}\left(\overline{\mathit{\text{z}}},t\right) is the along-wind dynamic force per unit length at \overline{\mathit{\text{z}}} and at time t, and \alpha is a nondimensional lateral stiffness ratio given by the following expression:

(39)\alpha =H\sqrt{\frac{G\text{}{A}_{\text{v}}}{E\text{}I}},

where G\text{}{A}_{\text{v}} is the overall effective shear stiffness of the building.

Assuming a linear elastic behavior, u\left(\overline{\mathit{\text{z}}},t\right) can be expressed in terms of a modal decomposition:

(40)u\left(\overline{\mathit{\text{z}}},t\right)=\sum _{i=1}^{\mathrm{\infty }}{\varphi }_i\left(\overline{\mathit{\text{z}}}\right)q_i\left(t\right),

where {\varphi }_{\xi }(\overline{\mathit{\text{z}}}) represents the ith mode shape and q_t(t) is the generalized displacement for the ith mode of vibration. According to Miranda & Taghavi [50], {\varphi }_s(\overline{\mathit{\text{z}}}) is given by

(41)\begin{array}{rl}{\varphi }_i\left(\overline{\mathit{\text{z}}}\right)=& \sin \left({\gamma }_i\overline{\mathit{\text{z}}}\right)-\frac{{\gamma }_i}{\sqrt{{\alpha }^2+{\gamma }_i^2}}\sin h\left(\overline{\mathit{\text{z}}}\sqrt{{\alpha }^2+{\gamma }_i{}^2}\right)\\ & -{\eta }_i\cos \left({\gamma }_i\overline{\mathit{\text{z}}}\right)+{\eta }_i\cos h\left(\overline{\mathit{\text{z}}}\sqrt{{\alpha }^2+{\gamma }_i^2}\right),\end{array}

where

(42){\eta }_i=\frac{{\gamma }_{\mathrm{i}}^2\sin {\gamma }_i+{\gamma }_i\sqrt{{\alpha }^2+{\gamma }_{\mathrm{i}}{}^2}\sin h\left(\sqrt{{\alpha }^2+{\gamma }_i^2}\right)}{{\gamma }_i^2\cos {\gamma }_i+\left({\alpha }^2+{\gamma }_i^2\right)\cos h\left(\sqrt{{\alpha }^2+{\gamma }_i{}^2}\right)},

and {\gamma }_i is the first root of the following characteristic equation:

(43)\begin{array}{rl}& 2+\left(2+\frac{{\alpha }^4}{{\gamma }_i^2\left({\alpha }^2+{\gamma }_i^2\right)}\right)\cos {\gamma }_i\cos h\left(\sqrt{{\alpha }^2+{\gamma }_i^2}\right)\\ +& \left(\frac{{\alpha }^2}{{\gamma }_i\sqrt{{\alpha }^2+{\gamma }_i^2}}\right)\sin {\gamma }_i\sin h\left(\sqrt{{\alpha }^2+{\gamma }_{\mathrm{i}}^2}\right)=0,\end{array}

On the other hand, the generalized displacement of the ith mode of vibration must be adapted to account the energy dissipation associated with FVDs, as follows

(44)\frac{\text{d}{\text{}}^2{q}_i\left(t\right)}{\text{d}t}+2\text{}{\xi }_i\text{}{w}_i\text{}\frac{\text{d}{q}_i\left(t\right)}{\text{d}t}+w_i^2\text{}{q}_i\left(t\right)=P_i\left(t\right),

where

(45){\xi }_{\text{}i}={\xi }_{\text{}s,j}+{\xi }_{\text{FVD}}^{\text{EE}},

(46)w_i=\sqrt{\frac{EI}{\overline{m}{H}^4}{\gamma }_i{}^2\left({\gamma }_i{}^2+{\alpha }_i{}^2\right)},

(47)P_i\left(t\right)=\frac{{\displaystyle {\int }_0^1\overline{f}\left(\overline{\mathit{\text{z}}},t\right){\varphi }_i\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}}{{\tilde{m}}_{\text{g},i}}=\frac{{\displaystyle {\int }_0^1\overline{f}\left(\overline{\mathit{\text{z}}}\text{,}t\right){\varphi }_i\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}}{\overline{m}{\displaystyle {\int }_0^1{\varphi }_i^2\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}},

where {\xi }_1 is the total damping ratio of the ith mode of vibration, {\xi }_{\mathrm{s},i} is the structural damping ratio (including aerodynamic damping) of the ith mode of vibration, {\xi }_{\text{FVD}}^{\text{EE}} is the equivalent energy damping ratio, that is, the added damping ratio provided by a set of linear FVDs that dissipate the same amount of energy as the nonlinear FVDs, and {\tilde{m}}_{\text{g},i} is the generalized mass for the ith mode of vibration.

The equation of motion of the generalized single-degree-of-freedom system can be solved by the state-space approach [51]. Accordingly, Eq. (44) can be rewritten as

(48)\left\{\begin{array}{c}{\displaystyle \frac{\text{d}{q}_i\left(t\right)}{\text{d}t}}\\ {\displaystyle \frac{{\text{d}}^2{q}_i\left(t\right)}{\text{d}t}}\end{array}\right\}=\left[\begin{array}{cc}0& 1\\ -w_i^2& -2{\xi }_i{w}_i\end{array}\right]\left\{\begin{array}{c}{q}_i\left(t\right)\\ {\displaystyle \frac{\text{d}{q}_i\left(t\right)}{\text{d}t}}\end{array}\right\}+\left\{\begin{array}{c}0\\ F_i^{\text{g}}\left(t\right)\end{array}\right\},

which leads to computing the dynamic response at time t_{v+1} as follows:

(49)\left\{\begin{array}{c}{q}_i\left(t_{v+1}\right)\\ {\displaystyle \frac{\text{d}{q}_i\left(t_{v+1}\right)}{\text{d}t}}\end{array}\right\}={\text{e}}^{\left[\begin{array}{cc}0& 1\\ -w_i^2& -2{\xi }_i{w}_i\end{array}\right]\mathrm{\Delta }t}\left\{\begin{array}{c}{q}_i\left(t_v\right)\\ {\displaystyle \frac{\text{d}{q}_i\left(t_v\right)}{\text{d}t}}\end{array}\right\}+{\text{e}}^{\left[\begin{array}{cc}0& 1\\ -w_i^2& -2{\xi }_i{w}_i\end{array}\right]\mathrm{\Delta }t}\left\{\begin{array}{c}0\\ P_i\left(t\right)\end{array}\right\},

where \mathrm{\Delta }t is the sampling interval of the along-wind loads.

3.1 Energy-dissipation equivalent linearization method

The manufactured FVDs exert a nonlinear force proportional to the relative velocity between the ends of the device [52], which can be obtained for the jth story using the following expression:

(50)f_{\text{FVD},j}\left(\mathrm{\Delta }{\dot{w}}_j\right)=c_{\text{FVD},j}{\left|\mathrm{\Delta }{\dot{w}}_j\left(t\right)\right|}^{{\beta }_j}\mathrm{s}\mathrm{g}\mathrm{n}\left[\mathrm{\Delta }{\dot{w}}_j\left(t\right)\right],

where

(51)\mathrm{\Delta }{\dot{w}}_j\left(t\right)=f_{\text{a},j}\text{}\mathrm{\Delta }{\dot{u}}_j\left(t\right),

(52)\mathrm{\Delta }{\dot{u}}_j\left(t\right)={\dot{u}}_j\left(t\right)-{\dot{u}}_{j-1}\left(t\right).

In Eqs. (50)–(52), j=1,2,3,\dots ,N, c_{\text{FVD,}j} denotes the damping coefficient of the jth FVD, {\beta }_j is the velocity exponent of the jth FVD, f_{\mathrm{a},j} is the amplification factor of the jth FVD based on the configuration of the bracing system (Fig. 5), {\dot{u}}_j\left(t\right) is the dynamic velocity of the building at jth story, \mathrm{\Delta }{\dot{w}}_j\left(t\right) is the relative velocity at the ends of the jth FVD, and \mathrm{\Delta }{\dot{u}}_j\left(t\right) is the relative velocity of the building at jth story. According to Eq. (40), a FVD shows nonlinear behavior when {\beta }_j\ne 1. Typically, {\beta }_j ranges between 0.1 and 1.0 for most practical applications, although values below 0.1 are also possible; for example, {\beta }_j=0 corresponds to an idealized friction damper. Additionally, values greater than 1.0 have been explored in research contexts, particularly in connection with dampers using shear-thickening fluids [53–56], which exhibit a superlinear force–velocity relationship. Although such devices are less common in structural engineering practice, acknowledging their existence provides a broader context for the interpretation of Eq. (40).

A lower exponent value {\beta }_j results in greater energy dissipation by FVDs, as noted by Ramírez et al. [57]. Lin and Chopra [58] introduced the concept of “equivalent energy” FVDs, considering harmonic motion with a displacement amplitude \mathrm{\Delta }w\left(t\right)=\mathrm{\Delta }w\sin {\left(w_{1}t\right)}_{\text{max}}, where w_1 is the fundamental angular frequency and \mathrm{\Delta }{w}_{\text{max}} represents the peak displacement. By equating the energy dissipation of both linear ({\beta }_j=1) and nonlinear ({\beta }_j\ne 1) FVDs, the hysteresis loops presented in Fig. 6 are obtained. Here f_{\text{FVD}}/f_{\text{FVD,max}}^{\text{L}} denotes the ratio of the nonlinear force of the FVD to the maximum amplitude of the corresponding equivalent linear force.

Figure 6(a) illustrates that as \beta decreases, the hysteresis cycle transitions from an elliptical to a rectangular shape, indicating increased energy dissipation. Similarly, Fig. 6(b) shows that for small values of \beta, the damper force changes more abruptly near zero velocity, which enhances the device’s vibration control performance. However, these abrupt changes in the force of a nonlinear FVD are characteristic of stiff differential equation systems, making the application of traditional numerical methods more challenging. Therefore, it is preferable to idealize the mechanical model in Fig. 5 as a system with equivalent linear FVDs that dissipate the same amount of energy as nonlinear FVDs.

The goal is to replace the original nonlinear system with a linear system, following the equivalent energy dissipation criterion of the FVDs. Therefore, the value of {\xi }_{\text{FVD}}^{\text{EE}} in Eq. (45) is provided by Ref. [50]:

(53){\xi }_{\text{FVD}}^{\text{EE}}=\frac{{\displaystyle \sum _{j=1}^{n_{\text{FVD}}}{\lambda }_j{c}_{\text{FVD,}j}{w}_{\text{l}}^{{\beta }_j-2}{\left|{\text{u}}_{\text{roof,l}}\right|}^{{\beta }_j-1}{\left|{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_j\right)-{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_{j-1}\right)\right|}^{{\beta }_j+1}{f}_{\text{a},j}{}^{{\beta }_j+1}}}{2\text{π}\overline{m}{\displaystyle {\int }_0^1{\varphi }_{\text{l}}^2\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}},

where

(54){\lambda }_j=\frac{2^{{\beta }_j+2}{\mathrm{\Gamma }}^2\left(1+{\displaystyle \frac{{\beta }_j}{2}}\right)}{\mathrm{\Gamma }\left(2+{\beta }_j\right)}.

where n_{\text{FVD}} is the total number of FVDs, w_1 is the fundamental angular frequency of the building, {\varphi }_1\left(\overline{\mathit{\text{z}}}\right) is the 1st mode shape, \mathit{\Gamma }(\cdot ) is the gamma function and u_{\text{roof},1} is the peak displacement at the rooftop of the building when only the first mode of vibration is considered, meaning it is the peak value of {\varphi }_1\left({\overline{\mathit{\text{z}}}}_{\text{N}}\right)q_1\left(t\right).

For practical construction purposes, it is convenient to assume that the damping coefficients of the nonlinear FVDs are equal at each story. According to Lin et al. [59], these are given by

(55)c_{\text{FVD,}j}=\frac{{\left(2\text{π}\right)}^{3-{\beta }_j}{\left|u_{\text{roof,l}}\right|}^{1-{\beta }_j}{\xi }_{\text{FVD}}^{\text{EE}}\overline{m}{\displaystyle {\int }_0^1{\varphi }_{\text{l}}^2\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}}{T_{\text{l}}^{2-{\beta }_j}{\lambda }_j{\displaystyle \sum _{j=1}^{n_{\text{FVD}}}{f}_{\text{a},j}^{1+{\beta }_j}{\left|{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_j\right)-{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_{j-1}\right)\right|}^{1+{\beta }_j}}},

where T_1 is the fundamental period of vibration of the building.

From Eq. (53), it can be seen that u_{\text{roof,1}} must be known beforehand to determine {\xi }_{\text{FVD}}^{\text{EE}}. However, {\xi }_{\text{FVD}}^{\text{EE}} is also required to compute u_{\text{roof,1}}. Consequently, finding {\xi }_{\text{FVD}}^{\text{EE}} requires an iterative approach, beginning with reasonable initial estimates (such as using u_{\text{roof,1}} for an initial assumption of {\xi }_{\text{FVD}}^{\text{EE}}=0) and repeating the process until convergence is reached. According to Hwang et al. [60], a target design value for the added damping ratio {\xi }_{\text{FVD}}^{\text{EE}} typically ranges between 5% and 15%, providing reasonable and cost-effective values while having a negligible impact on the fundamental period of vibration of the structure.

Linearization methods can serve as a useful tool, particularly during the preliminary design phase, and have been applied not only to structures under seismic loading but also under wind loading. Notably, Zeng et al. [61] demonstrated that the energy-dissipation equivalent linearization method can yield high accuracy when applied to buildings subjected to wind loads compared to traditional approaches commonly used in practice.

In the case of linear FVDs ({\beta }_j=1), the use of equivalent linearization methods is not required [62,63], as the damping coefficients can be directly determined using the formula provided in FEMA 356 [64], which, when adapted to the continuous model shown in Fig. 5, is given by

(56)c_{\text{FVD,}j}=\frac{4\text{π}{\xi }_{\text{FVD}}^{}\overline{m}{\displaystyle {\int }_0^1{\varphi }_{\text{l}}^2\left(\overline{\mathit{\text{z}}}\right)\text{d}\overline{\mathit{\text{z}}}}}{T{\displaystyle \sum _{j=1}^{n_{\text{FVD}}}{f}_{\text{a},j}^2{\left|{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_j\right)-{\varphi }_{\text{l}}\left({\overline{\mathit{\text{z}}}}_{j-1}\right)\right|}^2}},

where S_{\mathrm{F}\mathrm{V}\mathrm{D}} is the added damping ratio of the FVDs, which does not require the iterative process of Eq. (53).

For further validation, a brief comparison between the continuous model of Fig. 5, where the energy-dissipation equivalent linearization method is applied, and a six-degree-of-freedom lumped-mass shear building model analyzed via nonlinear step-by-step dynamic analysis is presented in the Appendix in Electronic Supplementary materials.

4 Probabilistic assessment

Based on wind codes, the criterion for determining when a building becomes susceptible to the turbulent component of wind varies. However, it is generally considered that buildings with a fundamental vibration period (T_1) greater than 1 s are more likely to be affected [6].

For this purpose, a 258 m tall building was chosen as the study model, with a structural system consisting of beams, columns, bracing, and shear walls at the core. It was considered that the building consists of 60 stories, each with a height of 4.3 m. For further details on the structural elements of the benchmark building, the reader is referred to the work of Lai et al. [65].

The 258 m benchmark building was selected to capture the full vertical distribution of hurricane winds, including low-frequency energy and turbulence effects. This height is representative of tall high-rises in hurricane-prone areas and within the practical range of existing skyscrapers, as very few exceed 300 m. For instance, in Fig. 2(c), corresponding to a Category 4 hurricane with a reference mean wind speed of 200 km/h, the maximum wind speed occurs at a height of approximately 441 m. For the upper limit of the fragility curves (Category 5 hurricane, 300 km/h), the building reaches close to the height of the maximum wind (382 m), capturing most of the energy in the super-gradient region. Figures 2 and 3 further illustrate significant differences in along-wind velocity profiles and power spectral densities between TC and NTS.

Based on Fig. 4, the building’s structural plan is elliptical (aerodynamic); however, for the purpose of evaluating the most critical case, the most unfavorable wind direction was considered, i.e., when the windward side is the long side of the ellipse (B=69.4\mathrm{m}, D/B=1.66 and C_{\mathrm{D}}=0.96). Accordingly, Table 3 shows the mechanical parameters for the benchmark building, assuming it is located in the city of Cancun, over terrain related to a coastal area exposed to the open sea \left({\mathit{\text{z}}}_0=0.001\mathrm{m}\right). Table 3 presents two case studies for the benchmark building: 1) the original model (C1) [65], where deformation is predominantly due to bending; 2) a modified model (C2), where effective shear stiffness is significant. In both cases, the fundamental period was kept constant to isolate the effect of the \alpha ratio on the response. Although the modified model (C2) is conceptual, similar variations in \alpha with minimal changes in period have been observed in real buildings, as reported by Huergo et al. [48].

It is to be expected that the benchmark building will exhibit a significant dynamic response associated with wind turbulence due to its large fundamental vibration period and small damping ratio ({\xi }_{\mathrm{s},i}=0.01), based on the values recommended by Tamura and Kareem [66]. For this reason, increasing the damping of the building represents a highly viable and less invasive solution from a structural perspective. In this study, an additional damping of {\xi }_{\text{FVD}}^{\text{EE}}=0.1 was chosen, which represents a realistic value for the design of this type of control devices. Consequently, Fig. 7 presents the values of c_{\text{FVD}} for FVDs with different values of \beta (computed using Eq. (45)), installed on each of the 60 stories of the benchmark building under the landfalling hurricane wind case (r = 10 km) shown in Fig. 2. The results indicate that for linear FVDs \left(\beta =1\right), the values of c_{\text{FVD}} remain constant regardless of wind intensity. In contrast, for nonlinear FVDs , the values of c_{\text{FVD}} increase as the wind intensity rises. Furthermore, for nonlinear FVDs under the same wind intensity, the target value of {\mathit{\xi }}_{\text{FVD}}^{\text{EE}} is achieved with lower values of c_{\text{FVD}}; that is, c_{\text{FVD}} decreases as \beta decreases, indicating that friction dampers \left(\beta =0\right) are significantly more efficient in controlling the structural response. Additionally, as \alpha decreases, c_{\text{FVD}} also decreases, suggesting that more cost-effective viscous fluid dampers can be designed for lateral load-resisting systems primarily governed by pure bending deformation.

To probabilistically assess the dynamic response of the benchmark building to both NTS and TC, a Monte Carlo simulation was performed following the flowchart shown in Fig. 8. According to this procedure, the reference wind speed \overline{U}\left({\mathit{\text{z}}}_r\right) was increased from 18 to 306 km/h in uniform increments, resulting in a total of 125 discrete wind speed values. For each wind intensity, along-wind loads were simulated 500 times, producing 62500 simulations per wind scenario. Additionally, three different radial distances to the hurricane eye were considered for the TC wind r = 50 km, r = 100 km, and r = 200 km. Including the NTS case, a total of four wind scenarios were analyzed, amounting to 250000 simulations. Moreover, four structural scenarios for the benchmark building were considered: C1 without FVDs, C1 with FVDs, C2 without FVDs, and C2 with FVDs, resulting in a total of 1000000 dynamic response simulations. In the probabilistic simulations, 60 vibration modes were included using Eq. (40) in accordance with the flowchart shown in Fig. 8. This ensures that higher-mode contributions are captured, which is particularly important in C1, where bending deformation is significant and higher modes, especially the second mode, significantly influence the building response.

As previously mentioned, the structural damping ratio and the target additional damping ratio used were {\xi }_{s,i}=0.01 and {\xi }_{\text{FVD,target}}^{\text{EE}}=0.1, respectively. Given the multi-story nature of the benchmark building, a single simulation can take several minutes on a standard computer; therefore, the DelftBlue supercomputer was employed for the computations.

In all cases involving FVDs, a single velocity exponent \beta =0.5 was adopted for the probabilistic assessment. This modeling decision was intentional: the chosen value represents a moderately nonlinear behavior, lying between the two common idealizations—linear viscous damping (\beta =1) and Coulomb-type friction damping (\beta =0). This intermediate value captures the nonlinear characteristics of practical FVDs without reducing the model to a purely linear or friction-based idealization. Moreover, the primary objective of this study is not to investigate the effects of varying the velocity exponent, but to assess the influence of wind load characterization, specifically, the contrast between NTS and TC, on the probabilistic structural performance of high-rise buildings, both with and without supplemental damping. Therefore, the velocity exponent was held constant throughout the analysis, and a parametric study on \beta is proposed as a valuable direction for future research.

To verify the reliability and statistical stability of the Monte Carlo simulations, a convergence analysis was performed for representative cases of the benchmark building. For clarity, only two key structural responses, peak rooftop displacement and RMS rooftop acceleration, under a representative reference wind speed of approximately 150 km/h in C1, both without and with FVDs, were considered. The analysis was carried out for two wind scenarios: a NTS and a TC at radial distances of r = 50 km. Figures 9 and 10 show the relative error of the mean and standard deviation as a function of the number of simulations. Thresholds of 1% for the mean and 3% for the standard deviation were chosen based on common practice in engineering Monte Carlo studies: the mean converges faster and requires a stricter criterion, while the standard deviation is more sensitive to data variability and extreme values. As shown in the figures, both quantities fall below these thresholds after approximately 300 simulations, whereas 500 simulations were performed for each wind intensity level, confirming that the results are statistically reliable for all considered dynamic response variables.

4.1 Bivariate analysis

Bivariate analysis and copula-based methods are essential in civil and structural engineering to accurately characterize the dependence between multiple variables during extreme events [67–70]. Incorporating such statistical dependence improves the reliability of risk assessments and vulnerability analyses, ultimately contributing to safer and more resilient structural designs against natural hazards and extreme loading scenarios.

Once the dynamic responses of the reference building were simulated through an incremental dynamic analysis, the statistical interpretation of the data set was performed (Fig. 8). This interpretation begins by analyzing the individual behavior of each output variable (rooftop displacement, rooftop velocity, rooftop acceleration, and peak interstory drift ratio) through marginal distribution fitting. Subsequently, the interaction between variables is evaluated in pairs, commonly known as bivariate analysis. For the bivariate analysis, copulas, semi-correlations, and the Cramér–von Mises criterion were used to gain a deeper understanding of the dependence between variables.

Copulas provide a powerful and flexible framework for modeling the dependence between random variables, separating the joint behavior from the individual marginal distributions. In the bivariate case, a copula is a joint distribution function with uniform marginals on the interval [0,1], meaning it is defined over the unit square [0,1]², which includes all pairs \left(u,v\right)where u=F_X\left(x\right) and v=F_Y\left(y\right) are the cumulative probabilities associated with the original variables X and Y. According to Sklar’s theorem [71], any bivariate joint distribution H_{XY}\left(x,y\right) can be expressed in terms of its marginal cumulative distribution functions (CDFs) F_X\left(x\right) and F_Y\left(x\right), and a copula C that captures their dependence structure:

(57)H_{XY}\left(x,y\right)=C\left\{F_X\left(x\right),F_Y\left(y\right)\right\}=C\left(u,v\right).

Several families of copulas are commonly used to model different types of dependence structures between variables. Among them, the Gaussian, Frank, Gumbel, and Clayton copulas are widely applied due to their distinct characteristics.

The Gaussian copula models symmetric-linear dependence and is based on the multivariate normal distribution. It is defined as

(58)C_{\text{Gauss}}\left(u,v;\rho \right)={\mathrm{\Phi }}_{\text{ρ}}\left\{{\mathrm{\Phi }}^{-1}\left(u\right),{\mathrm{\Phi }}^{-1}\left(v\right)\right\},

where {\mathrm{\Phi }}^{-1}(\cdot ) is the inverse of the standard cumulative function (CDF), and {\mathrm{\Phi }}_{\rho }\left\{\cdot ,\cdot \right\}is the bivariate normal CDF with correlation coefficient \rho \in \left(-1,1\right). This copula captures linear correlation but does not emphasize tail dependence, making suitable for symmetric and moderate dependencies.

The Frank copula is an Archimedean copula that captures moderate dependence evenly across the distribution range without focusing on extreme values. Its form is

(59)C_{\text{Frank}}\left(u,v;\theta \right)=-\frac{1}{\theta }\ln \left(1+\frac{\left(e^{-\theta u}-1\right)e^{-\theta v}-1}{e^{-\theta }-1}\right),

where \theta \in \mathbb{R}\mathrm{\setminus }\left\{0\right\} controls the strength and direction of the dependence. The Frank copula is symmetric and does not exhibit tail dependence, making it suitable for moderate, balanced relationships.

The Gumbel copula models stronger dependence in the joint upper tail, making it useful for extreme event analysis where high values tend to occur together. It is defined as

(60)C_{\text{Gumbel}}\left(u,v;\theta \right)=\exp \left\{-{\left[{\left(-\ln \text{}u\right)}^{\theta }+{\left(-\ln \text{}v\right)}^{\theta }\right]}^{1/\theta }\right\},

with parameter \theta \ge 1. As \theta increases, the strength of upper tail dependence grows. This copula is asymmetric and particularly useful for modeling extreme event co-occurrence.

The Clayton copula emphasizes lower tail dependence, capturing stronger dependence when both variables take low values simultaneously. It is given by

(61)C_{\text{Clayton}}\left(u,v;\theta \right)={\left(u^{-\theta }+v^{-\theta }-1\right)}^{-1/\theta },

where \theta >0. Larger \theta values imply stronger lower tail dependence. The Clayton copula is asymmetric and commonly used when joint low outcomes are of interest.

Accordingly, Table 4 and Fig. 11 presents the copula fitting results for each wind scenario, totaling 170 cases across 17 wind velocity intensities (in 18 km/h increments up to 306 km/h). These results indicate a strong tendency toward a Gaussian fit for both NTS wind and TC wind, regardless of the building’s lateral load-resisting system. Interestingly, for TC wind, Frank-type distributions emerge in the building without FVDs, while Gumbel-type distributions appear in the building with FVDs. Moreover, wind intensity leads to a diversification of copulas for data pairs across all case studies.

On the other hand, Fig. 12 shows the copula fitting results for each pair of dynamic responses, considering a total of 10 response pairs and 4 wind scenarios, resulting in 68 copula fittings (17 wind velocity intensities × 4 wind scenarios). The results in Fig. 12 exhibit a behavior very similar to that of Fig. 11, with most data pairs following a Gaussian fit. However, a notable exception is the data pair of peak displacement and peak drift in C1 (α = 2.95) without FVDs, where Gaussian fits are observed alongside a comparable occurrence of Frank- and Gumbel-type fits. On the other hand, when the building’s lateral load-resisting system is modified (C2, α = 10), this data pair returns to a Gaussian behavior. Likewise, in C1 (α = 2.95) with FVDs, the data pair u_{\text{peak}}\text{vs}{\text{Drift}}_{\text{peak}} varies only between the Gaussian and Frank fits, eliminating the tail dependence effect previously observed in the Gumbel fits. Furthermore, incorporating FVDs into the building enhances the likelihood of Gaussian fits for the data pairs, ensuring symmetric dependence on both low and high variable values while reducing the occurrence of tail-dependent fits.

For the bivariate analysis of dynamic responses, the correlation coefficient between paired data was examined. Spearman’s rank correlation coefficient was selected due to its ability to capture both linear and nonlinear monotonic relationships. It is defined as the Pearson correlation coefficient computed on the ranked data:

(62){\rho }_s=\frac{{\displaystyle \sum _{i=1}^n\left(R_{X_i}-{\overline{R}}_{{}_X}\right)\left(R_{Y_i}-{\overline{R}}_{{}_Y}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(R_{X_i}-{\overline{R}}_{{}_X}\right)}^2}}\sqrt{{\displaystyle \sum _{i=1}^n{\left(R_{Y_i}-{\overline{R}}_{{}_Y}\right)}^2}}},

where R_{X_i} and R_{Y_i} are the ranks of the two random variables X and Y, respectively, and {\overline{R}}_{{}_X}, {\overline{R}}_{{}_Y} are the corresponding mean ranks.

The analysis of the bivariate response of the benchmark building highlights three key observations. First, the presence of FVDs significantly increases the correlation between dynamic response data pairs, particularly those involving RMS acceleration. Second, the correlation of dynamic responses is generally higher under TC conditions compared to NTS. Third, correlations tend to increase with rising wind intensity, reflecting stronger coupling of structural responses at higher wind speeds. Figures 13–15 present the probability summaries for the most relevant data pairs supporting these conclusions, including copula fittings and Spearman correlation coefficients, for three wind intensities (\overline{U}\left({\mathit{\text{z}}}_r\right)=20\text{km/h}, \overline{U}\left({\mathit{\text{z}}}_r\right)=40\text{km/h} and \overline{U}\left({\mathit{\text{z}}}_r\right)=70\text{km/h}, etc.) and four wind scenarios: NTS, TC with r = 50 km (TC-50), TC with r = 100 km (TC-100), and TC with r = 200 km (TC-200).

From Figs. 13-15, there are certain notable situations. First, in the data pair of peak acceleration vs RMS acceleration, for the buildings without FVDs, the maximum correlation values are significantly lower than those presented in the other data pairs. Another interesting observation is that the correlation of the dynamic responses tends to increase as the wind intensity rises. Additionally, it is noteworthy that the dynamic response correlations are much higher for the TC case compared to the NTS.

Another observation is in C1 (α = 2.95) without FVDs, where all data pairs representing the dynamic response of RMS acceleration show correlation values lower than those in the other cases. A similar situation occurs in C2 (α = 10) without FVDs, though with slightly higher correlations. In contrast, for buildings with FVDs, the correlations of the data pairs increase significantly, with the effect being particularly noticeable for the data pairs involving RMS acceleration.

Building on the previous observations, when analyzing the wind intensity ranges at which these values occur, it is evident that in cases where shear lateral deformation predominates, the data pairs show high correlation values at lower wind intensities and low correlation values at medium to high wind intensities. In contrast, in cases where flexural lateral deformation predominates, the highest correlation values are observed at higher wind intensities, while the lowest values occur at lower wind intensities.

4.2 Fragility curves

Fragility curves are used to represent the probability of a structure or system exceeding a specific damage threshold under varying levels of external forces, such as wind. These curves are crucial in risk and reliability evaluation, as they provide essential information for assessing the vulnerability of structures.

To construct the fragility curves, the probability of exceeding a predefined damage threshold x_{\mathrm{c}} was evaluated as a function of the reference mean wind velocity \overline{U}\left({\mathit{\text{z}}}_r\right), measured at a standard height {\mathit{\text{z}}}_r (typically 10 m). For each wind intensity, the structural demand was compared against the threshold x_{\mathrm{c}}, and the exceedance probability was calculated as

(63)P_f\left(\overline{U}\left({\mathit{\text{z}}}_r\right)\right)=P\left(\text{Demand}>x_{\text{c}}|\overline{U}\left({\mathit{\text{z}}}_r\right)\right).

This exceedance probability was empirically estimated from the simulation results at each wind speed by computing the ratio of cases in which the demand surpassed the threshold. Then, different parametric CDFs, including Normal, Lognormal, Logistic, Weibull, and Exponential distributions, were fitted to the exceedance probabilities using Maximum Likelihood Estimation. The likelihood function provides a measure of how well a statistical model explains the observed data based on its parameters. To determine the most appropriate model, the Akaike Information Criterion (AIC) was used. The AIC is a widely accepted model selection criterion that balances model fit and complexity by penalizing models with a larger number of parameters. It is defined as

(64)AIC=2k-2\ln \hat{L},

where k is the number of parameters of the model and \hat{L} is the maximum value of the likehood function. The AIC evaluates the trade-off between the goodness of fit and model complexity. Models with lower AIC values are preferred, as they indicate a better balance between simplicity and fit.

In our analysis, a statistical model was fitted to each dynamic response at every wind intensity across all simulated scenarios. Table 5 presents the set of probability distributions that best fit each demand parameter, considering 17 wind velocity intensities (from 0 to 306 km/h in 18 km/h increments) under four wind scenarios: a NTS and three TC cases (with hurricane center distances of r = 50, 100, and 200 km). This results in a total of 68 fitted distributions per response parameter for the benchmark building.

In C1 (α = 2.95) without FVDs, the response distributions exhibit predominantly lognormal behavior. However, peak accelerations in approximately half of the responses deviate from this trend, with most following a normal distribution and a few fitting a logistic distribution. On the other hand, in C2 (α = 10) without FVDs, it is observed that peak accelerations tend to be more concentrated within a lognormal distribution.

The most representative observations occur in C1 (α = 2.95) with FVDs, where the frequency of normal distributions increases compared to the other cases. This results in a frequency similar to that of lognormal distributions, except for accelerations. In C2 (α = 10) with FVDs, a phenomenon similar to that in C1 with FVDs is observed, where the response results, except for the RMS acceleration, transition from a predominantly lognormal behavior to a more normal distribution of the building’s dynamic responses. In the uncertainty analysis of structural response in buildings subjected to turbulent wind, it is commonly accepted to assume lognormal distribution functions [16]. However, the results in Table 5 show that increasing structural damping through control devices can alter the tendency toward lognormal or normal distributions.

For the fragility curves, a finer discretization of the wind velocity was used by applying increments of approximately 0.68 m/s, resulting in 125 wind intensity values from 0 up to 306 km/h. This dense sampling was chosen to ensure that the fragility curves are smooth and well-defined across the full range of wind speeds considered, leading to a total of one million simulations across all scenarios.

For structures subjected to wind loads, wind design codes require the evaluation of two limit states: 1) ultimate limit state; 2) serviceability limit state. According to the Mexican wind code [6], the damage threshold x_{\mathrm{c}} for peak rooftop displacement is set at H/500, where H represents the total building height, while the threshold for peak inter-story drift ratio is specified as 0.002. For the serviceability limit state, it is necessary to ensure that the acceleration does not exceed a specified comfort level. Ferrareto et al. [72] suggest that peak acceleration is linked to the perception of motion, whereas RMS acceleration is a better parameter for evaluating comfort during sustained motion. In this study, the comfort threshold for RMS acceleration proposed by Irwin [73] is used, which is set at 0.04 m/s² for a building’s fundamental vibration period of T_1=5.17\mathrm{s}. Based on these considerations, fragility curves were constructed for peak rooftop displacement, peak inter-story drift ratio, and RMS rooftop acceleration, as these response parameters are directly associated with the ultimate and serviceability limit states defined in current wind design guidelines. These curves are presented in Figs. 16 to 18.

Based on the results of the fragility curves, benchmark buildings without FVDs have the highest probability of exceeding the maximum displacement threshold in a Category 4 hurricane when the wind is modeled as a NTS. However, when modeled as a TC, this threshold is exceeded in a Category 3 hurricane. On the other hand, it is observed that for a lateral resisting system with predominant lateral deformation in bending, high exceedance probabilities occur at lower wind intensities. Additionally, it is also observed that the probability of exceeding the damage threshold increases as the building gets closer to the eye of the hurricane. For the same case of fragility curves for peak rooftop displacements, it is observed that the use of FVDs shifts the critical zone of the fragility curves for TC from a Category 3 hurricane to a Category 4 hurricane. For the case of inter-story drift ratios, the critical zone is located between a Category 2 and Category 3 hurricane; however, the use of FVDs shifts these exceedance probabilities to the Category 4 and Category 5 hurricane range.

Additionally, the RMS acceleration fragility curves indicate that for α = 2.95, the highest probability of exceeding the comfort threshold occurs in Category 3 hurricanes, while for α = 10, this probability lies between Category 3 and Category 4 hurricanes. However, the use of FVDs was found to be highly effective in controlling RMS accelerations, shifting the critical zone to wind intensities of 300 km/h (Category 5 hurricane). This resulted in exceedance probabilities ranging from 0 to 70% for α = 2.95 and from 0 to 90% for α = 10.

5 Conclusions

This study evaluated the effects of turbulent wind from NTS and TC on the dynamic response of high-rise buildings, assessing the performance of FVDs in different lateral systems through probabilistic simulations and advanced analysis techniques. The findings support the development of wind-resilient high-rise design by improving understanding of response statistics, dependence patterns, and fragility under extreme winds. Key conclusions are as follows.

1) Dynamic responses generally follow a lognormal distribution, as assumed in many turbulence uncertainty analyses; however, when lateral deformation is dominated by flexure and FVDs are used, the distribution may shift toward a normal distribution, affecting serviceability limit state assumptions.

2) Bivariate analysis shows most data pairs follow a Gaussian copula, except in flexure-dominated deformation without FVDs, where Frank and Gumbel copulas fit peak displacement and drift. Adding shear stiffness or FVDs restores Gaussian behavior, reducing tail dependence and improving symmetry.

3) In buildings without FVDs, correlation between peak acceleration and RMS acceleration is lower than for other data pairs. Additionally, correlations involving RMS acceleration are generally lower than in other cases.

4) Dynamic response correlations increase with wind intensity and are higher for TC than NTS. Correlations in buildings with FVDs rise notably, especially for RMS acceleration pairs.

5) High correlations occur at low wind intensities in shear-dominated lateral deformation, and at high wind intensities in flexure-dominated cases.

6) At the ultimate limit state, modeling wind as a NTS in buildings without FVDs yields fragility curves similar to those with FVDs under TC winds. Thus, damper effectiveness may be overestimated if wind is simplified as non-tropical.

7) Increasing overall effective shear stiffness reduces fragility exceedance probability. This can be achieved via rigid connections, optimized interstory stiffness, stiffer cores, shear walls, and similar strategies.

8) FVDs effectively reduce displacements, interstory drifts, and accelerations from wind loads, with the greatest benefit in improving serviceability through acceleration reduction.

9) For linear FVDs, damping coefficients remain constant with wind intensity; for nonlinear FVDs, they increase as wind intensity rises. This suggests nonlinear FVDs may be more efficient, but this preliminary conclusion requires full probabilistic validation. Future research should analyze nonlinearity effects on fragility and response correlations.

These results reveal gaps in current wind regulations that simplify hurricane winds, risking underestimation of hazards. While some international standards consider dampers mainly for seismic design, Mexican regulations include dampers for seismic control but do not explicitly incorporate FVDs or other vibration control devices for wind loads. Our findings on response correlations, bivariate dependence modeled by copulas, and fragility curves highlight the importance of capturing complex statistical dependencies and accurate failure probabilities in design. Incorporating probabilistic hurricane modeling and realistic damper behavior into wind design regulations will improve accuracy and damper efficiency, ultimately enhancing structural resilience. This supports safer, more sustainable coastal cities better prepared for extreme winds.

References

[1] Holland G , Bruyère C L . Recent intense hurricane response to global climate change. Climate Dynamics, 2014, 42(3–4): 617–627

[2] Zhang J A , Rogers R F , Nolan D S , Marks F D . On the characteristic height scales of the hurricane boundary layer. Monthly Weather Review, 2011, 139(8): 2523–2535

[3] Powell M D , Houston S H . Surface wind fields of 1995 hurricanes Erin, opal, Luis, Marilyn, and Roxanne at landfall. Monthly Weather Review, 1998, 126(5): 1259–1273

[4] Fang G , Zhao L , Cao S , Ge Y , Pang W . A novel analytical model for wind field simulation under typhoon boundary layer considering multi-field correlation and height-dependency. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 175: 77–89

[5] Brown DKelly L. Hurricane Otis TC Update. Miami, Florida: National Hurricane Center, 2023

[6] Federal Electricity Commission (CFE). Manual for the Design of Civil Works: Wind Design, 2020

[7] Vickery P J , Wadhera D , Powell M D , Chen Y . A hurricane boundary layer and wind field model for use in engineering applications. Journal of Applied Meteorology and Climatology, 2009, 48(3): 381–405

[8] Snaiki R , Wu T . A semi-empirical model for mean wind velocity profile of landfalling hurricane boundary layers. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 180: 249–261

[9] Petersen E , Mortensen N , Landberg L , Højstrup J , Frank H . Wind power meteorology. Part I: climate and turbulence. Wind Energy, 1998, 1(1): 2–22

[10] Kaimal J C , Wyngaard J C , Izumi Y , Coté O R . Spectral characteristics of surface-layer turbulence. Quarterly Journal of the Royal Meteorological Society, 1972, 98(417): 563–589

[11] Von Kármán T . Progress in the statistical theory of turbulence. Proceedings of the National Academy of Sciences of the United States of America, 1948, 34(11): 530–539

[12] Schroeder J L , Smith D A . Hurricane Bonnie wind flow characteristics as determined from WEMITE. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91(6): 767–789

[13] Yu B , Chowdhury A G , Masters F J . Hurricane wind power spectra, cospectra, and integral length scales. Boundary-Layer Meteorology, 2008, 129(3): 411–430

[14] Jung S , Masters F . Characterization of open and suburban boundary layer wind turbulence in 2008 Hurricane Ike. Wind and Structures, 2013, 17(2): 135–162

[15] Li LXiao YKareem ASong LQin P. Modeling typhoon wind power spectra near sea surface based on measurements in the South China sea. Journal of Wind Engineering and Industrial Aerodynamics, 2012, 104–106: 565–576

[16] Kareem A . Effect of parametric uncertainties on wind excited structural response. Journal of Wind Engineering and Industrial Aerodynamics, 1988, 30(1–3): 233–241

[17] Chen J , Zeng X , Peng Y . Probabilistic analysis of wind-induced vibration mitigation of structures by fluid viscous dampers. Journal of Sound and Vibration, 2017, 409: 287–305

[18] Zhou Y , Wang D , Deng X . Optimum study of wind-induced vibration control of high-rise buildings with viscous dampers. Wind and Structures, 2008, 11(6): 497–512

[19] Sun X , Wu H , Wu Y , Su N . Wind-induced responses and control of a kilometer skyscraper with mass and viscous dampers. Journal of Building Engineering, 2021, 43: 102552

[20] Li M , Shi M , Fu W . A variable damping viscous damper for control of buildings under wind loading. Journal of Building Engineering, 2024, 97: 110960

[21] Li X , Li Q S . Monitoring structural performance of a supertall building during 14 tropical cyclones. Journal of Structural Engineering, 2018, 144(10): 04018176

[22] Martín del Campo J O , Pozos-Estrada A , Pozos-Estrada O . Development of fragility curves of land-based wind turbines with tuned mass dampers under cyclone and seismic loading. Wind Energy, 2020, 24(7): 737–753

[23] Hernández-Barrios H , Patlán-Manjarrez C M , Gómez-Martínez R . The passive control of nonlinear dynamic response using the force analogy method in a chimney located in a wind- and seismic-prone region in Mexico. Buildings, 2025, 15(3): 459

[24] National Weather Service. Hurricane statistics. 2024 (available at the website of National Weather Service)

[25] FEMA . Hazus Hurricane Model Technical Manual. Federal Emergency Management Agency, 2021

[26] Deaves D MHarris R I. A Mathematical Model of the Structure of Strong Winds. C.R.I.A. Report 76. 1981

[27] Fang G , Wei M , Zhao L , Xu K , Cao S , Ge Y . Site- and building height-dependent design extreme wind speed vertical profile of tropical cyclones. Journal of Building Engineering, 2022, 62(105322): 105322

[28] Martin del Campo J O. Numerical studies on the evaluation of the wind-induced dynamic response of wind turbine support structures. Dissertation for the Doctoral Degree. Mexico: National Autonomous University of Mexico, 2022

[29] Holland G J , Belanger J I , Fritz A . A revised model for radial profiles of hurricanes winds. Monthly Weather Review, 2010, 138(11): 4393–4401

[30] Vickery P J , Wadhera D . Statistical models of Holland pressure profile parameter and radius to maximum winds of hurricanes from flight-level pressure and H*Wind data. Journal of Applied Meteorology and Climatology, 2008, 47(9): 2497–2517

[31] Solari G . Gust buffeting, I: Peak wind velocity and equivalent pressure. Journal of Structural Engineering, 1993, 119(2): 365–373

[32] Davenport A G . The spectrum of horizontal gustiness near the ground in high winds. Quarterly Journal of the Royal Meteorological Society, 1961, 87(372): 194–211

[33] ESDU . Characteristic of Atmospheric Turbulence Near the Ground, Part II: Single Point Data for Strong Winds (Neutral Atmosphere). London: Engineering Sciences Data Unit, IHS Inc., 2001

[34] Huergo I F. Passive control of wind-induced vibrations in tall buildings with flexible base. Dissertation for the Doctoral Degree. Mexico: National Autonomous University of Mexico, 2021

[35] Huergo I F , Hernández-Barrios H , Gómez-Martínez R . Analytical simulation of 3D wind-induced vibrations of rectangular tall buildings in time domain. Shock and Vibration, 2022, 2022: 7283610

[36] Li L , Kareem A , Xiao Y , Song L , Zhou C . A comparative study of field measurements of the turbulence characteristics of typhoon and hurricane winds. Journal of Wind Engineering and Industrial Aerodynamics, 2015, 140: 49–66

[37] Amirinia G , Jung S . Along-wind buffeting responses of wind turbines subjected to hurricanes considering unsteady aerodynamics of the tower. Engineering Structures, 2017, 138: 337–350

[38] Amirinia G , Jung S . Along-wind response of high-rise buildings subjected to hurricane boundary layer winds. Journal of Structural Engineering, 2017, 143(9): 04017087

[39] Liu Z , Fang G , Hu X , Xu K , Zhao L , Ge Y . Stochastic power spectra models for typhoon and non-typhoon winds: A data-driven algorithm. Journal of Wind Engineering and Industrial Aerodynamics, 2022, 231(105214): 105214

[40] Fang G , Liu Z , Pang W , Zhao L , Xu K , Cao S , Ge Y . Probabilistic gust factor model of typhoon winds. Journal of Structural Engineering, 2024, 150(1): 04023205

[41] Dyrbye C DHansen S O. Wind Loads on Structures. Hoboken, NJ: John Wiley & Sons, 1997

[42] Vickery B J . Load fluctuations in turbulent flow. Journal of the Engineering Mechanics Division, 1968, 94(1): 31–46

[43] Shinozuka M , Yun C B , Seya H . Stochastic methods in wind engineering. Journal of Wind Engineering and Industrial Aerodynamics, 1990, 36(2): 829–843

[44] Gao Y , Wu Y , Li D , Liu H , Zhang N . An improved approximation for the spectral representation method in the simulation of spatially varying ground motions. Probabilistic Engineering Mechanics, 2012, 29: 7–15

[45] Davenport A G . The response of slender line-like structures to a gusty wind. Proceedings Institution of Civil Engineers, 1962, 23(3): 389–408

[46] Krenk S. Wind field coherence and dynamic wind forces. In: Proceedings of the IUTAM Symposium Held in Trondheim. Dordrecht: Springer Netherlands, 1996: 269–278

[47] Hansen S O , Krenk S . Dynamic along-wind response of simple structures. Journal of Wind Engineering and Industrial Aerodynamics, 1999, 82(1–3): 147–171

[48] Huergo I F , Hernández-Barrios H , Gómez-Martínez R . Along and across-wind vibration control of shear wall-frame buildings with flexible base by using passive dynamic absorbers. Wind and Structures, 2024, 38(1): 15–42

[49] Huergo IFHernández-Barrios HMorales-Nápoles O. Wind-induced vibration control of base-isolated high-rise buildings by using tuned mass damper inerter. Structure and Infrastructure Engineering, 2025: 1–26

[50] Miranda E , Taghavi S . Approximate floor acceleration demands in multistory buildings. I: Formulation. Journal of Structural Engineering, 2005, 131(2): 203–211

[51] Hart G CWong K. Structural Dynamics for Structural Engineers. Hoboken, NJ: John Wiley & Sons, 2000

[52] Seleemah AConstantinou M C. Investigation of Seismic Response of Buildings with Linear and Nonlinear Fluid Viscous Dampers. Report No. NCEER 970004. 1997

[53] Yeh F YChang K CChen T WYu C H. The dynamic performance of a shear thickening fluid viscous damper. Journal of the Chinese Institute of Engineers, 2014, 37(8): 983–994 (in Chinese)

[54] Zhou H , Yan L , Jiang W , Xuan S , Gong X . Shear thickening fluid–based energy-free damper: Design and dynamic characteristics. Journal of Intelligent Material Systems and Structures, 2014, 27(2): 208–220

[55] Wei M , Lin K , Guo Q , Sun H . Characterization and performance analysis of a shear thickening fluid damper. Measurement and Control, 2019, 52(1–2): 72–80

[56] Chen S , Fu X , Meng P , Cheng L , Wang L , Yuan J . Fluid flow modeling and experimental investigation on a shear thickening fluid damper. Buildings, 2024, 14(11): 3548

[57] Ramírez O MConstantinou M CKircher C AWhittaker A SJohnson M WGomez J DChrysostomou C Z. Development and Evaluation of Simplified Procedures for the Analysis and Design of Buildings with Passive Energy Dissipation Systems. Buffalo, NY: State University of New York at Buffalo, 2001

[58] Lin W H , Chopra A K . Earthquake response of elastic SDF systems with non-linear fluid viscous dampers. Earthquake Engineering & Structural Dynamics, 2002, 31(9): 1623–1642

[59] Lin T K , Hwang J S , Chen K H . Optimal distribution of damping coefficients for viscous dampers in buildings. International Journal of Structural Stability and Dynamics, 2017, 17(4): 1750054

[60] Hwang J S , Huang Y N , Yi S L , Ho S Y . Design formulations for supplemental viscous dampers to building structures. Journal of Structural Engineering, 2008, 134(1): 22–31

[61] Zeng X , Peng Y , Chen J . Serviceability-based damping optimization of randomly wind-excited high-rise buildings. Structural Design of Tall and Special Buildings, 2017, 26(11): e1371

[62] De Domenico D , Ricciardi G , Takewaki I . Design strategies of viscous dampers for seismic protection of building structures: A review. Soil Dynamics and Earthquake Engineering, 2019, 118: 144–165

[63] Huergo I F , Hernández-Barrios H , Patlán C M . A continuous-discrete approach for pre-design of flexible-base tall buildings with fluid viscous dampers. Soil Dynamics and Earthquake Engineering, 2020, 131: 106042

[64] FEMA ASCE . Prestandard and Commentary for the Seismic Rehabilitation of Buildings. Report FEMA-356. 2000

[65] Lai X , He Z , Li M . Strain energy-based uniformity spectrum for elastic deformation control of super high-rise buildings under near-fault pulse-like earthquake excitations. Soil Dynamics and Earthquake Engineering, 2024, 181: 108663

[66] Tamura GKareem A. Advanced Structural Wind Engineering. Berlin: Springer, 2013

[67] Mendoza-Lugo M A , Delgado-Hernández D J , Morales-Nápoles O . Reliability analysis of reinforced concrete vehicle bridges columns using non-parametric Bayesian networks. Engineering Structures, 2019, 188: 178–187

[68] Mendoza-Lugo M A , Nogal M , Morales-Nápoles O . Estimating bridge critically due to extreme traffic loads in highway networks. Engineering Structures, 2024, 300: 117172

[69] Mares-Nasarre P , Muscalus A , Haas K , Morales-Nápoles O . The probabilistic dependence of ship-induced waves is preserved spatially and temporally in the Savannah River (USA). Scientific Reports, 2024, 14(1): 28154

[70] Paprotny D , Hart C M P , Morales-Nápoles O . Evolution of flood protection levels and flood vulnerability in Europe since 1950 estimated with vine-copula models. Natural Hazards, 2024, 121(5): 6155–6184

[71] Sklar A . Multivariate distribution functions and their margins. Publications of the Institute of Statistics of the University of Paris, 1959, 8: 229–231

[72] Ferrareto J , Mazzilli C , França R . Wind-induced motion on tall buildings: A comfort criteria overview. Journal of Wind Engineering and Industrial Aerodynamics, 2015, 142: 26–42

[73] Irwin A W. Motion in tall buildings. Second Century of the Skyscraper: Council on Tall Buildings and Urban Habitat. Boston, MA: Springer US, 1986: 759–778

RIGHTS & PERMISSIONS

Higher Education Press