Feasibility of a nonlinear spectral approach for peak floor acceleration of multi-story bilinear hysteretic buildings

To protect the acceleration-sensitive non-structural components in a multi-story building, a proper seismic design of them against the peak floor acceleration (PFA) is needed. The PFA, therefore, should be estimated in advance. To efficiently estimate the PFA of multi-story buildings at the inelastic stage, a nonlinear modal combination approach is provided in this study. First, for the single-story building (single-degree-of-freedom system), a PFA spectral ratio δ, defined as the ratio between the inelastic PFA and the elastic spectral amplitude, is proposed to measure the nonlinear reductions in PFA. A series of time history analysis (THA) revealed that there is a dependable relationship between the PFA spectral ratio δ, the strength reduction factor R, the post-yield stiffness ratio α, and the period T of the structure. The scatter of δ is notably smaller than the ductility demand, showing the good feasibility of a δ-R-T-α relationship for predicting δ. Second, for the multi-story building (multi-degree-of-freedom system), the floor acceleration is decoupled in the modal space, the contribution of each mode to the elastic PFA is multiplied by δ to account for the nonlinear reduction in PFA at the inelastic stage. Here δ can be determined by an appropriate δ-R-T-α relationship, while R and α are determined via a modal pushover analysis. The modified modal contributions are combined via the complete-quadratic combination (CQC) rule, forming the inelastic CQC approach. The PFAs of some 2-, 5-, and 8-story building structures with different degrees of nonlinearity are computed via the THA and the inelastic CQC approach. Results demonstrate the satisfactory accuracy of the latter. Notably, the effects of nonlinearities associated with the higher-order modes on the nonlinear PFA are well considered in the inelastic CQC, thus the un-acceptable over-estimations of the inelastic PFA, which occurs in conventional methods, are avoided.

The ith and jth modal participation factor α Post-yield stiffness ratio α 0 Post-yield stiffness ratio of each story α i α For the ith modal oscillator δ PFA spectral ratio δ CoV CoV of δ δ (R, T, α) δ Determined by T, R, and α δ i (R i , T i , α i ) δ For the ith modal oscillator ζ Damping ratio ζ i , ζ j Damping ratio of the ith and jth mode ι Influence vector μ Ductility demand μ CoV CoV of μ ρ ij Cross-modal correlation coefficient The ith and jth mode shape vector ϕ i,k , ϕ j,k The kth element of ϕ i and ϕ j ω i , ω j Circular frequency of the ith and jth mode

