APPLICATION OF AEM IN PROGRESSIVE COLLAPSE DYNAMICS ANALYSIS OF R.C. STRUCTURES

The Finite Element Method (FEM) and the other numerical strategies are viably actualized in linear and non-linear analysis of structures. Recently, a new displacement based on Applied Element Method (AEM) has been developed. It is applicable for static and dynamic for both linear and non-linear analysis of framed and continuum structures. In AEM, the structural member is partitioned into virtual elements connected through normal and shear springs representing stresses and strains of certain portion of structure. FEM assumes the material as continuous and can indicate highly stressed region of structure, however it is difficult to model separation of element unless crack location is known. The main advantage of AEM is that it can track the structural collapse behavior going through all phases of the application of loads.


INTRODUCTION
Progressive collapse became a hot point for research because of the several collapses that have occurred since the beginning of the twentieth century.For example, progressive collapse that happened in Ronan Point apartments tower -London in 1968, 2000 Commonwealth Avenue building -Boston in 1971, Murrah Federal Office building -Oklahoma in 1995 and World Trade Center -New York in 2001.
Therefore, many international structural codes and standards started to consider the progressive collapse resistance, such as the General Services Administration (GSA) [1] published in 2003, and the Unified Facilities Criteria (UFC) [2] by Department of Defense (DoD), USA published in 2005 and 2009.Precautions can be taken in the new design of structures to limit the effect of the local failure and resist progressive collapse.According to DoD guidelines, two general approaches are used for reducing the possibility of progressive collapse; Direct and Indirect Design.For the Indirect Design approach, the structure resistance to progressive collapse is considered implicitly through the provision of minimum levels of strength, continuity and ductility.Direct Design incorporates explicit consideration of resistance progressive collapse through two methods.One is the Alternate Path Method, which requires that the structure be capable of bridging over a missing structural element, with the resulting extent of damage being localized.The other method is the Specific Local Resistance Method, which seeks to provide sufficient strength to resist a specific threat.
Simulation of progressive collapse behaviour is not an easy job using available numerical techniques.Therefore, the use of highly efficient numerical modelling procedures for progressive collapse of reinforced concrete (R.C.) structures is needed.The available numerical methods for structural analysis can be classified into two categories.In the first category, model is based on continuum material equations.The Finite Element Method (FEM) is typical example of this category.The mathematical model of the structure is modified to account for reduced resistance of yielding components.However, to perform collapse far exceeding their elastic limit is a difficult task to many numerical methods.There are several limitations in adopting FEM in high non-linear case where crack has initiated and element is not detached from the structure, smeared crack approach is followed.However, smeared crack approach cannot be adopted in zones where separation occurs between adjacent structural elements.Therefore, FEM assumes the material as continuous and is able to indicate highly stressed region of structure but it is difficult to model separation of element unless crack location is known.While, Discrete Crack Methods assume that the location and direction of crack propagation are predefined before the analysis.With this group of the methods, analysis of concrete structures can be performed at most before collapse.The second category of methods uses the discrete element techniques, like the Rigid Body and Spring Model, RBSM [3] and Extended Distinct Element Method, EDEM [4].The main advantage of these methods is that they can simulate the cracking process with relatively simple technique compared to the FEM.On the other hand, the main disadvantage is that crack propagation depends mainly on the element shape, size and arrangement.Analysis using the RBSM could not be performed up to complete collapse of the structure.Whereas, the EDEM can follow the structural behaviour from zero loading and up to complete collapse of the structure.
The Applied Element Method (AEM) is simple in modelling and programming and has high accuracy of the results with relatively short CPU time.Also, the response of the structure can be followed all the way to collapse involving the elastic stage, crack initiation and propagation, reinforcement yielding, large deformations, element failure and separation, rigid body motion of falling elements, collision between elements, and impact forces resulting from falling debris.This method can be easily applied to a wide range of applications; both for small and large deformation ranges under static or dynamic loading conditions, and for linear or non-linear materials.In this paper, AEM is used to perform the analysis.
There are many available researches that used AEM.Helmy et al. [5], used Extreme Loading for Structures (ELS) software which is based on AEM and following GSA and UFC guidelines to study the progressive collapse resistance of 10-storey reinforced concrete building.Helmy et al., [6] used ELS to investigate the effect of slab modelling considering slab reinforcement and slab thickness on progressive collapse.Lupoae and Bucur [7 and 8], presented a controlled implosion for structure collapse behaviour using ELS.Lupoae et al. [9], discussed the effect of infill walls provided on periphery of frame with and without openings on progressive collapse using ELS.Salem et al. [10], used AEM to propose an optimum and economic way of design to prevent progressive collapse of multi-storey reinforced concrete structure due to column removal.Khalil [11], studied a progressive collapse potential of 4-storey moment resisting steel building using AEM by following alternate load path method specified in UFC guidelines and made a comparison between results obtained by AEM and FEM.Dessousky [12], discussed a special algorithm to model the interface between the blocks and added to an Applied Element-Based solver.These algorithms predicted the strength and stiffness at interfaces when cracks opened and closed.Tagel-Din and Meguro [13], gave a simulation of collapse processes of scaled reinforced concrete structure by using AEM and compared with results obtained by shake table experiments.Raparla et al. [14], considered four bare frames designed as per Indian Standards for studying performance of building up to collapse using AEM.All bare frames were subjected to Northridge earthquake ground motion (freq.1-4 Hz).Vikas et al. [15], performed linear static analysis of portal 2-D frame subjected to lateral loading to study the effect of mesh size and spring numbers and compared the results with FEM.Ismail, [16] studied the effect of the element size, spring distribution and shear stiffness of the spring on the accuracy of the AEM in obtaining the buckling load of a simply supported rectangular plate and compared the results with the theoretical solution.

