Modeling of RC Frame Buildings for Progressive Collapse Analysis

: The progressive collapse analysis of reinforced concrete (RC) moment-frame buildings under extreme loads is discussed from the perspective of modeling issues. A threat-independent approach or the alternate path method forms the basis of the simulations wherein the extreme event is modeled via column removal scenarios. Using a prototype RC frame building, issues and considerations in constitutive modeling of materials, options in modeling the structural elements and speciﬁcation of gravity loads are discussed with the goal of achieving consistent models that can be used in collapse scenarios involving successive loss of load-bearing columns at the lowest level of the building. The role of the ﬂoor slabs in mobilizing catenary action and inﬂuencing the progressive collapse response is also highlighted. Finally, an energy-based approach for identifying the proximity to collapse of regular multi-story buildings is proposed.


Introduction
The systematic development of numerical models so as to gain an in-depth understanding of the behavior of concrete frame buildings subjected to extreme loading conditions, such as the sequential removal of load carrying vertical elements, is presented in this paper. This preliminary study is part of a more comprehensive research directed towards the formulation of robustness indices to assess the resistance of reinforced concrete (RC) structures to disproportionate collapse. Assessing the probability of building failure during either an ordinary or exceptional event deserves particular attention because of its relevance on the safety of human communities and the consequent economic impact on society. These studies also encourage the development of new structural design procedures and assessment criteria to prevent failure or minimize damage due to unexpected extreme events.
Within the above framework, most of the research effort during the last decade has been devoted to the study of disproportionate collapse of multi-story buildings. Progressive or disproportionate collapse occurs when a structure has its load pattern or boundary conditions altered in a manner such that some structural elements are loaded beyond their capacity and fail (Krauthammer et al. 2003). As well-documented in El-Tawil et al. (2014), several studies have been carried out in this field especially after the partial collapse of Ronan Point tower in England in 1968 (Pearson and Delatte 2005), with a significant increase in published papers in the field in the last 10 years. Early studies, including the papers by Lewicki and Olesen (1974), Arora et al. (1980), Gross and McGuire (1983), and McConnel and Kelly (1983), while reviewing alternative design methods to prevent progressive collapse and proposing the first computational approach for performing simulations, already highlight the need for unified analytical procedures for evaluating the resistance of structures to progressive collapse. In subsequent years, aided by the increasing availability of powerful computational tools, a large number of numerical studies have been carried out on different types of structures and with different modeling approaches, leading to a wide database of results and a range of research findings.
The simulation models used in previous studies differ from each other in many ways and can be classified into distinct groups based on the model features: (1) the type of analysis: linear or nonlinear-that can significantly affect the response of the structure as documented in Marjanishvili and Agnew (2006), where a comparison between different analyses shows that the predicted responses can vary significantly when performing static/dynamic and linear/nonlinear analyses, (2) the typology of elements: continuum, beam-column elements or a combination of both have been successfully used for modeling local and global phenomena of progressive collapse: examples of micro-models can be found in Khandelwal and El-Tawil (2007), Sasani and Kropelnicki (2008), Kwasniewski (2010), and Bao et al. (2008), while examples of macro-models are those utilized in Kaewkulchai and Williamson (2004), Bao et al. (2008), and Bao and Kunnath (2010) and an example of the use of hybrid models is reported in the work of Alashker et al. (2011), (3) the dimension of the model: planar or three dimensional-that is crucial in capturing spatial effects: most of the published work were conducted on two-dimensional structures (Khandelwal and El-Tawil 2007;Bao et al. 2008;Kim et al. 2009) and very few on three-dimensional models (Ruth et al. 2006;Alashker et al. 2011) making the comparison of results very challenging, (4) the floor system, modeled using a collection of beam-column, shell or brick elements , which plays a key role in determining the response of three-dimensional structures, as shown by Yu et al. (2010), Alashker et al. (2011) and, (5) the loads applied on the structure, which can include only the self-weight or the self-weight and a fraction of the design dead/live loads.
Despite the availability of a large number of numerical studies, there still remains the need for procedures and numerical modeling guidelines to carry out simple yet reliable progressive collapse studies of buildings. The development and use of unified modeling criteria also derives from the necessity of having consistent models so that the results of different simulations can be effectively compared. In the long term this could lead to the development of a database of uniform data, useful to define and validate robustness indices for different typologies of structures. In this context, this study aims at identifying basic modeling features that have to be properly taken into account when analyzing the large-deformation behavior of RC buildings in response to extreme loading conditions. The simulations are based on the alternate path method (APM) wherein the extreme event itself is not simulated but the consequence of the event, i.e., the removal (due to failure) of a critical element is considered. Sensitivity analyses are carried out to assess the adequacy of element and section discretization and the efficacy of alternate options for modeling the floor slab are investigated. An energy-based approach to define proximity to partial or total collapse of the structure is also proposed.

