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

Abstract

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 algorithm, 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.

Keywords:
Multibody dynamics; positional finite element method; geometric nonlinear analysis; Lagrange multipliers; progressive collapse

1 INTRODUCTION

One very impressive effect that structural engineers are frequently concerned about is the progressive collapse of structures. Progressive collapse is characterized by the failure of the total or great amount of the structure due to an initial localized problem, like the rupture or instability of one single or few structural members, which shoots the trigger for failure propagation. Many researchers addressed this subject along the years, starting in the 1970s, as one can mention important works of (Ellingwood and Leyendecker, 1978) and (Leyendecker and Ellingwood, 1977) that contributed to starting a discussion on design against progressive collapse. However, the precise consideration of dynamic effects during load redistribution seems to be a more recent subject, explored by authors like (Kaewkulchai and Williamson, 2004), (Hakuno and Meguro, 1993) and (Botez et al., 2015).

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 and Coda, 2013 and 2014).

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 in (Ibrahimbegović and Mamouri, 2000 and 2002), (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.

2 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.

Figure 1
Frame element mapping.

The mappings fm0 and fm1, from the auxiliary non-dimensional configuration (I) respectively to the initial (II) and final configurations (III), have its components written as follows:

(1)

and

(2)

Where X^m _ji and x^m _ji 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 (g⁰) is normal to the neutral surface, and is not necessarily normal to the current neutral surface (g¹), allowing distortion due to shear stresses. The generalized vector is unconstrained, and so, allows deformation in the height direction (see Fig. 2).

Figure 2
Deformation kinematics.

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:

(3)

and

(4)

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 g⁰ and g¹ components can be written as:

(5)

and

(6)

where e ⁰ _ij are the j-th nodal values for the i component of the unity vector normal to the initial configuration middle surface, e ¹ _ij are the j-th nodal values for the i component of the generalized vector (not necessarily normal or unity) at current configuration, h ₀ is the initial cross section height and a_j 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 e⁰ and e¹ can be written as:

(7)

and

(8)

where θ⁰ and θ¹ 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:

(9)

where X is the vector of nodal values at initial configuration.

Considering the gradient matrixes A⁰ and A¹ of initial and current mappings (3) and (4), with it components given by:

(10)

and

(11)

where the index “, j” indicates derivatives regarding ξ₁ if j = 1, or ξ₂ if j = 2. The deformation gradient matrix A is expressed as:

(12)

After evaluating the gradient A, the Green strain tensor is given by:

(13)

where the variables C_ij 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:

(14)

Such constitutive law relates second Piola-Kirchhoff stress tensor and Green strain tensor according to:

(15)

where the elastic tensor D is given by:

(16)

with G and v 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):

(17)

The finite element adopted is an isoparametric curved line with arbitrary number of nodes. The shape functions are Lagrange polynomials of degree n-1 where n is the number of nodes of the considered element. Each node has 5 nodal parameters: 2 middle surface position vector components x_i with i = 1 or 2, the current angle to obtain the non-normal generalized position vector component e1, the current thickness 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.

2.1 Time Integration

As mentioned before in section 1, we use the Newmark β algorithm for numerical time integration. The time marching process is summarized as:

(18)

and

(19)

where x_s 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).

2.2 Nonlinear System Solution

From preceding developments, one may write the equilibrium equation at the instant t_{S +1} using the stationary potential energy principle as:

(20)

where П is the total potential energy functional, U_e 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:.

(21)

where _{S +1} are the nodal values for the prescribed positions at S+1 and

_S+1 are the prescribed values for traction.

From Newmark β method, (18) and (19), equation (20) becomes:

(22)

where the vectors Q_s and R_s represent the dynamic contribution from the past and are expressed by:

(23)

and

(24)

The components of the gradient matrix of (22) (∇ℱ), are second derivatives of the energy functional regarding the unknowns, expressed by:

(25)

From (25), one can use the Newton-Raphson iterative process so solve (22), with x0_s +1 = xs as the first guess for the current x_s +1, leading to the iterative process:

(26)

and

(27)

where l is the iteration number. The iterative process is interrupted when reached the admissible user prescribed error.

3 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:

(28)

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 (x_i = x_j ).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:

(29)

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:

(30)

where H is the hessian matrix for the structure without the Langrange multiplier enforced constraints, given by:

(31)

and the matrix L has its terms given by:

(32)

As all the elements are connected by Lagrange multipliers, matrix H has the form:

(33)

Where H₁, H₁, ..., H_n 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):

(34)

combined to an iterative procedure to enforce the subsidiary condition.

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:

(35)

with

(36)

and

(37)

where m is the Uzawa iteration number and ρ is a scaling factor which can be used for convergence acceleration.

4 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.

4.1 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. The beam has Young modulus E = 210 GPa, density ρ = 7860 kg/cm³, cross section area A = 23.91 cm², cross section moment of inertia I = 47.619 cm⁴ and length. L = 10m

Figure 3
Cantilever beam.

In order to check Lagrange multipliers values for small displacements, we apply a load F = 5 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.

Figure 4
Internal forces vs. Lagrange multipliers.

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 (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.

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

Moving forward, we consider the dynamic behavior of the cantilever considering the force F as a 50N constant impact load applied at instant t = 0. At the segments CD and BC are released respectively at instants t = 5.4 s and t = 6.0 s by enforcing the Lagrange multipliers related to its connections to zero, according to figure 6.

Figure 6
Cantilever segments detachment.

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.

Figure 7
Vertical displacement vs. time.

From Fig. 7 one can observe, as expected, that the full beam vibrates initially predominantly according to its first mode free vibration natural frequency f = 0.409Kz 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 (f = 1.63Kz). In figure 8 we present snapshots at each 0.25 s, from t = 0 to t = 8s.

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

4.1 Six Storey Building Progressive Collapse

We consider a 6 pavement building structure, discretized by 84 elements with 4 nodes, according to figure 9.

Figure 9
Geometry and discretization for the 6 pavement building.

The entire structure is made of AISC W 410x60 (metric system) I-shaped steel bars with Young modulus E = 205GPa, density ρ = 7860kg/m³, yield stress fy = 250MPa.The geometric properties are moment of inertia I = 2.17∙10^-4m⁴, cross section area A = 7.62∙10^-3m².

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 M_Rd = 103.18kNm, design resisting shear force V_Rd = 374.85kN and a design axial resisting force N_Rd = 607.9kN. Considering the design soliciting bending moment M_Sd , design soliciting axial force N_Sd and design soliciting shear force V_Sd , we consider the structural element to fail if any of the conditions bellow occurs:

1. Bending rupture:

(38)

1. Shearing rupture:

(39)

1. Normal rupture:

(40)

1. Combined bending/normal rupture:

(41)

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 ρ = 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 g = 9.81m?s².

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:

(42)

The impact is enforced by introducing new lagrange multipliers when a given node crosses the ground surface according to:

(43)

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.

Figure 10
Internal forces for static analysis.

4.2.1 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 (t = 0s). 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 11
Static internal forces redistribution.

Figure 12 shows the internal forces for the dynamic analysis at instant t = 0.025s 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.

Figure 12
Internal 0.025 s after column rupture.

Figure 13
Progressive collapse snapshots for the column rupture problem.

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.

4.2.2 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:

(44)

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.

Figure 14
Horizontal displacement at the top of the building.

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.

Figure 15
Progressive collapse snapshots for the base movement problem.

5 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 of Lagrange Multipliers, the technique can be combined with interactive methods and easily parallelized.

REFERENCES

• NBR 8800:2008 (ABNT) (2008) Norma técnica - Projeto de estruturas de aço e de estruturas mistas de aço e concreto de edifícios.
• Betsch, P. and Steinmann, P. (2003) Constrained dynamics of geometrically exact beams. Computational Mechanics, 31(1):49-59.
• Bonet, J., Wood, R. D., Mahaney, J., Heywood, P. (2000) Finite element analysis of air supported membrane structures. Computer Methods in Applied Mechanics and Engineering, 190(5-7):579 - 595.
• Botez, M., Bredean, L., Ioani, A. M. (2015) Improving the accuracy of progressive collapse risk assessment: Efficiency and contribution of supplementary progressive collapse resisting mechanisms. Computers & Structures.
• Carrazedo, R., and Coda, H. B. (2010) Alternative positional FEM applied to thermomechanical impact of truss structures. Finite Elements in Analysis and Design, 46(11):1008 - 1016.
• Coda, H.B. and Greco, M. (2004) A simple FEM formulation for large deflection 2d frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering, 193(33-35):3541 - 3557.
• Coda, H.B. and Paccola, R.R. (2011) A FEM procedure based on positions and unconstrained vectors applied to non-linear dynamic of 3d frames. Finite Elements in Analysis and Design, 47(4):319 - 333.
• Coda, H.B. and Paccola, R.R. (2014) A total-Lagrangian position-based FEM appied to physical and geometrical nonlinear dynamics of plane frames including semi-rigid connections and progressive collapse, Finite Elements in Analysis and Design, 91, 1-15, 2014.
• Ellingwood, B.R. and Leyendecker, E.V. (1978) Approaches for design against progressive collapse. Journal of the Structural Division/ASCE, 104(ST3):413 - 423.
• Ghosh, S. and Roy, D. (2009) A frame-invariant scheme for the geometrically exact beam using rotation vector parametrization. Computational Mechanics, 44(1):103-118.
• Hakuno, M. and Meguro, K. (1993) Simulation of concrete-frame collapse due to dynamic loading. Journal of Engineering Mechanics, 119(9):1709 - 1723.
• Hughes, T.J.T. (1976) Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics. Computers & Structures, 6(4):313 - 324.
• Ibrahimbegović, A. and Mamouri, S. (2000) On rigid components and joint constraints in nonlinear dynamics of flexible multibody systems employing 3d geometrically exact beam model. Computer Methods in Applied Mechanics and Engineering, 188(4):805 - 831 {IVth} World Congress on Computational Mechanics. (II). Optimum.
• Ibrahimbegović, A. and Mamouri, S. (2002) Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations. Computer Methods in Applied Mechanics and Engineering, 191(37-38):4241 - 4258.
• Ibrahimbegović, A., Taylor, R.L., Lim, H. (2003) Non-linear dynamics of flexible multibody systems. Computers & Structures, 81(12):1113 - 1132.
• Jelenić, G. and Crisfield, M.A.(2001) Dynamic analysis of 3d beams with joints in presence of large rotations. Computer Methods in Applied Mechanics and Engineering, 190(32-33):4195 - 4230.
• Kaewkulchai, G. and Williamson, E.B. (2004) Beam element formulation and solution procedure for dynamic progressive collapse analysis. Computers & Structures, 82(7-8):639 - 651.
• Leyendecker, E.V. and Ellingwood, B.R. (1977) Design Methods for Reducing the Risk of Progressive Collapse in Buildings: NBS Building Science Series. 98. U.S. Department of Commerce / National Bureau of Standards, Washington.
• Mattiasson, K. (1981) Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals. International Journal for Numerical Methods in Engineering, 17(1):145-153.
• Ogden, R.W. (1984) Non-linear elastic deformations. Engineering Analysis, 1(2):119 - 119.
• Reis, M.C.J. and Coda, H.B. (2014) Physical and Geometrical non-linear analysis of plane frames considering elastoplastic semi-rigid connections by the positional FEM, Latin American Journal of Solid and Structures, 11 (7), 1163-1189.
• Romero, I. (2008) Formulation and performance of variational integrators for rotating bodies. Computational Mechanics, 42(6):825-836.
• Romero, I. and Armero, F.(2002) An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy-momentum conserving scheme in dynamics. International Journal for Numerical Methods in Engineering, 54(12):1683-1716.
• Sanches, R.A.K. and Coda, H.B. (2013) Unconstrained vector nonlinear dynamic shell formulation applied to fluid structure interaction. Computer Methods in Applied Mechanics and Engineering, 259:177 - 196.
• Sanches, R.A.K. and Coda, H.B. (2014) On fluid-shell coupling using an arbitrary Lagrangian-Eulerian fluid solver coupled to a positional lagrangian shell solver. Applied Mathematical Modelling, 38(14):3401 - 3418.
• Schiehlen, W., Guse, N., Seifried, R. (2006) Multibody dynamics in computational mechanics and engineering applications. Computer Methods in Applied Mechanics and Engineering, 195(41-43):5509 - 5522. John H. Argyris Memorial Issue. Part {II}.
• Ahmed A. Shabana, A.A. (2013) Dynamics of multibody systems. Cambridge University Press, Cambridge, UK, 4th edition.
• Uzawa, H. (1958) Iterative methods for concave programming. In K. J. Arrow, L. Hurwicz, and H. Uzawa, editors, Studies in linear and nonlinear programming. Stanford University Press.

Publication Dates

• Publication in this collection
  Jan 2017

History

• Received
  25 May 2016
• Reviewed
  29 July 2016
• Accepted
  01 Nov 2016