Flexible Multibody Dynamics Finite Element Formulation Applied to Structural Progressive Collapse Analysis

This paper presents a two-dimensional frame finite element methodology to deal with flexible multi-body dynamic systems and applies it to building progressive collapse analysis. The proposed methodology employs a frame element with Timoshenko kinematics and the dynamic governing equation is solved based on the stationary potential energy theorem written regarding nodal positions and generalized vectors components instead of displacements and rotations. The bodies are discretized by lose finite elements, which are assembled by Lagrange multipliers in order to make possible dynamical detachment. Due to the absence of rotation, the time integration is carried by classical Newmark algo-rithm, which reveals to be stable to the position based formulation. The accuracy of the proposed formulation is verified by simple examples and its capabilities regarding progressive collapse analysis is demonstrated in a more complete building analysis.

Several terrorist attacks in the end of 20th and start of 21st century attempted against public or business buildings, featuring the World Trade Center collapse in New York in 2001 as one of the most economically and emotionally affecting terrorist action.Nowadays buildings in populated areas are in general considered potential targets for terrorist attacks, demanded more interest in the design of buildings to resist to explosions avoiding progressive collapse.
The problems considered here are initially a continuum media or composed by several continuous parts linked in a statically stable way, however, due to rupture or to some applied action over the structure, the structural members are converted into a set of separated bodies.This situation requires a multibody simulation.
Multibody dynamics studies are related to dynamic analysis of interacting bodies.Recent studies on multi-body dynamics are dedicated to several fields, from structural engineering to biomechanics, and include the search for real time simulation, study of contact and impact problems, extension to electronics and mechatronics, dynamic strength analysis, optimization of design and control devices (see e.g.Shiehlen et al. (2006) and Shabana (2013)).
In this paper, we develop a finite element strategy based able to model flexible multi-body which can be applied to simulate progressive collapse problems.This strategy consists in assembling finite elements using Lagrange multipliers, so that the connections can be dynamically broken by enforcing the corresponding Lagrange multipliers to zero.
One can consider this formulation, as described here, conservative, as structures are considered to fail in the elastic domain, with rupture criteria based on standard design equations leading to a fragile behavior of the analyzed structure.However, it is possible to extend the proposed formulation to consider ductile or combined ruins introducing plastic joints, Coda and Paccola (2014) and Reis and Coda (2014).
To perform this analysis we employ a frame finite element geometric nonlinear dynamics formulation, which has the following characteristics:  (i) Absence of large rotation description naturally resulting in an energy conserving time integration algorithm;  (ii) Natural achievement of high-order curved elements;  (iii) Adequate kinematics for thin or thick frame elements developing large displacements.Such characteristics are achieved by using positions and generalized vectors as nodal parameters instead of displacements and rotations.This formulation, motivated by the work of Bonet et al. (2000), has being successfully applied in various applications as one can see, among others, in Coda and Greco (2004), Coda and Paccola (2011), Carrazedo and Coda (2010), Sanches andCoda (2013 and2014).
Energy-conserving is one of the most controversial subjects related to bars and shells nonlinear dynamic analysis.This controversy occurs because the most employed formulation (co-rotational description) uses finite rotations as degrees of freedom.Finite rotations are objective only when small increments are adopted.Moreover, co-rotational formulations result in a non-constant mass matrix, which forbids the use of well-established time integration procedures for linear analysis.In this sense, special time integrators for co-rotational formulations were developed, as it can be seen Latin American Journal of Solids and Structures 14 (2017) 52-71 in Ibrahimbegović andMamouri (2000 and2002), Ibrahimbegović et al. (2003), Romero (2008), Ghosh and Roy (2009), Betsch and Steinmann (2003), Jelenić and Crisfield (2001) and Romero and Armero (2002).Most of those studies state that the Newmark time integration procedure cannot be used in nonlinear dynamics as a whole, however, finite rotation based methods are not the unique strategy to solve nonlinear dynamics.
The total Lagrangian frame formulation based on positions and unconstrained vectors employed in this work, avoiding the use of large rotation schemes, enables the use of Newmark time integration procedure in nonlinear dynamic multibody problems with large displacements and rigid body rotations.In this case the formulation with constant mass matrix is equivalent to the one studied by Hughes (1976), where one can find a good description of stability and energy conserving properties for the average acceleration time integration (equivalent to the Newmark version employed here).In Sanches and Coda (2013) the authors give a proof of the linear and angular momentum conserving property and perform stability and energy conservation tests for the Newmark  integration procedure including large rotations but keeping small strains for the position based formulation.
In section 2 we describe the employed frame finite element as well as the nonlinear dynamics solution procedure, in section 3 we develop the procedure for dynamical splitting of structural members with Lagrange multipliers and discuss the system solving procedure.In section 4 we test this methodology with a simple cantilever example and then we apply it to multi-storey building progressive collapse problem.