OBJECTIVE
The objective of this paper is to verify the capability of AEM in simulating the non-linear dynamic analysis for progressive collapse of the experimental work conducted by Park et al., [17].The 1/5 scaled 3 and 5 stories reinforced concrete structures are modelled using ELS software that follows the AEM and compares the numerical results with the experimental ones.Moreover, a study has been made to investigate the effect of contribution of the R.C. slab on the overall collapse mechanism.

APPLIED ELEMENT METHOD OVERVIEW
The Applied Element Method (AEM), was developed by Tagel-Din and Meguro [13].This method combines the advantages of Finite Element Method (FEM) and the Discrete Element Method (DEM).AEM discretizes the domain into a grid of rigid finite elements with six degrees of freedom, three translations and three rotations that represent the rigid body motion of the element located in the geometric centre of the element as shown in Figure 1.The connection of the elements is established through a mesh of springs on the contact faces of the elements, Meguro and Tagel-Din [18].Two elements shown in Figure 1-b are assumed to be connected by one normal and two shear springs located at contact points, which are distributed around the elements edges.Each group of springs completely represents stresses and deformations of a certain representative volume of the element as illustrated in Figure 1-b.The spring stiffness is determined as shown in Equations 1 and 2.
Where, Kn is the stiffness of normal spring; Ks is the stiffness of the shear spring; d is the distance between the springs; t is the element thickness; a is the length of the representative area; and E and G are the Young's modulus and shear modulus of the material, respectively.To consider the effects of collision, it is necessary to check the collision between elements during analysis.To simplify the problem, element shape is assumed as circle during collision.This assumption is acceptable if the element size is relatively small.Even in case of relatively large elements used, it may be reasonable because in the deformation range of collision, the sharp corners of elements are broken due to the stress concentration and the edge of the elements becomes rounded shape.Based on these assumptions, only distance between the centres of the elements is calculated to check collision, Tagel-Din [19].Contact and collision forces are considered in the analysis through collision springs, which are added at the locations where collision occurs.These effects can be considered easily with high accuracy and without large increase of the CPU time.
Figure 2 shows the constitutive models adopted in AEM.For modelling of concrete under compression, Maekawa [20] and Okamura and Maekawa [21], model is adopted.In this model, the initial Young's modulus, the fracture parameter and the compressive plastic strain are introduced to define the envelope for compressive stresses and strains.Therefore, unloading and reloading can be conveniently described.The tangent modulus is calculated according to the strain at the spring location.After peak stresses, spring stiffness is assumed as a minimum value (1% of initial value) to avoid negative stiffness.This results in difference between calculated stress and stress corresponds to the spring strain.These residual stresses are redistributed by applying the redistributed force values in the reverse direction in the next loading step.
For concrete springs subjected to tension, spring stiffness is assumed as the initial stiffness until reaching the cracking point.After cracking, stiffness of springs subjected to tension is set to be zero.The residual stresses are then redistributed in the next loading step by applying the redistributed force values in the reverse direction.For concrete springs, the relationship between shear stress and shear strain is assumed to remain linear till cracking of concrete.Then, the shear stresses drop down as shown in Figure 2 Fig. 2 -Tension, compression and shear models for concrete [20 & 21] For reinforcement springs, the model presented by Ristic et al. [22], is used as shown in Figure 3.The tangent stiffness of reinforcement is calculated based on the strain from the reinforcement spring, loading status (either loading or unloading) and the previous history of steel spring which controls the Bauschinger's effect.The main advantage of this model is that it can consider easily the effects of partial unloading and Baushinger's effect without any additional complications to the analysis.After reaching 10% tensile strain, it is assumed that the reinforcement bar is cut.The force carried by the reinforcement bar is redistributed, when it reaches the failure criterion, Tagel-Din and Meguro [23], by applying the redistributed force to the corresponding elements in the reverse direction.