Prototype Model, Materials and Elements
The study was conducted on prototype concrete momentresisting building frames of varying height with and without the incorporation of floor slabs, using explicit time integration in LS-DYNA (Hallquist 2007), a general-purpose finite element code. Primary beams and columns were modeled using Hughes-Liu elements with multiple integration points along the length and cross section integration by means of fiber discretization at integration points, while slabs were represented as layered shells with smeared reinforcement. Both material and geometrical nonlinearities, including damage and fracture, were considered. Center-to-center dimensions were used to define beam and column lengths and joint shear deformations were ignored because separate studies carried out by the authors indicated that the overall displacement responses were not influenced by incorporating special joint elements for the moment frame configurations considered in this study whose behavior is controlled by flexure. Considering joint deformations will be important for the case of shear-critical members (Jeong and Kim 2014).
A prototype building, having plan dimensions 45 m 9 30 m and story height 4.6 m for the 1st floor and 3.7 m for the remaining floors ( Fig. 1), was designed to meet the requirements of seismic design category C for a site in Atlanta, as per ASCE-7 (2010) with detailing conforming to ACI 318 specifications (2014). Design details are summarized in Table 1. The building is characterized by two missing lines of columns-along B and E-to maximize floor space for practical considerations. Transverse reinforcement in the beams consists of single closed ties whereas for columns the closed ties are augmented with a central cross-tie in each direction (as required by ACI 318).
Nonlinearity of steel and concrete in tension and compression was modelled with isotropic elastic-plastic model (material model 124 in LS-DYNA), by defining constitutive laws expressed through effective stress versus plastic strain curves. The basic properties of the adopted materials are: unconfined compressive strength of concrete = 27.6 MPa, yield strength of reinforcement = 413.8 MPa, Young's modulus of concrete = 24.8 GPa, and for steel = 200 GPa. The peak stress and corresponding strain due to confinement as well as the ultimate stress and strain are obtained using the Mander et al. (1988) model. The Scott et al. (1982) model is used to represent the shape of the stress-strain curve up to peak compressive stress based on the material Concrete02 in OpenSees (2015) and a linear representation is used for the post-peak response. Overall three different types of confined concrete and one unconfined concrete were defined, as indicated in Table 2. The stress-strain curve for concrete is shown in Fig. 2a and that for steel in Fig. 2b. The assigned stress-strain curves are different in compression and tension, as shown in Fig. 2. The specified yield stress was the same in tension and compression. The stress-strain relationship of steel in compression accounts for buckling by following the model proposed by Gomes and Appleton (1997). The complete stress-strain curves in tension and compression, respectively were obtained using the ReinforcingSteel material in OpenSees (2015). It was observed that though confinement has an effect on the strength and deformation capacity of RC members, failure at large deformations is more influenced by the ductility of the steel reinforcement.  Shear reinforcement 2: /12 @152 mm 2: /12 @76 mm 2: /12 @178 mm 2: /12 @76 mm 3: /12 @203 mm 3: /12 @178 mm Shear reinforcement 2: /12 @152 mm 2: /12 @76 mm 2: /12 @178 mm 2: /12 @76 mm 3: /12 @203 mm 3: /12 @178 mm Shear reinforcement 2: /12 @152 mm 2: /12 @76 mm 2: /12 @152 mm 2: /12 @76 mm 3: /12 @203 mm 3: /12 @152 mm

