CASCO: a simulator of load paths in 2D frames during progressive collapse

Modern structural design software can simulate complex collapse dynamics, but the main physical processes driving collapse propagation are often hidden among structure-specific details. As a result, it is still unclear which structural geometries and material properties should be preferred when approaching the design of a damage-tolerant structure. This manuscript presents a new approach to explore the relationships between structural geometry, local mechanical properties, and collapse propagation. The insight comes from a unique ability to trace the evolution of load paths during collapse, achieved by combining energy conservation with local mechanisms of plastic failure and a few simplifying assumptions. The method is implemented in a new simulator of collapse of 2D frames, called CASCO and programmed in MATLAB. Simulation results for reinforced concrete frames predict collapse loads and mechanisms in agreement with fully non-linear, dynamic simulations, while also providing a graphical description of the evolving structural topology during collapse. A first application of CASCO to mechanically homogeneous and heterogeneous frames, indicates certain evolutions in number and density of load paths during collapse that may be targetted to improve collapse resistance.


Introduction
Resistance to progressive collapse is the ability to withstand local damage without triggering a chain of failures with disproportionately severe consequences. The Alternative Load Path Method (ALPM) is an established framework to analyse progressive collapse [9,10]. The ALPM considers a model structure initially at equilibrium under characteristic loads, and then removes of one or more elements, usually columns, to represent accidental damage. The concept behind the ALPM is that each load has its own path, relying on internal forces, to reach the external supports. Accidental damage removes some paths, forcing the loads to alternative paths. Such redistributions might amplify the internal forces and trigger new failure events, prompting further redefinitions of paths and redistributions of loads. This process can cause avalanches of failures and compromise the stability of large portions of the structure. Going from this intuitive concept to actually tracing load paths is not straightforward.
The current tools for numerical simulations can realistically describe structural collapse using dynamic and nonlinear analyses. These simulations are typcially based on Finite Element Analysis (FEA) [15,18,26] or on the Discrete Element Method [22], and can also include stochastic analyses [35]. Detailed simulations, however, produce complex distributions of internal forces from which it is generally difficult to trace load paths, as illustrated in Fig. 1. Detailed simulations are the state-ofthe-art in structural engineering, and they are certainly needed when designing a structure professionally ow when planning its demolition. There is however scope for more conceptual, simplified analyses when training structural engineers in progressive collapse analyses or when approaching the first stages of a design. The value of conceptual analyses is to provide an insight into which structural features may lead to a better intrinsic ability to resist local damage; for example, which geometric ratios (e.g. beam lengths to column heights), materials properties, or combinations of strong and weak structural elements and connections would optimise damage tolerance. In this persepctive, a graph-based topological approach has been proposed in the early 1990's, which analyses the connectivity map of truss structures to identify critical sequences of element failures and provide a measure of vulnerability [1,33,34]. Conceptually similar approaches have been proposed since then, e.g. measuring resilience via measures or srtuctural robustness [6] or combinatorial analyses of collapse extents [5]. The challenge for such approaches is to account for loads redistributions during collapse [8].
Another class of conceptual models uses structural mechanics to provide analytical expressions, highlighting how certain design variables affect progressive collapse resistance. These approaches have been applied both to structural sub-components [14,30,36] and to full frames and structures [3,17,24]. They often exploit energy conservation to account for dynamics, consider either individual structural elements, subsets of them, or whole frames, and can consider various levels of detail in the mechanical description of structural elements and connections. The challenge for such approaches is to go beyond the first failure event after the initial accidental damage, and also account for damage propagation due to impacts between structural element. A few simulation approaches combine simplified analyses of local failures with Finite Element analyses of 2D frames, incorporating impacts between falling structural elements via kinematic rules and energy conservation [11,32]. None of these approaches, however, explicitly identifies and quantifies load paths.
This manuscript presents a new algorithm to trace load paths in 2D frames and follow their evolution during collapse. The algorithm is named Conceptual Analysis of Structural COllapse, or CASCO. It is implemented in MATLAB, open source, and is attached to this manuscript [20] Sect. 2 explains the new methodology, starting with the construction of initial 2D frame geometries and the parameters required as inputs. The section then explains how load paths are computed and how local plastic failures are sampled using an energy-based criterion. Energy conservation, propagation, transmission by falling debries, and dissipation, are central in CASCO, consistently with their recognized importance in progressive collapse analysis [9,21,31]. Section 2 ends with explaining how a sudden initial damage is applied. Section 3 presents new results obtained using CASCO. The section first discusses the validation of the new method against state-of-the-art, fully dynamic, nonlinear simulations of collapse of 2D frames made of reinforced concrete. After that, new results compare the response to local damage of "homogeneous" structures, with same mechanical properties everywhere, and "heterogeneous" structures, with disorder in the spatial distributions of beams and columns strength. The relationships between topology evolution and progressive collapse resistance emerging from the simulations, prompt the discussion of potential targets for structural optimisation.

Methods
This section describes the CASCO simulator, whose algorithm is in Fig. 2.