Dynamic Equation Formulation in of Large Displacement
In the AEM, the equation of motion under any dynamic loading is described in Equation 3, [13 & 24]: Where, [M] is the mass matrix, [C] is the damping matrix, and [K] is the matrix of global nonlinear stiffness, ∆̈, ∆ ̇and ∆u are the incremental acceleration, velocity, and displacement vectors, respectively, ∆f(t) is the incremental applied load vector, RG is the residual load vector due to geometric changes and, RM is the residual load vector due to cracking or incompatibility between spring stress and the corresponding strain.The effects of material nonlinearity are considered in RM and [K], [24].
After applying a small load increment, the structure geometry is modified and hence, incompatibility between external loads and other forces occurs.The main difference between the AEM and the conventional methods is that the geometrical stiffness matrix is omitted and its effects were replaced by adding the geometrical changes effects as an additional load vector RG as shown in Equation 4.
It should be emphasized that this technique can be applied in both static and dynamic loading conditions.In case of static loading condition, the mass and damping matrices are set equal to zero.The main limitation in static analysis is that separation of elements is prohibited during analysis as it makes the stiffness matrix singular.On the other hand, analysing structures subjected to dynamic loading condition enables us to follow both geometrical changes of the structure and the rigid body motion during failure.As the deformations are assumed to be small in each load increment, small time increment should be used.In the following subsections, the matrices in Equation 3 are presented.

Formulation of stiffness matrix
The stiffness matrix components corresponding to each degree of freedom are determined by assuming a unit displacement in the studied direction and by determining forces at the centroid of each element.The element stiffness matrix size is (6x6).The matrix in Equation 5shows the components of the upper left quarter of the stiffness matrix [24].All used notations in this matrix are shown in Figure 4.It is clear that the stiffness matrix depends on the contact spring stiffness and the spring location.(5) The stiffness matrix in Equation 5 is for only one pair of contact springs.However, the global stiffness matrix is determined by summing up the stiffness matrices of individual pair of springs around each element.Consequently, the developed stiffness matrix is an average stiffness matrix for the element according to the stress situation around the element.This technique can be used in load and displacement control cases.

