Abstract
This study investigates the behavior and the load-bearing mechanism of a typical flat slab with rectangular panels in several scenarios including the removal of a corner, penultimate, and internal columns. The scenarios are rather similar to those used in the conventional evaluation of the progressive collapse potential; however, application of the uniformly distributed loading over panels adjacent to the removed columns was not limited to twice the value of the initial load. Thus, load-deflection curves were drawn up to the point in which a great number of longitudinal slab bars ruptured. Introducing 5 stages on each curve, finite element outputs on concrete cracking pattern and rebar stress state were presented. A significant increase in the stresses along the diagonals of the slab panels accompanied by bar ruptures around columns adjacent to the removed column proved contribution of an important load-bearing mechanism in addition to the behavior called “quasiframe action.” Consecutive rupture of bars showed formation of a zipper-type collapse mode as well as a great tendency to transfer load share of missing column mainly along shorter direction of slab panels. Moreover, the findings indicated that the slab damaged zone could exceed the panels under uniform overloading.
1. Introduction
Progressive collapse is defined as gradual spread of local damage from one element to another until collapse of the whole or a part of a structure which is disproportionate to the initial damage [1]. This event can be triggered in frames, space structures, and domes [2–4]. Progressive collapse can occur in a variety of forms such as domino (e.g., a series of adjacent buildings overturn on each other), pancake (for example, upper floors fall down on the ones below), instability (e.g., structural elements buckle consecutively), zipper-type (for example, suspension bridge cables tear away one after another), or combined collapse modes [5]. Although the behavior and propagation pattern differs from one case to another, many structural collapse incidents, including the twin towers of the World Trade Center (WTC), Pipers Row Car Park, and Sampoong Department Store, were considered to mainly include pancake collapse mode [6–8]. Concerns about the possible consequences of such disasters in very commonly used flat slab system led the researchers to evaluate their behavior after column loss events [9]. These notional scenarios might represent the impact of a wide range of unusual loading conditions [10, 11]. Recently, Vieira et al. [12] with employing the special techniques conducted four quasistatic experiments on half-scaled specimens considering a central column loss scenario in a perimetral portal frame.
There are a vast number of studies on the progressive collapse potential of RC frames and subassemblages [13, 14]. Fascetti et al. [15] proposed a new procedure derived from nonlinear static and dynamic analyses for comparing the relative robustness of RC frame buildings against progressive collapse. Genikomsou et al. [16] conducted nonlinear finite element analyses of reinforced concrete slab-column connections under static and pseudodynamic loadings to investigate their failure modes in terms of ultimate load and cracking patterns. Du et al. [17] investigated the performance of a RC spatial frame structure and found that short-span beams fail earlier than long-span beams on asymmetrical structures subjected to progressive collapse. Yu et al. [18] showed that progressive collapse resistance of RC beam-slab assemblies subjected to perimeter middle column removal scenarios was affected by loading positions in testing. The effects of supporting column stiffness and restraint due to building continuity on the behavior of beam-slab assemblies in rather similar loading events were also investigated by Yu et al. [19]. Through the proper test setup and instrumentation, Lim et al. [20] elucidated interactions between the columns and beams in RC substructures under loss of exterior and corner columns.
Although experimental study on the overstrength capacity of slabs has a relatively long history [21, 22], ease of access to updated laboratory equipment has launched a new wave of research, especially in recent years. For instance, Muttoni [23] investigated the effect of slab reinforcements on preventing progressive collapse, and the scaled models of Qian and Li [24] provided information on the overall shape of load-displacement curves and slab cracking pattern. Valuable 1/3 scale experimental studies of Ma et al. [25, 26] evaluated the alternate load path and bearing mechanisms of RC flat plate substructures subjected to corner column, edge-column, and edge-interior-column removal scenarios. However, limitations and costs of conducting full-scale laboratory tests have aroused the need for other alternatives. Latest development of analytical tools has paved the way for a more accurate assessment of the structural behavior in support failure scenarios [27]. For example, Bredean and Botez [28] examined the performance of slab structural element and its resistance mechanisms in addition to beam’s arch action which is commonly considered in progressive collapse analysis. Pham et al. [29] also studied the behavior of beam and beam-slab subassemblages with different loadings and showed that removal of a penultimate column could be even more critical than a corner column loss. Russell et al. [30] evaluated the effect of some design parameters on the nonlinear behavior of flat slabs considering different column removal locations. However, as the panels connected to the removed column were only overloaded up to initial uniform loading on the entire floor, no limit states such as extensive failure of slab bars were reported in that study.
Summarizing past studies suggests that full-scale models of typical slab floors have been less studied. Some limitations in numerical modeling have also kept the flat slab behavior due to abnormal events little known to the engineering community. In other words, the pattern of stress redistribution and damage propagation in concrete and slab bars after exposure to large deformations corresponding to support removal have not received much attention. Additionally, the effect of this damage on load-displacement curves and floor failure modes has not been well investigated yet. For this purpose, three separate finite element (FE) models of a typical floor of reinforced concrete (RC) flat slab consisting of similar rectangular panels are developed in Abaqus software [31]. All the models study the effects of increasing uniformly distributed loading on the slab behavior, but the location of overloaded panels differs from one to another, depending on the missing column location (interior, corner, or penultimate). Completing the verification process, load-deflection curves are presented and for better comparison of the results, five stations (starting from yielding of first series of slab bars and ending by chain rupture of a significant number of them) are considered on each curve. Then, after removing mentioned columns, cracking pattern and stress condition of the slab bars were compared. The load redistribution mechanisms and failure mode are also discussed.
2. Numerical Modeling
2.1. Design Assumptions and Geometry
According to Figure 1, a RC slab floor of 200 mm thickness with repetitive rectangular panels and two reinforcing layers was designed in accordance with ACI318-14 [32]. Each layer contains continuous longitudinal bars of 12 mm in diameter at a center distance of 250 mm. The slab represents the floor of a typical residential building and is supported by 500 mm square RC columns. The slab overhang length is 100 mm in all edges. The design dead load includes self-weight of the structural floor in addition to an overburden of 1 kN per square meter corresponding to architectural flooring. The live load was considered supposing that lightweight partitions are used. It has been attempted to comply with the minimum reinforcement requirements; therefore, providing some additional bars along the slab edge at the upper reinforcing layer of slab and in the vicinity of exterior columns (excluding corner columns) was inevitable. It is also notable that a larger number of additional bars were embedded along both directions of the slab around interior columns C2 and B2.

