Performance Assessment of an Energy–Based Approximation Method for the Dynamic Capacity of RC Frames Subjected to Sudden Column Removal Scenarios

: The alternative load path method is widely used to assess the progressive collapse performance of reinforced concrete structures. As an alternative to an accurate non–linear dynamic analysis, an energy–based method (EBM) can also be adopted to approximately calculate the dynamic load–bearing capacity curve or the dynamic resistance based on a static capacity curve. However, dynamic effects cannot be explicitly taken into account in the EBM. The model uncertainty associated with the use of the EBM for evaluating the dynamic ultimate capacity of structural frames has not yet been quantiﬁed. Knowledge of this model uncertainty is however necessary when applying EBM as part of reliability calculations, for example, in relation to structural robustness quantiﬁcation. Hence, this article focuses on the evaluation of the performance of the EBM and the quantiﬁcation of its model uncertainty in the context of reliability–based assessments of progressive or disproportionate collapse. The inﬂuences of damping effects and different column removal scenarios are investigated. As a result, it is found that damping effects have a limited inﬂuence on the performance of the EBM. In the case of an external column removal scenario, the performance of the EBM is lower as the response is not a single deformation mode according to the results in the frequency domain. However, a good performance is found in the case of an internal column removal scenario in which the assumption of a single deformation mode is found to be sufﬁciently adequate. Probabilistic models for the model uncertainties related to the use of the EBM compared to direct dynamic analyses are proposed in relation to both the resistances and the associated displacements. Overall, the EBM shows to be an adequate approximation, resulting in a small bias and small standard deviation for its associated model uncertainty.


Introduction
Progressive collapse triggered by extreme events is often accompanied by catastrophic losses in terms of human lives and the economy, e.g., as experienced in the 1995 Alfred P. Murrah Federal Building collapse [1][2][3][4][5]. Therefore, the robust design of structures has gained importance during the last decades [6,7]. Both in experimental tests [8][9][10][11] and some design codes, e.g., DoD [12], the threat-independent alternative load path (ALP) method is widely used to evaluate the performance of reinforced concrete (RC) building structures against disproportionate or progressive collapse through the notional removal of one or more load-bearing structural members.
Considering that tests in this regard are very time-consuming and costly, numerical simulation approaches are usually adopted to study the influences of different parameters on the performance against progressive collapse in detail [6,[13][14][15]. Generally, dynamic effects should be considered in the ALP method. A dynamic non-linear time history analysis (NTHA) is currently the most accurate procedure to simulate the sudden column the literature [19]. The model uncertainty information of the EBM is crucial to adequately quantify the failure probability of a system when the EBM is adopted.
In light of the above-mentioned state-of-the-art, this paper hence focuses on evaluating the performance of the EBM and quantifying the model uncertainty in relation to the use of the EBM in case of the assessment of RC frames subjected to different sudden column removals. The layout of the paper and research approach is as follows, see Figure 1. Section 2 presents the methodology of the EBM and the procedures on how to quantify the model uncertainty of the EBM is illustrated in Section 3. The numerical modelling techniques are validated based on an experimental result in Section 4. In Section 5, the EBM is then validated by comparing with the results of the NTHA and both the influences of damping effects and different column removal scenarios are studied. Furthermore, stochastic analyses are carried out and model uncertainties to be considered when using the EBM are calculated in Sections 6 and 7, respectively. Concluding remarks are addressed in Section 8.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 23 quantitative assessment of the model uncertainty associated to the EBM becomes important when the EBM is applied to quantify the reliability or robustness of a RC building structure following a sudden column removal scenario. Such information regarding RC frames is currently lacking in the literature [19]. The model uncertainty information of the EBM is crucial to adequately quantify the failure probability of a system when the EBM is adopted.
In light of the above-mentioned state-of-the-art, this paper hence focuses on evaluating the performance of the EBM and quantifying the model uncertainty in relation to the use of the EBM in case of the assessment of RC frames subjected to different sudden column removals. The layout of the paper and research approach is as follows, see Figure  1. Section 2 presents the methodology of the EBM and the procedures on how to quantify the model uncertainty of the EBM is illustrated in Section 3. The numerical modelling techniques are validated based on an experimental result in Section 4. In Section 5, the EBM is then validated by comparing with the results of the NTHA and both the influences of damping effects and different column removal scenarios are studied. Furthermore, stochastic analyses are carried out and model uncertainties to be considered when using the EBM are calculated in Sections 6 and 7, respectively. Concluding remarks are addressed in Section 8.

Energy-Based Method
The total energy is considered constant throughout the dynamic deformation of a structure ( Figure 2a) following a sudden column removal scenario, conditional on the assumption that the energy dissipated by other sources, e.g., heat, are neglected [18,19,23]. Three components of energy are taken into account when balancing the energy: kinetic energy (KE) originating from the moving sub-structure, strain energy (SE) due to the structural deformation, and release of potential energy (PE) due to the downward movement of the structure. The released potential energy resulting from the unbalanced gravity loads for the remaining structural system leads to dynamic motions and an increase of kinetic energy. In contrast, elastic and inelastic strain energy accumulated by the deformed structure dissipate the kinetic energy, counteracting the motions. If a structure can redistribute the unbalanced loads, the potential energy can be absorbed by the strain energy. Consequently, the kinetic energy can attain a value of zero (i.e., KE = PE -SE = 0).
Based on the aforementioned energy balance during a typical column removal scenario, the EBM can be easily interpreted and applied. Generally, the strain energy stored in the deformed system can be calculated through a non-linear static pushdown analysis [18]. In that case, the strain energy is the hatched area under the static capacity curve up to a displacement, namely us (Figure 2b). On the other hand, the released potential energy generated by the unbalanced loads following a sudden column removal is equal to the hatched rectangular area for a constant load level, namely Qd. If the kinetic energy is assumed to be zero (at the displacement level, us = ud), the potential energy is therefore equal to the strain energy (PE = SE). Note that only the situation in which no collapse

Energy-Based Method
The total energy is considered constant throughout the dynamic deformation of a structure ( Figure 2a) following a sudden column removal scenario, conditional on the assumption that the energy dissipated by other sources, e.g., heat, are neglected [18,19,23]. Three components of energy are taken into account when balancing the energy: kinetic energy (KE) originating from the moving sub-structure, strain energy (SE) due to the structural deformation, and release of potential energy (PE) due to the downward movement of the structure. The released potential energy resulting from the unbalanced gravity loads for the remaining structural system leads to dynamic motions and an increase of kinetic energy. In contrast, elastic and inelastic strain energy accumulated by the deformed structure dissipate the kinetic energy, counteracting the motions. If a structure can redistribute the unbalanced loads, the potential energy can be absorbed by the strain energy. Consequently, the kinetic energy can attain a value of zero (i.e., KE = PE − SE = 0).
Based on the aforementioned energy balance during a typical column removal scenario, the EBM can be easily interpreted and applied. Generally, the strain energy stored in the deformed system can be calculated through a non-linear static pushdown analysis [18]. In that case, the strain energy is the hatched area under the static capacity curve up to a displacement, namely u s (Figure 2b). On the other hand, the released potential energy generated by the unbalanced loads following a sudden column removal is equal to the hatched rectangular area for a constant load level, namely Q d . If the kinetic energy is assumed to be zero (at the displacement level, u s = u d ), the potential energy is therefore equal to the strain energy (PE = SE). Note that only the situation in which no collapse occurs is considered. Eventually, the predicted dynamic load Q d can be calculated as follows [18,19]: where Q d is the load in the dynamic capacity curve, Q s is the load in the static capacity curve, u d is the peak dynamic deflection, and u s is the static deflection corresponding to the load Q s .
occurs is considered. Eventually, the predicted dynamic load Qd can be calculated as follows [18,19]: where Qd is the load in the dynamic capacity curve, Qs is the load in the static capacity curve, ud is the peak dynamic deflection, and us is the static deflection corresponding to the load Qs. The EBM approach is a compromise between accuracy and complexity. The following simplifications are incorporated: • The moving sub-structure, subjected to a column failure, is assumed to behave like a single degree of freedom (SDOF) system [17,18,23]. The response is controlled by a single deformation mode and the mode keeps constant during the dynamic response. Therefore, the energy of the whole system can be linked to the energy of a point, i.e., every point in the system reaches its maximum displacement response at a same time. However, this never happens in a real structure as the existing infinite number of deformation modes will reach their maximum response at different moments. Consequently, the stored strain energy counted by the EBM is overestimated alongside its calculated deflections [23]. For a structural frame system subjected to an exterior column removal scenario, a non-single deformation response may occur due to complex load redistribution mechanisms (e.g., overturn), such as the results by [4,13,46] in relation to the exterior or side column removal scenarios.
• All the energy introduced into a system by the loads is switched into pure strain energy. The EBM neglects the energy dissipated by damping or other mechanisms. Therefore, the maximum deflection response will be overestimated. Moreover, it is still a controversial issue regarding how to model the damping mechanism for a sudden column removal scenario [19,47,48], e.g., viscous damping or Coulomb damping [23]. For instance, the studies have identified drawbacks associated with the use of Rayleigh damping based on initial stiffness, in which the initial stiffness may introduce in the unwanted artificial viscous forces [14,49]. This is because the stiffness term involved in the Rayleigh damping should also be updated accordingly when a system responds in the inelastic stage [48].
• The strain energy storage capacities of a system for a given displacement in static and dynamic situations are different [50]. The EBM cannot take this into account. However, the influence of the strain rate is low according to both previous experimental results [10] and numerical studies [19,23], as the maximum strain rates occur only in a small area and during a short time duration for sudden column removals.

Quantification of the Model Uncertainty Associated to the Application of EBM
As discussed before, an approximate result is obtained using the EBM. The EBM curve completely depends on the static pushdown curve (see Figure 2b). Considering it is The EBM approach is a compromise between accuracy and complexity. The following simplifications are incorporated:

•
The moving sub-structure, subjected to a column failure, is assumed to behave like a single degree of freedom (SDOF) system [17,18,23]. The response is controlled by a single deformation mode and the mode keeps constant during the dynamic response. Therefore, the energy of the whole system can be linked to the energy of a point, i.e., every point in the system reaches its maximum displacement response at a same time. However, this never happens in a real structure as the existing infinite number of deformation modes will reach their maximum response at different moments. Consequently, the stored strain energy counted by the EBM is overestimated alongside its calculated deflections [23]. For a structural frame system subjected to an exterior column removal scenario, a non-single deformation response may occur due to complex load redistribution mechanisms (e.g., overturn), such as the results by [4,13,46] in relation to the exterior or side column removal scenarios. • All the energy introduced into a system by the loads is switched into pure strain energy. The EBM neglects the energy dissipated by damping or other mechanisms. Therefore, the maximum deflection response will be overestimated. Moreover, it is still a controversial issue regarding how to model the damping mechanism for a sudden column removal scenario [19,47,48], e.g., viscous damping or Coulomb damping [23]. For instance, the studies have identified drawbacks associated with the use of Rayleigh damping based on initial stiffness, in which the initial stiffness may introduce in the unwanted artificial viscous forces [14,49]. This is because the stiffness term involved in the Rayleigh damping should also be updated accordingly when a system responds in the inelastic stage [48].

•
The strain energy storage capacities of a system for a given displacement in static and dynamic situations are different [50]. The EBM cannot take this into account. However, the influence of the strain rate is low according to both previous experimental results [10] and numerical studies [19,23], as the maximum strain rates occur only in a small area and during a short time duration for sudden column removals.

Quantification of the Model Uncertainty Associated to the Application of EBM
As discussed before, an approximate result is obtained using the EBM. The EBM curve completely depends on the static pushdown curve (see Figure 2b). Considering it is an approximate approach, it is therefore important to quantitatively assess the performance of the EBM, i.e., quantifying its model uncertainty through comparison to the more accurate dynamic analysis results. To evaluate this model uncertainty, the following three major steps are applied (see also the flowchart for the model uncertainty quantification of the EBM shown in Figure 3): (a) Initially, the appropriate random variables and a column removal scenario are selected. (b) Subsequently, both the non-linear static and dynamic analyses are carried out. In case of the former, the pushdown curve is subsequently used to derive the EBM curve and the corresponding dynamic resistance for every realization (the left branch in Figure 3). In case of the latter, incremental dynamic analyses (IDA) are executed to accurately determine the dynamic resistance for every realization (the right branch in Figure 3). (c) Eventually, the model deviation for the EBM comparing to the IDA are calculated for both the resistances and the corresponding displacements: where K EBM is the ratio of EBM/IDA. R EBM is the ultimate load-bearing capacity according to the EBM, while R IDA is the ultimate capacity according to the IDA. D EBM is the maximum displacement calculated according to the EBM, while D IDA is the maximum displacement calculated from the IDA.
an approximate approach, it is therefore important to quantitatively assess the performance of the EBM, i.e., quantifying its model uncertainty through comparison to the more accurate dynamic analysis results. To evaluate this model uncertainty, the following three major steps are applied (see also the flowchart for the model uncertainty quantification of the EBM shown in Figure 3): a) Initially, the appropriate random variables and a column removal scenario are selected. b) Subsequently, both the non-linear static and dynamic analyses are carried out. In case of the former, the pushdown curve is subsequently used to derive the EBM curve and the corresponding dynamic resistance for every realization (the left branch in Figure 3). In case of the latter, incremental dynamic analyses (IDA) are executed to accurately determine the dynamic resistance for every realization (the right branch in Figure 3). c) Eventually, the model deviation for the EBM comparing to the IDA are calculated for both the resistances and the corresponding displacements: where KEBM is the ratio of EBM/IDA. REBM is the ultimate load-bearing capacity according to the EBM, while RIDA is the ultimate capacity according to the IDA. DEBM is the maximum displacement calculated according to the EBM, while DIDA is the maximum displacement calculated from the IDA.