Localization of Inelastic Behavior
One of the considerations in specifying constitutive model properties, particularly for the Liu-Hughes beam-column elements used in the current simulation, is the likelihood of localization of inelastic behavior. Such localization can lead to a non-objective response during softening (Bao et al. 2014). Progressive collapse simulations involve large strains with failure occurring during the post-peak phase; hence care must be taken to adjust post-peak parameters to obtain the desired stress-strain response. Shown in Fig. 3 is the localized response if the target stress-strain curve is input without accounting for mesh size (in the case of the reinforcing steel bars, mesh size refers to the length of the discretized element). After adjusting the input stress-strain curve based on the chosen mesh size, the desired stressstrain curve is achieved as shown in the figure labeled ''Nonlocal response''. This is not truly a non-local response (resulting from a non-local model) but an indirect approach that modifies the input constitutive model to avoid localization. Strain estimates in the localized zones from such an approach will be inaccurate; however, the overall displacement estimates (the primary parameter of interest in the present collapse simulations) are reliable.

Sensitivity and Convergence Studies
Prior to finalizing the simulation model, it is essential to gain confidence in the level of discretization of elements as well as the choice of integration points in the cross-section. This is achieved by carrying out sensitivity studies with varying mesh sizes and integration points and examining the distribution of deformations and forces in critical elements. The objective of the simulations is to achieve a reasonable balance between accuracy and computational efficiency. In structural modeling for progressive collapse analysis, the model accuracy has to be verified for both linear and nonlinear responses, when large deformations involving material and geometrical nonlinearities are expected. As stated earlier, the simulations in this study involve column-removal scenarios, consistent with APM, to examine the collapse resistance of concrete frame structures following an extreme event. Consequently, the focus of the modeling is related to the nonlinear response of a frame model to sequential column removal.
Since an explicit time integration method is used in all the simulations, damping across a frequency range needs to be defined. For the presented case studies the damping ratio is set equal to 0.05 and the frequency range is set to cover all relevant frequencies of the structure. The studies in this section were carried out on a single-story building because the validity of the sensitivity study is not affected by the height of the building: such as the localization of inelastic behavior, the element discretization and the cross section integration schemes. As for the modeling of gravity loads and the analysis of alternative grid beam models presented later in the paper, the one-story building case actually represents an upper bound condition-since additional stories would only increase the gravity loads supported by the firststory columns.