(a)

(b)
2.2. Loading Protocol and Column Removal Scenarios
The structural damage scenario is defined by three separate analyses. As seen in Figure 1, these scenarios are entitled as RP, RI, and RC. The common initial letter R denotes the rectangular geometry of slab panels, whereas the letters I, P, and C, respectively, demonstrate that the model lacks a corner, a penultimate, or an internal column.
The loading process in each of three scenarios consists of two main steps. In the first step, a uniformly downward distributed load, , is applied to all slab upper surfaces of the column-removed structure. This load is ramped linearly, from zero to its full value during 1.5 seconds:where LL and DL are live and dead loads on the slab, respectively. The coefficients of this gravity load combination have been selected based on the recommendations of design guidelines against progressive collapse for floor areas away from removed column [33]. It is notable that DL includes the self-structural weight plus 1 kN/m2 corresponding to the flooring weight, while LL is considered equal to 2.5 kN/m2.
In the second step, an increasing distributed overload is incrementally applied to the panels immediately adjacent to the removed column and the behavior of the concrete slab and steel bars during the loading process is evaluated. Note that the application of uniformly distributed overload is not limited to the value of initial loading and continues up to the point in which a great number of slab bars fail. Applying proper nonstructural mass to the concrete elements and using a time period of 6 seconds, quasistatic loading rate is ensured (see Section 3.2).
3. Specifications of the FE Model
In the present study, the response of concrete slabs considering a support failure is simulated in Abaqus/Explicit. The explicit solution strategy is inherently based on dynamic equilibrium. Hence, it needs to reasonable computational memory and time as huge stiffness matrices are not required. Moreover, it is capable of modeling highly nonlinear effects without the common convergence problems which often emerge in implicit formulations.
3.1. Simulation of Material Behavior
The compression behavior after the initial linear elastic stage (approximately equivalent to 30% of the 28-day concrete compressive strength, i.e., ) was generated using equation (2) [34] and is shown in Figure 2.where is the ratio of compressive strain to concrete cracking strain and is a stress calculated as a function of compressive strain over relevant range. k is a numerical constant which is considered equal to 2.15 for grades C25 and C30 concrete [30].