Formulation of mass matrix
The mass of the rigid element is lumped at the element centre.As the size of the AEM element is relatively small, the effect of the lumped mass is similar to that of a distributed mass.The mass matrix of a square element of thickness t is determined in Equation ( 6) [24]: Where, D is the element size, ρ is the material density, M1 and M2 are the element mass in X and Y directions, and M3 is the element mass relative to the moment of inertia about the axis which passes through the element centre.The mass matrix is diagonal positive definite, so the diagonal elements of the mass matrix should be also positive.The addition of the mass matrix to the stiffness matrix in dynamic analysis prevents the occurrence of a singular stiffness matrix, especially where element separation occurs.Element separation is prohibited in static analysis, due to the absence of the mass matrix.

Formulation of damping matrix
In the non-linear response stage of R.C. structures, internal damping can arise due to the following [24]: 1. Concrete cracking.
2. Energy dissipation due to the loading and unloading of compression springs.
3. Unloading of reinforcement after yielding.4. Energy dissipation due to crack opening and closure. 5. Friction between elements during contact.6. Unloading factors during contact.
The sources of internal damping forces mentioned above, are automatically considered in the non-linear response of structures using the AEM.In the elastic range, these effects from internal damping forces are neglected.Thus, an external damping matrix, [C], is required at this stage to account for internal damping [25].The damping matrix is mass proportional and depends on the first mode of vibration, as in Equation 7.
Where, ξ is the damping ratio, and ω1 is the first natural frequency.

Solving the Nonlinear Dynamic Equation
The non-linear Equation 3 is implicitly solved utilizing a step-by step integration Newmark-Beta approach [26 and 27].The equilibrium equations represent a linear system of equations for each step.The solution of the equilibrium equations is commonly solved using Cholesky upperlower decomposition.Once elements are separated, the stiffness matrix becomes singular.However, the stability of the overall system of equilibrium equations is kept because of the existence of the mass matrix.Separated elements may collide with other elements.In that case, new springs are generated at the contact points of the collided elements.

VALIDATION OF AEM USING ELS
In this paper, Extreme Loading for Structure software (ELS) which is based on the AEM is used [25].In order to validate the AEM models, the analytical results are compared with existing experimental data obtained from previously performed tests.The details of the comparisons for each case are described in this section.The effects of failure criteria on the collapse process are investigated, since these are anticipated to be one of the key parameters that can control the collapse phenomenon.Experimental test conducted by Park et al. [17] at Chonbuk National University is modelled using ELS (V.2.3) to validate the AEM.The R.C. structures for 3-story and 5-story models are illustrated in Figure 5.The reinforcement details of the structural components and the measured material properties are presented in Tables 1 and 2. Fig. 5 -Model dimension, cross section of members, and arrangements of steel reinforcement (mm) [17] Tab. 1 -The steel reinforcement details [17] Type Sudden failure columns

Modelling and Simulating
As in the experiment, the numerical simulation of a sudden failure of the three floors structure, the two main columns labelled C2 at first-storey are performed in a dynamic after the self-weight of the structural components was completely loaded.Also, for the five floors structure, the numerical simulation of a sudden failure of the two main columns labelled C2 at first-storey and secondstorey are performed in a dynamic after the self-weight of the structural components was completely loaded (refer to the marked columns in Figure 5).However, the sensitivity analysis for mesh, springs, time steps and calibration of the constitutive models should be firstly done.

Sensitivity analysis
A mesh sensitivity analysis is performed for the experimental works to study the effect of the element size and the number of connecting springs between elements.The experiments are analysed using three increasingly smaller-sized elements (approximately cubic shapes with edges of 100, 50 and 25 mm).For each different element size, two models are considered using 5 and 10 connecting springs for each pair of adjacent element faces.The convergence is considered to be achieved when the changes in the load-displacement result from one analysis to the next are too small to be visually noticeable.For all experiments, the estimates for load-displacement achieved for a mesh with edge of 50 mm and 5 connecting springs between element faces and therefore, this combination is considered an appropriate mesh and used for all tests reported in this paper.
A similar analysis is performed to calibrate the time step on model.Time effects are continuous from start till end of analysis.As the nonlinear dynamic phenomenon is very complicated to be solved using exact solutions, approximate numerical solution for the dynamic equation is adopted.The numerical solution is based on assumption of a small-time step that can follow the structural behaviour.The selection of the time step is very important.Having too short time step will result in very long analysis time.Whereas, using a large time step will result a short analysis time as well as a less accurate analysis.The numerical solution may fail to converge if the selected time step is large.Time -displacement curves are obtained using 0.10, 0.05, 0.01, 0.005, 0.001, 0.0005, and 0.0001 seconds loading steps and no noticeable differences were obtained between the resulting capacity curves.So, a time step of 0.001 second is used.

