Reinforced Concrete Beam under Support Removal—Parametric Analysis

This paper investigates the behaviour of a reinforced concrete beam under a support removal. A detailed parametric analysis is carried out, covering the effect of support removal rate on dynamic response. The linear elastic and nonlinear inelastic responses are computed and studied in detail. Critical parameters during the structural response are identified. In order to determine the ultimate load, the vertical pushover analysis is performed. The key parameters driving the beam response are assumed as random variables, and respective reliability study makes it possible to check the overall uncertainty of the dynamic response. In particular, the response spectrum measuring the effect of support removal rate has been computed. It has been demonstrated that the critical vertical response occurs when the time of support removal is up to to 17% of the first natural period. The vertical pushover analysis results in obtaining capacity curves and showed the order in which two plastic hinges occur for various load patterns. Finally, the reliability-based sensitivity analysis indicates the geometric cross-section cover and height are the most sensitive parameters of the beam response.


Introduction
The recent (24 June 2021) partial collapse of the Champlain Towers South condominium located at the Miami suburbs of Surfside gave rise to increasing attention of better understanding of progressive collapse of buildings. The research of this problem has a long tradition starting from the Ronan Point apartment building partial collapse in London 1968 [1]. The problem came back to the spotlight of engineers after the collapse of the Alfred P. Murrah Federal Building (Oklahoma City, 1995) (e.g., Osteraas [2] and Kazemi-Moghaddam and Sasani [3]) or World Trade Center collapse in 2001 (e.g., Bazant and Verdure [4]).
The progressive collapse phenomenon happens when a local failure causes successive structural damages leading to the partial or total collapse of a building. The problem of progressive collapse of buildings often leads to human life losses enhancing fear in the society and generating substantial economic costs. Therefore, further progress in understanding this phenomenon at various stages (in designing, erecting and maintaining buildings) is of great importance to mitigate its negative consequences. An extensive literature review on research and practice in the field of progressive collapse can be found in Mohamed [5], Nair [6], Kokot [7], Kokot and Solomos [8] and Adam et al. [9].
Since the reinforced concrete composite is one of the most popular materials used in structures, this paper is devoted to the behaviour of a reinforced concrete beam structure under a support removal. Nonlinear modelling of reinforced concrete structural members has already been analysed by the author of this paper using concentrated plastic hinges [10,11]. However, this study focuses on the distributed plasticity using the so-called fibre modelling of reinforced concrete cross-sections. Such more detailed, parametric analyses of structural response under support removal can be beneficial before planning complicated 3D finite element analyses of buildings in progressive collapse.
The nonlinear inelastic modelling of reinforced concrete is based on the guidelines recommended in Bao et al. [51] and Alashker et al. [52] and is validated by reproducing the results of the vertical pushover curves in Bao et al. [51] where the authors reported the behaviour of reinforced concrete subassemblies modelled with micro-and macro-elements. The micro-elements are referred to continuum type finite elements (the model composed of such elements was considered as a high fidelity model), whereas the macro-elements are the beam or frame finite elements using various formulations (displacement-based or force-based elements; see e.g., Neuenhofer and Filippou [53], Calabrese [54]). According to Alashker et al. [52], in order to perform progressive collapse analyses effectively (in reasonable time and computer memory usage, preferably on personal computers), the modelling simplifications are inevitable; however, the degree of simplification must be controlled. The authors of this paper proved that a properly calibrated model using macroelements can be regarded as reliable as far as accuracy of the material and structural behaviours are concerned.
Although there are many research studies on more complex structures, this paper focuses on a reinforced concrete two-span beam under a support removal. Such beams can represent, for example, footbridges or regular bridges.
The three objectives of this paper are: (1) the examination of the impact of support removal rate on the dynamic response, (2) the observation of the nonlinear behaviour of reinforced concrete beams under increasing deformation up to collapse (through the vertical pushover analysis using various load patterns), and (3) reliability-based sensitivity studies leading to the identification of the most important geometric and material parameters under vertical pushover loading conditions.
In the author's opinion the original/novel contributions and findings of this paper are: 1. Establishing a characteristic shape of a response spectrum curve.

2.
Comparison of various load patterns in the vertical pushover analyses and finding the percentages of the reduced ultimate capacities for concentrated load patterns compared to the uniformly distributed loading.