The tensile concrete behavior was assumed to be linearly elastic up to the point of cracking, but was followed by a nonlinear softening stage in accordance with equation (3) [35].where , , , and are the concrete tensile strength (about 10% of its compressive strength), the varying tensile strain, the initial modulus of elasticity, and the strain corresponding to the occurrence of tensile cracks, respectively. It is worth noting that specifications of grade C30 concrete were used in equations (2) and (3).
Using concrete damaged plasticity (CDP) model [36], the behavior of concrete after cracking could be estimated from the zone of plastic strain on the slab. For this purpose, indices of tensile damage and compressive damage were used to consider the reduction of concrete stiffness after cracking. According to equations (4) and (5), these indices depend on the level of tensile or compressive stresses applied to the material and are equal to zero and 1.0 for the undamaged and fully damaged concrete, respectively [27]:
A number of parameters should be also given in order to implement the CDP model in Abaqus. The parameters were selected on the basis of Abaqus user manual [31], Russell [27], and Hafezolghorani et al. [37] and were used in this study (see Table 1) to convert the uniaxial stress-strain relationship for compression and tension into the two-dimensional yield surface.
Modeling the behavior of steel material in this study was primarily based on the nominal stress-strain curve of grade 60 steel with a yield stress of 420 MPa. Although the erosion phenomenon (deletion of largely deformed elements) was not directly simulated, rupture of bars was taken into account. To simulate bar rupture at the ultimate strain, a ductile damage model was used, which made it possible to determine sequence of bar failures. Fracture strain, strain rate, stress triaxiality, and fracture energy were needed for this option. Fracture strain was set as 0.119. The minimum and maximum strain rate values were 0.0001 and 10, respectively. Considering that the bars were modeled using truss elements, the stress triaxiality limit values were defined by +0.3333 and −0.3333. The fracture energy was set to 0, so the bar rupture occurred instantaneously when the fracture strain was reached [28].
3.2. Specifications of the Full-Scale FE Models
Figure 3 depicts geometry of the full-scale FE models built for this research. The length and boundary conditions of columns are selected based on the assumption that the contraflexure point forms at the midheight of upper-story columns. The lower surfaces of the lower columns are restrained against both rotation and translation, while the top surfaces of upper columns are restrained in a way that rotational degree of freedom remains free (see Figure 1). Additionally, the interfaces of columns and the slab are faced to tie bond, whereas slab edges are exposed to no restraints.

Concrete with 28-day compressive and tensile strength of 30 and 3 MPa, respectively, was modeled through 8-node cubic elements with reduced integration points (C3D8R). Steel reinforcements, which are considered from grade 60 steel with Young’s module of 200 GPa, Poisson’s ratio of 0.3, and yield and ultimate stress of 420 and 560 MPa, respectively, were modeled using T3D2 truss elements.
Clearly, attention should be paid to the mesh size so that the equilibrium equations would be solved in a reasonable time with acceptable accuracy. A coarse mesh does not show the local decrease of stiffness, due to yielding or cracking, and keeps some parts of the model undamaged or too stiff. Alternatively, in situations with localized areas of high tensile stress, very fine mesh sizes result in narrower crack bands and inaccurate distribution of stresses and strains. It is seen that fine cubic elements (of dimension 25 mm) were used for the columns. Relatively similar elements were used to mesh slab panels and slab-column intersections, but their thickness was approximately reduced to 6.67 mm. Such a fine mesh size for elements has been chosen based on the sensitivity analysis and the recommendations of previous studies in the field of slabs under column loss scenarios. Moreover, it has been confirmed that application of the same mesh size in the slabs brings about acceptable accuracy as well as reducing the possibility of hour-glassing phenomenon [27]. Note that it took on average 16 days for the supercomputer of University of Kashan available at the High-speed Processing Center (HPC) to perform each of the mentioned analyses (6 processing cores were assigned to simultaneously do the task).
The total number of cubic elements used for slabs and columns in each model was 7809120 and 792000, respectively, and 36160 elements were applied in modeling of bars. Considering that it was too time-consuming to analyze this FE model with such a large number of cubic elements in a commonly used computer, the supercomputer of University of Kashan available at the High-speed Processing Center (HPC) was assigned to perform the analyses. The analysis process was also made more efficient by increasing the mass matrix. The stable time increment is inversely proportional to the square root of the mass terms; however, mass scale was set in a manner that the ratio of kinetic energy to internal energy remained under control (less than 10%).
3.3. Validation of the Developed FE Model
The FE model of this research was validated with the findings of an existing experimental program, where Russell et al. [30] evaluated the behavior of a scaled two-panel flat slab supported by 5 columns (lacking a corner column). The slab was under a uniform downward increasing load on its entire top surface and the pattern of crack propagation on slab bottom surface at the end of the loading process was reported. Figure 4 validates capability of the developed FE model in predicting crack pattern of the support-removed slabs as there is a good agreement between numerical results and experimental data. In fact, in both cases, the major right-panel and left-panel cracks occurred along and inclined to the direction of longitudinal bars, respectively.