FRAME FINITE ELEMENT FOR NONLINEAR DYNAMICS
The unconstrained vector frame element employed here is designed to model large rotation problems without lack of objectivity and conserving energy in general nonlinear dynamic applications with small strains.
The employed methodology is based on the stationary potential energy theorem written regarding nodal positions and generalized unconstrained vectors instead of displacements and rotations, inspired by Bonet et al. (2000).As mentioned before, this characteristic avoids the use of finite rotation approximations, what makes it very simple especially for 3D cases.The frame formulation is total Lagrangian, and due to its unconstrained vector mapping, it presents constant mass matrix and therefore, it is possible to apply the Newmak  integrator as a momentum conserving algorithm, like in a geometric linear problems with physical non-linearity (see Sanches and Coda (2013) for more details on the position based formulation and Newmark  integrator).
Frame structures consist of solids with one dimension much larger than the other two.Therefore, the frame kinematics can be written considering the middle surface mapping as a reference, as depicted in Fig. 1, and then employ kinematics hypothesis to map the whole element.
Latin American Journal of Solids and Structures 14 (2017) 52-71 The mappings m0 f and m1 f , from the auxiliary non-dimensional configuration (I) respectively to the initial (II) and final configurations (III), have its components written as follows: and Where m ji X and m ji x are the j -th nodal value of middle surface position in i direction respectively for initial and current configuration.
In order to map completely the solid domain, the Timoshenko-Reissner kinematics is adopted.We consider a generalized vector which, in the initial configuration ( 0 g ) is normal to the neutral surface, and is not necessarily normal to the current neutral surface ( 1 g ), allowing distortion due to shear stresses.The generalized vector is unconstrained, and so, allows deformation in the height direction (see Fig. 2).In this sense, the Cartesian components of the mapping of any point in the whole element from non-dimensional to initial and to current configurations are given respectively by: and The changes in the cross section height direction can be better represented by introducing a nodal enrichment to consider linear strain rate in this direction, so that the vectors 0 g and 1 g components can be written as: and where 0 ij e are the j -th nodal values for the i component of the unity vector normal to the initial configuration middle surface, 1 ij e are the j -th nodal values for the i component of the generalized vector (not necessarily normal or unity) at current configuration, 0 h is the initial cross section height and j a is the j-th nodal value of strain rate along thickness (see Fig. 2).
Considering two-dimensional frame elements, at any point of coordinate  the vectors 0 e and 1 e can be written as: and where 0  and 1  are the tangent angles to the neutral surface respectively at initial and current configuration.
Finally, the mapping from initial to current configurations is represented by: where X is the vector of nodal values at initial configuration.
Considering the gradient matrixes 0 A and 1 A of initial and current mappings (3) and ( 4), with it components given by: Latin American Journal of Solids and Structures 14 (2017) 52-71   0 and where the index " , j " indicates derivatives regarding 1  if After evaluating the gradient A , the Green strain tensor is given by: where the variables ij C and ij  are the right Cauchy-Green stretch tensor and the Kronecker delta, respectively.The following quadratic strain energy per unit of initial volume is obtained considering Saint-Venant Kirchhoff constitutive law: Such constitutive law relates second Piola-Kirchhoff stress tensor and Green strain tensor according to: where the elastic tensor D is given by: with G and  being respectively the transverse modulus of elasticity and the Poisson ratio.
The true stress tensor (Cauchy stress) is achieved directly from the Second Piolla-Kirchhoff stress following (see Ogden (1984) for details): The finite element adopted is an isoparametric curved line with arbitrary number of nodes.The shape functions are Lagrange polynomials of degree 1 n  where n is the number of nodes of the considered element.Each node has 5 nodal parameters: 2 middle surface position vector components i x with = 1 i or 2 , the current angle to obtain the non-normal generalized position vector component 1 e , the current thickness 1 h and the linear strain rate along thickness a .In several problems with small strains and or small Poisson ration, the user may chose to use only 3 nodal parameters (current nodal positions and tangent angle) saving computational effort.