Calibration of the constitutive models
Due to the lack of information for steel stress-strain curves used in the experiments, the material characterization under tensile tests were simulated numerically and used to calibrate the parameters defining the Ristic constitutive model of steel [22].The mechanical properties presented in experimental were used, considering a Young's modulus, a shear modulus presented in each experimental and a post yield stiffness factor of 0.01.
All concrete material characterization is used to calibrate parameters defining Maekawa [20], constitutive model for concrete, especially normal and shear contact stiffness ratio and also loading ratio.Many groups are examined and it is found that, the best group giving good simulation with experiment is [0.0001,0.000001,10].Also, damping will be calibrated as 0%, 1%, 2%, 4%, and 5% damping ratios are examined.It is found that the best ratio which gives good simulation with experiment is 0%.Another important parameter is the separation strain.This parameter defines the strain value in the springs located between two neighbouring elements where the elements are considered to be physically separated.According to the ELS modelling manual [25], for reinforced concrete elements, the separation strain should be higher than the ultimate tensile strain of the rebar.Mesh division and reinforcement details similar to experimental works are shown in Figure 6.

RESULTS AND DISCUSSION
Figures 7a and 7b show the collapse shapes for the three-story structure at every 200 ms obtained from experiment and ELS model performed by Park et al. [28] and the author.From these figures, it can be observed that the collision with the ground occurs from time 800 ms to 1200 ms.Also, no plastic hinges were formulated until the building completely falls on the land at 1200 ms.The first cut off steel between column of second floor and first floor beam occurred at 1400 ms.The comparison shows that the results obtained from present model are similar to those obtained either experiment or analytically by Park et al. [28].

Fig. 7b -Experimental and analytical collapse shapes for 3-storey building
Similarly, Figures 8a and 8b show the collapse shapes for the five-story structure at every 200 ms obtained from experiment and ELS model performed by Park et al. [28] and the author.It can be observed that, the collision with the ground occurs from time 1000 ms to 1600 ms.Also, no plastic hinges formulated until the building completely falls on the land at 1000 ms and the first cut off steel between column of second floor and first floor beam occurred at 1400 ms.Again, the comparison shows that the results obtained from present model are similar to those obtained either experimentally or analytically by Park et al. [28].

EVALUATION OF CONTRIBUTION OF FLOOR SLABS
Numerical models are investigated where slabs do not exist for both three and five-story buildings, as depicted in Figures 7 to 10.In these models the beam weights are increased to account for slab dead loads.Hence, the density of the concrete beams is taken as 40 kN/m 3 instead of 25 kN/m 3 .From Figures 7 to 10, it can be noticed that the numerical ELS results with and without considering concrete slabs are matching the experimental results until the structure completely collision with the ground.After that, the collapse of the model without considering slab is different from both the experiment and the model considering slab because of the plastic hinge formulation in different places.The floor slab connects all beams and acts as a rigid diaphragm.On the other hand, every beam in the floor behaves like an individual element when slabs do not exist.Also, the contribution of slabs decreases the beams deflection and straining actions.Therefore, it is

Fig. 1 -
Fig. 1 -Modelling in AEM[18] -b.The level of drop of shear stresses depends on the aggregate interlock and friction at the crack surface.d i n g R e l o a d i n g U n l o a d i n g

Fig. 6 -
Fig. 6 -ELS Models for a) three-story with slab model, b) three-story without slab model, c) fivestory with slab model, and d) five-story without slab model