(a)

(b)
To ensure the accuracy of the FE model in predicting load-bearing capacity, one can also monitor the axial force of the three supports adjacent to the removed column location during the overloading process and compare them with laboratory data. As can be seen in Figure 5, the results agree well with the experimental curves. The observed difference might be attributed to limitations of the numerical model in simulating sudden drop in the concrete stiffness after crack development. The sensitivity analysis approved that an element size of 25 × 25 × 6.67 mm could be appropriate for the slab elements. The ratio of each slab span to maximum element size was 2100/25 = 84 in the validation model. As the same mesh size has been also used for the full-scale numerical models, the ratio increased up to 4000/25 = 160 and subsequently the accuracy of main FE models was approved.

4. Results
In this section, the results of the separate analyses RP, RI, and RC are compared from three aspects of load-deflection curves, the stress state of the bars up to slab near-collapse state, and concrete crack pattern at a moment corresponding to first bar yielding.
4.1. Load-Deflection Curve
The changes in the amount of uniformly distributed load applied to unit area of panels adjacent to the column removal location as a function of the slab maximum deflection (normalized relative to its thickness) are briefly referred to as “load-deflection curve” and are shown in Figure 6 for the three models RP, RI, and RC. The maximum displacement of the slab occurs at the location of removed columns. To allow for a robust interpretation of the results, five hypothetical stations are considered on each curve. The curve shape, however, follows a relatively similar trend in all 3 cases. There is an elastic region with a constant stiffness at the beginning of each curve. By increasing the overload induced deflection in the panel (s), cracks grow in the slab and the stiffness (the slope of the curve) gradually decreases until the first series of slab bars around the columns adjacent to the removal location yield (station A). By continuing the loading process, not only do the cracks propagate extensively and open considerably, but also more bars experience strains beyond their yield limit. This results in a drastic change in the slope of the curve. The stiffness reaches zero and becomes negative before station C due to the failure of one or more slab bars (simultaneously or consecutively) at station B. In fact, panel(s) facing increasing downward overload cannot carry existing applied loads anymore. Growth of deformations up to hypothetical station D is also associated with the continued softening behavior in the curve as more bars reach failure state. Finally, where the chain failure occurs in a significant number of slab longitudinal bars or the entire reinforcements around one or more columns rupture, very large deformations occur, and slab collapse seems to be almost imminent (station E). The normalized peak load (NPL), as well as those corresponding to the stations A to E, is presented in Table 2. By comparing NPL of the three investigated models, it can be seen that the slab was able to withstand the lowest and highest amount of overload in the models lacking corner and interior columns, respectively, meaning that RC was the most critical case. It should be noted that the local fluctuations in the curve path can be somehow attributed to the opposite effect of strain hardening of some bars accompanied by rupturing of others.