Time Integration
As mentioned before in section 1, we use the Newmark  algorithm for numerical time integration.The time marching process is summarized as: and where s x is the unknown nodal values vector at instant t  s and the over dots represent time derivatives.Sanches and Coda (2013) proved that for the positional total Lagrangian description, the Newmark  with = 1/2  and = 1/4  presents momentum and energy conservative properties, and so we use these constants values, making the algorithm equivalent to the trapezoidal rule, for which stability and energy conserving properties for constant mass matrix nonlinear problems are studied by Hughes (1976).

Nonlinear System Solution
From preceding developments, one may write the equilibrium equation at the instant 1 S t  using the stationary potential energy principle as: where  is the total potential energy functional, e U is the total strain energy, obtained by integrating ( 14) over the domain, F is the external forces vector, C is the viscous damping matrix and M is the mass matrix.Neumann boundary conditions are taken into account when computing F at each instant S+1, as well as the Dirichlet boundary conditions are imposed in the approximated solution as standard for finite elements so that:. where x are the nodal values for the prescribed positions at S+1 and 1 S  t are the prescribed values for traction.
From Newmark  method, ( 18) and ( 19), equation (20) becomes: where the vectors s Q and s R represent the dynamic contribution from the past and are expressed by: and The components of the gradient matrix of (22) (  ), are second derivatives of the energy functional regarding the unknowns, expressed by: From ( 25), one can use the Newton-Raphson iterative process so solve ( 22), with 0 x as the first guess for the current x , leading to the iterative process: and where l is the iteration number.The iterative process is interrupted when reached the admissible user prescribed error.

UNDAMAGED STRUCTURE ASSEMBLING WITH LAGRANGE MULTIPLIERS
In spite of the fact that the introduction of Lagrange multipliers increases the number of degrees of freedom to the problem of the continuous structure, it makes the system solver easy to be parallelized with iterative methods such as the Usawa method.
The problem is discretized with all elements unconnected, i.e., the elements do not share the same node at intersections, instead there are multiple nodes in the same position.Continuity across elements is enforced by adding subsidiary conditions with Lagrange multipliers according to: where i  is the i -th Lagrange multiplier employed to make the j -th degree of freedom equal to the k -th degree of freedom ( = i j x x ).One can observe that equation (28) represents the potential energy introduced to the system by the Lagrange multipliers.Adding (28) to the total potential energy, we have a new energy functional, written in the index notation as: Latin American Journal of Solids and Structures 14 (2017) 52-71 It is important to notice that subsidiary conditions can be used to make position or generalized vector components (for so, tangent angle and thickness strain rate) equal for coincident nodes.The subsidiary condition can be canceled based on a given rupture criteria by simply enforcing the associated Lagrange Multipliers to zero.The Lagrange multipliers in this case are the vertical and horizontal internal forces components when respectively used to link horizontal positions and to link vertical positions and bending moment when used to link tangent angle, as one can see from the first example.
Appling the same procedure described in section 2.2 to the new functional, including the Lagrange multipliers parameters, the resulting linear system to be solved at each Newton method iteration is given by: where H is the hessian matrix for the structure without the Langrange multiplier enforced constraints, given by: and the matrix L has its terms given by: As all the elements are connected by Lagrange multipliers, matrix H has the form: Where 1 H , 1 H , ..., n H are the n local Hessian matrixes of the n finite elements.
Special care is need when choosing the solver for equation ( 30), as it is a not definite system, containing zeros over the diagonal.Nowadays there are several efficient libraries for solving kind of linear systems available in any programming language.
On the other hand, one can notice that for a dynamic problem H is never singular due to the presence of the mass matrix terms.As the local matrixes are small, it is easy to evaluate its inverse.It makes easy to parallelize the code by solving the system (notice that we ommited the time step and Newton-Raphson iteration indexes from equations ( 34)-(37) for clarity): combined to an iterative procedure to enforce the subsidiary condition.
Latin American Journal of Solids and Structures 14 (2017) 52-71 The system (34) is composed by n independent small systems, and so can be easily distributed over different CPU's.One iterative procedure for the enforcement of the subsidiary condition is given by the Uzawa method, introduced by Uzawa (1958), following the equations bellow: and ( 1) where m is the Uzawa iteration number and  is a scaling factor which can be used for conver- gence acceleration.