Finite Element Modelling
The modelling strategy for RC structures using OpenSees software package [51] is first introduced. Regarding the beam and column components, the force-based fiber beam-column element in OpenSees is used. Shear deformation and shear failure are not accounted for as the element is based on the Euler-Bernoulli beam theory. Nevertheless, an accurate prediction of the responses for slender flexure-dominated (large span-todepth ratio) elements is possible, as the shear deformation is not significant in this situation [52,53]. The cross-section of the fiber element is discretized into several fibers which

Finite Element Modelling
The modelling strategy for RC structures using OpenSees software package [51] is first introduced. Regarding the beam and column components, the force-based fiber beamcolumn element in OpenSees is used. Shear deformation and shear failure are not accounted for as the element is based on the Euler-Bernoulli beam theory. Nevertheless, an accurate prediction of the responses for slender flexure-dominated (large span-to-depth ratio) elements is possible, as the shear deformation is not significant in this situation [52,53]. The cross-section of the fiber element is discretized into several fibers which are subjected to a uniaxial stress state. Different stress-strain relationships can be assigned to the different fibers ( Figure 4a). Finally, the mechanical behavior of the section is obtained through integration over the entire section. Additionally, co-rotational transformation is adopted to consider geometrical non-linearity [51].
For beam-to-column connections, the Joint2D element in OpenSees is applied, which is idealized as a parallelogram-shaped shear panel with adjacent elements connected to its midpoints ( Figure 4b). The Joint2D element consists of five rotational springs to simulate the shear behavior of the joint panel (central spring) and the moment-rotation behavior of the four sections at beam or column interfaces (interface connections). For the shear panel, the envelope of shear stress-shear strain relation τ-γ behavior is determined by the modified compression-field theory (MCFT) [54]. Eventually, the moment-rotation relation derived from the shear stress-shear strain relation is simplified into a multilinear relationship and then assigned to the central spring by the uniaxial Pinching 4 material model (Figure 4c) in OpenSees. Stiffness and strength are assumed to deteriorate due to the imposed loading history, for which the parameters for cyclic behavior are determined according to the recommendations by Lowes and Altoontash [55]. The four beam-column interface connections are assumed to be rigid [3,29,46], i.e., no moment-rotation relations are assigned to the four beam end and column end springs.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 23 are subjected to a uniaxial stress state. Different stress-strain relationships can be assigned to the different fibers ( Figure 4a). Finally, the mechanical behavior of the section is obtained through integration over the entire section. Additionally, co-rotational transformation is adopted to consider geometrical non-linearity [51].
For beam-to-column connections, the Joint2D element in OpenSees is applied, which is idealized as a parallelogram-shaped shear panel with adjacent elements connected to its midpoints (Figure 4b). The Joint2D element consists of five rotational springs to simulate the shear behavior of the joint panel (central spring) and the moment-rotation behavior of the four sections at beam or column interfaces (interface connections). For the shear panel, the envelope of shear stress-shear strain relation τ-γ behavior is determined by the modified compression-field theory (MCFT) [54]. Eventually, the moment-rotation relation derived from the shear stress-shear strain relation is simplified into a multilinear relationship and then assigned to the central spring by the uniaxial Pinching 4 material model (Figure 4c) in OpenSees. Stiffness and strength are assumed to deteriorate due to the imposed loading history, for which the parameters for cyclic behavior are determined according to the recommendations by Lowes and Altoontash [55]. The four beam-column interface connections are assumed to be rigid [3,29,46], i.e., no moment-rotation relations are assigned to the four beam end and column end springs.