4.2. Stress State of Slab Bars
The stress state of slab bars is discussed in the following to determine in what direction and how the loads on the floor are redistributed after column removal. The findings can also provide valuable information on the potential failure mode of the slab and provide a basis for better design of structures against progressive collapse. The study, however, excludes modeling of behaviors such as concrete spalling at the upper surface of the slab due to ripping out of negative bars adjacent to the column.
4.2.1. Model RC
After removing corner column A1, the stresses of the bars located at the panel A1-B2 (rectangular panel with an assumed diagonal passing through the two columns A1, B2) undergo significant changes. Table 3 summarizes the stress state of the slab bars around critical areas of the model RC, i.e., in the vicinity of columns A2 and B1. Except for some slab bars which reach yield limit near the two mentioned columns, the stress of other bars remains within the elastic range at station A. Increasing the applied overloads, the strain rises in the bars until at station B; seven slab bars (consisting of four add bars and three continuous bars) at the location of column A2 experience failure almost simultaneously. In station C, the first failure of slab bars around column B1 is observed, while more bars around column A2 enter inelastic phase and failure of bars even progresses into the lower reinforcement layer of slab. The slab continuous bars fail one after another at the end of stations D and E. The graphic view of the stress state of slab bars in Figure 7 provides evidence on initiation of zipper-type collapse pattern in the slab bars because the reinforcement failure penetrates inward the panel at the end of station E.

(a)

(b)
Redistribution of loads to adjacent columns located in the same frames formerly including the removed element through bending in adjoining beams is a particular case of frame action and is one of the most important load-bearing mechanisms of moment frames against progressive collapse [38]. According to Figure 8, a rather similar mechanism (called “quasiframe action”) has also occurred to a large extent in the studied flat slab because the slab bars around the columns A2 and B1 have entered inelastic deformation phase after removing the corner column A1. However, the priority of yielding and rupture of slab bars around column A2 over those located in the vicinity of column B1 shows that the mentioned mechanism mainly tends to develop along shorter spans of the flat slab.

(a)

(b)
4.2.2. Model RP
The pattern of change in stress of the slab bars at different stages of the model RP is largely similar to what was previously stated for the model RC, for the quasiframe action mechanism also plays a key role against the progressive collapse hazard in this case. Transferring load share of the removed element along the diagonals of the slab panels is certainly another mechanism which contributes to load-bearing; as shown in Figure 8(a), considerable stress around the columns A2 and C2 which are not in the same axis with missing column B1 is an evidence to this fact. Moreover, rupture of bars (corresponding to station B) not only starts from the upper layer of slab bars in the vicinity of the corner column A1, but also continues up to the point in which all bars around this column fail (Figure 9). Thus, it is not far from expectations that the occurrence of a punching shear around column A1 precedes other locations. It is worth noting that the observed progressive failure of bars at station E compared to station D could be regarded as a symptom of zipper-type collapse mode.

A glance at Figures 8 and 9 also demonstrates that most of continuous rebars of severely deflected slab panels along axis 1 have already yielded. Such an event is a symptom of tension membrane action (TMA) which could be regarded as an additional load-bearing mechanism.
4.2.3. Model RI
The stress state of the flat slab bars in the vicinity of various columns at station E of the model RI is demonstrated in Figure 10. In this case, yielding of slab bars which bridge over the tributary area of removed column B2 provides evidence of TMA along both directions. In addition to diagonal transfer of loads to columns A2 and C2, compressive membrane action (CMA) also could be distinguished in model RI as NPL values go far beyond the predictions of the yield line theory. It is noticeable that axis C can provide significant lateral restraint for parts of the deformed slab panels. Even if the support conditions of slabs are such that lateral movement at the edges can occur freely, a ring of compression supporting tensile membrane forces in the inner region of the slab might form [39].

Upon reaching stage E in the load-deflection curve, all slab bars around the columns B1 and B3 fail. The coincidence of bar failures with the crushing of concrete slab around the mentioned columns can be considered as an indication of the possibility for punching shear. As the slab is completely detached from each vertical load-bearing element, its load share is inevitably transferred to other slab areas where their bars themselves are on the verge of failure, so the collapse of the floor is imminent. Since the experience of collapse events such as WTC towers has clearly shown that the impact load induced by collision of falling floor(s) might be generally intolerable for the lower ones [6], the vertical propagation of damage toward the foundation is conceivable. It is also notable that, among the columns adjacent to the removed column B2, the more critical condition occurred around columns B1 and B3 (and not A2 and C2). These exterior columns not only are located in frame 2 which includes the missing column B2, but also lie along shorter direction of the panels facing the overload. Hence, slab bar detailing especially at discontinuous edges of the exterior slab panels along their long direction and in the vicinity of the columns is found to be of crucial importance to mitigate the progressive collapse risk.
4.3. Concrete Cracking Pattern
Abaqus software output in the field of concrete plastic strain can be used to predict cracking pattern, as presented in Figure 11, at a stage corresponding to the first bar(s) yielding (station A). A brief look at the reported pattern for the models RC and RP reveals that the positive cracks (located at the bottom of the slab) of overloaded panels tend to stretch vertically into the neighboring panels. These cracks propagate along the y-direction because, in this case, a shorter path and subsequently less energy are needed to form the ultimate failure mechanism. For instance, in model RC, the crack entered the panels A2-B3 instead of the panels B1-C2, which in turn emphasizes the possibility of exceeding slab collapsed area beyond overburdened panels. Though the positive cracks in the model RC occur only on one side of the removed column and do not show forking. Unlike the two mentioned models, positive cracks around the removed column form only in one direction (x) of model RI.