NUMERICAL EXAMPLES
In this section, we firstly use a simple cantilever beam to validate the proposed methodology and then apply it to the study of building progressive collapse under two circumstances: 1) considering the subtle failure of a column, what may happen due to car accidents or due to explosions and 2) considering a dynamic horizontal base movement, like during earthquakes.

Cantilever Beam
We chose this example to verify the proposed method because dynamic linear analytic solution and static geometric nonlinear Bernoulli-Euler analytic solution are available to be used as a reference.
The beam is considered to be formed by 8 elements with 4 nodes each (cubic Lagrange polynomial shape functions), connected by Lagrange multipliers, according to figure 3 In order to check Lagrange multipliers values for small displacements, we apply a load 5 F  N. The Lagrange multipliers for the tangent angle are compared to the bending moment absolute value (in order to take into account signal convention) and the Lagrange multipliers for vertical positions are compared to the shear force absolute value in Fig. 4, revealing to be coincident as expected.
Latin American Journal of Solids and Structures 14 (2017) 52-71 We firstly perform a static analysis of the entire beam in order to compare the proposed numerical formulation to the analytic solution considering Bernoulli-Euler kinematics, which should be very close to the present numerical solution as it is a very slender beam.
Non-dimensional vertical and horizontal displacements ( / v L  and / u L  ) vs. non-dimensional force ( 2 /( ) PL EI ) are compared to the exact solution numerically obtained via elliptic integrals given by Mattiasson (1981) in Fig. 5.One may observe that the obtained results are coincident to the analytic solution, attesting that the proposed model is adequate for static analysis of unbroken structures.Moving forward, we consider the dynamic behavior of the cantilever considering the force F as a 50N constant impact load applied at instant = 0 t .At the segments CD and BC are released respectively at instants = 5.4 t s and = 6.0 t s by enforcing the Lagrange multipliers related to its connections to zero, according to figure 6.
Figure 7 shows the vertical displacement for points D , C and B vs. time and compares to the linear analytic solution for point D vertical displacement for the unbroken beam, considering only the first vibration mode.
Latin American Journal of Solids and Structures 14 (2017) 52-71  From Fig. 7 one can observe, as expected, that the full beam vibrates initially predominantly according to its first mode free vibration natural frequency = 0.409Hz f and after detaching the two segments, the remaining clamped part vibrates according to the first mode free vibration natural frequency for a 5m long cantilever ( = 1.63Hz f ).In figure 8 we present snapshots at each 0.25 s, from = 0 t to = 8s t .We consider a 6 pavement building structure, discretized by 84 elements with 4 nodes, according to figure 9. .The adopted resistance criteria for the structural elements are based on the limit state theory according to the Brazilian standard for steel structure ABNT NBR 8800:2008(2008).During the analysis if one given element is considered unsafe according the mentioned criteria, then the closest extremity node to the rupture point is detached from the structure.
The chosen cross section will present a design resisting bending moment = 103.1) Bending rupture: 2) Shearing rupture: 3) Normal rupture: 4) Combined bending/normal rupture: For the dynamic analysis we consider a distributed mass of 59.89 kg/m for the columns (only the structural mass considering steel density of 3 = 7860kg/m


) and 1600 kg/m for the beams (considering structural and static loading masses), i.e., we consider the loading masses to be firmly attached to the beams.The considered gravity acceleration is 2 = 9.81 m/s g .
We carried two analyses, in the first a column is considered to be suddenly broken and in the second the building is subjected to a dynamic harmonic base movement.For both cases, the initial condition is the one of static equilibrium for the weight and static loading.This condition is achieved by an initial static nonlinear analysis, resulting the internal forces of figure 10.
For both cases, no contact is considered among elements after rupture, however we consider friction-less impact with the ground.In order avoid instabilities due to flat impact, the ground is considered to be a rigid parabolic line described by the equation: The impact is enforced by introducing new lagrange multipliers when a given node crosses the ground surface according to: where u is the vector that goes from the position of the impacting node to the crossing point at ground surface and n is the ground surface normal vector.