Material Models
The uniaxial plastic-damage model in OpenSees is adopted for concrete fibers, i.e., the ConcreteD material [46,56]. The confinement effect caused by the stirrups on the compressive behavior is taken into account adopting the Mander model [57] (curve 1 vs. curve 3 in Figure 5a). The tension stiffening effect is taken into account using the model by Stevens et al. [58], considering the effect of longitudinal reinforcement (curve 2 vs. curve 4 in Figure 5a).

Material Models
The uniaxial plastic-damage model in OpenSees is adopted for concrete fibers, i.e., the ConcreteD material [46,56]. The confinement effect caused by the stirrups on the compressive behavior is taken into account adopting the Mander model [57] (curve 1 vs. curve 3 in Figure 5a). The tension stiffening effect is taken into account using the model by Stevens et al. [58], considering the effect of longitudinal reinforcement (curve 2 vs. curve 4 in Figure 5a). are subjected to a uniaxial stress state. Different stress-strain relationships can be assigned to the different fibers ( Figure 4a). Finally, the mechanical behavior of the section is obtained through integration over the entire section. Additionally, co-rotational transformation is adopted to consider geometrical non-linearity [51].
For beam-to-column connections, the Joint2D element in OpenSees is applied, which is idealized as a parallelogram-shaped shear panel with adjacent elements connected to its midpoints (Figure 4b). The Joint2D element consists of five rotational springs to simulate the shear behavior of the joint panel (central spring) and the moment-rotation behavior of the four sections at beam or column interfaces (interface connections). For the shear panel, the envelope of shear stress-shear strain relation τ-γ behavior is determined by the modified compression-field theory (MCFT) [54]. Eventually, the moment-rotation relation derived from the shear stress-shear strain relation is simplified into a multilinear relationship and then assigned to the central spring by the uniaxial Pinching 4 material model (Figure 4c) in OpenSees. Stiffness and strength are assumed to deteriorate due to the imposed loading history, for which the parameters for cyclic behavior are determined according to the recommendations by Lowes and Altoontash [55]. The four beam-column interface connections are assumed to be rigid [3,29,46], i.e., no moment-rotation relations are assigned to the four beam end and column end springs.