The findings indicate that the negative cracks on slab surface, especially in the model RI, were mainly limited to the area of overloaded panels and developed mostly along slab continuous edge (adjacent to panels C2-D3 and C1-D2) which was significantly constrained to rotation. However, these cracks at the slab discontinuous edges have shown a strong tendency to move around the edge columns, confirming the possibility of occurrence of punching shear failure. Note that the end panels C2-D3 and C1-D2 are not shown in Figure 10 as they do not suffer significant cracks.
5. Conclusion
In this study, the behavior and the loading mechanism of a RC flat slab with repetitive rectangular panels was investigated in separate corner, penultimate, and internal columns loss scenarios. However, the applied load over the panels adjacent to the removal location was not limited to twice the value of initial load. After validation of FE models, load-deflection curves corresponding to increasing uniform overload at the location of panels immediately adjacent to the removed support were extracted until failure occurred in a significant number of slab bars. Introducing five stations on each curve, Abaqus software output in the field of concrete cracking pattern and stress state of slab bars were presented. Some of the most important results are as follows:(i)The pattern of overall changes in the load-deflection curves is rather similar in all the analysis cases. There is an elastic region with a constant stiffness at the beginning of the curve which is followed by a gradual decrease in stiffness. As the loading process continues, due to the failure of one or more slab bars (simultaneously or consecutively) the curve slope first reaches zero and then a general softening phase is seen as a large number of bars fail.(ii)The corner column loss is the most critical scenario because the slab is able to withstand the least amount of uniform overloading. In this case, the normalized peak load is about 85% and 89% of the values corresponding to the models RI and RC, respectively.(iii)Cracking pattern at first bar yielding of models lacking corner column (RC) and penultimate column (RP) indicates that the positive cracks of overloaded panels perpendicularly enter into adjacent panels and propagate along shorter span of slab panels. Such a path requires the least energy to form a collapse mechanism and also proves that the slab damaged (or possibly collapsed) zone can exceed overloaded panels. On the other hand, the growth of negative cracks at the edges of the flat slab around the exterior columns indicated the possibility for a punching shear mechanism.(iv)The results from model RP not only confirmed the possibility of punching shear phenomenon around the corner columns of the structure, but also, along with model RC, demonstrated that the load share of removed elements can be significantly transferred along the diagonals of the slab panels. Such a load redistribution mechanism collaborates with “quasiframe action” (the transfer of loads mainly to adjacent columns which are located in the same axis with the removed column) in enhancing structural resistance against progressive collapse. Some evidence on the formation of the mechanisms of TMA and CMA was also seen in model RI.(v)At the end of the model corresponding to an internal column loss, all the slab bars around the penultimate exterior columns B1 and B3 failed. These columns provide the shortest path for redistribution of loads. It is therefore conceivable that, to counter the vertical expansion of damage area (caused by partial collapse of upper floors) in flat slab structures, one can pay special attention to slab bar detailing at discontinuous edges of the exterior slab panels along their long direction and in the vicinity of the columns.(vi)Comparing the state of bars in the last 3 stations of load-deflection curves suggested the possibility of a zipper-type collapse pattern in slab bars.
As a suggestion for future work, one can develop analytical formulas for simulating different areas of load-displacement curves of the support-removed slab floors. The effect of drop panels on the behavior of slab floors subjected to overloading could also be studied.
Data Availability
Required data have been provided in the manuscript. Additional information can be available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.