3.
Identification of important input parameters (according to the so called importance measures) under the vertical downward loading as it is the case of progressive collapse.
The paper begins with Section 2 which introduces the computational model of the analysed beam. Section 3 is devoted to the determination of the response spectrum of the analysed beam for wide range of support removal rates. Section 4 deals with the vertical pushover analysis enabling us the assessment of the ultimate load and understanding the failure mechanism. Section 5 addresses the sensitivity studies using the the so called important measures of the first-order reliability method. Finally, Section 6 summarises and concludes this paper.
It should be emphasised that all models used in this paper have been created and analyses have been performed in the OpenSees free open-source software framework (OpenSees is a free software framework that enables users to create finite element applications for simulating the behaviour of structures under static and dynamic loading. It is an open source project started by Frank McKenna and other researchers joined to its developmentdeveloped by other researchers [55][56][57]. OpenSees users can create application using Tcl [58,59] and Python [60] scripting programming languages.

Computational Model of the Analysed Structure
Consider a two-span beam which a computational model is shown in Figure 1a. The beam has three supports: the left fixed support (node A), the inner roller one (node B) and the right roller support (node C). The lengths of the spans are L 1 = 6 m and L 2 = 4 m. The beam is made of concrete with the initial modulus of elasticity E c = 10.6 GPa, and the cross-section is rectangular for the whole beam. Assuming that the inner support is to be destroyed due to an unspecified accidental loading, the beam becomes a single span beam of length L = 10 m, when the reaction R B vanishes to zero (see Figure 1b).
Computational model of (a) a two-span beam with the inner support removal scenario and (b) an equivalent single span beam with reaction R B to be removed.
In the linear elastic analysis, the OpenSees 'elasticBeamColumn' Euler-Bernoulli finite element is used, whereas for the nonlinear inelastic analyses, the 'forcedBeamColumn' element is chosen. This force-based element is based on the formulation proposed and developed in Neuenhofer and Filippou [61] and Neuenhofer and Filippou [53], where the plasticity is controlled in five Gauss points located along the beam, according to the Gauss-Lobatto scheme, and the reinforced concrete cross-sections are modelled using the so called fibre modelling. To each fibre of discretised cross-section, a constitutive material relationship is assigned for both steel and concrete. The steel constitutive model is based on the Menegotto-Pinto formulation [62][63][64][65] (OpenSees 'Steel01' uniaxial material model; see Figure 2a for the backbone of the steel stress-strain relationship). The concrete material model uses OpenSees 'Concrete01' implementation, which neglects the concrete tensile capacity while in compression, is characterised by the initial parabolic stress-strain curve up to the ultimate strength f c0 (e.g., Hognestad [66]), followed by linear descending branch, down to the residual stress f cu . The backbone of the concrete stress-strain relationship is shown in Figure 2b. For beams, the cross-section is usually considered unconfined with zero residual stress. The inclination of the descending line depends on the amount and distance of transverse reinforcement. The theory behind this concrete material follows the works of Kent and Park [21], Hognestad [66], Karsan and Jirsa [67] and Scott et al. [22]. In the subsequent sections, in the nonlinear inelastic analyses, the following material parameters have been used: f c0 = 10.6 MPa, ε c0 = 0.002, f cu = 0.2 f c0 , f y = 500 MPa and E s = 200 GPa, and the steel strain-hardening ratio b (ratio between post-yield tangent and elastic tangent) is assumed as 0.001. It should also be noted that in order to account for large displacements, the so called corotational formulation (e.g., De Souza [68], Crisfield [69]) is used in the nonlinear static and dynamic analyses.
The cross-section is designed for the uniformly distributed loading of 5 kN/m, resulting in the rectangular concrete section: 0.2 m wide and 0.4 m high. The longitudinal reinforcement represents 3φ16 bars in the bottom and 3φ16 in the top (symmetry about the horizontal neutral axis). The cross-section is modelled and analysed in OpenSees and its moment-curvature characteristics are shown in Figure 2c. The first steep increasing segment (meaning the increase both in concrete stress of the compressed zone and in axial force of the tensile rebars) ends with the onset of yielding of the tensile rebars at the bending moment of 125 kNm. Then, in the second segment, the concrete stress further increases up to the compressive strength f c0 , and at the same time, the compressed concrete zone moves up the cross-section until the strain of tensile rebars exceeds the fracture limit, which happens at the ultimate curvature of about 0.32 m −1 . Due to symmetry of the cross-section, the moment-curvature curve is antisymmetric, and only the right side of the plot is shown. The typical shape of the moment-curvature curve agrees with experimental tests for rectangular reinforced concrete beams (see, e.g., Park and Paulay [70], Espion and Halleux [71], Srikanth et al. [72]). In dynamic analyses, it is valuable to perform the modal analysis to get information on dynamic properties of the structure, and therefore, first, the eigenvalue problem is solved for the initial two-span beam. The first two natural periods are equal to T 1 = 0.054 s and T 2 = 0.032 s, and their associated mode shapes are presented in Figure 3a.  Since the inner support is to be removed, it is even more important to get the dynamic characteristics of a single span beam (without the inner support). In this case, the first natural period is equal to T 1 = 0.168 s (more than three times longer than for the initial model), and the second one is T 2 = 0.052 s (1.6 times longer than previously). Obviously, the natural periods are greater in the second case due to smaller stiffness of the resulting one span fix-roller support beam. The first two mode shapes are illustrated in Figure 3b.
It should be noted that although for beams, the mode shapes and vibrations are only in the vertical direction, for frames, often the dominating motion (first mode shape) is in the horizontal direction. Therefore, the modal analysis should mainly identify the modal characteristics in the vertical (gravitational) direction when analysing the progressive collapse problem.

The Impact of the Sudden Support Removal
Assuming that the computational model of the structure is adapted to the formulation of the finite element method, for the nonlinear inelastic model, the following matrix equation of motion needs to be iteratively solved where M is the mass matrix,Ü is a vector of accelerations, P r (U,U) is the vector of nodal resisting forces (which can be linearised using the Newton-Raphson method), and P(t) is the vector of external time-dependent loads. When the model is linear elastic, the above equations can be simplified to where C and K are the damping and stiffness matrices, whereasU and U are the vectors of velocities and displacements, respectively. Equations (1) and (2) can be numerically solved using direct time integration methods such as Newmark family or Hilber-Hughes-Taylor (HHT) methods (e.g., Newmark [73] and Hilber et al. [74]). In this study, the Newmark method is used as implemented in the OpenSees finite element method framework. In Sections 3.1 and 3.2, the dynamic behaviour of the two-span beam under sudden (shortest possible duration) and increasing durations of support removal is investigated.

The Sudden Support Removal
First, let us begin with the analysis of the most adverse dynamic response. To this end, assume that due to an unspecified cause, the inner support is suddenly destroyed. Therefore, the support removal time t r is equal to the time interval dt = 0.001 s used in the Newmark integration scheme when solving the system of the equations of motion. It is also assumed the so called Rayleigh damping is in the form of αM + βK where the parameters are proportional to the mass and stiffness matrices. Given 5 % damping ratio for the first two modes and using the formula it results in α = 2.3 and β = 0.000162.
In any finite element method software which offers a module for the solution of linear or nonlinear dynamic problems, the sudden support removal can be modelled in two stages. In the first stage, the support to be removed is included in the model, and after applying the gravity loads, the reaction values need to be recorded. In this case, the value of R B at support B must be noted. In the second stage, in the other modified model, the support (not present) is replaced by the reaction force to obtain a model equivalent to the initial one. Then, a decreasing in time force (of the initial value equal to the support reaction and the downward direction) is applied. In this way, we can control the time in which the reaction is cancelled and hence simulate instantaneous support removal.
The following procedure is implemented, and the analysis is performed in OpenSees software framework. To this end, a uniformly distributed load is statically applied within 1 s. Then, at time 1.001 s, the reaction R C is cancelled by counter-reaction −R C . The beam begins to deform downward, reaching the maximum deflection 0.072595 s at time 1.110 s, following by the free damped vibration stage as illustrated in Figure 4. The vertical displacements at two characteristic nodes are shown: midspan node at distance x = 5 m from the left fixed support and node B (x = 6 m) of the removed support (see Figure 1). Due to different types of supports (the left support is fixed, and the right one is a roller support) a slightly larger displacement is developed at node B.

The Response Spectrum to Various Support Removal Rates
To analyse the impact of the support removal time on the dynamic response (here the vertical displacement) at node B of the removed support, a set of analyses have been run for removal times t r changing from 0.001 s to 6 s. Figure 6 shows how the maximum vertical displacement u 6y at node B changes with the increasing values of t r . On the right-hand side in this plot, there is also information on the dynamic factor in reference to the static displacement (that is for the case when the support is removed within t r = 6 s). The values collected in Figure 6 are extracted from individual 75 time histories shown in Figure 7. The plot in Figure 6 can also be regarded as a response spectrum of a structure undergoing the support removal for a wide range of removal durations. To have some insight into the development of the first crucial part of vibrations, an enlarged version of Figure  Note that the maximum vertical displacement of u 6y is obtained for t r containing in the range of 0.001 − 0.01 s. For t r = 0.001 s, the maximum vertical displacement equals u 6y = 0.0726 m and for t r = 0.01 s, u 6y = 0.0724 m with only 0.2% relative difference. For t r = 0.052 s, u 6y = 0.069 m, the relative difference increases to 5%. It means that sudden support removal can be considered for up to t r = 0.052 s, which is 30 % of the first natural period (T 1 = 0.168 s). With further increase in t r , we can observe a fast decrease in the maximum vertical displacement down to the point of (t r = 0.168 s, u 6y = 0.040 m). The ordinate of this point coincides with the fundamental period of beam after the support removal (see Figure 3b). Following further the curve in Figure 6, there are six noticeable cycles where alternately the maximum vertical displacement increases and decreases.
At first glance, it seems counter-intuitive that for greater values of t r , the maximum vertical displacement is greater than for a smaller t r . However a closer look at the associated velocity values helps understand this phenomenon. The related time instants where the vertical displacement takes the maximum value are shown in Figure 9 with dots whereas asterisks depict the times when the cancellation of the support reaction ended. In this way, the velocities indicate that (a) the maximum vertical displacement occurs when the respective velocity is zero (dots in Figure 9), (b) when t r exceeds the period T 1 = 0.168 s, the second cycle of vibration begins in the downward direction and hence delays approaching the value of zero velocity. The length of these cycles is of about 0.17 s and coincides with the first natural period. Eventually with longer support removal times (t r > 2.5 s), the maximum vertical displacement stabilises and converges to the static displacement of 0.039 m.
It should also be noted that the maximum dynamic factor reaches the value of 1.8 (see Figure 6). Inner support removal -linear analysis A similar study can be carried out when the material properties for concrete and steel are nonlinear inelastic. Let us start with the individual time-histories of the vertical displacement at node B versus the increasing support removal times t r as shown in Figure 10. For very short removal times (t r < 0.06 s), the beam response is highly inelastic, and after the first waves of vibrations, a significant plastic deformation is noticeable. Moreover, for all individual time-histories, the duration time of the free vibration is shorter than in the case of linear analysis, which can be explained by the fact that, in the nonlinear case, the tangent stiffness changes in different phases of deformation, and this fact can be regarded as an additional damping factor.
The maximum values of the individual 75 time-histories shown in Figure 10 are collected and presented in the form of the inelastic response spectrum in Figure 11 where the maximum vertical displacement u 6y at node B is plotted versus the increasing values of t r . On the right of this plot, we can see the dynamic factor scale in reference to the static displacement. The maximum dynamic factor reaches the value of 1.91 for t r = 0.005 s and quickly decreases to the values of 1.1-1.2 for t r = 0.20 s; then, the maximum displacement of u max 6y slowly approaches the static displacement. It is interesting that the characteristic cycles visible in the linear elastic response spectrum (see Figure 6) are not present here, mainly, because the natural period changes in time as the nonlinear inelastic vibrations proceed.
It should be emphasised that in the nonlinear case, the structural response depends on the excess of capacity with which a structure is designed against a particular support or other load-bearing element removal scenario. Would there be no redundancy in the capacity, the structure would not survive the first downward movement stage leading to partial or total collapse.
It could be argued that the present study contains the systematic assessment of the influence of support removal rates on the dynamic response. Existing publications (e.g., Russell et al. [75,76] and Adam et al. [77]) usually consider selected column removal durations for specific assumed fast-dynamics failure events (e.g., blast), but simulations or experimental tests of slower-dynamics element failures such as impacts or material deterioration require longer load-bearing element removal duration. In this light, this analysis gives further insight into the impact of support/column removal rates on the structural dynamic behaviour which, in turn, can be helpful in planning expensive experimental tests and/or designing an appropriate support/column removal technique.

The Vertical Pushover Analysis
Now, let us perform a vertical pushover analysis of the two-span beam to evaluate its ultimate capacity against the inner support removal. It should be noted that the initial system is twice statically indeterminate which means that formation of the second plastic hinge transforms the system into a mechanism, and the structure can be regarded as collapsed. From the numerical point of view (and the finite element method) the onset of a mechanism can be indicated by the lack of numerical convergence.
According to the theory of pushover analysis in the seismic engineering, the reference loads and controlled displacement are considered in the horizontal direction. Here, in the progressive collapse problems, the pushover analysis can be adapted for the vertical direction when the support is removed from the model and the vertical deformation is induced until reaching the ultimate capacity of the system. The reference vertical loads are applied to the model in the displacement control strategy to reach a prescribed target displacement at a chosen node and vertical degree of freedom. In many experimental tests (e.g., Lew et al. [78,79]), one concentrated force is applied to the structure mainly to simulate the inertia forces generated by sudden removal of a load-bearing element. One concentrated force is selected due to simplicity reasons; however, in realistic progressive collapse events in smaller frames, a uniformly distributed gravity load can better represent the real external force distribution. Therefore, in the following considerations, four cases of applied loading distribution are used. The assumption is that in all cases, the sum of vertical loads is equal to the reference force P = 1 kN. This is the basic simplified case shown in Figure 12a. In the second and third cases, the reference loading is represented by two P/2 and four P/4 forces, respectively, as illustrated in Figure 12b,c. The last case is the realistic load pattern in the form of the uniformly distributed load q = P/L (see Figure 12d). The main outcome of the vertical pushover analysis is the capacity curve which represents the total vertical loading (hence the sum of vertical support reactions) as a function of the increasing vertical displacement u 6y at the point of the removed support. Figure 13 shows that when a single vertical force "1P" is applied, the maximum vertical load the beam can withstand is the smallest value of V = 77.6 kN among the four cases. The first drop at u 6y = 0.436 m is caused by the formation of the first plastic hinge at point B. This fact is visible in Figure 14 where at the same time instance, the moment capacity drops to zero. However, there still exists some residual vertical capacity of about (V = 22.4 kN) up to u 6y = 0.703 m, where the second plastic hinge at point A causes the moment M A drop to zero (see Figure 15) and eventually transform the beam into a mechanism. The order of plastic hinges, in which they develop, is associated with the bending moment distribution presented in Figure 16 (which shows the bending moment diagram for a unit force). It can be seen that the maximum bending moment is located at point B and is equal to 1.728, compared to the smaller value of 1.68 at point A. The same diagram indicates that for other load pattern cases: "2P", "4P", and "q", the maximum bending moment is located at point A, while the bending moment at point B is smaller. The relevant values for M A are: 1.8, 1.5, 1.25 and for M B : 1.28, 0.8, 0.7, respectively. This observation will imply the order of plastic hinges in which they occur in the "2P", "4P" and "q" cases of the vertical pushover analysis. Moreover, Figure 13 shows that when the load is distributed over a wider span, then the vertical capacity increases: for "2P": V max = 91.8 kN; for "4P": V max = 122.1 kN and for "q": V max = 152.4 kN. On the other hand, the more vertical load is distributed, the earlier the first hinge develops (u 6y : 0.216 m, 0.194 m and 0.187 m for the "2P", "4P" and "q" load patterns, respectively); however, this difference is not significant. Then, upon reaching the maximum vertical capacity, the beam in these three cases is transformed into a simply supported beam, which in turn modifies the bending moment distribution as presented in Figure 17. Now, the next plastic hinge is expected to occur at the point of the maximum bending moment which is the midspan of the simply supported beam (point D). In fact, Figure 18 shows that at u 6y = 0.832 m, the moment M D drops to zero, indicating the formation of the plastic hinge at point D. Figure 13 shows that up to this point, the beam for the "q" case, still demonstrates significant vertical load capacity of 108.2 kN (which is greater than the maximum capacities for the "1P" and "2P" load patterns. When it comes to the "2P" and "4P" cases, the computations have been stopped due to lack of convergence, which means that the last successful pseudo-time step takes place right before the occurrence of the second plastic hinge at point B and the resulting instability caused by the mechanism preventing convergence. Looking at Figure 17 can lead to the conclusion that the lack of convergence can also be caused by simultaneous formation of multiple plastic hinges because for the "2P" and "4P" cases, the maximum bending moment is not located at one particular point but spread over the distance of two meters (x ∈ 4 − 6 m ). This situation seems to be similar to the buckling bifurcation or snap-through issues encountered in frame or shell structures. From the physical viewpoint, it is interesting to compare the obtained capacity curves to the vertical pushover curves reported in many studies on the widespread subassembly (a two-span beam with two side columns and the removal of the central column) (e.g., Bao et al. [51], Alashker et al. [52], Lew et al. [78,79] Bao et al. [80]). This subassembly is considered as an extracted fragment of larger multi-span, multi-storey plane frames and is useful as a specimen for experimental tests. It can be summarised (according to Bao et al. [51]) that the capacity curve in such a subassembly is characterised by three stages: (1) the ascending stage responsible for the concrete arch effect, (2) the descending one indicating the development of a plastic hinge and finally (3) the second ascending stage due to the so called catenary action which ends with the collapse. In the first stage, the axial forces in the beams are in compression because the side columns resist the outward movement. However, in real frames, where many bays and storeys are involved, the outward movement of columns can be negligible. From the reinforced concrete cross section viewpoint, in this stage, the concrete cracking in the tensile zone is developed, and the steel bars exhibit increase in both tensile and compressive forces. The descending stage is responsible for further increase of concrete compressive stress up to crushing in the extreme fibres, and the cracks in the tensile zone continue to grow. The side columns gradually move back towards the initial position. In the second ascending stage, the catenary action develops because of (1) further growth of cracks, (2) change of rebar forces in the initial compressive zone from compressive to tension and (3) inward movement of side columns. As showed in Bao et al. [51], Alashker et al. [52] and Bao et al. [80], the realistic physical phenomena described above were both represented in both micro and macro finite element models, which proves that the macro finite elements are effective in reflecting the physical behaviour of reinforced concrete beams or columns under progressive collapse loading conditions. Finally, the collapse can happen due to, for example, the pull-out of rebars from the anchorage zones or the fracture of rebars (depending on the detailing). The numerical simulations and the described phenomena have been confirmed in full scale experimental tests ( e.g., Lew et al. [78,79] and Bao et al. [80], Gu et al. [81]).
In the light of the above, the ascending stage in Figure 13 is characterised by segments of lines with decreasing inclination because each of the two critical cross sections separately enters the hardening part of the moment-curvature relationship (see Figure 2). The arch effect (reflected by the compression in beams) is not present, here, because the boundary condition (the right roller support) does not impose any lateral resistance. In fact, in this stage, the axial force in beam (near midspan) is close to zero, and the horizontal displacement of the roller support exhibits a small negligible value to the right. Then the descending stage is abrupt, signifying that one of the critical sections reached the ultimate curvature due to the tensile rebar fracture (see Figure 2c), and the resulting capacity depends only on the other critical cross section. Finally, in the second ascending stage, the concrete cracking increases, and the catenary action is reflected by large vertical displacement and the slowly increasing horizontal displacement to the left. Looking at the forces in beams and horizontal displacement of the roller support reveals that in the re-ascending stage of Figure 13, the axial force in beam takes relatively small values of 0.3-0.8 kN, and the horizontal displacement of the roller support increases to 0.12 m when the vertical displacement u 6y reaches 0.83 m (the second plastic hinge and rebar fracture). It should be noted that without using the corotational formulation taking into account the geometric nonlinearity, both the arch effect and catenary action would be unidentifiable.
The insight from this two-span beam example can be helpful in understanding more complex frame structures, where yielding and plastic hinges are expected at several locations up to the point where the structure becomes unstable.
It should also be noted that in the many experimental tests (e.g., [78,79,81]), the increasing deformation is carried out by pushdown techniques concentrated at one node (the location where the removed column is connected to the remaining structure). Therefore, the current study sheds some light on potential underestimation in the ultimate capacity when performing the vertical pushover test. For practical reasons, it can be difficult to control deformation using actuators in many directions; hence, if only one node is loaded, the underestimation due to simplified load pattern should be taken into account.
Since this study focuses on beams, it is of interest to perform similar investigations pertaining to various configurations of frame structures (including other models of frame subassemblies) by applying realistic distributed load patterns.

Sensitivity Study
During the modelling stage, the values various material and geometric parameters (e.g., concrete compressive strength, concrete strain at peak stress, tensile concrete strength, concrete softening properties, steel yield and ultimate stress, uncertainty in final cross-section width, height and cover executed on building site, etc.) are assumed. These assumptions are due to limited knowledge or due to the inherent randomness of concrete or steel mechanical properties. Small changes in those parameters can produce significant unforeseen changes in local or global structural response. Therefore, it is valuable to perform sensitivity analysis to gain insights into the importance of parameters used in this study.
The calculation of sensitivity of a structural response with respect to uncertain properties indicates those parameters which have a significant influence on the structural behaviour. The sensitivity can be determined calculating the gradient of the structural response using the finite difference method or through direct differentiation with regards to parameters, which are included in the computational model of the structure and the applied loads. It is of utmost importance to identify parameters when their small change causes meaningful changes in structural response.
The structural response obtained using the finite element method depends highly on the assumptions of the model properties imposed in the analysis. These properties, specified by the parameters, describe the material, the finite elements it is made of, and the applied loads to the structure as well as the location of nodes. The parameters which describe the structural system can be collected in a vector x, and then, the equilibrium equation takes the form Both the vector of nodal forces P r and the vector of external forces P f depend on the parameter, explicitly. Additionally, the vector of nodal forces can depend implicitly on the parameters through the vector of nodal displacements U.
There are basically two approaches to calculate the gradient ∂U/∂x of nodal displacements with respect to changes in a parameter x. The first approach is based on the direct differentiation method, where the equations of the system are differentiated analytically [82][83][84] whereas the other approach utilises the finite difference method, in which the analyses are repeated with a small change (perturbation) of the parameter x + ∆x. Thus the gradient of the response can be evaluated as the difference between the response with and without the parameter perturbation divided by the parameter change It seems obvious that the direct differentiation method is more effective than the calculation based on the finite difference approach because the gradient of the response is obtained during a single analysis. Furthermore, the round-off errors are avoided, and no additional analyses with changes in parameter values are necessary. The detailed derivation of all required quantities in force-based element formulation both for linear and nonlinear geometry can be found in Scott [83], Scott et al. [84], Haukaas and Scott [85] and Scott [86].
Although the calculations of response sensitivities using the direct differentiation method or finite difference method can help understand what are the parameters raw sensitivities to a specified fixed parameter changes, for example, ±10% change, parameters have different accuracy in estimation which in turn leads to probabilistic approach, where they can be treated as random variable variables with a specified probability distribution function and its related mean and standard deviation characteristics. Therefore, in the following section, another approach, based on a particular reliability method, is briefly presented, and the notion of importance measures is presented after [87][88][89][90].

Sensitivity Based on the Theory of Reliability
Apart from the raw sensitivity methods, the theory of reliability ( [89,[91][92][93][94][95][96]) revealed other approaches to identifying which input parameters play a significant role in the structural response. To present the important measures which help identify the parameters contributing the most to the structural response, a brief introduction to reliability theory is needed. The considerations are narrowed to the class of time-invariant reliability problems.
Due to the well established literature notation of the standard normal random variable as a lowercase u, exceptionally in this section, a lowercase vector u stands for the vector of standard normal random variables and not the vector of element displacements as in earlier sections. However, still in this section, the uppercase U refers to the global displacement vector.
The structural reliability problem in the case of a single limit-state function (the so called component problem as opposed to the series, parallel and general system problems) can be formulated as n-fold integral where p f is an unknown probability, g is a limit-state function (also called a performance function) which determines the state of failure (unsafe state) for which the probability is to be found, and f (x) is the joint probability distribution function of random variables which constitute an n-dimensional vector x. In the reliability analysis using the finite element method, the limit-state function can be defined with respect to the response quantities, obtained from the finite element analysis. In engineering practical applications, the closedform solution of Equation (5) is very difficult to obtain, and therefore, the approximate solution can be found using various methods such as the First-Order Reliability Method (FORM), the Second-Order Reliability Method (SORM) or Monte-Carlo techniques. The thorough descriptions of these methods can be found in Ditlevsen and Madsen [93], Der Kiureghian [97] and Chang [98]. In particular, the FORM method seems to be appealing because it requires only a few numbers of limit-state function evaluations during the finite element method analyses. Additionally, the FORM method provides valuable information on the significance of particular random variables which represent material or geometric parameters of the computational model. In the FORM method, after transforming random variables x to the uncorrelated standard normal random variables u = u(x) and the limit-state function g to g u in the standard normal space, the probability in Equation (5) can be approximated as where ϕ(u) is the standard normal joint probability density function and the boundary of the limit-state function g u = 0 is approximated in the transformed standard normal space. When the limit-state function is nonlinear, the desirable approximation point is the point on the surface g u = 0, whose location is the nearest to the origin of the standard normal space. Such a point is called the design point u * , and results from the solution of the following minimisation problem u * = arg min{ u subject to g u = 0} (7) where 'arg min' is the argument of a function at which the function is minimised.
When g u (u * ) = 0, the linearised limit state function can be expressed as (8) where ∇g u (u) is the gradient row vector, α = −∇g u (u * )/ ∇g u (u * ) is the normalised negative gradient vector at point u * , and β = αu is the reliability index. Gradients of the limit-state function ∇g u = ∂g u /∂u are often used in the most effective algorithms which enable us to solve this optimisation problem. The differentiation of the limit-state function with respect to u can be expressed as where ∂x/∂u is the Jacobian of the probability transformation, ∂U/∂x stands for the gradients of the global structural response vector U with respect to parameters x, and finally ∂g u /∂U is the derivative of the limit-state function with respect to the response vector U.
The exchange of information between the reliability algorithm and the module of the finite element software consists in updating the finite element model by the realisation of random variables x and returning the global response vector U and its derivative ∂U/∂x to evaluate the limit-state function. Note that this technique has already been implemented and available in the OpenSees software framework Haukaas and Der Kiureghian [89], Scott and Haukaas [99].
After finding u * , the probability p f in the FORM method can be approximated using the following formula where Φ stands for the standard normal cumulative distribution function, and β is the reliability index (see, e.g., Haukaas and Scott [85], Ditlevsen and Madsen [93]). In civil engineering structure applications, typically, the marginal distribution of random variables and their correlation coefficients are relatively straightforward to identify (see, e.g., Benjamin and Cornell [100]). The most effective probability transformation of non-normal random variables x into the standard normal random variables u is the Nataf transformation [90,91,101,102]. The Nataf model can be characterised by two desirable features: (1) the model can take into account broad span of correlation values between the random variables x, and (2) the transformation to standard normal space is relatively straightforward regardless of the order of random variables.
In the Nataf transformation, the joint normal distribution f (x) can be calculated from the marginal probability distribution functions f (x i ) and the matrix of correlation coefficients R as follows where ϕ n (z n ) is an n-dimensional joint probability distribution function with zero mean value, unit variance and correlation matrixR. The term ϕ(u i ) is the standard normal probability distribution. The reciprocal relationship (e.g., Der Kiureghian [90], Liu and Der Kiureghian [91], Gu [96]) of a random variable x i by an arbitrary probability distribution function to a random variable z i (characterised by the standard normal distribution) can be expressed as The terms F i (x i ) and Φ(u i ) stand for the cumulative distribution function of variables (of distribution f (x i )) and standard normal probability density function ϕ(u i ).R =ρ ij whileρ ij and ρ ij are related together by the following expression whereρ ij are the correlation coefficients between z i and z j , and ϕ 2 (z i , z j ,ρ ij ) is the normal distribution function of a two-dimensional random variable with zero mean value and unit variance. In the iterative process, the coefficientsρ ij can be calculated from the prescribed values of ρ ij [90,91].

Importance Measures of Parameters
The FORM formulation has an interesting feature because it provides the so called importance and sensitivity measures which basically tell us how important are the individual random variables compared to the other random variables. To have a broader context of how these importance measures result from the FORM method, a brief derivation (after Hohenbichler and Rackwitz [87], Bjerager and Krenk [88], Haukaas and Der Kiureghian [89], Der Kiureghian [90]) is presented below.
Assuming that g u1 (u) = ∇g u (β − αu) is the linearised limit-state function, noticing that the mean value of u equals zero, and its covariance is a unit matrix, the mean and variance of g u1 (u) can be determined from where α = ∇g u / ∇g u is a unit vector, and β = µ g u1 /σ g u1 = αu * is the reliability index of the linearised problem. Moreover, Equation (15) shows that the components α 2 i are a function of the contribution of random variable u i to the total variance of the linearised limit-state function. It should be noted that the more significant the contribution, the more important is the influence of the random variable u i . In this way, the elements of α give relative measures of importance for a standard normal variable u i . In addition, analysing the expanded form of g u1 (u) = ∇g u (β − α 1 u 1 − . . . − α n u n ), it should be noted that the positive value indicates that the random variable has a loading nature while the negative value points to a nature of resistance.
If the main random variables in the reliability problem are statistically independent, there exist reciprocal correspondence between the initial random variables x i and standard normal random variables u i . Thus, the importance order and nature of random variables (loading or resistance) are similar to corresponding u i and can be obtained with respect to the vector α. If, however, the random variables x are statistically dependent, then there is no reciprocal correspondence. In such a situation, the vector α does not give information on the relative importance of the initial random variables. Therefore, to determine the measures of importance of the random variables, a linearised transformation is defined at the design point u * wherex can be regarded as a counterpart of x at the design point and differs inconsiderably from the vector x while the covariance matrix of u can be written asΣ = J −1 u,x (J −1 u,x ) T . Using Equation (16) in the formula for g u1 (u) and taking into account that β = αu * leads to g u1 (u) = − ∇g u αJ u,x (x − x * ). Hence, the variance of g u1 (u) can be expressed as whereD = diag[σ i ] is a diagonal matrix of standard deviations of the vectorx (e.g., Der Kiureghian [90]). The first component in the above expression includes the contributions to the variance g u1 (u) which result from the individual element variancesx, while the second component represents the contributions resulting from the covariances of random variable pairs. Finally, a vector representing the relative importance of the random variables x can be defined which is usually normalised to a unit vector. It can be proved that when the random variables are statistically independent, vector γ takes the form of vector α. Again, the formulas and the algorithm described in this section have been implemented in the OpenSees software framework [89,99] which is used in the calculations of the importance measures for the parameters of the analysed two-span beam.

Results of the Sensitivity Analysis
This section presents the results of sensitivity of input parameters on the nonlinear static response when applying to the vertical pushover analysis. Selected material and geometric parameters are treated as random variables with specified probability distribution, mean value, standard deviation or coefficient of variation. These probability properties are established based on recommended values found in literature (e.g., [85,89,103]). However, it should be noted that the assumed probability properties can differ and be updated if more accurate information is available.
The main concrete parameters subjected to probabilistic analysis are: compressive strength f c0 , ultimate compressive strength f cu and their associated strains ε c0 and ε cu . In this example, it is assumed that the mean values of these quantities are their nominal values. Therefore, for the compressive strength f c0 , the coefficient of variation equals cov = 0.1 which results in the standard deviation σ = 1.06 MPa. Similarly, it is assumed that for f cu , ε c0 and ε cu , the coefficient of variation is equal to 0.1. The coefficients of variation of the reinforcement bars f y and elastic modulus E s for steel are equal to 0.05 and 0.085, respectively. The area of each steel rebar is characterised by the coefficient of variation 0.01. The second group of uncertainty is related to both cross-section geometry and nodal coordinates. Uncertainty related to the concrete cross-section width and depth (height) is assumed in the standard deviation equal to 0.01 m which for the section width b and height h results in the coefficients of variation equal to 0.05 and 0.025, respectively. The standard deviation of the cover thickness is equal to 0.01 m which results in a relatively large coefficient of variation of 0.25. It is also assumed that all nodes have standard deviations 0.01 m from their nominal x and y coordinates. The summary of the uncertainty parameters and their probability distributions are given in Table 1. In this study, the overall number of material and geometric parameters mapped to random variables equals 28.
As described in Section 5.2, the first-order reliability method has a useful side-effect which allows us to investigate the importance of individual parameters in a particular loading scenario. In this Section, the first-order reliability method is applied to the vertical pushover analysis. The problem is formulated as follows: what is the probability to exceed a specified threshold value of vertical displacement at node B when several material and structural parameters are treated as random variables (see Table 1). The limit-state function is specified as: g = 0.072 − U B . However, unlike the regular reliability analysis using FORM, in this study, the probability of failure is not as interesting as the importance measures.
First it is assumed that all random variables are uncorrelated and, after running the first-order reliability analysis, the following importance measures in the form of the components of the vector γ (see Equation (18)) are presented in the 4th column of Table 2.  If assuming that the random variables for the two bay elements are correlated with the correlation coefficients specified in Table 1, then the ranking of parameters can be found in the 5th column of Table 2. From this table, we can observe that the correlation does not change the order of the ranking and only slightly changes the values of γ and corresponding percentage contribution (γ 2 × 100%). Therefore, looking at columns 4 and 5, we can find that the most sensitive parameters are the cover thickness in element 1 (1st position: 43-44%) and element 2 (3rd position: 15-16%). This is mainly due to a relatively large coefficient of variation. Moreover, very high in the ranking is the depth h of element 1 (2nd position: 18%) and element 2 (5th position: 6%) because this parameter is part of flexural rigidity (I = bh 3 /12). Next in the ranking, we can find the steel yield stress f y for element 1 (4th position: 8%) and element 2 (8th position: 1%), followed by the modulus of elasticity E s of element 1 (6th position: 4%) and element 2 (7th position: 2%). It is interesting that the concrete compressive strength f c is at 9th and 10th positions (less than 1%) in the ranking. The positive values of γ mean that the parameter plays the role of resistance, while the negative, the role of loading.
On the other hand, we assume that only one random variable is assigned to both elements 1 and 2, and the obtained results are presented in Table 3. This table summarises the averaged importance of parameters regardless of their location in the structure. The most important parameter for the whole structure is the concrete cover thickness followed by the concrete section depth h, the steel yield stress f y , the elastic modulus of steel E s and the compressive concrete strength f c0 with relative contributions of 59%, 24.4%, 7.6%, 6.6% and 1.2%, respectively. While all other parameters (including the concrete strain ε c0 , the concrete section width b, the bottom and top reinforcement areas A s , the X and Y coordinates and the concrete ultimate stress and strain) have a contribution of less than 2%. It is interesting to note that similar investigations on the importance measures for steel and reinforced concrete structures have been published in Haukaas and Der Kiureghian [89] and Haukaas and Scott [85]; however, they deal with the classic horizontal pushover analysis as opposed to the vertical one performed in this study. For instance, in Haukaas and Der Kiureghian [89], the most important parameters were horizontal loads, compressive strength and imperfections in horizontal nodal coordinates, whereas in Haukaas and Scott [85], the highest importance measures are related to concrete cross-section heights, the yield stress of rebar steel and a cross-section concrete cover. Nevertheless, direct comparison between these studies cannot be made mainly due to different model assumptions and probability characteristics. The author's research plans include a comparison of importance measures of a nonlinear inelastic frame model under horizontal (seismic) and vertical (progressive collapse) pushover analyses.

Conclusions
A parametric analysis of the effects of a sudden support removal of a two-span reinforced concrete beam has been presented. Linear and nonlinear dynamic responses are analysed in detail and supported by vertical pushover computations and reliability-based sensitivity studies.
In the linear dynamic analysis, from the study of the influence of support removal rate, it can be concluded that the time of support removal has an important role when it is up to 30% of the first natural period of the beam without the inner support. In this case, the support removal can be regarded as "sudden", and it still generates substantial dynamic effects (dynamic factor of about 1.8). Moreover, the plot of the maximum vertical displacement versus the time of the support removal (response spectrum) presents a characteristic pattern (Figure 6) where the value of maximum vertical displacement is equal to the static vertical displacement when the time of support removal is the multiple of the first natural period. The maximum dynamic factor for the shortest time of support removal reaches the value of 1.8 when (5%) damping is applied in dynamic calculations.
When the material nonlinearity is included in the analysis, the characteristic pattern in the linear response spectrum is less pronounced and can be explained by the following facts: (a) in the nonlinear analysis, the stiffness contributing to the mode shape is calculated according to the current tangent stiffness and thus is responsible for the constant change in natural frequencies, which in turn desynchronise the contribution of mode shapes to the total structural dynamic response; (b) in the nonlinear cyclic analysis, at each deformation and loading cycle, the energy is dissipated providing additional damping to the assumed viscous damping and further influences the current natural frequencies and mode shapes (see, e.g., Rodrigues et al. [104]).
The vertical pushover analysis has been adopted from the horizontal pushover analysis from earthquake engineering. Four cases (one, two, and four concentrated loads and uniformly distributed load) of reference vertical load schemes have been considered to obtain the capacity curves of the vertical loading versus the vertical displacement at the point of the removed support. Moreover, the bending moments have been monitored at three critical nodes (the left fixed support, the midspan of the resulting one span beam, and the node of the removed support). It has been observed that the largest vertical capacity is obtained for the uniformly distributed load (100 %), while for "4P", "2P" and "1P", the capacity is reduced to 80%, 60% and 51%, respectively. It is also noted that the failure mechanism depends on the vertical load distribution, and for the case "1P", the first plastic hinge occurs at the point of removed support, while for "2P", "4P" and "q" cases it occurs at the left fixed support which is related to the bending moment distributions.
It is also noted that comparing the current study to the studies of widespread subassemblies, ( e.g., Bao et al. [51], Alashker et al. [52], Lew et al. [78,79], Bao et al. [80]) representing an extracted fragment larger frames; in the case of two-span beam with one the roller support, (1) no arch effect is observed, (2) plastic hinges develop in discrete points equal to the number of the degree static indeterminacy and (3) the catenary action is only reflected by large vertical deformation and horizontal movement of the roller support.
A sensitivity analysis based on the importance measures of the first-order reliability method has been performed, and the result of all material and geometric parameters of the model has been presented using specified probability properties (probability distribution, mean and standard deviation). First, the sensitivity analysis results of the correlated and uncorrelated random variables are similar. Second, the first span parameters have higher sensitivity than the parameters of the second span; when only one parameter is specified for both spans, the most sensitive parameters when modelling and analysing the two-span beam of this paper are the geometrical dimensions of the cross-section, followed by the steel mechanical properties, while other parameters play a minor role in the sensitivity analysis. Note also that these first four group of parameters represent more than 97% of the total contribution of all parameters.
The assumption of the boundary conditions can be further investigated by changing the clamped support to a pinned one and adding rotational spring with its rotational stiffness factor k ϕ . This additional parameter can be treated as random and can be incorporated in the sensitivity analysis which can be a subject of further investigations.
The results of this paper can help in planning experimental and numerical analyses of more complicated systems susceptible to progressive collapse.
It should also be emphasised that the results presented in this study, containing hundreds of finite element simulations, can be compared to the analyses using more sophisticated continuum-type finite elements (micro-modelling), which, however, will pose other modelling challenges such as computer memory and analysis time consumption. Moreover, a validation experimental test would be of interest; however, depending on the number of analysed support/column removal durations, it can face challenging, prohibitive costs or alternatively only selected representative cases could be validated.

Data Availability Statement:
The data presented in this study are available on request.