Material Models
The uniaxial plastic-damage model in OpenSees is adopted for concrete fibers, i.e., the ConcreteD material [46,56]. The confinement effect caused by the stirrups on the compressive behavior is taken into account adopting the Mander model [57] (curve 1 vs. curve 3 in Figure 5a). The tension stiffening effect is taken into account using the model by Stevens et al. [58], considering the effect of longitudinal reinforcement (curve 2 vs. curve 4 in Figure 5a).  The uniaxial Giuffre-Menegotto-Pinto steel material model, i.e., the Steel02 material in OpenSees, is adopted for reinforcing bars [59,60]. The typical hysteretic behavior is presented in Figure 5b, for which the recommended parameters (R0 = 20, cR1 = 0.925, and cR2 = 0.15 when applying the Steel02 materials) for the cyclic behavior are adopted [51].
Material failure may occur when a structure is undergoing progressive collapse, e.g., reinforcement fracture and concrete crushing. The min-max material model in OpenSees can be adopted to model the material failure [46,51]. If the strain of the material exceeds the predefined minimum and maximum values, the material is assumed to have failed and both zero stress and stiffness are returned. Specifically, the thresholds for reinforcement are ε min = −0.01 (buckling) and ε max = the ultimate strain (fracture); for unconfined concrete, these are ε min = −0.0035 (crushing) and ε max = 0.001 (tensile failure); and for confined concrete, these are ε min = −0.03 (crushing) and ε max = 0.001 (tensile failure) [46,[60][61][62].

Validation of the Modelling Techniques
To validate the previously described modelling techniques, a one-third-scale threestorey four-bay RC frame tested by Yi et al. [63] is employed. The configuration of the frame is shown in Figure 6. The middle column at the ground floor was replaced by a hydraulic jack and a concentrated load was imposed on the top of the column at the third floor. Subsequently, the hydraulic jack in the first storey was gradually unloaded to simulate the progressive loss of the middle column, in which the vertical displacement was recorded. More details can be found in the article from Yi et al. [63]. Regarding the material mechanical properties, parameters obtained from the test are adopted [63].
Material failure may occur when a structure is undergoing progressive collapse, e.g., reinforcement fracture and concrete crushing. The min-max material model in OpenSees can be adopted to model the material failure [46,51]. If the strain of the material exceeds the predefined minimum and maximum values, the material is assumed to have failed and both zero stress and stiffness are returned. Specifically, the thresholds for reinforcement are ɛmin = -0.01 (buckling) and ɛmax = the ultimate strain (fracture); for unconfined concrete, these are ɛmin = -0.0035 (crushing) and ɛmax = 0.001 (tensile failure); and for confined concrete, these are ɛmin = -0.03 (crushing) and ɛmax = 0.001 (tensile failure) [46,[60][61][62].

Validation of the Modelling Techniques
To validate the previously described modelling techniques, a one-third-scale threestorey four-bay RC frame tested by Yi et al. [63] is employed. The configuration of the frame is shown in Figure 6. The middle column at the ground floor was replaced by a hydraulic jack and a concentrated load was imposed on the top of the column at the third floor. Subsequently, the hydraulic jack in the first storey was gradually unloaded to simulate the progressive loss of the middle column, in which the vertical displacement was recorded. More details can be found in the article from Yi et al. [63]. Regarding the material mechanical properties, parameters obtained from the test are adopted [63]. The difference between the imposed load P and the axial load N in the middle column, i.e., P -N, against the displacement is shown in Figure 7a, i.e., the pushdown curve. A good agreement is found between the test and the simulated pushdown curve, in which both the beam mechanism and the catenary action (CA) are well captured. The ultimate load-bearing capacities are 106.1 kN and 107.3 kN for the test and simulated results, respectively. In addition, the horizontal displacements at the beam-column joints in the first floor ('3-1′, '3-2′, '3-3′, and '3-4′ in Figure 6) are plotted and compared to the test results (positive values correspond with a movement towards the central removed column), as in Figure 7b. It can be seen that a good performance is obtained for the previously mentioned numerical modelling approach.
Subsequently, the predicted dynamic capacity curve using the EBM is calculated and plotted in Figure 7a (curve 'EBM'). The ultimate dynamic capacity is 87.2 kN. To validate the EBM, a series of NTHA, i.e., the incremental dynamic analysis (IDA) technique, are carried out to determine the dynamic capacity curve by progressively increasing the imposed concentrate loads from a small value to the ultimate load-bearing capacity. The result is plotted in Figure 7a (curve 'IDA'). Two situations are considered to study the influence of damping effects (Rayleigh damping is adopted, in which the involved stiffness matrix is updated in each step), i.e., the damping ratios of DR = 0% and DR = 5%. The EBM curve is almost identical to the IDA curve in the case of DR = 0% (IDA-0%), in which The difference between the imposed load P and the axial load N in the middle column, i.e., P-N, against the displacement is shown in Figure 7a, i.e., the pushdown curve. A good agreement is found between the test and the simulated pushdown curve, in which both the beam mechanism and the catenary action (CA) are well captured. The ultimate loadbearing capacities are 106.1 kN and 107.3 kN for the test and simulated results, respectively. In addition, the horizontal displacements at the beam-column joints in the first floor ('3-1 , '3-2 , '3-3 , and '3-4 in Figure 6) are plotted and compared to the test results (positive values correspond with a movement towards the central removed column), as in Figure 7b. It can be seen that a good performance is obtained for the previously mentioned numerical modelling approach.
Subsequently, the predicted dynamic capacity curve using the EBM is calculated and plotted in Figure 7a (curve 'EBM'). The ultimate dynamic capacity is 87.2 kN. To validate the EBM, a series of NTHA, i.e., the incremental dynamic analysis (IDA) technique, are carried out to determine the dynamic capacity curve by progressively increasing the imposed concentrate loads from a small value to the ultimate load-bearing capacity. The result is plotted in Figure 7a (curve 'IDA'). Two situations are considered to study the influence of damping effects (Rayleigh damping is adopted, in which the involved stiffness matrix is updated in each step), i.e., the damping ratios of DR = 0% and DR = 5%. The EBM curve is almost identical to the IDA curve in the case of DR = 0% (IDA-0%), in which the resistance is 87.0 kN. In the case of DR = 5% (IDA-5%), the IDA curve is almost identical to the EBM curve in the small displacement stages, while some differences are found at the larger displacement levels. The ultimate capacity is 93.0 kN for that case. Overall, a good performance for the EBM can be observed. the resistance is 87.0 kN. In the case of DR = 5% (IDA-5%), the IDA curve is almost identical to the EBM curve in the small displacement stages, while some differences are found at the larger displacement levels. The ultimate capacity is 93.0 kN for that case. Overall, a good performance for the EBM can be observed.

Description of the Structural Model
A 5-storey and 4-bay RC frame is designed according to the Eurocodes [64,65]. The frame presented in Figure 8a is assumed to be a part of an office building. The height of the first storey is 4.5 m, while the height of the other four storeys is 3.6 m each. The span for each bay is 6.0 m and the distance between two adjacent frames is also 6.0 m. The cross-sections and reinforcement layout of beams and columns are shown in Figure 8b Concrete of type C20/25 [64] is used, i.e., the characteristic cylinder compressive strength is fck = 20 MPa. According to fib Model Code 2010 [61], the mean compressive strength of concrete is assumed to be fc = fck + 8 = 28 MPa, while the mean tensile strength is fct = 0.3 (fck) 2/3 = 2.2 MPa. The mechanical properties of the concrete are summarized in Table 1.
The characteristic yield stress and tensile strength of the reinforcing steel are fyk = 500 MPa and fuk = 575 MPa (ductility class C), respectively [64]. The mean yield stress of steel

Description of the Structural Model
A 5-storey and 4-bay RC frame is designed according to the Eurocodes [64,65]. The frame presented in Figure 8a is assumed to be a part of an office building. The height of the first storey is 4.5 m, while the height of the other four storeys is 3.6 m each. The span for each bay is 6.0 m and the distance between two adjacent frames is also 6.0 m. The cross-sections and reinforcement layout of beams and columns are shown in Figure 8b. The dimensions of the beams and columns are 250 mm × 500 mm and 500 mm × 500 mm, respectively. The concrete cover is 30 mm in all elements. Dead loads of both floors and roof are assumed to be DL = 5.0 kN/m 2 (characteristic value), while the live loads on both the floors and roof are LL = 3.0 kN/m 2 (characteristic value) [65]. The columns on the ground floor are labeled as A, B, C, D, and E from left to right (Figure 8a). the resistance is 87.0 kN. In the case of DR = 5% (IDA-5%), the IDA curve is almost identical to the EBM curve in the small displacement stages, while some differences are found at the larger displacement levels. The ultimate capacity is 93.0 kN for that case. Overall, a good performance for the EBM can be observed.

Description of the Structural Model
A 5-storey and 4-bay RC frame is designed according to the Eurocodes [64,65]. The frame presented in Figure 8a is assumed to be a part of an office building. The height of the first storey is 4.5 m, while the height of the other four storeys is 3.6 m each. The span for each bay is 6.0 m and the distance between two adjacent frames is also 6.0 m. The cross-sections and reinforcement layout of beams and columns are shown in Figure 8b Concrete of type C20/25 [64] is used, i.e., the characteristic cylinder compressive strength is fck = 20 MPa. According to fib Model Code 2010 [61], the mean compressive strength of concrete is assumed to be fc = fck + 8 = 28 MPa, while the mean tensile strength is fct = 0.3 (fck) 2/3 = 2.2 MPa. The mechanical properties of the concrete are summarized in Table 1.
The characteristic yield stress and tensile strength of the reinforcing steel are fyk = 500 MPa and fuk = 575 MPa (ductility class C), respectively [64]. The mean yield stress of steel Concrete of type C20/25 [64] is used, i.e., the characteristic cylinder compressive strength is f ck = 20 MPa. According to fib Model Code 2010 [61], the mean compressive strength of concrete is assumed to be f c = f ck + 8 = 28 MPa, while the mean tensile strength is f ct = 0.3 (f ck ) 2/3 = 2.2 MPa. The mechanical properties of the concrete are summarized in Table 1.
The characteristic yield stress and tensile strength of the reinforcing steel are f yk = 500 MPa and f uk = 575 MPa (ductility class C), respectively [64]. The mean yield stress of steel is assumed to be f y = f yk + 2σ 1 = 560 MPa, where σ 1 = 30 MPa is the standard deviation [66]. The mean tensile strength is assumed to be f u = f uk + 2σ 2 = 655 MPa, where σ 2 = 40 MPa is the standard deviation [66]. The ultimate strain of the reinforcement is assumed to be 12% [44]. The Young's modulus E s = 205 GPa is adopted. The mechanical properties of the reinforcement are summarized in Table 1. A 2D finite element (FE) model of the RC frame is built in OpenSees using the previously introduced modelling techniques (Section 4). On the basis of the FE model and the mean values of the previously mentioned material properties, deterministic analyses of both incremental dynamic analyses (IDA) and static pushdown analyses are carried out. Three cases of column removal scenarios are investigated: Case A: loss of the column A; Case B: loss of the column B; and Case C: loss of the column C.

Dynamic Non-Linear Time History Analysis (NTHA)
To determine the dynamic capacity curve, an IDA is carried out by successively increasing the vertical imposed loads on the structure by P i+1 = P i + ∆P, where ∆P is the load increment. A NTHA is applied for each load level, fully taking into account the dynamic effects. This gives the most accurate results, but one NTHA can only obtain the dynamic response for one imposed load case for a structural system and thus a series of NTHA (i.e., an IDA) is needed to determine the dynamic capacity curve. As the dynamic structural response becomes highly non-linear when approaching the ultimate load-bearing capacity, the load increment ∆P of the IDA is gradually reduced. To accurately determine the resistance, a final load increment resolution of 0.05 kN/m is applied near the load for which failure occurs. Note that shear failure cannot be reproduced with the current model.
For every load level, one NTHA is executed, in which the uniform loads on all the beams are first applied, followed by the instantaneous removal of column A, B, or C (a removal time duration of 0.001 s). The explicit Kolay-Ricles-alpha algorithm in OpenSees is employed to execute the dynamic column removal analysis, for which the time step is set as 0.001 s [46]. The dynamic response of the first 4 s for each NTHA is recorded at the top of column A, B, or C (at the control point). Two different values for the damping ratio are investigated in the dynamic analyses (DR = 0% and DR = 5%, where 5% is a representative value for RC structures and Rayleigh damping is adopted) [16,67]. As the EBM cannot take into account any dynamic effects, DR = 0% in IDA is therefore investigated for a comparison purpose. Actually, the following results show that almost identical results are obtained between EBM and IDA if DR = 0%.
The time-history displacement responses at the control points from the IDA when DR = 5% are presented in Figure 9a (Figure 10a), and Case C (Figure 11a), respectively. A summary of the results is presented in Table 2.  (Figure 9a), Case B (Figure 10a), and Case C (Figure 11a), respectively. A summary of the results is presented in Table 2.   The time-history displacement responses from the IDA when DR = 0% are presented in Figures 9b, 10b, and 11b for Cases A, B, and C, respectively. The system oscillates   The time-history displacement responses from the IDA when DR = 0% are presented in Figure 9b, Figure 10b, and Figure 11b for Cases A, B, and C, respectively. The system oscillates around the equilibrium position and the oscillation does not decay as the damping ratio is set as 0%. A summary of the resistances is presented in Table 2, in which the resistance in the case of DR = 0% is slightly lower than in the case of DR = 5%.
From Figure 9a,b it can be seen that the oscillations are irregular in Case A. The reason for this can be found in Figure 9c,d, and e that show the Fourier amplitude spectra (FAS) of the time-history displacement responses using the Fourier transform technique [68][69][70]. It can be seen that there are two main frequency components in the frequency domain as shown in Figure 9c,d, and e. Therefore, this configuration cannot be assumed as a single deformation mode or a SDOF dynamic process, which may give rise to an error from this assumption as discussed in Section 1. The oscillations may, on one hand, disappear in the time-history displacement when DR = 5% (Figure 9a), as the damping effects dissipate the kinetic energy, especially under the higher loading levels in which the oscillations decay fast. On the other hand, there is little oscillation before the displacement response reaches its first reversal point, in which a longer duration is needed to reach the first reversal point for a larger load level (Figure 9a,b). The two aspects result in a wide bond in the frequency domain, as can be seen in Figure 9c,d. However, when only the responses once the system starts oscillating (after the first reversal point as shown in Figure 9b in the case of DR = 0%) are considered, two peaks in the frequency domain (dominant frequencies) can be observed clearly, as can be seen in Figure 9e. around the equilibrium position and the oscillation does not decay as the damping ratio is set as 0%. A summary of the resistances is presented in Table 2, in which the resistance in the case of DR = 0% is slightly lower than in the case of DR = 5%. From Figure 9a,b it can be seen that the oscillations are irregular in Case A. The reason for this can be found in Figure 9c,d, and e that show the Fourier amplitude spectra (FAS) of the time-history displacement responses using the Fourier transform technique [68][69][70]. It can be seen that there are two main frequency components in the frequency domain as shown in Figure 9c,d, and e. Therefore, this configuration cannot be assumed as a single deformation mode or a SDOF dynamic process, which may give rise to an error from this assumption as discussed in Section 1. The oscillations may, on one hand, disappear in the time-history displacement when DR = 5% (Figure 9a), as the damping effects dissipate the kinetic energy, especially under the higher loading levels in which the oscillations decay fast. On the other hand, there is little oscillation before the displacement response reaches its first reversal point, in which a longer duration is needed to reach the first reversal point for a larger load level (Figure 9a,b). The two aspects result in a wide bond in the frequency domain, as can be seen in Figure 9c Moreover, the frequencies associated with the two peaks (Figure 9c-e) are found to be around the first and second natural frequencies obtained from a linear elastic modal analysis, as can be seen in Figure 12a,b. This means the two main frequency components around the equilibrium position and the oscillation does not decay as the damping ratio is set as 0%. A summary of the resistances is presented in Table 2, in which the resistance in the case of DR = 0% is slightly lower than in the case of DR = 5%. From Figure 9a,b it can be seen that the oscillations are irregular in Case A. The reason for this can be found in Figure 9c,d, and e that show the Fourier amplitude spectra (FAS) of the time-history displacement responses using the Fourier transform technique [68][69][70]. It can be seen that there are two main frequency components in the frequency domain as shown in Figure 9c,d, and e. Therefore, this configuration cannot be assumed as a single deformation mode or a SDOF dynamic process, which may give rise to an error from this assumption as discussed in Section 1. The oscillations may, on one hand, disappear in the time-history displacement when DR = 5% (Figure 9a), as the damping effects dissipate the kinetic energy, especially under the higher loading levels in which the oscillations decay fast. On the other hand, there is little oscillation before the displacement response reaches its first reversal point, in which a longer duration is needed to reach the first reversal point for a larger load level (Figure 9a,b). The two aspects result in a wide bond in the frequency domain, as can be seen in Figure 9c Moreover, the frequencies associated with the two peaks (Figure 9c-e) are found to be around the first and second natural frequencies obtained from a linear elastic modal analysis, as can be seen in Figure 12a,b. This means the two main frequency components Moreover, the frequencies associated with the two peaks (Figure 9c-e) are found to be around the first and second natural frequencies obtained from a linear elastic modal analysis, as can be seen in Figure 12a,b. This means the two main frequency components in the responses may be linked to the first two modes (note that the former is obtained from the non-linear dynamic analyses, while the latter is based on the linear modal analyses), i.e., evidently the response cannot be assumed as a single deformation mode.
Regarding the time-history displacement responses for Case B and Case C, the oscillations are more regular as only one peak (i.e., a single deformation mode) is found in the corresponding FAS (calculated once the system starts oscillating from the first reversal point and DR = 0%), as shown in Figures 10c and 11c, respectively. Moreover, the frequency peak (Figure 10c) is around the third natural frequency (obtained from a modal analysis) for Case B, as can be seen in Figure 13c. This means the third mode may be associated with the response as both its model shape and the displacement response in Case C are in the vertical direction. Note that the first two modes are mainly in the horizontal direction, although the second natural frequency is also close to the dominant frequency. Moreover the other natural frequencies are much higher than the dominant frequency. A similar phenomenon can be observed in Case C, i.e., a single deformation mode (Figure 11c) associated with the third mode (Figure 14c). in the responses may be linked to the first two modes (note that the former is obtained from the non-linear dynamic analyses, while the latter is based on the linear modal analyses), i.e., evidently the response cannot be assumed as a single deformation mode. Regarding the time-history displacement responses for Case B and Case C, the oscillations are more regular as only one peak (i.e., a single deformation mode) is found in the corresponding FAS (calculated once the system starts oscillating from the first reversal point and DR = 0%), as shown in Figures 10c and 11c, respectively. Moreover, the frequency peak (Figure 10c) is around the third natural frequency (obtained from a modal analysis) for Case B, as can be seen in Figure 13c. This means the third mode may be associated with the response as both its model shape and the displacement response in Case C are in the vertical direction. Note that the first two modes are mainly in the horizontal direction, although the second natural frequency is also close to the dominant frequency. Moreover the other natural frequencies are much higher than the dominant frequency. A similar phenomenon can be observed in Case C, i.e., a single deformation mode ( Figure  11c) associated with the third mode (Figure 14c).    In order to construct the IDA curves (i.e., the envelope obtained from the series of NTHA), the vertical displacement peaks obtained in the NTHA (Figures 9a,b, 10a,b, and  11a,b) against the corresponding imposed loads are presented in Figure 15 (IDA curve) for the three cases and for both DR = 5% and DR = 0%. in the responses may be linked to the first two modes (note that the former is obtained from the non-linear dynamic analyses, while the latter is based on the linear modal analyses), i.e., evidently the response cannot be assumed as a single deformation mode. Regarding the time-history displacement responses for Case B and Case C, the oscillations are more regular as only one peak (i.e., a single deformation mode) is found in the corresponding FAS (calculated once the system starts oscillating from the first reversal point and DR = 0%), as shown in Figures 10c and 11c, respectively. Moreover, the frequency peak (Figure 10c) is around the third natural frequency (obtained from a modal analysis) for Case B, as can be seen in Figure 13c. This means the third mode may be associated with the response as both its model shape and the displacement response in Case C are in the vertical direction. Note that the first two modes are mainly in the horizontal direction, although the second natural frequency is also close to the dominant frequency. Moreover the other natural frequencies are much higher than the dominant frequency. A similar phenomenon can be observed in Case C, i.e., a single deformation mode ( Figure  11c) associated with the third mode (Figure 14c).    In order to construct the IDA curves (i.e., the envelope obtained from the series of NTHA), the vertical displacement peaks obtained in the NTHA (Figures 9a,b, 10a,b, and  11a,b) against the corresponding imposed loads are presented in Figure 15 (IDA curve) for the three cases and for both DR = 5% and DR = 0%. in the responses may be linked to the first two modes (note that the former is obtained from the non-linear dynamic analyses, while the latter is based on the linear modal analyses), i.e., evidently the response cannot be assumed as a single deformation mode. Regarding the time-history displacement responses for Case B and Case C, the oscillations are more regular as only one peak (i.e., a single deformation mode) is found in the corresponding FAS (calculated once the system starts oscillating from the first reversal point and DR = 0%), as shown in Figures 10c and 11c, respectively. Moreover, the frequency peak (Figure 10c) is around the third natural frequency (obtained from a modal analysis) for Case B, as can be seen in Figure 13c. This means the third mode may be associated with the response as both its model shape and the displacement response in Case C are in the vertical direction. Note that the first two modes are mainly in the horizontal direction, although the second natural frequency is also close to the dominant frequency. Moreover the other natural frequencies are much higher than the dominant frequency. A similar phenomenon can be observed in Case C, i.e., a single deformation mode ( Figure  11c) associated with the third mode (Figure 14c).    In order to construct the IDA curves (i.e., the envelope obtained from the series of NTHA), the vertical displacement peaks obtained in the NTHA (Figures 9a,b, 10a,b, and  11a,b) against the corresponding imposed loads are presented in Figure 15 (IDA curve) for the three cases and for both DR = 5% and DR = 0%. In order to construct the IDA curves (i.e., the envelope obtained from the series of NTHA), the vertical displacement peaks obtained in the NTHA (Figures 9a,b, 10a,b, and  11a,b) against the corresponding imposed loads are presented in Figure 15 (IDA curve) for the three cases and for both DR = 5% and DR = 0%.

Non-linear Static Analysis
Before applying the EBM, one pushdown analysis is carried out for each column removal case. Uniform loads are gradually increased on all beams. A displacement-controlled analysis is executed, in which the displacement at the top of column A, B, or C is controlled, depending on the column removal scenario that is considered. The displacements at the control points against the imposed loads are recorded, resulting in the pushdown curves in Figure 15a Moreover, the displacements at the first failure (concrete crushing at beam ends in the first floor at the first load peaks in the pushdown curves) according to the static analyses are almost identical to those of the IDA results for Cases A, B, and C. The first failure in Case C corresponds with 37.9 kN/m (279.2 mm) for the IDA. The ultimate load-bearing capacity is, however, higher at 38.3 kN/m (341.4 mm). In Case C, the maximum displacement following the IDA (corresponding to the dynamic ultimate load-bearing capacity) is closer to the displacement at the second load peak in the pushdown curve (when failure occurs at the beam ends in the second floor), as can be seen in Figure 15c. This explains why the ultimate load-bearing capacities are obtained at different displacement levels when comparing the static result to the dynamic results in Case C. For cases A and B, the static and dynamic ultimate load-bearing capacities are obtained at almost the same displacement levels.
For all three cases, the second load peak in the pushdown curve is lower than the first load peak (failure only in the first floor). Consequently, only the first load peak is used here as it relates to the ultimate load-bearing capacity and the failure mode (first failure) is identical to the IDA. Furthermore, the first failure is a signal of the start of the system failure.

Dynamic Amplification Factor
Typically, the DAF is defined as the ratio of the maximum dynamic displacement (ud) to the static displacement (us) for an elastic SDOF system under an imposed loading (Qd) [71]. Similarly, a displacement-based DAFD can also be defined with regard to the IDA and pushdown curves (Figure 16a): where ud is the dynamic displacement in the dynamic capacity curve and us is the static displacement under the same imposed load in the static capacity curve. Note that DAFD

Non-Linear Static Analysis
Before applying the EBM, one pushdown analysis is carried out for each column removal case. Uniform loads are gradually increased on all beams. A displacementcontrolled analysis is executed, in which the displacement at the top of column A, B, or C is controlled, depending on the column removal scenario that is considered. The displacements at the control points against the imposed loads are recorded, resulting in the pushdown curves in Figure 15a Moreover, the displacements at the first failure (concrete crushing at beam ends in the first floor at the first load peaks in the pushdown curves) according to the static analyses are almost identical to those of the IDA results for Cases A, B, and C. The first failure in Case C corresponds with 37.9 kN/m (279.2 mm) for the IDA. The ultimate load-bearing capacity is, however, higher at 38.3 kN/m (341.4 mm). In Case C, the maximum displacement following the IDA (corresponding to the dynamic ultimate load-bearing capacity) is closer to the displacement at the second load peak in the pushdown curve (when failure occurs at the beam ends in the second floor), as can be seen in Figure 15c. This explains why the ultimate load-bearing capacities are obtained at different displacement levels when comparing the static result to the dynamic results in Case C. For cases A and B, the static and dynamic ultimate load-bearing capacities are obtained at almost the same displacement levels.
For all three cases, the second load peak in the pushdown curve is lower than the first load peak (failure only in the first floor). Consequently, only the first load peak is used here as it relates to the ultimate load-bearing capacity and the failure mode (first failure) is identical to the IDA. Furthermore, the first failure is a signal of the start of the system failure.

Dynamic Amplification Factor
Typically, the DAF is defined as the ratio of the maximum dynamic displacement (u d ) to the static displacement (u s ) for an elastic SDOF system under an imposed loading (Q d ) [71]. Similarly, a displacement-based DAF D can also be defined with regard to the IDA and pushdown curves (Figure 16a): where u d is the dynamic displacement in the dynamic capacity curve and u s is the static displacement under the same imposed load in the static capacity curve. Note that DAF D can only be calculated up to the ultimate dynamic response, even if a higher portion still exists in the static capacity curve (Figure 16a).
A force-based DAF F [71,72] can also be calculated by dividing the static load-bearing capacity by the dynamic load-bearing capacity at the same displacement (Figure 16a): where Q s is the load in the static capacity curve and Q d is the load in the dynamic capacity curve. Based on the previous pushdown and IDA curves, the displacement-based DAF D is calculated and presented in Figure 16b. Regarding Case A, the DAF D increases from 1.74 to 4.11 when the imposed load increases. In terms of Case B and Case C, the DAF D increases from 2.00 to 4.93 and 4.95 when the imposed load increases, respectively. It is worth noting that the results are obtained using a resolution of 0.05 kN/m for the imposed load increment in the IDA. A similar result is obtained for a bilinear SDOF model with different non-linear parameters as studied by Tsai [71], in which the DAF D may increase from 2.00 to infinity as the applied loading increases. The values lower than 2.00 in Case A in Figure 16b can be attributed to the fact that the response in this case is not a single deformation mode as mentioned before. Additionally, yield points (YP), i.e., the moment when the reinforcement starts yielding (see Table 1), in the curves are marked (Figure 16b), for which DAF D is 2.10, 2.28, and 2.32 for the three cases. can only be calculated up to the ultimate dynamic response, even if a higher portion still exists in the static capacity curve (Figure 16a). A force-based DAFF [71,72] can also be calculated by dividing the static load-bearing capacity by the dynamic load-bearing capacity at the same displacement (Figure 16a): where Qs is the load in the static capacity curve and Qd is the load in the dynamic capacity curve. Based on the previous pushdown and IDA curves, the displacement-based DAFD is calculated and presented in Figure 16b. Regarding Case A, the DAFD increases from 1.74 to 4.11 when the imposed load increases. In terms of Case B and Case C, the DAFD increases from 2.00 to 4.93 and 4.95 when the imposed load increases, respectively. It is worth noting that the results are obtained using a resolution of 0.05 kN/m for the imposed load increment in the IDA. A similar result is obtained for a bilinear SDOF model with different non-linear parameters as studied by Tsai [71], in which the DAFD may increase from 2.00 to infinity as the applied loading increases. The values lower than 2.00 in Case A in Figure 16b can be attributed to the fact that the response in this case is not a single deformation mode as mentioned before. Additionally, yield points (YP), i.e., the moment when the reinforcement starts yielding (see Table 1), in the curves are marked ( Figure  16b), for which DAFD is 2.10, 2.28, and 2.32 for the three cases. The force-based DAFF is calculated and shown in Figure 16c. For Case A, the DAFF decreases from 2.23 (when the vertical imposed load is 4.0 kN/m in the IDA curve) to 1.10 (at the displacement where first failure occurs) when the displacement increases. In terms of Cases B and C, the DAFF decreases from 2.02 (at 4.0 kN/m) to 1.08 and 1.04 (at the displacement where first failure occurs), respectively. A similar result is obtained for a 5storey RC frame subjected to a column removal scenario as studied by Brunesi and Nascimbene [72], in which the DAFF decreased from almost 2.70 to 1.11 with the increasing vertical displacement. Again, yield points (YP) in the curves are marked (Figure 16c), for which DAFF is 1.65, 1.73, and 1.72 for the three cases.

Energy-Based Method
The EBM is subsequently applied to calculate the dynamic load-bearing capacity curve using Equation (1). The results for the different cases are presented in Figure 15. Generally, the EBM is observed to have a good performance in relation to predicting the maximum dynamic response in comparison with the more accurate IDA results. Regarding the IDA curves with DR = 5%, a slight deviation is observed for the EBM curve as the The force-based DAF F is calculated and shown in Figure 16c. For Case A, the DAF F decreases from 2.23 (when the vertical imposed load is 4.0 kN/m in the IDA curve) to 1.10 (at the displacement where first failure occurs) when the displacement increases. In terms of Cases B and C, the DAF F decreases from 2.02 (at 4.0 kN/m) to 1.08 and 1.04 (at the displacement where first failure occurs), respectively. A similar result is obtained for a 5-storey RC frame subjected to a column removal scenario as studied by Brunesi and Nascimbene [72], in which the DAF F decreased from almost 2.70 to 1.11 with the increasing vertical displacement. Again, yield points (YP) in the curves are marked (Figure 16c), for which DAF F is 1.65, 1.73, and 1.72 for the three cases.

Energy-Based Method
The EBM is subsequently applied to calculate the dynamic load-bearing capacity curve using Equation (1). The results for the different cases are presented in Figure 15. Generally, the EBM is observed to have a good performance in relation to predicting the maximum dynamic response in comparison with the more accurate IDA results. Regarding the IDA curves with DR = 5%, a slight deviation is observed for the EBM curve as the damping effect is considered in the IDA, which is not possible to account for in the case of the EBM. However, the overall influence remains rather small. The deviation in Case A is slightly larger, which can be attributed to the fact that, as found before, this case is not a SDOF dynamic process. Nevertheless, an excellent performance is observed for all three cases if DR = 0%, for which the EBM curves are almost identical to the IDA curves. A summary of the results can be found in Table 2.

Probabilistic Models of Random Variables
The RC frame is consecutively also used to investigate the influence of uncertainties and to quantify the model uncertainty related to the use of the EBM. In order to calculate the stochastic response for the three cases (i.e., Cases A, B, and C), eight parameters are selected as random variables, which are presented in Table 3. c is the concrete cover, which is modelled as a bounded Beta distribution [73], with a mean value equal to a specified value of 30 mm, a standard deviation of 5 mm, a lower bound of 0 mm, and an upper bound of three times that of the mean value, i.e., 90 mm. f c is the concrete compressive strength, while ε c1 is the peak compressive strain [44,61,66,73]. For modelling the concrete tensile strength, Y 2,j is employed to reflect variations due to factors being not well accounted for by the concrete compressive strength, i.e., f ct = 0.3(f ck ) 2/3 Y 2,j [66]. E s , f y , f u , and ε u are the modulus of elasticity, the yield stress, the tensile strength, and the ultimate strain of the reinforcement, respectively [44,66,73].

Stochastic Analysis
Based on the probabilistic models presented in Table 3, 10 000 Latin Hypercube Samples (LHS) are generated. Considering that the standard LHS may bring undesired spurious correlations into the sample scheme, correlation reduced Latin Hypercube Sampling (CLHS) is used to avoid this unwanted effect [74]. It is worth noting that all the parameters are assumed to be independent, except for the yield stress f y and the tensile strength f u of the reinforcement, for which a correlation coefficient of 0.86 is considered [66]. This is an assumption adopted in many studies in relation to the consideration of material uncertainties for RC structures [44,[75][76][77][78].
Subsequently, both a pushdown analysis and an IDA are carried out for every sample for every column removal case, i.e., Cases A, B, and C. The loading schemes are identical to the corresponding deterministic analyses in the previous section. In order to reduce the computational demand in the stochastic dynamic analyses, a resolution level of 0.5 kN/m (compared to 0.05 kN/m used for the previous deterministic dynamic analysis) is used for the IDA when determining the ultimate load-bearing capacity. The damping ratio DR = 5% is considered in the stochastic dynamic analysis.
The histograms of the ultimate load-bearing capacity from both static and dynamic analyses are shown in Figure 17 and a summary of the results is presented in Table 4. Note that the resistances obtained through the pushdown analyses are considered to correspond to first failure as discussed before. Cases A, B, and C are presented in Figure 17a-c, respectively, in which both the results from the pushdown analyses and IDA are shown in order to compare with each other. Lognormal distributions are found to fit the histograms well for the three cases. The resistance from the pushdown analysis is found to be larger than that of the IDA, as also already evidenced in the previous sections. Mean values and standard deviation (St.D.) for the stochastic pushdown and IDA analyses are summarized in Table 4. Furthermore, it was found that the resistance in the case of the external column loss scenario (Case A) is significantly lower than those for the internal column loss scenarios (Case B and Case C), as it is more difficult to redistribute the unbalanced load in the former case.
that the resistances obtained through the pushdown analyses are considered to correspond to first failure as discussed before. Cases A, B, and C are presented in Figure 17ac, respectively, in which both the results from the pushdown analyses and IDA are shown in order to compare with each other. Lognormal distributions are found to fit the histograms well for the three cases. The resistance from the pushdown analysis is found to be larger than that of the IDA, as also already evidenced in the previous sections. Mean values and standard deviation (St.D.) for the stochastic pushdown and IDA analyses are summarized in Table 4. Furthermore, it was found that the resistance in the case of the external column loss scenario (Case A) is significantly lower than those for the internal column loss scenarios (Case B and Case C), as it is more difficult to redistribute the unbalanced load in the former case.  --1.06% --1.34% -Next, the EBM is adopted to calculate the dynamic capacity curves from the pushdown curves. It is important to indicate that here, as well, the EBM curve is calculated based on the pushdown curve up to the first failure (i.e., the first load peak as it is the ultimate load-bearing capacity) in the frame. The predicted dynamic load-bearing capacity curves are similar to those of the deterministic analyses. The histograms and PDFs of the ultimate load-bearing capacities from EBM are presented in Figure 18 for the three cases, for which lognormal distributions are found to be appropriate to fit the histograms. Only a small discrepancy between the EBM result and the IDA result is observed for Case A (Figure 18a), while little difference is observed between the EBM results and the IDA results for both Cases B and C (Figure 18b,c). The mean resistances from the EBM and their deviations compared to the IDA results are summarized in Table 4. The maximum deviation is -3.17% in Case A and much smaller deviations are found for the other two cases, which shows that the EBM has a good performance to predict the dynamic resistance, in addition to when uncertainties of variables are taken into account.
In order to better compare EBM and IDA results, the resistances and corresponding displacements are compared in Figures 19 and 20, respectively. Figure 19a-c show the comparisons of the resistances between the EBM and IDA for Cases A, B, and C, respectively. It can be seen that the data points are located along the line y = x, which means that  Next, the EBM is adopted to calculate the dynamic capacity curves from the pushdown curves. It is important to indicate that here, as well, the EBM curve is calculated based on the pushdown curve up to the first failure (i.e., the first load peak as it is the ultimate load-bearing capacity) in the frame. The predicted dynamic load-bearing capacity curves are similar to those of the deterministic analyses. The histograms and PDFs of the ultimate load-bearing capacities from EBM are presented in Figure 18 for the three cases, for which lognormal distributions are found to be appropriate to fit the histograms. Only a small discrepancy between the EBM result and the IDA result is observed for Case A (Figure 18a), while little difference is observed between the EBM results and the IDA results for both Cases B and C (Figure 18b,c). The mean resistances from the EBM and their deviations compared to the IDA results are summarized in Table 4. The maximum deviation is −3.17% in Case A and much smaller deviations are found for the other two cases, which shows that the EBM has a good performance to predict the dynamic resistance, in addition to when uncertainties of variables are taken into account.
In order to better compare EBM and IDA results, the resistances and corresponding displacements are compared in Figures 19 and 20, respectively. Figure 19a-c show the comparisons of the resistances between the EBM and IDA for Cases A, B, and C, respectively. It can be seen that the data points are located along the line y = x, which means that similar results are obtained between the EBM and IDA. Moreover, coefficients of determination (R 2 ) are calculated to be 0.67, 0.93, and 0.92 for the three cases, respectively. Apparently, the deviation in Case A is larger (R 2 = 0.67) as the dynamic response is not a SDOF oscillation, as discussed before. In general, the IDA results are found to be slightly larger as the assumptions (e.g., damping effect) discussed in Section 1 will result in a higher resistance compared to when using the EBM. This means the EBM provides a conservative result for resistance assessment. similar results are obtained between the EBM and IDA. Moreover, coefficients of determination (R 2 ) are calculated to be 0.67, 0.93, and 0.92 for the three cases, respectively. Apparently, the deviation in Case A is larger (R 2 = 0.67) as the dynamic response is not a SDOF oscillation, as discussed before. In general, the IDA results are found to be slightly larger as the assumptions (e.g., damping effect) discussed in Section 1 will result in a higher resistance compared to when using the EBM. This means the EBM provides a conservative result for resistance assessment.   similar results are obtained between the EBM and IDA. Moreover, coefficients of determination (R 2 ) are calculated to be 0.67, 0.93, and 0.92 for the three cases, respectively. Apparently, the deviation in Case A is larger (R 2 = 0.67) as the dynamic response is not a SDOF oscillation, as discussed before. In general, the IDA results are found to be slightly larger as the assumptions (e.g., damping effect) discussed in Section 1 will result in a higher resistance compared to when using the EBM. This means the EBM provides a conservative result for resistance assessment.   similar results are obtained between the EBM and IDA. Moreover, coefficients of determination (R 2 ) are calculated to be 0.67, 0.93, and 0.92 for the three cases, respectively. Apparently, the deviation in Case A is larger (R 2 = 0.67) as the dynamic response is not a SDOF oscillation, as discussed before. In general, the IDA results are found to be slightly larger as the assumptions (e.g., damping effect) discussed in Section 1 will result in a higher resistance compared to when using the EBM. This means the EBM provides a conservative result for resistance assessment.    20a-c show the comparisons of the corresponding displacements between the EBM and IDA for Cases A, B, and C, respectively. Additionally, in this case, the data points are located along the line y = x, which indicates a good agreement between the EBM and IDA. However, a larger variation is observed for the displacement responses as the response varies fast in the highly non-linear stage (see the load-displacement curves in Figure 15). This can also be verified from the coefficients of determination (R 2 ), which are calculated to be 0.77, 0.74, and 0.80 for the three cases, respectively.

Model Uncertainty Quantification of the EBM vs. IDA Analyses
To quantitatively evaluate the performance of the EBM, the ratios K EBM of the resistances of the EBM to those from the incremental dynamic analyses (IDA) are calculated according to Equation (2), as well as for the ratios of the displacements. The mean and standard deviation of the ratios K EBM for Cases A, B, and C separately and as a total for the three cases considered are shown in Table 5. Histograms and PDFs of the ratios with regard to the resistances are presented in Figure 21.  Figure 21d presents the ratios of the ultimate capacities of all three cases as a whole (case 'all' in Table 5), for which a lognormal distribution LN(0.98, 0.02) is fitted. and IDA. However, a larger variation is observed for the displacement responses as the response varies fast in the highly non-linear stage (see the load-displacement curves in Figure 15). This can also be verified from the coefficients of determination (R 2 ), which are calculated to be 0.77, 0.74, and 0.80 for the three cases, respectively.

Model Uncertainty Quantification of the EBM vs. IDA Analyses
To quantitatively evaluate the performance of the EBM, the ratios KEBM of the resistances of the EBM to those from the incremental dynamic analyses (IDA) are calculated according to Equation (2), as well as for the ratios of the displacements. The mean and standard deviation of the ratios KEBM for Cases A, B, and C separately and as a total for the three cases considered are shown in Table 5. Histograms and PDFs of the ratios with regard to the resistances are presented in Figure 21. Figure 21a-c show ratios in terms of the resistances for Cases A, B, and C, respectively, in which lognormal distributions LN(0.97, 0.01), LN(0.99, 0.01), and LN(0.99, 0.01) are found to represent the histograms. Figure 21d presents the ratios of the ultimate capacities of all three cases as a whole (case 'all' in Table 5), for which a lognormal distribution LN(0.98, 0.02) is fitted. Histograms and PDFs of the ratios in terms of the corresponding displacements for the three cases are presented in Figure 22.  According to these results, a great performance has been found for the model uncertainties involved in the application of the EBM instead of the more complex IDA, especially for the prediction of the resistance, i.e., resulting in a small bias and small standard deviation. The model uncertainty models can be adopted in a probabilistic design situation, for instance, as follows. 1) Perform stochastic static non-linear analyses with specified random variables and Histograms and PDFs of the ratios in terms of the corresponding displacements for the three cases are presented in Figure 22.  Histograms and PDFs of the ratios in terms of the corresponding displacements for the three cases are presented in Figure 22.  According to these results, a great performance has been found for the model uncertainties involved in the application of the EBM instead of the more complex IDA, especially for the prediction of the resistance, i.e., resulting in a small bias and small standard deviation. The model uncertainty models can be adopted in a probabilistic design situation, for instance, as follows. 1) Perform stochastic static non-linear analyses with specified random variables and column removal scenarios. According to these results, a great performance has been found for the model uncertainties involved in the application of the EBM instead of the more complex IDA, especially for the prediction of the resistance, i.e., resulting in a small bias and small standard deviation. The model uncertainty models can be adopted in a probabilistic design situation, for instance, as follows.
(1) Perform stochastic static non-linear analyses with specified random variables and column removal scenarios. (2) Apply the EBM (Equation (1)) to calculate dynamic ultimate capacities R according to the static capacity curves from step 1. (3) Evaluate the failure probability P f through the following limit state function Z: where K R and K L are the resistance and load model uncertainties, respectively, [66]; L is the load effect; and K EBM is the additional model uncertainty relating to the application of the EBM instead of the IDA (Table 5).
Although the performance of the EBM is verified and its model uncertainty is calculated herein, it must be emphasized that the results are based on a particular RC frame and that further investigations might be considered in a similar manner as illustrated in this contribution before coming to a more generally applicable proposal for the model uncertainty involved.

Conclusions
In this article, the effectiveness of an energy-based method for the prediction of the dynamic resistance of RC frames was verified through comparing the approximation to more complex non-linear dynamic analyses. First, the performance of the EBM was studied. The influences of damping effects and different column removal scenarios on the performance of the load-bearing capacity prediction were investigated. It was found that the response is not a single deformation mode in the case of an external column removal scenario, which may influence the performance of the EBM and the dynamic amplification factor. Model uncertainties related to the application of the EBM instead of incremental dynamic analyses (IDA) are quantified. The detailed findings are summarized as follows: Although no dynamic effects such as damping can be explicitly considered in the EBM, a good performance was found for the EBM to predict the dynamic capacity curve through comparison with results of the IDA. Almost identical results were obtained between the EBM and IDA curve if a damping ratio of 0% was considered, which indicated that the EBM can effectively calculate the maximum dynamic responses. Only a small difference was found between EBM and IDA results if the damping ratio is taken as 5%, which demonstrates that the damping effects have only a limited influence. Moreover, the EBM provides a conservative result for the resistance assessment.
The EBM assumes that a structure subjected to a column removal scenario deforms in a single deformation mode. An external column removal scenario may influence the performance of the EBM and the dynamic amplification factor as its dynamic response may not be represented by a single deformation mode according to the responses in the frequency domain. However, the internal column removal scenarios were found to behave similar to SDOF dynamic systems, for which the EBM overall has a very good performance. A better performance of the EBM was observed in an internal column loss case rather than for an external column loss case.
In order to allow for the use of the EBM as part of probabilistic failure analyses and robustness assessments, probabilistic models have been proposed for the model uncertainty of the EBM compared to the IDA. For the different removal scenarios considered as a whole, a lognormal distribution with mean of 0.98 and a standard deviation of 0.02, i.e., LN(0.98,0.02), was fitted for the ratio of the EBM resistance to the IDA resistance, while a lognormal distribution of LN(1.07,0.08) was obtained for the ratio of the EBM to IDA displacements. These values are of particular interest to be included as an additional model uncertainty in probabilistic analyses in case the EBM is applied for the quantification of alternative load path resistances considering dynamic behavior.
The values of the model uncertainties associated with the resistances indicate a small bias and small standard deviation. This means that the EBM has a good accuracy in calculating the dynamic resistances. In contrast, a larger model uncertainty was found for the computation of the corresponding displacements, as the mean values were slightly larger than one and the standard deviations were considerably larger.