Input, classification, hierarchy, and load paths
CASCO considers 2D frames made of vertical columns and horizontal beams, with moment resisting connections and fixed to the ground. Bays and columns are discretised into an arbitrary number of beam elements (BEs) and column elements (CEs): see example in Fig. 3. For discussion in the rest of this manuscript, each BE and CE is followed by a unique numerical label, see for example BE71 in in Fig. 3.
Before applying any load, BEs and CEs can be subtracted by the full frame in Fig. 3, to create a more complex initial geometry, such as that in Fig. 4. This subtraction procedure can also involve the disconnection of certain BEs or CEs, to simulate the presence of structural joints. An example  Geometry of an initial frame, before assigning properties to beam and column elements (BE, CE). All images from CASCO are visualised using OVITO [29] is shown in Fig. 4, where a cross section between BE71 to BE75 is disconnected, turning the beams to the right and to the left of said cross section into cantilevers.
The materials and cross sections of the BEs and CEs are assumed to be rigid until yielding, and then perfectly plastic until a threshold failure strain. Therefore only 8 mechanical properties are required for each BEs and CEs: the tensile and compressive axial force at yielding, N p , the tensile and compressive displacements at failure, u (assigning displacement rather than strain at failure implies assuming localisation within a plastic zone), the positive and negative yield (or plastic) moments, M p (with positive moments assumed to compress the top of the BEs' sections), and the positive and negative ultimate plastic rotations p of a localised plastic hinge.
CASCO samples possible formation of plastic hinges at the left and right end cross sections of each BE: these hinges, being part of the same BE, have identical plastic properties. However, by assigning different plastic properties to different BEs, one can model complex distributions of cross sectional properties within a structure.
After defining the initial geometry, each BE and CE is assigned a uniformly distributed load per unit length and an initial input of kinetic energy (units of energy), both positive if pointing downward or to the right for BEs and CEs respectively. The initial kinetic energy may be used to simulate vehicle impacts or earthquakes. Some BEs can be flagged "transfer beam" (TB). CASCO's default assumption is that BEs transfer their load to the nearest column, which approximates the limit of infinitely flexible beams in bending. TBs cover the opposite case, of infinitely rigid beams in bending. Results for structures with TBs will not be presented here, so CASCO's handing of TBs is described in the Supplemental Material only.
In addition to initial loads and energy inputs, CASCO can add accidental events by removing or disconnecting some BEs or CEs or by suddenly applying additional loads or kinetic energy. In this way one can mimic accidental events such as impacts, explosions, or even an earthquake.
Initial intact structure After reading the input file, CASCO considers the structure before the accidental damage to simulate any collapse that might already take place. First, CASCO removes any CE that cannot contribute to load transfer; in Fig. 4, these are the CEs on the left side of BE32, the CE100 element, and CE112. Also the CEs including and above CE13 are removed, because the broken section in the column interrupts any load path there. This clean-up leads to the structure in Fig. 5. The mass of the removed CEs is considered to be negligibly small compared to the vertical loads, hence their weight and free-fall energy are not transmitted to elements below. Instead, any unsupported BE would not be removed, and the algorithm in Fig. 2 accounts for its free fall at step 8.
Classifying beams and columns CEs and BEs are classified as different types of columns and beams to later define a hierarchy of elements guiding loads from each BE and CE to the ground. Beams are classified first. They are defined as groups of adjacent and connected BEs. Missing BEs or broken cross sections terminate a beam. Fig. 5 shows the classified beams in the initial structure.
Load bearing columns, LBCs, and transfer beams, TBs, are classified next. Starting from the left end (bay 1), CASCO classifies all the columns at floor 1 as LBCs, because being connected to the ground they certainly provide potential load paths. Then the beams at floor 1 are considered, classifying as TBs ant connected set of adjacent transfer BEs (assigned in the input file) supported by at least two LBCs (see Supplemental Material for more details).
CASCO then moves to floor 2, classifying as LBCs the columns supported by an LBC or a TB underneath. The TBs at floor 2 are then classified with the same criteria as for floor 1, and the procedure is repeated for all the Defining the hierarchy of beams Each beam and TB is assigned a hierarchy level, later used to trace load paths. In CASCO, the loads preferentially travel to more fundamental levels of the hierarchy, the most fundamental one being the ground (hierarchy level 0). The criteria to trace load paths are given in the next subsection; here the focus is on constructing the hierarchy. From now on, only structures without TBs will be considered, the reason being that we will only presents results for structures without TBs here. Hierarchy and load path definition in the presence of TBs are discussed in the Supplemental Material.
Beams connected to one LBC at least are classified as "supported beams", SBs. For example, Beams 3, 4, 8, 9 and 10 in Fig. 5 are SBs. The others beams have higher hierarchy level: Beams 5, 6, and 7, have level 1, B h1 , being connected to SBs via columns that are not LBCs (the same would apply if they were hanging from TBs via non-LBCs). Beam 2 has level 2, i.e. B h2 , because it hangs from Beam5 that is a B h1 .
Computing load paths This section discusses how load paths are traced for structures featuring only supported beams, SBs, and lower-hierarchy beams connected to them, B (see the Supplemental Material for load paths in the presence of transfer beams).
The first step to compute load paths is to scan all beams and determine whether each node (beam-column connection) on them is a "sender" or a "receiver". Within each beam, the loads travel to the closest sender node, which directs them to its corresponding receiver node, either on another beam with more fundamental hierarchy level, or on the ground.
All nodes on the ground are potential receivers. For the SBs, sender nodes are those connected vertically to the ground below. SBs do not pass loads to other SBs: this is why Fig. 6 there is no load path from SB2 and SB3 to SB1 despite there is a column between them, and instead all their loads go directly to the ground as a preferential path. For each other beam, B, all the nodes are considered, searching underneath and above them for receiver nodes either on an SB or on another B with more fundamental hierarchy level and connected to it on a vertical line; if several possible receiver nodes are found for a given node on a B, only those on the Bs or SBs with most fundamental hierarchy level are classified as receivers. As an example, in Fig. 6, B1 transfers load to SB1 on one side and to both SB2 and SB3 (because they have same hierarchy level) on the other. B2 and B4 pass their loads respectively to SB1 and SB2, but not to B3 which  Fig. 4 after removing CEs that cannot contribute to load paths and after classifying load bearing columns (LBCs, in black), transfer beams (TBs, in black), other beams (in grey), and other columns that may still provide load paths (also in grey). Notice that Beams 5 and 9, originally flagged as transfer beams in Fig. 4, are now considered as normal, nontransfer, beams. The reason is discussed in the text of the manuscript Fig. 6 a Example of load paths in a framed structure featuring only supported beams and lower hierarchy beams connected to them. b Underlying topology (a) (b) has a less fundamental hierarchy level. B3 instead passes load to both B4 and B2 because they have identical and more fundamental hierarchy level. The rules presented here to identify sender and receiver nodes assume that structural elements with more fundamental hierarchy level i are infinitely stiffer than elements with less fundamental level i + 1 . So, in terms of stiffness, the hierarchy is: ground >> TB h1 >> TB h2 >> TB hn >> SB >> B h1 >> B hn . This is of course an approximation, as in reality all elements have finite stiffness and all contribute to the redistribution of load. This approximation implies that the code cannot compete in precision with Finite or Discrete Element simulations, when trying to capture complicated collapse modes. However, the approximation is conceptually useful because it leads to unique and non-overlapping load paths. The load path can be represented as in Fig. 6a and their evolution during collapse (viz. the creation of alternative load paths) can be traced too. The uniqueness of load paths also enables a representation of the underlying topology, e.g. in Fig. 6b, which will be used later to discuss the evolution of load paths during collapse.

Sampling collapse mechanisms
Having defined the load paths, one must decide the possible mechanisms of local collapse to be sampled. The number of such mechanisms is vast, and they can be described at various levels of detail. Since the principal scope of this manuscript is to present the CASCO algorithm, complexity in the description of local mechanisms is kept to a minimum. Five possible mechanisms of local failure are sufficient to capture the main patterns of collapse in 2D frames: (i) compressive or tensile failure of columns; (ii) formation of one plastic hinge causing collapse of cantilever beam elements; (iii) beam bending failure via two plastic hinges; (iv) three point bending failure of beams via three plastic hinges; (v) four point bending failure of beams via four hinges. Some simplifying assumptions will be made to describe these mechanisms; the quality of the resulting simulations will be assessed in the Results section.
For each mechanism, CASCO computes energies to use in an energy-based failure criterion for identifying mechanisms that are thermodynamically possible. Enegy criteria are often used in collapse-resistant design, as they allow simplified dynamic analyses using internal forces from static analyses [9]. According to the criterion, a mechanism occurs if: E k is the kinetic energy in all the beam and column elements (BEs and CEs) in the mechanism. U f is the failure energy of the mechanism, e.g. for several plastic hinges to form and rotate until a failure angle p . W q is the work of loads and forces involved in the mechanisms, from the initial undeformed state (having assumed rigid-plasticity) to collapse (e.g. the p rotations of certain plastic hinges). E d is the energy U f dissipated by previous failure events, redistributed between the BEs and CEs involved in the previous collapse mechanisms and propagated through the load paths (details in a later section).
Column failure in compression or tension The columns are assumed to fail in ideally plastic compression or tension, only due to axial force. Other mechanisms involving bending or buckling are not considered, but can be added in future versions of the simulator. Buckling, in particular, is the main collapse mode for steel structures; adding it to CASCO would require calculating the energy to activate and realise the local buckling mechanisms, for which there are methods available in the literature [4,31]. Figure 7a shows how load paths are used to calculate axial forces and kinetic energy through each column. Col3 and Col4 share equal parts of the distributed load between them, developing tensile forces of N = 10 kN that are passed as concentrated forces to points C and D. Col2 gathers all the load and forces between points B and D, and half of the load to the left of B, becoming compressed by N = −30 kN (besides shear and bending, not shown because of the assumption that they do not contribute to column failure here). Col1 is instead compressed by N = −10 kN . Figure 7a also depicts an input of kinetic energy per unit length, e.g. due to an impact from above. This is also transmitted via the load paths, thus Col2 takes 2 kJ of kinetic energy favouring compression, and Col1 takes 1 kJ.
The energy criterion (Eq. 1) for axial forces and energy in each column is: E ax k,i is the kinetic energy passing axially through column i: positive for tension and negative for compression. N p,i and u,i are the plastic axial force and ultimate axial displacement for column i; different values for tension t and compression c are assigned by the input file. N i is the axial force in column i, from the load paths. N p , u and N are positive when tensile, negative when compressive. E d,i is the previously dissipated energy transmitted through column i: its calculation will be discussed later, but it always opposes the kinetic energy, hence it is positive both in tension and in compression.
One-hinge plastic failure of beams Beams are assumed to be rigid-plastic in bending, failing due to formation of plastic hinges. For cantilever elements, one hinge is sufficient to cause failure. For example, for the A-B segment in , there is an externally applied E k that goes into the energy criterion in Eq. 1. The failure energy U f has two contributions. First, the energy dissipated by the plastic hinge, M − p ⋅ − p , with negative plastic moment M p and rotation at failure p as per convention in CASCO, because they set the top of the cross section in tension. The second contribution is due to catenary action, for which CASCO assumes an axial force N in the beam equal to to the smallest axial force at tensile yielding among all the BEs involved the mechanisms. N generates work by the elongation of these BEs, which is linked to their rotation at failure. Thus, for a cantilever of length L, the catenary energy is N 2 p L . Additional work W q comes from the load q times the vertical displacements at failure of all the BEs in the mechanism. For Fig. 7b, this leads to W q = qL CASCO samples cantilever mechanisms also assuming upward movements, thus positive moment M + p and rotation + p in the plastic hinge. Fro example, if such upward motion is sample for the mechanism in Fig. 7b, E k and W q would be negative and the energy criterion in Eq. (1) would not be satisfied.
For the position of the plastic hinge, still referring to the example in Fig. 7b, CASCO samples all possible positions on all right-end and left-end cross sections of each BE between points A and B. All mechanisms from hinge positions that satisfy the energy criterion are recorded as possible candidates for the next failure. One-hinge mechanisms instead cannot occur between points B and C or between D and E, because CASCO assumes: (1) that sender nodes (B and C) cannot move, because linked to elements with more fundamental hierarchical levels, which are considered as rigid (in Fig. 7b, the ground); (2) that nodes connected to columns cannot rotate (a cantilever mechanism between C and D would require D and E to rotate, which is not allowed). Other simplified simulation [11,32] and fully dynamic ones would allow the rigid rotation of the structural cluster to the right of C, should a plastic hinge form between C and D. CASCO does not consider cluster rotations for two reasons. First, to keep this manuscript brief and simple, as cluster rotation analysis would require rather complex algorithms, whose importance is not central here. Second, cluster rotation involves significant rotational inertia, making the mechanism slower than vertical displacements from shear propagation. Sudden dynamic responses to accidental events are more likely to generate two hinges near C and D and cause vertical motion of all elements to the right of C, rather than cluster rotation.
Two-hinge plastic failure of beams Figure 7c shows some examples of two-hinge mechanisms sampled by CASCO, where one plastic hinge only rotates and the other one rotates and moves vertically along with the rest of the structure attached to it. The plastic rotation at failure is the smallest between those of the two plastic Fig. 7 a Load paths for axial forces and energy. Beam bending collapse mechanisms involving b one, c two, d three, and e four plastic hinges hinges (one positive and one negative). If receiver nodes translate due to a two-hinge mechanism, for example the mechanism with hinges between C-D in Fig. 7c, the forces received by the nodes generate work W q by the vertical displacements of the nodes themselves.
CASCO samples all possible hinge positions in each beam. For example, between A and B in Fig. 7c, pairs of hinges are placed in all left and right cross sections of all BEs. For each pair of hinge positions, CASCO evaluates the energy criterion in Eq. (1) with analogous calculations as for the one-hinge mechanism (including unlikely mechanisms causing upward translations).
Two-hinge mechanisms where the segment between the hinges includes a node connected to a column are not sampled, due to the aforementioned assumption that such nodes cannot rotate. CASCO also ignores twohinge mechanisms causing vertical movement of sender nodes, whose motion would imply deformation of hierarchically more fundamental hence rigid elements.
Three-hinge plastic failure of beams Three-hinge mechanisms are only sampled between two neighbouring columns, e.g. in B-C, C-D, D-E and F-G in Fig. 7c, not between A-B because there is no column in A and onehinge mechanisms would prevail there. Figure 7d shows an example of three-hinge mechanism. All the involved BEs translate and rotate, hence nodes connected to a column cannot be involved. CASCO samples all the possible triplets of plastic hinge positions on the left and right cross sections of all the BEs between B and C. The same sampling is repeated for all BEs and sections in A-B, C-D, and in E-F. Energy contributions are computed with the same procedure as for the one hinge mechanism. The maximum plastic rotation is the minimum one from the hinges involved ( p positive or negative depending on the directions of rotation). The three plastic rotations are mutually dependent, so once the maximum rotation is identified, the other two follow by kinematics. The total elongation of all BEs involved in the mechanisms, used for the energy of catenary action, is L 1−2 ( + 1 ) , where L 1−2 is the distance between the first plastic hinge (in B in Fig. 7d) and the middle hinge in the mechanism.
Four-hinge plastic failure of beams Four-hinge mechanisms are only sampled between two consecutive sender columns, e.g. between nodes A-D and E-F in 7e. Groups of four hinges are placed in all possible combinations of left and right sections of all BEs between the two sender columns. The mechanisms causes the two lateral parts of the involved segment to rotate, and the middle segment to translate. CASCO makes two simplifying assumptions: (1) the middle segment does not rotate, but this assumption might be removed in future implementations; (2) the nodes connected to columns, e.g. B and C in Fig. 7e, can only belong to the middle segment of a mechanism, as they cannot rotate. Figure 7e shows two examples of four-hinge mechanisms. One of those, with first and last hinges near B and C, shows that the length of the middle segment can tend to zero, meaning that the two central hinges form on adjacent cross sections on two neighbouring BEs. This is sufficient to enable the displacement of the node C without rotation. A three-hinge mechanism in the same position would not be possible, as the plastic hinge would be on one or the other side of node C, and thus the node connected to the column should rotate. The elongation for the catenary action is computed as for the the three hinge mechanism, with and 1 now being the plastic rotations of the lateral hinges.

Realising a collapse mechanism
Out of all the sampled mechanisms, CASCO realises to the one that satisfies Eq. (1) by the largest margin. If this is a column failure in tension or compression, the CEs of the column are removed from the structure. A failing column is necessarily linking some sender nodes to some receiver nodes on a vertical line. Each sender node pertains to a beam, whose motion will be favoured by the loss of the column but which will also benefit from the energy dissipated by the column failure, U f , which is thus passed to the beams of those sender nodes as an E d term. Actually, all the beams and columns whose load paths rely on the failed column should receive a fraction of the column's U f as E d . However, identifying all these dependencies would be computationally intensive, and deciding the fraction of U f to assign as E d to each BE and CE would be arbitrary. Therefore, for simplicity, CASCO only assigns E d to beams that are immediately supported by the failed column. All the BEs of these beams take an equal share of E d . The mass of the columns is considered as negligible compared to external loads and weight of the beams, therefore the additional load and kinetic energy from the failed CEs is neglected.
If the mechanism to realise is a beam bending failure, only the cross sections reaching their ultimate plastic rotation p are broken. If a mechanism involves multiple hinges, only one will certainly reach its p : the others may or may not depending on kinematics and on their values of p . The failure energy U f is the sum of all M p p only from the broken sections. An equal share of it is assigned as E d to the BEs involved in the mechanism.
Fall of debris After realising a mechanism, CASCO reclassifies beams and columns, updates the hierarchy of beams and the load paths, and samples again all the possible collapse mechanisms to realise the next one. This iterative process causes the progressive removal of CEs and rupture of BEs, changing the topology of the structure. At some point, entire clusters may end up disconnected from the rest of the structure and from the ground. The elements pertaining to disconnected clusters are not tested for collapse mechanisms any more, because they should undergo free fall. CASCO carries out free falls only when no collapse mechanism satisfies Eq. (1) any more. The rationale is that wave propagation causing failures is much faster than free fall, hence the latter is realised later.
The CEs in a free-falling cluster are removed and disregarded, due to their negligible mass. Fallen BEs instead transfer their loads, kinetic energy E k , and dissipated energy E d , to the first BE underneath them on a vertical line, if that pertains to a beam that is not part of a collapsed cluster itself. The transferred kinetic energy is also increased by the potential energy of gravity that is lost during the fall (load on the BE times height of free fall before impacting a beam below). After load and energy are transmitted by free fall, CASCO restarts the sampling of possible collapse mechanisms.
In terms of topology, redistributions by free fall do not follow the previously defined load paths. Figure 8 shows the topology for the structure in Fig. 6 highlighting also the paths that of beams free fall. Such presence of two categories of links, governed by different rules of load and energy transmission, makes the topology of the system unique to the process of structural collapse.

Accidental event, end of simulation, and output
CASCO first samples the collapse of an intact structure, without any accidental damage or load/energy input. The simulation of the intact structure is over when no mechanism satisfies Eq. (1) and no cluster undergoes free fall any more. Then CASCO applies the accidental event, assigned by the user, and starts a new simulation of progressive collapse for the damaged structure. When the conditions of no failure nor free fall are satisfied again, the simulation is over.
In its current version, CASCO produces three output files, all in text format. One, with default name dump.dat, contains the spatial configuration of BEs and CEs, plus information about their geometry and type (in the sense of classification and hierarchy). The format of this file is the same as for ellipsoidal elements in the software LAMMPS [27]. CASCO saves one configuration after each collapse mechanisms, free fall, or accidental event. The configurations can be visualised with many standard programs: for example, the images in this manuscript have been obtained reading the dump file with the open-source visualiser OVITO [29].
A second output file, hierarchy.dat, prints one row at each step of collapse, counting the number of beams in each hierarchy level and separating the hierarchy groups of transfer beams, supported beams, and other beams.
The last output file, thermo.dat, contains the history of certain relevant quantities at each step of collapse. These are: the type of collapse mechanism triggered at the generic step; the id of the beam affected by the mechanism; the leftmost and rightmost BEs involved in the collapse mechanism; the number of sections broken by the mechanism; the energy dissipated by the mechanism; the total potential energy of the structure (i.e. all the loads on all BEs times their height from the ground); the total energy dissipated by all previous failure mechanisms; and the total kinetic energy in the system (obtained as the sum of E k on all BEs). These quantities enable a quantitative description of the collapse process, as shown in the following section.

CASCO has been validated against various examples with analytical solution, shown in the Supplemental Material.
Here, instead, CASCO is validated against previous, fully non-linear, dynamic simulations of collapse of 2D frames made of reinforced concrete [23]. The simulations from the literature were themselves validated against experiments of dynamic collapse of reinforced concrete beams, mimicking sudden column removal of accidental origin [22]. Direct experimental validation of CASCO would be desirable, but the literature lacks systematic experimental campaigns of progressive collapse of reinforced concrete frames, including energy transfer by falling debris, which is where the focus and originality of CASCO lie. CASCO is then used to explore how collapse resistance and propagation are affected by mechanical heterogeneity, viz. structural elements with random strength within certain ranges. The results are discussed in usual terms of collapse loads and mechanisms, but also looking at evolutions of load paths and topology, thanks to the new insight from CASCO.

Progressive collapse of homogeneous frames after columns removal
Three 2D frames are considered, with same overall size but different number of beams and columns: see Fig. 9. An  Fig. 6: load (solid) and free fall (dashed) paths initial damage destroys the shaded area in Fig. 9, which is the same for all frames. The collapse loads and mechanisms are already known from fully non-linear, dynamic, Discrete Element simulations [23], which considered two scenarios: one where the columns are strong, inducing collapse by beam bending, and the other where columns are weak, inducing collapse by column crushing. These global collapse modes are respectively called "bending" and "pancake". The geometric and mechanical parameters used in CASCO are the same as in [23]. The beams have rectangular cross section, with height h = L∕10 and depth d = 2 3 h . The columns have square cross section with edges b = d = H∕10 . The tensile axial force at yielding N + p is assigned based on assumptions on reinforcement content [23], leading to N + p = s f y A , where A is the area of the cross section, s is the fraction of area occupied by steel reinforcement, and f y is the yield strength of steel. Here f y = 440 MPa and s = 0.0029 and 0.0226 respectively for beams and columns [23]. The axial force at compressive yielding is N − MPa for simulations wit strong columns favouring bending collapse, and f c = −0.35 for simulations favouring pancake collapse. The axial strain at tensile and compressive failure are 0.1 and −0.007 , same for beams and columns. The plastic moments, same for positive and negative bending, are M + p = −M − p = s Af y h , with = 1∕2 for the beams and 3/8 for the columns, having to do with the arrangement of reinforcement bars in the cross sections as per [23] (h to be replaced by b for the columns). The ultimate plastic rotation is 0.2, same for beams and columns, and for positive and negative bending. The load is uniform on all beams: its value is not assigned but is sampled to find the collapse loads. Figure 10 shows the loads causing dynamic collapse after damagem comparing results from CASCO with those from nonlinear dynamic simulations in [23]. The results are similar for pancake collapse, as shown in Fig. 10b. For the bending collapse mode instead, in Fig. 10 am CASCO underestimates the collapse loads when catenary actions are neglected (as expected, given the importance of catenary action for progressive collapse resistance, proven both experimentally [12] and by modelling [2,7,16]). On the other hand, when considering catenary actions, CASCO predicts higher collapse loads than those from the dynamic simulations. The reason is that the dynamic simulations assumed that the ultimate plastic rotation p decreases when a beam is under tension, whereas p in CASCO does not. Overall, Fig. 10 shows that despite its approximations, CASCO predicts collapse modes and loads that agree with more detailed simulations.
The bending collapse of the frame with 5 bays is detailed in Fig. 11, to show the additional insight provided by CASCO. The snapshots show an expected progression of four-point bending mechanisms. Initially, the mechanisms on floors 2 and 3 develop asymmetrically, with the cross sections on the right failing at mid-bay rather than next to a beam-column node. This asymmetry would not emerge in dynamic simulations. CASCO predicts asymmetry because it all the allowed mechanisms without  (1) by a larger margin than their symmetric counterparts, hence be chosen by CASCO as most probable. For the case in Fig. 11, the energy contribution favouring the asymmetric mechanism is due to catenary action: simulations without catenary actions indeed led to symmetric mechanisms only. CASCO could be forced to respect symmetry, but real structures are never perfectly symmetric, so sampling also asymmetric mechanisms is more cautious.
CASCO enables a quantitative analysis of the evolutions of structural topology and load paths during collapse. For Fig. 11 Topology evolution during the collapse of the 5-bay frame. Red dots on the frames are failed cross sections. SB are supported beams connected to the ground. B are beams supported by other beams or unsupported. The first number following SB or B is the floor level; the second one is the position on that level from left to right. In the graphs, free fall paths are dashed and multiple solid load paths indicate that multiple columns provide them the frame with 5 bays, these are shown in Fig. 11 and are quantified in Fig. 12. At collapse step 0 the structure is intact. At step 1 the initial damage is applied. At step 6 the structure is only one failure event away from a large cluster-collapse event. At step 7 the central part of the structure is disconnected and ready to fall. Fig. 12a shows the that the number of beams, viz. sets of connected beam elements BEs, increases while damage propagates. Incipient collapse features a sharp increase of collapsed beams and a corresponding decrease of not-collapsed beams (collapsed beams have no load path left to transfer their loads to the ground). Figure 12b shows the evolution of links and load paths numbers during collapse. The total number of links and the number of inactive free fall links (terminology explained in the figure's caption) increase until step 6 and then decrease at incipient collapse between steps 6 and 7. These trends reflect the analogous evolution in the number of not-collapsed beams during collapse propagation, in Fig. 12a. Also the number of "non free fall" (NFF) links initially increases and then decreases during collapse propagation, but their decrease starts at step 5, thus before incipient collapse. This indicates that progressive collapse resistance might rely on a structure's ability to generate new NFF links while failures propagate. The increase of NFF links means that columns that were not part of load paths just after the initial damage get involved later during collapse progression. In Fig. 11 these are the columns just above the initially damaged area, which start to act as vertical ties under tension only when collapse progression leaves the beams in the centre of the structure hanging from the beams above. These new active links are quantified by the "Beam-SB NFF" curve in Fig. 12b, which also peaks at step 5. Therefore, the increase of NFF links is related to the transformation of some directly supported beams SB into beams with lower hierarchical level, hanging from other SBs. Enabling a structure to create new lower hierarchy levels during collapse, might be a conceptual strategy to improve collapse resistance.
Another quantity that peaks at step 5, before incipient collapse,is the number of load paths in Fig. 12b. The paths are computed combinatorially: for example, at step 4 in Fig. 11, beam B2,3 is linked to SB4 by two links, and SB4 is linked to the ground by 4 links, giving 2x4=8 load paths. Out of the three types of links peaking at step 5, only the load paths retain the peak when expressed as a density, i.e. divided by the number of not-collapsed beams. This is shown in Fig. 12c, where the density of NFF links does not peak at all and the density of Beam-SB links peaks immediately after the initial damage. Therefore, a structure's ability to increase its number and density of load paths during collapse might be a target for optimisation towards collapse resistance.

Progressive collapse of heterogeneous frames after columns removal
This section considers the same 5-bay frame as in the previous section, but adding heterogeneity in the bending and axial strengths of individual BEs and CEs. Heterogeneity is introduced assigning plastic moments M p and axial forces at failure N p , both positive and negative, with same average value as in the previous section but adding or subtracting M p or N p chosen randomly and uniformly within a range, ± M p and ± N p . With values between 0 and 1, is an uncertainty parameter measuring the mechanical heterogeneity. Figure 13a shows that uncertainty has a negative impact on strength. = 0 corresponds to the collapse load q = 19.5 kN∕m of the homogeneous frame with catenary action in Fig. 10. At low , different realisations Fig. 12 Evolution of topological quantities during collapse. a Number of beams: collapsed beams are linked to the ground only by free fall links, SBs by one or more columns. b Number of links, excluding active free fall links that pertain to collapsed beams. "Beams" indicate non-collapsed beams that are not SBs, thus sup-ported by other beams. NFF are "Non free fall" links provided by columns. "Inactive free fall links" involve beams that are not collapsed. "Total" links are the NFF plus the inactive free fall links, excluding active free fall links. c Link densities: same as in b divided by number of non-collapsed beams of randomly distributed strength lead to similar collapse loads. Instead, at high , the collapse load can vary significantly, e.g. between 5 and 7.5 kN/m for = 0.9. Figure 13b, e show that moderate uncertainty, = 0.5 , has little effect on collapse extent and dissipated energy, with all repetitions experiencing no collapse propagation for loads below 12 and 14 kN/m, and same propagation as the homogeneous frame for greater loads. Also at high disorder, = 0.9 in Fig. 13c, most frames go directly from no collapse propagation to similar results as for the homogeneous frame. However there are exceptions. One is that some frames display a final collapse extent that is greater than that of the homogeneous frame. This is due to some columns on the left and right of the initial damage failing during collapse propagation. This indicates that strong columns are key to limiting horizontal propagation of collapse. Another exception is the frame corresponding to the the bold curves in Fig. 13. This frame features some proportionality between load intensity and collapse extent, which is a desirable feature in collapse-resistant design [9].
The dissipated energy in Fig. 13f shows similar trends as the collapse extent, which simply indicates that wider collapses cause more energy dissipation. One may wonder whether the dissipated energy per unit of collapse extent affects the extent of collapse itself. Figure 13d shows that this is not the case, because the same value of specific dissipated energy can correspond to broad or limited collapse extents. On the other hand Fig. 13d also shows that the possible outcomes are clustered around two regions. One region is centred around point (1,1), representing a similar mechanism as for the homogoenous frames, with some variability due to the disorder . The other region entails similar specific dissipated energy, but less collapse: the points in these region are the same points in bold in Fig. 13c when the load is q < 7 kN/m. Figure 13 raises the question of what makes the frame in bold develop proportionality between load intensity and collapse extent. The first thing to consider is the qualitative progression of damage, shown in Figs. 14 and 15, respectively for loads just below and just above the threshold load q ≈ 7 kNm in Fig. 13c. When q = 6 kN/m , collapse is different from that of the homogeneous frame, in that its final extent is very limited, it starts on the third floor and not the second floor, and it propagates downward and not upward, driven by redistribution of internal forces and falling debris. However, despite these differences, Fig. 13 Impact of uncertainty on collapse load, extent, and dissipated energy of the 5 bay frame in Fig. 9. Each curve is a different realisation of heterogeneous frames with same average parameters but different spatial distributions of strength. Bold curves highlight one realisation whose behaviour is particularly interesting. (a) Collapse load reduced by mechanical disorder. (b,c) Collapse extent after damage, viz. metres of collapsed beams, as functions of the applied load, for frames with = 0.5 and 0.9. The collapse extent is divided by 60 metres of beams collapsing in the homogeneous frame. (e,f ) Energy dissipated by all failed cross sections after damage, divided by 1.043 MJ dissipated by the homogeneous frame when collapsing. (d) Relationship between energy dissipated per unit length of collapsed beam, and collapse extent: two regions emerge, indicating limited and broad collapse propagation the progression of collapse via four-hinge mechanisms is qualitatively similar to the first steps of the homogeneous frame in Fig. 11, and also the evolution of topology is very similar to the homogeneous case (hence not shown here). Interestingly, the final configuration of the partially collapsed heterogeneous frame in Fig. 14 corresponds to the configuration of the homogeneous frame at its collapse step 5 in Fig. 11. Step 5, for the homogeneous frame, was just before incipient collapse, and was when the number of NFF links and load paths started to decrease in Fig. 12. This again suggests that the ability to multiply links and load paths might underlie progressive collapse resistance, whereas a decreasing number of those might a driver for collapse. Figure 15 shows the collapse progression for the same frame as in Fig. 14, but under a larger load q = 7 kN/m. The first steps are analogous to those in Fig. 14, but a difference emerges in step 5, when the tensile failure of a column at the top floor triggers the four-hinge bending failure of the third floor beam. This failure is followed by tensile failure of another column at the top floor (step 7), resulting in a wider final collapse, although still smaller than the homogeneous frame in Fig. 11. From the topological analysis, step 4 in 16.a coincides with step 5 of the homogeneous frame in Fig. 11, hence all the quantities in 16b, c and d at Collapse sequence for the 5 bay frame with = 0.9 and load q = 7 kN/m step 4 also coincide with the same quantities at step 5 in Fig. 12. After step 4 in the heterogeneous frame, the number of NFF links and load paths in Fig. 16c stop increasing and start to decrease, as well as their densities in Fig. 16d. This is similar to what was observed at step 5 in the collapse of the homogeneous frame, suggesting again that the ability of the structure to multiply links and load paths may be important for progressive collapse resistance. This however does not explain why the heterogeneous frame undergoes less collapse than the homogeneous one.
A feature that might correlate to the final collapse extent, instead, is the rate at which the number of beams and links change after the critical step 4 of collapse (step 5 for the homogeneous frame). Figure 16b shows that the rate of increase in number of beams is smaller for the heterogeneous frame than for the homogeneous one (cf. Fig. 12a). Furthermore, the number of load paths in Fig. 16.c decreases linearly after peaking at step 4, whereas in the homogeneous frame the loss of NFF links accelerates after step 5 (see Fig. 12b). Similarly, also the loss of total links and NFF is more gradual in the heterogeneous frames, which experiences less collapse. Finally, the density of load paths in Fig. 16d shows a gradual decrease after step 4, suggesting a horizontal asymptote towards step 7. Similarly, the density of NFF links decreases more gradually compared to its counterpart for the homogeneous frame in Fig. 12a, and the density in total number of links even increases between steps 6 and 7. Overall, these results indicate that the rate of topological changes during collapse progression might contribute to determining the final extent of collapse, thus providing a possible target for structural optimisation.

Conclusion
The new algorithm presented here, CASCO, analyses the progressive collapse of 2D frames made of reinforced concrete, by quantifying the evolutions of load paths and other topological features. Despite introducing simplifying assumptions and considering dynamics only via energy conservation, CASCO can predict collapse loads and mechanisms that are comparable to those obtained from fully nonlinear dynamic simulations. Applications to the collapse of mechanically homogeneous and heterogeneous frames indicate that structures whose number of load paths increases during damage propagation, are more likely to avoid collapse initiation after initial damage. If collapse starts, structures that lose fewer load paths after each collapse propagation step are likely to experience a more limited final extent of collapse. Creating structures with such abilities might thus provide targets for optimisation towards progressive collapse resistance. The arising challenge is now to understand which geometries and material properties may lead to desirable evolutions of the number of load paths. CASCO can be exploited in this sense, by testing various strcutural design solutions and determining how these affect the evolution of load paths and the progression of collapse (for now only in 2D frames, but with scope to extend the algorithm to 3D structures too). Furthermore, CASCO's topological interpretation of structures provides a starting point to create more abstract models of structural collapse, for example network models. These would differ from existing network models, in that they should capture some peculiar topological features that emerged from the simulations in this manuscript, such as the coexistence of different types of beams and links, and the possibility for their number to decrease as well as increase during failure propagation. Such conceptual models would help making structural collapse amenable to analysis using the tools of statistical mechanics [19,28], for a more fundamental understanding of the phenomenon. In this sense, CASCO is a first effort to develop a minimal model of structural collapse that can help understand the main drivers of progressive collapse, akin to how the Ising [25] and Fibre Bundle models [13] have scaffolded the rationalisation of phase transitions and fracture in materials.

Compliance with ethical standards
Conflict of interest The author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.