Introduction
To mitigate the earthquake-induced damage of the non-structural components (NSCs) in a multi-story building, the connection between the NSCs and the primary structure of the building should be appropriately designed (Anajafi et al. 2020;Villaverde 2006;Kazantzi et al. 2020). In doing so, the strength demand of the NSC-structure connections, which generally equals to the seismic force on a NSC, should be well predicted. To obtain the seismic force on a NSC, a floor response spectrum (FRS) is needed (Shooshtari et al. 2010). While in constructing the FRS, the peak floor acceleration (PFA) should be estimated as a priority. This is because the PFA defines the amplitude of FRS at T = 0 s (T is period), which in most cases shows the seismic force on a NSC with a short natural period (Singh et al. 2006). Many studies on the elastic PFA have been conducted, and some of them are incorporated in the FRS model (Lucchini et al. 2017;Jiang et al. 2015;Li et al. 2015;Miranda and Taghavi 2005;Pinkawa et al. 2014). The response spectrum method for the elastic PFA has been well developed by Moschen and Adam (2017), Moschen et al. (2016), Taghavi and Miranda (2009), and Pozzi and Der Kiureghian (2015). Their proposals involved combination of modal accelerations via the complete-quadratic combination (CQC) rule. Specially, Moschen and Adam (2017) and Moschen et al. (2016) addressed some approximations and simplifications related to the modal combination, while Pozzi and Der Kiureghian (2015) processed the rigid contribution of the truncated higher-order modes.
For the inelastic structures, on the other hand, predicting PFA remains a challenge. Due to the structural nonlinearity, the linear relationship between the PFA and the excitation intensity (usually measured by the peak ground acceleration, PGA) no longer exists, and the nonlinear time history analysis (THA) is needed to obtain the inelastic PFA. The nonlinear THA on RC structures conducted by Surana et al. (2017) and Petrone et al. (2015Petrone et al. ( , 2016 revealed that the ratio between PFA and PGA decreased as the value of R increased. Here R is the structural strength reduction factor, i.e., R = F e /F y , whereas F e and F y is the elastic seismic force demand and the structural yield force, respectively. Vukobratović et al. (2021) provided vertical distributions of the inelastic PFA. Their results were obtained from three shake-table tests on RC building structures. According to the results, the PFA estimations given by the prevailing codes are conservative. Berto et al. (2020) found that the PFA of RC structures was significantly mitigated as the nonlinearity of buildings became remarkable. Lepage et al. (2012) used shaking table tests and nonlinear THA results to generate an empirical prediction equation for the PFA of RC structures, in which an empirical response modification factor was adopted to account for the structural nonlinearity. In Lepage et al. (2012), the predicted inelastic PFA increased linearly along the height of the building, which was an over-simplification. Ray-Chaudhuri and Hutchinson (2011) conducted THA on eight steel moment resisting frames, demonstrating that the PFA could be either moderately reduced or significantly amplified due to the structural nonlinearity. Based on parametrical analyses on the ductile frame models, Medina et al. (2006), Sankaranarayanan and Medina (2007), Wieser et al. (2013), Lucchini et al. (2014), Uma et al. (2010), and Flores et al. (2015 manifested the dependency of PFA on the structural nonlinearity as well. Meanwhile, Miranda and Taghavi (2009) and Taghavi and Miranda (2012) considered different types of lateral force-resisting systems and revealed that, in addition to nonlinearity, the structural dynamic properties were quite influential to the inelastic PFA.
The preceding reports on the inelastic PFA were mostly based on the nonlinear THA of specific structural models. Therefore, their conclusions may exhibit certain limitations and may not be suitable for the general application. Due to the lack of a generalized model or approach, some over-simplified methods are recommended in the design codes, such as ASCE 7-10 (2010) and Eurocode 8 (2004) wherein the inelastic PFA are assumed to have the same pattern with the elastic scenario, which is, apparently, inaccurate.
To address the problem, efforts have been made by some forerunners. Rodriguez et al. (2002) proposed a "first mode reduction (FMR)" method to account for the nonlinear reduction in PFA. However, the FMR ignored the nonlinearities associated with the higher-order modes, thus the outcome of FMR was quite conservative. This indicates that, although it may be feasible to assume the higher-modes to be elastic in estimating the story drift (Chopra et al. 2004), it is not appropriate to do so in predicting the PFA. Moschen et al. (2013) suggested that the inelastic PFA could be quantified via a pushover analysis, by which the PFA reduction and the structural ductility could be connected. Still, their numerical example showed some large errors because the nonlinearities of the higher-order modes were omitted. Politopoulos and Feau (2007) focused on the nonlinear single-degree-of-freedom (SDF) structures, and proposed a modified equivalent oscillator to account for the structural nonlinear effect. Sullivan et al. (2013) and Calvi (2014) analysed the influence of a nonlinear SDF structure on the amplitude of FRS. The shortcomings of Politopoulos and Feau (2007), Sullivan et al. (2013), andCalvi (2014) are: 1) they ignored the effect of structural post-yield stiffness, 2) the multi-modal contribution to the inelastic PFA was not sufficiently incorporated. Based on the extensive parametrical analyses, Fajfar (2015, 2016) quantified the influence of structural ductility and hysteresis behaviour on the FRS of the SDF structures, and combined the SDF results with the pushover method for the multi-degree-of-freedom (MDF) structures. However, Fajfar (2015, 2016) mainly considered the ideal elasto-plastic system and system with a decreasing stiffness.
In view of the above-mentioned problems, the inelastic PFA is further studied in this paper. The objective of this study is to propose a generalized spectrum method for estimating the PFA of the bilinear hysteretic MDF systems with a positive post-yield stiffness. Such systems mainly represent the multi-story buildings structured in ductile momentresisting steel frames, as shown in Fig. 1a. For clarity, the hysteresis of each story of the MDF system can be simplified as a bilinear curve, as shown in Fig. 1b.
In validating the feasibility of the proposal, two aspects are specifically discussed: (1) the effect of a positive post-yield stiffness on the PFA of the SDF systems, (2) the influence of nonlinearities associated with the higher-order modes on the PFA of the MDF systems. Firstly, for the SDF systems (representing single-story structures), the feasibility of predicting the inelastic PFA based on R, T, and α is studied. Here α is the post-yield stiffness ratio, i.e., α = k py /k e , where k py and k e is the post-yield stiffness and the elastic stiffness of the system, respectively. Secondly, for the MDF system (multistory structure), an inelastic CQC approach is proposed and validated. In the inelastic CQC approach, the structural PFA associated with each natural mode is adjusted based on a δ-R-T-α relationship and then combined via the conventional CQC rule. Extensive numerical computations are conducted to prove the good feasibility of the proposal for the practical application.

PFA of SDF system
For the SDF system, the nonlinear reduction in PFA (quantified by δ) is related to the dynamic property of the system (represented by period T) and the degree of nonlinearity (measured by strength reduction factor R and post-yield stiffness ratio α). Given this observation, one would like to know if it is possible to get a dependable relationship between R, α, T, and δ. For brevity, we denote such a relationship as where δ(T,R,α) represents the nonlinear reduction in PFA for a pre-set T, R, and α of a SDF system, which is called the "PFA spectral ratio" hereafter, S a (T) is the elastic pseudoacceleration spectral amplitude at period T, while PFA NL (T,R,α) is the peak floor acceleration (total acceleration) of a nonlinear SDF system with period T, strength reduction factor R, and post-yield stiffness ratio α.
To prove the existence of a δ-R-T-α relationship, or in other words, to validate the feasibility of developing an Eq. (1), a series of nonlinear THA is conducted subsequently. Two groups of hazard-consistent earthquake records are adopted in the computation, as instructed in the following context. For clarity, it is worth noting that the target of this work is to prove the feasibility of a δ-R-T-α relationship, not to propose any specific δ(T,R,α) formula for any hazard scenario. In other words, the following context is to prove the predictability of δ, not to generate the expression of Eq. (1). (1)

Yield point (dy, Fy)
Inter-story shear force, F Inter-story drift, d

Earthquake records
The earthquake records, which were selected from certain consistent seismic hazard scenarios, were downloaded from the Japan NIED KiK-net. The first collection was from the Kyushu area of Japan, containing 92 records. The second was from the Osaka-Nagoya-Kanazawa area and consists of 100 records. The hazard scenario (including location of hypocentre, moment magnitude, JMA intensity) are shown in Fig. 2. The 3%-damped pseudo-acceleration spectra of the records are plotted in the figure as well. In general, the predominate period of the spectra for the Kyushu record collection (set A) is longer than that for the Osaka-Nagoya-Kanazawa collection (set B) because the soil of the latter is relatively firm.

Parametrical analyses
A wide range of T, R, and α were considered in computing the values of PFA NL (T,R,α) and S a (T) in Eq.
(1). The period T ranged from 0.1 to 2.0 s, the strength reduction factor R varied between 1.0 (elastic) and 5.0, while the post-yield stiffness ratio α spanned from 0.05 to 0.80. In the computation, the damping ratio ζ was set to 3% (a common value recommended by many design codes for steel structures). Additional computations revealed that the δ-R-T-α relationship is not quite influenced by the value of ζ in the low-damping range (ζ ≤ 5%). Therefore, although the properties of the δ-R-T-α relationship is discussed based on ζ = 0.03, the results are deemed suitable for the other low-damping situations. A portion of the numerical results are plotted in Fig. 3 to show the effects of R and α on the δ-T spectrum. In Fig. 3a, the δ-T spectra under different values of R are plotted (with α = 0.1). Each spectrum is the average of the δ-T curves under the many earthquake excitations from the same record collection. As expected, a larger R, which means a severer nonlinearity, gives a more pronounced reduction in the PFA (compared with the spectral acceleration S a ). In Fig. 3b, the δ-T spectra with different α (under R = 2.5) are plotted. As the value of α increases, the nonlinear reduction in PFA decreases, and the value of δ gets larger. Compared with the effect of R, the effect of α on δ is less dramatic. Nonetheless, one still needs to consider the α-effect on δ because such an effect is not ignorable. For example, for a steel dual system (with a fundamental period T 1 equals to, say, 0.5 s) showing an  Fig. 3b-2. If δ = 0.45 (corresponding to α = 0.1) is used to compute the PFA, there would be an error as large as 30%. In general, there is a notable decrease in the δ-T curve at T ≤ 0.3 s, and a smooth increase in the value of δ after T surpasses the short period range, i.e., T ≥ 0.3 s in this study. For an elastic system with a very short period, its PFA is almost equal to the PGA, once it starts yielding (F e becomes larger than F y ), it would exhibit a prominent ductility demand because its elastic spectral displacement is tiny. Such a strong nonlinear behavior indicates a notable reduction in its acceleration response. On the other hand, for a system with a very long period, e.g., T approaching to ∞, its PFA would be near zero, no matter what its degree of nonlinearity is. In this case, the value of δ should get close to unit. Given such an inference, it is reasonable to obtain a smooth increase in δ as T gets longer because the value of δ would get closer to 1.0 as T increases.
As a supplementary result, the inelastic spectra (in the form of μ-T curves) are shown in Fig. 3c-1 and Fig. 3c-2, corresponding to the two groups of earthquake excitation records, respectively. In specific, the μ-T curves are computed with α = 0.1. For the other values of Comparing Fig. 3a with Fig. 3b, the δ-T curves under the two earthquake hazard scenarios are different. This proves the dependency of the δ-R-T-α relationship on the seismological property. In this regard, the authors of this paper are not proposing any specific δ-R-T-α formulae because such proposals can hardly be universal. Instead, only the feasibility of a δ-R-T-α formula for the practical application is validated in this study (see the following Sect. 2.3). Based on the validation, one can develop his/her own δ-R-T-α formula for any hazard scenarios or seismological situations to be encountered.
Additionally, it is worth noting that the prediction of δ is a little bit more complicated than the prediction of μ (ductility demand). This is because the α is explicitly incorporated in Eq. (1), while for μ, generally the effect of α can be ignored and an R-μ-T relationship is sufficient, especially in the medium-to-long period range where the peak inelastic displacement follows the "equal displacement rule" (Veletsos and Newmark 1960).

Feasibility of a δ-R-T-α model
In this section, the predictability of δ is compared with that of μ (ductility ratio), and the feasibility of a δ-R-T-α relationship is validated. In Fig. 4, the coefficients of variation (CoV) of δ and μ are plotted for the earthquake record collection A. Specifically, shown in Fig. 4a-c are the CoV of δ and μ (denoted by δ CoV and μ CoV , respectively) under R = 2.0, 3.0, and 4.0 with α = 0.1, while in Fig. 4d-f are those under α = 0.10, 0.25, and 0.40 with R = 2.5. From Fig. 4, the δ CoV is smaller than the μ CoV over the entire period range (0 ≤ T ≤ 3.0 s) for all the values of R and α. Thus, the δ exhibits a smaller dispersion than μ under the earthquake uncertainty. From this observation, δ is more predictable than μ. Therefore, it is feasible to predict the value of δ using a relevant δ-R-T-α relationship. Similar with Fig. 4, the results for earthquake record collection B are plotted in Fig. 5. From Fig. 5, the value of δ CoV is smaller than μ CoV in most situations. Some exceptions occurred in Fig. 5e and f whereas δ CoV is slightly larger than μ CoV at T ≥ 2.0 s. In such cases, the predictability of δ and μ are basically alike. Given the results shown in Figs. 4 and 5, the δ exhibits a better or similar predictability with μ. Since the μ can be estimated by a R-μ-T relationship, one can certainly establish a δ-R-T relationship to Fig. 4 δ CoV -T and μ CoV -T curves for earthquake record collection A estimate the value of δ. In special, as mentioned previously, the effect of α needs to be considered in predicting δ, thus the established δ-R-T relationship is in fact a δ-R-T-α formula. For completeness, a δ-R-T-α formula obtained by fitting the δ-R-T-α curves generated by realistic earthquake records is given in the Appendix. Note that such a δ-R-T-α formula is an empirical δ-R-T-α relationship, which reveals the general trend of the δ-R-T-α formulations to be developed in the future for many different seismic scenarios.

PFA of MDF system
Given a δ-R-T-α formula, one can predict the δ for the nonlinear SDF system based on the definition of δ, as shown in Eq. (1). For the inelastic MDF systems, one can use a nonlinear modal combination process to estimate the PFA. In doing so, the elastic PFA associated with each mode is revised based on an appropriate δ, while the value of δ is obtained from the δ-R-T-α relationship. Following this process, an inelastic CQC (ICQC) approach is proposed in this section. Since the ICQC approach is derived based on the elastic CQC approach, the elastic CQC approach for PFA is briefly revisited before the ICQC approach is introduced.

Revisit the CQC approach for PFA of elastic MDF systems
For an elastic MDF structure, the dynamic equilibrium is: where M, C, and K is the mass, damping, and stiffness matrix, respectively, , ̇ , and ̈ is the structural displacement, velocity, and acceleration relative to the ground, ι is an influence vector (Chopra 2011), and a g means ground acceleration. The relative acceleration response ̈ can be expressed in the modal space, as where n is the number of degree-of-freedom, Γ i and ϕ i is the ith modal participation factor and the ith mode shape vector, respectively, s i is the acceleration response of the ith modal oscillator, which is controlled by: where s i and ̇s i are the modal displacement and velocity of the ith mode, ζ i and ω i is the damping ratio and circular frequency of the ith mode, respectively.
Note that ̈ is the relative acceleration responses, the total (absolute) acceleration ̈ t of the structure should be calculated by incorporating the ground acceleration: where s t i is the total modal acceleration of the ith mode. From Eq. (5), the peak total acceleration can be predicted by combining the Γ i is t i of each mode. According to Pozzi and Der Kiureghian (2015) and Chopra (2011), the peak total acceleration s t i of an oscillator can be approximated by the pseudo-acceleration: where S a (T i , ζ i ) is the ζ i -damped pseudo-acceleration at T i . Therefore, the peak Γ i is t i can be replaced by Γ i i S a T i , i in the combination. Based on this approximation, the elastic CQC approach is derived as ( In Eq. (7), PFA k is the PFA of the kth story in a n-story structure, ρ ij is the cross-modal correlation coefficient between the ith and the jth modes, ϕ i,k and ϕ j,k is the kth element of ϕ i and ϕ j , respectively, Γ j is the modal participation factor of the jth mode, while S a (T j , ζ j ) is the ζ j -damped pseudo-acceleration at T j . In this study, the value of ρ ij is computed by: where β ij = T j /T i . Note that Eq. (8) is from Kiureghian (1981) and is widely used nowadays.

ICQC approach for PFA of inelastic MDF system
By modifying the elastic modal PFAs and combining them via the CQC rule, the ICQC approach can be established. In doing so, the value of Γ i is t i in Eq. (5) is revised firstly. In the revision, it is assumed that: (1) the mode shape Γ i i of the inelastic MDF structure remains the same with its elastic counterpart; (2) the amplitude of s t i is mitigated due to the structural nonlinearity. To be clear, assumption (1) is not theoretically rigorous because the mode shapes of an inelastic system differ from its elastic ones. Given the complexities in computing the nonlinear modes, using the elastic mode shapes in the subsequent process is simple and straightforward. The mitigated s t i is computed by multiplying the elastic S a with the δ defined by a proper δ-R-T-α relationship. By doing so, the contribution of the ith mode to the total acceleration becomes δ of δ for the ith modal oscillator, R i , T i , and α i is the strength reduction factor, period, and post-yield stiffness ratio of the ith mode, respectively. Notably, the values of R i and α i can be obtained via a modal pushover analysis procedure (Chopra and Goel 2002). The above processing leads us to the ICQC approach for the inelastic PFA, which is written as: In Eq. (9), δ j (R j ,T j ,α j ) is the value of δ for the jth modal oscillator, R j , T j , and α j are the strength reduction factor, period, and post-yield stiffness ratio of the jth modal oscillator. Comparing Eq. (9) with Eq. (7), the spectral amplitude S a (T i , ζ i ) in Eq. (7) is replaced by δ(R i ,T i ,α i )S a (T i , ζ i ) in Eq. (9) to account for the effect of structural nonlinearity, while the mode shape, the modal participation factor, and the cross-modal correlation coefficient remain unchanged. In Eq. (9), the values of δ i (R i ,T i ,α i ) and δ j (R j ,T j ,α j ) are to be determined via an appropriate δ-R-T-α relationship.

Illustration via an example of a 2-story structure
An example is presented in this section to illustrate the application of the ICQC approach. The structure has 2 stories, and the mass lumped on each story is m = 100 t. The interstory stiffness ratio between the 1st and the 2nd story is k 1 : k 2 = 3: 2, and the fundamental period of the structure is T 1 = 0.2 s. Damping of the structure is simulated by the Rayleigh damping model with a modal damping ratio of ζ 1 = ζ 2 = 0.03. The structure is subjected to an earthquake record shown in Fig. 6a. The pseudo-acceleration spectrum of the excitation is plotted in Fig. 6b. For this excitation, there is S a (T 1 ) = 1.283 g and S a (T 2 ) = 0.402 g at T 1 = 0.2 s and T 2 = 0.082 s, respectively. The mode shapes (Γ i ϕ i,k ) of the structure are shown in Fig. 6c.

Modal pushover analysis
The degree of structural nonlinearity is quantified via a modal pushover analysis. There are several variants of the modal pushover procedure in the literature, in this study, the one proposed by Xiang et al. (2017Xiang et al. ( , 2018) is adopted. Subjecting the structure to the static modal pushover load MS a (T 1 ,ζ 1 )Γ 1 ϕ 1 and MS a (T 2 ,ζ 2 )Γ 2 ϕ 2 for mode 1 and mode 2, respectively, the modal pushover curves are obtained (by plotting the equivalent pushover force F eq against the equivalent modal displacement d eq ), as shown in Fig. 6d and e. From Fig. 6d, the 1st mode exhibits an obvious nonlinearity as its modal oscillator shows a strength reduction factor of R 1 = 1269 kN / 571 kN = 2.22 and a post yield stiffness ratio of α 1 = 0.10. From Fig. 6e, the 2nd mode remains elastic.

Determine the PFA spectral ratio δ
According to R i , α i , and T i , the PFA spectral ratio δ i for each modal oscillator can be obtained. In the engineering practice, one can use a δ-R-T-α relationship to determine the value of δ i (R i ,T i ,α i ). While in this example, the THA of the accordant nonlinear SDF system is conducted to get the value of δ i . For the 1st mode, an oscillator with R 1 = 2.22, α 1 = 0.10, and T 1 = 0.2 s is subjected to the earthquake record plotted in Fig. 6a. The total acceleration Fig. 6f. From Fig. 6f, the peak absolute value of s t is 0.611 g. Dividing this value by S a (T 1 ) = 1.283 g, it gives δ 1 = 0.476. Analogously, for the 2nd mode, there is δ 2 = 1.0 because its modal oscillator keeps elastic throughout the response history.

Calculate the PFA via ICQC
Substituting δ 1 = 0.476, S a (T 1 ) = 1.283 g, δ 2 = 1.0, and S a (T 2 ) = 0.402 g in Eq. (9), the PFAs predicted by the ICQC approach are obtained. For story 1 and story 2, there is PFA 1 = 0.401 g and PFA 2 = 0.732 g, respectively. A comparison between the ICQC results and the exact results (obtained from nonlinear THA) are given in Fig. 6g and (h). As expected, the ICQC results are quite accurate (errors are -10.7% and 1.8% for story 1 and story 2).

Validation of the ICQC approach
In this section, the ICQC approach is further validated via extensive numerical computations. In addition to the 2-story structure adopted above, a 5-and an 8-story structure are subjected to 300 randomly-selected earthquake records. The record collection consists of three parts, each includes 100 records under a specific site class (site class B, C, and D according to the Vs30-based category in ASCE 7-10 (2010), where Vs30 is the timeaveraged shear wave velocity of the top 30 m of the soil). These records were downloaded from the Japan NIED KiK-net. The Vs30 of each site, the magnitude of earthquake, and the PGA of each record are marked in Fig. 7. The 3%-damped pseudo acceleration spectra of the records are plotted in the figure as well.

Structural models
The fundamental period of the 2-, 5-, and 8-story model is T 1 = 0.2 s, 0.5 s, and 0.8 s, respectively. The lumped masses on all the stories are identical (denoted by m hereafter). The structures are designed following three steps: (1) The elastic seismic force demand is determined via the equivalent lateral force procedure (ASCE 7-10, 2010), whereas the base shear demand is calculated as F b = nmS r a T 1 , 1 . Here n equals to 2, 5, and 8 for the three models, respectively, S r a T 1 , 1 is the ζ 1 -damped spectral acceleration at T 1 for the rth earthquake record.
(2) The inter-story drift vs. story shear relationship is bilinear (representing ductile steel structures or structures equipped with metallic dampers). The yield strength F yi of the ith story is calculated as: dividing the elastic seismic force demand by a design strength reduction factor R 0 . For example, for the 1st story, there is F y1 = F b /R 0 , as detailed in Fig. 8. (3) The inter-story stiffness k i is in proportion to the inter-story shear force determined by the equivalent lateral force procedure. For clarity, the height-wise distributions of k i are given in Fig. 8 as well.
For each model, the post-yield stiffness ratios of all stories (denoted by α 0 hereafter) are identical. In the computation, several values of R 0 (from 1.5 to 4.0) and α 0 (from 0.1 to 0.5) are considered.
The natural modes (Γ i ϕ i,k ) of the 2-story structure have been shown in Fig. 6c, and the accordant natural periods are T 1 = 0.2 s and T 2 = 0.082 s. The natural modes and periods of the 5-and 8-story structures are shown in Fig. 9a and b. Damping of the structures are simulated by the Rayleigh damping model. A modal damping ratio of 0.03 is assigned to the 1st and 2nd modes, the 1st and 3rd modes, and the 1st and 4th modes for the 2-, 5-, and 8-story structure, respectively. For clarity, the values of ζ i for the 5-and 8-story structures are marked in Fig. 9c. In the ICQC computation, the modal damping ratios shown in Fig. 9c are used in determining S a (T i , ζ i ). This is because the value of damping significantly affects the value of S a (T i , ζ i ) in the low-damping range (Xiang and Huang 2019), and such an effect may further influence the value of PFA in a notable way (Pinkawa et al. 2014).

PFA computation
Firstly, the structural PFA is computed using the ICQC approach. In the ICQC, the modal pushover analyses are conducted to obtain the values of R i and α i for the modal oscillators. Usually, the pushover curve for a MDF structure may exhibit a multistage yielding process. In such cases, the multi-linear pushover curve is equivalently transformed to a bilinear curve based on the "equal energy" criterion demonstrated in Fig. 10, and the values of R i and α i for the modal oscillator are determined from the bilinearized curve.
Based on R i , T i , and α i , the value of δ i (R i ,T i ,α i ) for each modal oscillator is determined. In the engineering practice, an δ-R-T-α model can be used for computing δ i . In this section, however, the THA on the nonlinear modal oscillator is conducted to obtain the value of δ i under each earthquake excitation (denoting δ i by r i for the rth excitation). By doing so, the fitting errors of the δ-R-T-α model are fully avoided. Putting the values of r i and S a (T i , ζ i ) in Eq. (9), the nonlinear PFA is obtained. For clarity, it is denoted by PFA ICQC hereafter. On the other hand, the PFAs are computed by the directly conducting nonlinear THA on the MDF structures, and the results are denoted by PFA THA .

Accuracy of the ICQC approach
The accuracy of ICQC is checked via a parameter R ICQC THA defined as: Shown in Fig. 11 are the values of R ICQC THA for the 2-, 5-, 8-story structures under all the earthquake excitations. The data in this figure corresponds to R 0 = 2.5 and α 0 = 0.1. In the figure, the results for the different site classes are plotted separately. Each individual data point represents a value of R ICQC THA under a specific earthquake excitation. The mean (ave.) and the standard deviation (std.) of R ICQC THA on each site class are marked as well. From Fig. 11, the value of R ICQC THA is generally within 0.5 ~ 1.5. For the 2-story structure, the mean of R ICQC THA is almost identical to 1.0 on site class B and C, as shown in Fig. 11a-1 and b-1, respectively. In Fig. 11c-1, the mean of R ICQC THA is less than 1.0 for the 1st story on site class D, showing a roughly 20% underestimation of PFA by the ICQC approach. For the 5-story structure, the ICQC approach generates excellent results for the 2nd, 3rd, and 4th stories on all site classes, as shown in Fig. 11a-2, b-2, and c-2. For the 1st story, the value of R ICQC THA decreases as the site becomes softer (Vs30 gets smaller). This may be induced by the varied modal contributions to the PFA on the different site conditions. For the 5th story, there is about a 30% over-estimation of PFA by the ICQC approach on all site classes. The results of the 8-story structure are quite similar with those of the 5-story model, as shown in Fig. 11a-3, b-3, and c-3. Generally, the dispersion of R ICQC THA decreases on a softer site (with smaller Vs30). In all, despite of the moderate errors in some situations, i.e., for specific stories on certain sites, the PFA ICQC is close to its accurate counterpart PFA THA . This demonstrates the good accuracy of the ICQC approach. Note that the data in Fig. 11 corresponds to R 0 = 2.5 and α 0 = 0.1. To be comprehensive, some other R 0 -α 0 combinations are presented as well. Shown in Fig. 12 are the outcomes of the 5-story structure with R 0 = 1.5 and 4.0, and α 0 = 0.25 and 0.40. Comparing the data in Fig. 11a-2, b-2, and c-2 with that in Fig. 12, the mean of R ICQC THA does not change much when the values of R 0 and α 0 alter. Although not very pronounced, the dispersion of R ICQC THA increases as R 0 gets larger and α 0 gets smaller. This indicates an enlarged uncertainty of the accuracy of ICQC when the structure exhibits a severer nonlinearity. In general, from the data in Fig. 12, the ICQC approach performs well under the different R 0 -α 0 combinations. The results in Figs. 11 and 12 demonstrate that the ICQC approach is suitable for application.

Effect of nonlinearity of higher-order modes
It is worthwhile to check the effect of nonlinearity of the higher-order modes on the PFA. In some preceding studies (Rodriguez et al. 2002;Moschen et al. 2013), the higher-order modes are (implicitly) assumed elastic, and the nonlinear reduction in PFA (or total acceleration s t i ) is ignored. Consequently, some unacceptable over-estimations of PFA are generated. While in ICQC, the nonlinearities of all the modes are considered by incorporating δ(R i ,T i ,α i ) in Eq. (9), and this could lead to a better estimation of the nonlinear PFA.
For demonstration, the degree of modal nonlinearity for the 5-story structure with R 0 = 2.5 and α 0 = 0.1 are presented in Fig. 13. The strength reduction factor R i of each modal oscillator is marked. Here the values of R i are obtained from the modal pushover curves, while these curves are derived by subjecting the structure to the modal pushover force S r a T i , i Γ i i . For example, for the 1st mode, R 1 is determined under the pushover force S r a T 1 , 1 Γ 1 1 . Since the design strength reduction factor R 0 corresponds to the seismic force F b = 5mS r a T 1 , 1 , there is a stable relationship between R 0 (global measure of structural strength) and R 1 (measure of the 1st modal oscillator) because both parameters are associated with S a (T 1 ,ζ 1 ). In this example, in accordance with R 0 = 2.5, there is R 1 = 2.05. For the other modes (i ≥ 2), the value of R i varies as S a (T i ,ζ i )/S a (T 1 ,ζ 1 ) changes. From Fig. 13, the nonlinearities of the higher-order modes become severer as the site gets firmer. This is because the value of S a (T i ,ζ i )/S a (T 1 ,ζ 1 ) (i ≥ 2) becomes larger on a firm site. For example, the predominant period of the response spectra on site class B is within 0.1 s < T < 0.2 s, which corresponds to the natural period of the 2nd mode. This makes the value of S a (T 2 ,ζ 2 )/S a (T 1 ,ζ 1 ) larger than 5.0 in most cases, leading to a severer 2nd modal nonlinearity. Consequently, the value of R 2 is mostly larger than R 1 in this situation, and the contribution of the 2nd mode to the PFA can be significantly mitigated by the modal nonlinearity.
The effect of the higher-order modal nonlinearity on the PFA is quantified in this section. In doing so, the nonlinear PFA is calculated via a simplified version of ICQC (denoted by ICQC# hereafter). In the ICQC#, the higher-order modes are assumed to be elastic and the values of δ i (R i ,T i ,α i ) and δ j (R j ,T j ,α j ) in Eq. (11) are set to 1.0 at i ≥ 2 or j ≥ 2. Denoting the accordant results by PFA ICQC# , the accuracy of ICQC# is measured by: For consistency, the results of the 5-story model with R 0 = 2.5 and α 0 = 0.1 are presented and discussed. The values of R ICQC# THA and R ICQC THA are marked in Fig. 14a and b, respectively. From Fig. 14a, the value of R ICQC# THA is close to 1.0 on site class D, yet it is significantly larger than 1.0 on site class B. This shows that the higher-order modal nonlinearity is quite influential on the firm sites. Specially, the error of ICQC# for the 5th story is un-acceptably severe on site class B. Also, there is a large dispersion of R ICQC# THA , which means a notable uncertainty in the accuracy of PFA ICQC# .
Comparatively, the results of ICQC are much better. From Fig. 14b, the PFA given by the ICQC approach, which considers the higher-order modal nonlinearities, is generally identical with the THA results. The site-dependency of R ICQC THA is much less than that of R ICQC# THA . In addition, the dispersion of the results is well mitigated in ICQC. From Fig. 14, it is needed to incorporate the higher-order modal nonlinearities in computing the PFA. Therefore, the ICQC approach formulated in Eq. (9) is deemed better than those approaches that ignore such an issue, e.g., the FMR method in Rodriguez et al. (2002).

Summary and conclusion
Proposed and validated in this paper is a nonlinear spectral approach for the earthquakeinduced inelastic PFA of the MDF bilinear hysteretic systems. Such systems represent the multi-story buildings structured with ductile steel frames. The main findings are: (1) For a bilinear hysteretic single-degree-of-freedom system, a PFA spectral ratio δ could be adopted to quantify the nonlinear reduction in PFA. The δ is defined by δ = PFA NL (R, α, T)/S a (T), wherein PFA NL (R, α, T) is the nonlinear PFA of a bilinear hysteretic oscillator with strength reduction factor R, post-yield stiffness ratio α, and vibrating period T, while S a (T) is the elastic spectral acceleration at period T. Numerical computations reveal that there is a dependable relationship between δ, R, α, and T. The predicted δ based on the δ-R-T-α relationship exhibits a small dispersion under the earthquake uncertainty. Thus, it is feasible to estimate the value of δ via a δ-R-T-α relationship for the SDF systems.
(2) For a bilinear hysteretic multi-degree-of-freedom system, a nonlinear spectral approach (inelastic CQC approach) is proposed to estimate its inelastic PFA. In the ICQC approach, the contribution of each mode to the elastic PFA is modified to account for the nonlinear effect. The modification could be conducted based on a proper δ-R-T-α relationship, while the values of R and α are determined by a modal pushover analysis procedure. The modified modal contributions to the PFA are then combined via the complete-quadratic combination rule. The good accuracy of the inelastic CQC approach is validated by the extensive numerical computations involving some 2-, 5-and 8-story structural models with different degrees of nonlinearity.
(3) The effect of nonlinearities of the higher-order modes on the PFA can be significant, especially on the firm sites. The inelastic CQC approach can appropriately account for such effects, thus it avoids over-estimations of the inelastic PFA. This is an advantage of the inelastic CQC approach over the conventional methods.
The δ-R-T-α relationship defined by Eq. (12) and (13) exhibits a good accuracy. For demonstration, the mean values of δ from the nonlinear RHA are plotted against those given by Eq. (12) and (13), as shown in Fig. 15. The RHA results are marked by black dots, while the simulations are shown using surfaces in blue-brown. For brevity, only the results at T = 0.1 s, 0.2 s, 0.5 s, 1.0 s, 1.5 s, and 2.0 s are given here. From Fig. 15, the simulation well captures the RHA results for T ≥ 0.2 s, while for the very short period (T ≤ 0.1 s), some minor discrepancies can be found. Lastly, it is emphasized that the δ-R-T-α formula given here is an empirical δ-R-T-α relationship established based on limited earthquake