Element Discretization
In the first set of simulations, beams were discretized into 8, 10 and 12 elements of equal length. Since beams are deemed to be critical for the evaluation of the structural behavior in response to the sudden loss of a vertical loadcarrying element, element discretization in beams were distributed into two different ways, resulting in six different cases to compare, as displayed in Fig. 4. Another consideration in determining an appropriate discretization was motivated by the fact that plastic hinge zones are crucial for the activation of catenary effect (Yi et al. 2008(Yi et al. , 2014. Two initial options were considered: in the first case 8, 10 and 12 elements were uniformly subdivided into equal lengths along the beam span-Option A, in the next case elements were concentrated in the hinge zones (25 % of beam span) leaving just one node at the mid-span-Option A*. For each case a nonlinear dynamic analysis is performed Fig. 3 Modified stress-strain input to avoid localization: a reinforcing steel model, and b concrete model. first under gravity loads (self-weight) followed by the removal of two columns. Figure 5 shows the deflection of the continuous beam segments A1-3 under gravity loads for both discretization options. Option A resulted in generally similar shapes while Option A* expectedly led to different profiles at the midspan of each beam. A closer view of the deflected shape in the support regions (potential plastic hinge locations) revealed that Option A* provides more consistent shape for all discretizations. However, the need to better represent the mid-span profile led to a third option (A?): starting with Option A comprising 8 elements two additional nodes were added at the mid-span to generate 10 elements evenly distributed in each zone (see Fig. 4).
The deflection profile resulting from this option is shown in all the plots of Fig. 5 and labeled as 10?. The accuracy of the adopted solution, when compared with Options A and A*, was also verified when considering larger deformations. Figure 6 shows the deflection profile of beams A1-3 after removing columns A1 and A2. It is seen that Option A? (with only 10 elements) provides as accurate a response as attained with 12 elements. Next, considering a deformed profile in the short direction of the building, Fig. 7 shows the deflection profile of beams D2-C2-A2 following removal of columns A1 and A2. This confirms the previous observation with respect to beam deflection in the long direction that Option A? provides the optimal discretization to achieve an acceptable response. The adopted discretization can capture nearly the same deflection profile obtained with 12 elements in addition to saving about 30 % in computational time.

Cross-Section Integration
Once the optimal number and discretization of the beams were defined, three different options to model the crosssection of the beams were compared: the cases of 23, 42 and 52 integration points to model the reinforcing bars and cover and core concrete were considered. The accuracy of each option in predicting the distribution of stress/strain in the cross section was evaluated by comparing the variation of the axial force in a critical beam following the removal of two columns. Note that the axial force is computed by considering the average stress on the cross-section. The simulation consists of first applying the gravity load followed by removing column A1 at time t = 2 s and column A2 at t = 8 s. Figure 8 shows the time history of the axial force in beams A2-A3. Interesting fluctuations in the axial force in the beam are observed-the beam transitions from compression to tension after the removal of column A1 and then reverses back to compression after the removal of the second column. Although all three options give comparable responses after the first column removal, the section with 23 integration points overestimates the axial force following the second column removal. Since the results obtained with 42 and 52 integration points are nearly identical, the option with 42 integration points was adopted, representing an acceptable balance between solution accuracy and computational efficiency.

Modeling of Floor Slabs
A separate study was then conducted on slabs for defining the optimal mesh refinement for shells, modeled as fournode layered shell elements in LS-DYNA. Each of the 10 layers in the cross section-203 mm thick-is assigned a specific mechanical property, representing concrete and steel: concrete properties are those reported in Table 2, while  for steel are the same adopted for beams and columns. Reinforcement is modeled through a smeared area corresponding to 13 mm dia @203 mm reinforcing bars at the top and bottom (in both directions). For the sensitivity analysis, mesh sizes of 300 9 300, 500 9 500 and 700 9 700 (mm 2 ) were investigated. Deflected shapes and axial force history along the mid-line of the slabs under gravity loads and after two column removals-columns A1 and A2-were assumed as the basic loading scenarios to compare the accuracy/acceptability of the three different mesh sizes. Figure 9a shows the deflected shape along the mid-line of slab panel A1-B1-B3-A3 in the long direction under gravity loads. The comparison at points of maximum deflection suggests that 500 9 500 mesh size gives a solution fairly close to 300 9 300, in addition to savings in computational time. This observation is confirmed after examining the deflected shape following the removal of columns A1 and A2, as shown in Fig. 9b. Figure 10 presents the results of the same analysis for slab panel A1-C1-C2-A2.
The three shell refinements were also compared in terms of axial force variation in the slab. Figures 11a and 11b show axial force-time history at the mid-line of slab panel A1-A2-B1-B2 in long and short direction, respectively, under gravity load between 0 and 2 s and then after removing      column A1 at 2.0 s and after removing column A2 at 5.0 s. The difference in the three discretizations becomes significant after the second column removal, especially in short direction, where the 700 9 700 mesh size underestimates the force while the other two solutions give nearly the same values. This difference is confirmed in plots (a) and (b) of Fig. 11, showing the variation of the peak axial force in the exterior slab panels in the long and short direction: A1-B1-B6-A6 and A1-F1-F2-A2, respectively-after columns A1 and A2 are removed. The resultant axial force, calculated along the mid-line of the slab, follows the expected trend when moving from the panel next to the removed column to those further away and confirms that the mesh refinement of 700 9 700 underestimates the magnitude of the axial force by about 12 % on average. Based on these findings and those reported previously on slab deflections, the 500 9 500 mesh was adopted, having the added merit of saving approximately 40 % in total computational time.

Modeling of Gravity Loads
In some progressive collapse simulations reported in the literature, floors slabs are not explicitly modeled to save computational time and modeling effort. Loads from the slabs are therefore applied on the beams. However, when simulating column removal scenarios the load re-distribution that occurs after each removal alters the way loads are transferred to the beams. A study was hence undertaken to examine load distributions that simulate load transfer from slabs to beams for different scenarios.
Structure self-weight can typically be specified in most structural analysis software through body force loads after both the mass density of the adopted materials and the geometry of the structure are defined. This feature allows the automatic evaluation of the actual self-weight when slabs are included in the model. When slabs are not modeled, the selfweight of slabs (a significant source of gravity loading on buildings) and additional live loads should be distributed on beams. The nature of this distribution is not straight-forward and depends on the configuration of the building in plan, the stiffness of structural elements and the sequence of column removals.
In this study three different models were compared. In all cases, beams and columns are modeled with Hughes-Liu elements, having characteristics derived from the previously described sensitivity study: 10 elements for each beams based on discretization A?, 5 evenly distributed elements for columns and fiber sections with 42 integration points. The first model represents the ''exact'' solution incorporating two-way 203 mm thick slab, modeled with 500 mm 9 500 mm shells connected to beams with shell edge-to-surface contacts. The second model is a frame model, where the self-weight of the slabs is evenly distributed on beams, representing the solution often adopted in structural design. The third model is also a frame without slabs, but a modified load distribution is used, derived to obtain the same beam deflection and the same magnitude of shear forces calculated in the so-called ''exact'' model, resulting in the definition of a specific load distribution for each beam. The contribution of the slab in terms of stiffness is not considered in the frame models. Figure 12a and 12b show the deflection of beams A1-3 under gravity loads and after the loss of column A1, respectively, for the three models mentioned above.
It is observed that under gravity loads the frame model with uniform loads applied on the beams considerably overestimates the deflection of beams, whereas that with the modified distribution, as expected, gives the same deflected profile as the baseline ''exact'' model. Figure 12b shows the deformed shape of the same beams after removing column A1. It can be seen that the profile representing the frame model with uniform loads is missing because, right after column A1 is removed, beams A1-2 collapses under the imposed gravity load. Next, comparing the vertical displacement of node A1 using the frame model with modified load distribution it is seen that the predicted displacement is highly inaccurate. This is attributed to two factors: first, the distribution of loads based on a target deflection (shape and magnitude) and shear distribution that result from gravity loads leads to a concentration of loads in the plastic hinge zones and consequently to very large displacements after column removal, secondly, unlike models which incorporate slabs, frame models are unable to adequately account for alternate load paths and force redistribution after sudden column removals.

Grid Beam Models
Efficient slab models are not commonly available in structural analysis software thereby prompting engineers and researchers to seek alternative schemes using available beam models. Previous studies (Nurhuda and Lie 2004;Tian et al. 2012) have identified some configurations of beam grids suitable for the analysis of RC flat-plate structures. Sasani et al. (2011) used a two-directional grid of beams (similar to Fig. 13a) with modified torsional properties and an effective width to capture nonlinear effects in the slab region and reported satisfactory results.
To verify the efficiency and accuracy of alternative modeling methods for simulating the contribution of slabs in progressive collapse analysis, the nonlinear response of some simplified grid beam models under gravity loads and column removal scenarios was investigated. The validity of four different grid beam configurations, as depicted in Fig. 13, was verified by comparing their responses with the so-called ''exact'' solution (the response obtained with a model that simulates the slab using layered shell elements with nonlinear material properties). Additionally, for each configuration a separate parametric study was also conducted, wherein finer discretizations and modified stiffness properties were considered. The validation included comparisons of deflection (shape and magnitude), shear and axial forces in the beams and axial forces in the columns, under both gravity loads and column removal scenarios, when large deformations and pronounced nonlinear behavior are expected.
The first analyzed configuration has two-way grid beams, as shown in Fig. 13a: the cases of 1, 3, 5, and 7 evenly distributed girders for each panel were analyzed. In all cases grid beams were assigned the same depth-203 mm-and the same reinforcement ratio as the slab. When compared to the ''exact'' model that incorporates nonlinear slabs, the twoway grid beam model produces acceptable results when comparing axial forces under gravity loads, however the errors in peak displacements range from 10 to 25 % (in the region of the missing column lines) and even larger errors are noted in predicting shear forces in the beams. Increased refinement of the grids and modifying the stiffness properties of the equivalent beams do result in improved predictions but the overall effort in identifying a suitable grid and stiffness modifiers was deemed too cumbersome to merit further investigation.
Both the second and the third configurations, Figs. 13b and 13c, represent an attempt in redistributing loads on the main beams based on the column configuration and expected load transfer and both these configurations were evaluated for 3, 5 and 7 grid beams. The second configuration is that of one-way grid beams, running in the short direction. When tested under gravity loads, this configuration with 7 grid beams in each panel gave acceptable approximations of both resultant forces and displacements. However, when evaluated for the scenario involving a column removal, the errors were significant (with displacements being overestimated by a factor of 3 or more). Likewise, the third configuration, wherein one-way grid beams run in the short direction in the exterior panels and in the long direction in the interior panels, was also unsatisfactory in predicting either displacements or internal forces to an acceptable degree. The fourth configuration, depicted in Fig. 13d, basically represents a modification of the first configuration (a), where elastic diagonals were added at the corners of the panels. The overall arrangement of the grid beams was conceived with the idea of reproducing in each panel the distribution of bending moments and the formation of yield lines observed in slab panels. A separate parametric study was conducted to calibrate the number of grid beams and the stiffness of diagonals to produce acceptable agreement with the ''exact'' model under large deformations. Though this configuration provided the best estimates of displacement, the model was incapable of reproducing correct internal forces.
This part of the study demonstrated the difficulty in simulating true slab effects using a grid of beams. While calibrating the stiffness properties of the beams to produce acceptable results under gravity loads appeared feasible, the task of recalibrating the models for each successive column removal makes this approach unreliable and unfeasible.

Energy-Based Approach to Assess Proximity to Collapse
This final section presents an energy-based approach to assess the proximity of the building to progressive collapse. Starting from the method proposed by Dusenberry and Hamburger (2006), in which the potential of a structure to survive disproportionate collapse is determined by an indicator proportional to the kinetic energy (KE), a more comprehensive approach is developed in this study based on a careful consideration of the energy contributions of the different structural components to the redistribution of forces following each column removal. A ''collapse index'' I c is derived that tracks the variation in the system energy following each sudden removal of a critical element up to the collapse of the structure. The collapse index is used in parallel with another damage measure, the displacement ratio of the model that is defined as the ratio of the peak downward displacement of any removed column to the floor height. It is postulated that a collapse condition is imminent when the displacement ratio exceeds 0.6. Though the numerical solution may yet converge to an equilibrium condition, it is evident that the physical damage from a displacement that exceeds half the story height is likely to be severe and irreparable.
The proposed approach recognizes the fact that the system KE and strain energy (SE) are reliable indicators of the motion of a system caused by the sudden loss of bearing elements and also of the energy absorbed by the system as new load paths are established and equilibrium is being restored, if possible. It is also well-recognized that the variation of these two quantities is not uniform across the entire structure, but mostly involves the local area where the structural elements are removed. For this reason a separate evaluation and a subsequent comparison of global and local energy change is assumed as an indicator of how much the local failure of the structure affects the stability of the overall system. Consequently, the collapse index I c is defined as follows: In the above expression, j global and j local represent the ratio of the KE to the SE calculated over the entire system and over the local portion involved in the collapse, respectively. In the present study, these energy quantities are directly computed in LS-DYNA. However, the computation of SE and KE can be accomplished using well-known relationships as follows: The following notations are used: L = length, N, V, M, T = axial force, shear, moment and torsion, respectively, E, G, A, A*, I, J = Young's modulus, shear modulus, crosssectional area, shear area, moment of inertia and polar moment of inertia, respectively, {r}, {e} = vectors of stresses and corresponding strains, nt, nr = number of translational and rotational degrees-of-freedom, m j , I h = translational mass and rotational inertia, _ u; _ h = translational and rotational velocities, respectively. When following the APM, I c indicates the proximity to collapse of the structural system during a sequence of column removals. Ideally, at the ultimate condition I c = 1 though it was necessary to introduce the displacement ratio as an additional criterion to assess the severity of the damage state.
The proposed method was applied to three buildings (1-, 3and 6-story), and both the models with and without floor slabs were analyzed. A strain-based failure criterion was also assigned in which the rupture strain of reinforcing steel was specified as 0.15. However, as noted previously, this failure strain is not a reliable measure since a modified post-peak slope was specified for the reinforcing steel material to avoid localization. Yet, it does provide an indication that the ultimate stress in the bar has been exceeded. Figures 14a and 14b show the variation of I c and the corresponding displacement ratios for three different column removal sequences (paths A-C) for the single-story building for both the bare frame model and the full model with slabs. Likewise, Figs. 15a and 15b present the change in the collapse index with successive column removals for the 3-and 6-story models. Based on the results shown in these plots, the following facts emerge: (1) The collapse index, I c , shows an increasing trend with each subsequent column removal hence validating the normalized energy-based approach for defining the onset of progressive collapse. (2) Incipient ''collapse'' is initiated when the displacement ratio exceeds 0.6. This corresponds to an I c value close to unity for the 1-story building and approaches a ratio of approximately 0.8 for the 3-and 6-story buildings. This damage state was also accompanied by simultaneous rupture of reinforcing steel bars in critical locations where the largest deformations are localized.
(3) In all the analyzed cases, the bare frame model (without slabs) approaches the collapse state with fewer column removals than the model incorporating floor slabshighlighting the role of slabs in finding alternate load paths after column removal and contributing to progressive collapse resistance of the system. (4) An assessment of the plots shown in Figs. 14 and 15 indicate that progressive collapse is ''sequence-dependent'' and that there exists a critical sequence of column removals that results in the minimum number of removed columns before collapse is initiated.
In the case of the model incorporating floor slabs, path A corresponds to the sequential removal of the following columns: A1-A2-B1-A3, path B corresponds to removal of C2-C1-B1-A1-A2 and path C refers to A3-A2-A1-B1. For the bar frame models, the removal sequences were as follows: path A: A1-A2-B1, path B: C2-C1-B1-A1, path C: A3-A2-A1. For the 1-story frame model, paths A and C require only three column removals before collapse is initiated whereas path B involves the removal of four columns. Similar observations are evident from the response of the full model incorporating slab elements (Fig. 14a) though an additional column (compared to the frame model) needs to be removed to initiate failure. The results shown in Fig. 15a and b are for the models in which the floor slabs are included. Results for the case of the bare frame are not shown here and trends similar to that for the single-story building were observed. The proposed method, therefore, in addition to providing a measure of the reserve capacity of the structure, can also be used to identify critical load-bearing elements in the structure that provide the greatest resistance to progressive collapse.

Conclusions
The results presented in this study offer useful guidelines on modeling and simulation of progressive collapse of RC frame structures within the context of the APM. Using a simple case study of a RC frame structure, various issues in modeling of materials and elements are presented with a view to providing practical insights into progressive collapse simulations. The importance of adequately modeling the floor slab is highlighted and the often employed approach of using a grid of beams to model the slab system is shown to be impracticable for simulating collapse within the APM framework that involves the removal of column elements.  Finally, a new collapse index is proposed that can be used to assess both the damage state and the reserve capacity of the system and can also serve as a means to identify critical load-bearing elements in the structure that provide the greatest resistance to progressive collapse.
One of the objectives of the study was to consider simplified approaches for progressive collapse analysis. Consequently, though an advanced finite element software was used in the simulations presented in this study, beams and columns were modeled as line elements with cross-section integration (available in both open-source and commercial software used in engineering practice) and slabs were modeled using layered shells with smeared reinforcement. Some structural analysis software may not provide convenient options for sudden column removal. If dynamic analysis is feasible, the application of short-duration pulse loads corresponding to the axial force in the removed columns is a viable option. Another alternative is to carry out a static analysis in the absence of the removed columns but to magnify the applied loads to consider dynamic effects. But both these approaches become more involved and complex for multiple column removals.
Another critical aspect in the development of a simulation model is its validation using experimental data. The availability of appropriate experimental data to validate the response to sudden column removals in a building is extremely limited. A number of recent tests carried out by researchers in China (Yi et al. 2008(Yi et al. ,2014Xiao et al. 2015) provide an initial starting point for such validations but additional experimental data on large-scale experiments is still needed to calibrate and assess the adequacy of the proposed modeling schemes.