Column Rupture
In this example we consider the second column from the left to the right to be cut in the middle in the initial instant ( = 0 s t ).Considering only static effects, the loading finds alternative ways and the building remains stable, as one can see from the static analysis in Fig. 11.However, when considering dynamic effects, the effect of mass acceleration leads to a sequence of ruptures.Figure 12 shows the internal forces for the dynamic analysis at instant = 0.025 s t and figure 13 shows snapshots of the progressive collapse.In Fig. 13 (c), considering the large displacement scale, one can observe that at t=0.095 s, most of the beam-column joints have already failed.Following the column on the left (Fig. 13 (d)) also fail and the others structural elements break when touching the ground.It is important to remember that we are considering frictionless impact with the ground, so that the elements are pushed up again.
From these results, a first solution to be tested can be strengthening the beams from the first pavement, as ruptures start propagating from them before any column rupture.

Base Movement
In order to show that this methodology can also be applied to earthquake induced structural collapse, we consider a horizontal base movement imposed to building with initial condition of static equilibrium, according to the equation: with time in seconds and displacement in meters.From Fig. 14, which shows the horizontal displacements for the center of the top along time in the first 1.2 s (before ruptures starting), one can observe that the vibration amplitude increases until the building collapses.
Again we observe a sequence of ruptures represented by the snapshots in Fig. 15 and we can conclude the methodology is also robust and efficient to deal with such problems.The beam-column joint is the first place to present rupture again and following, with less lateral support, the second and the forth columns also present failure as well as most of the remaining beams.

CONCLUSION
A 2D frame finite element methodology to deal with flexible multibody systems dynamics is presented and applied to building progressive collapse analysis.Initially we describe the large displacement position based frame dynamics formulation, which is based on a total Lagrangian description, with no displacements or rotations degrees of freedom and is suitable for simulate short bars, considering searing deformations, being at same time stable and robust.Following we develop one strategy based on Lagrange multipliers to dynamically separate the elements, allowing to simulate dynamical rupture.Finally we test the proposed methodology by numerical examples and confirm its robustness and efficiency, including multi-storey building collapse.It is important to mention that, in spite of introducing more degrees of freedom to the initial structure by the introduction Latin American Journal of Solids and Structures 14 (2017) 52-71 of Lagrange Multipliers, the technique can be combined with interactive methods and easily parallelized.

Figure 5 :
Figure 5: Displacements of the cantilever beam obtained with static analysis.

Figure 8 :
Figure 8: Snapshots at each 0.25 s for the cantilever dynamic problem.

Figure 9 :
Figure 9: Geometry and discretization for the 6 pavement building.
ering the design soliciting bending moment Sd M , design soliciting axial force Sd N and design soliciting shear force Sd V , we consider the structural element to fail if any of the conditions bellow occurs:

Figure 10 :
Figure 10: Internal forces for static analysis.

Figure 13 :
Figure 13: Progressive collapse snapshots for the column rupture problem.

Figure 14 :
Figure 14: Horizontal displacement at the top of the building.

Figure 15 :
Figure 15: Progressive collapse snapshots for the base movement problem.