ECCOMAS Thematic Conference on Computational Methods
in Structural Dynamics and Earthquake Engineering
M. Papadrakakis, D.C. Charmpis, N.D. Lagaros, Y. Tsompanakis (eds.)
Rethymno, Crete, Greece, 13-16 June 2007

A DYNAMIC FINITE ELEMENT CODE FOR ANALYZING COLLAPSEBEHAVIORS OF FRAMED STRUCTURES UNDER SEISMIC LOADS
D. Isobe1 and N. Katahira2
1 Dept. of Eng. Mech. and Energy, Univ. of Tsukuba

1-1-1 Tennodai Tsukuba-shi, Ibaraki 305-8573, Japane-mail: isobe@kz.tsukuba.ac.jp
2 College of Eng. Systems, Univ. of Tsukuba

1-1-1 Tennodai Tsukuba-shi, Ibaraki 305-8573, Japane-mail: e0311273@edu.esys.tsukuba.ac.jp

Keywords: Framed Structures, Seismic Load, Collapse Behavior, ASI-Gauss Technique, Finite
Element Method.
Abstract. A new dynamic finite element code to simulate collapse behaviors of structures
under excitation at fixed points is developed using the ASI-Gauss technique. The technique
is a modified version of the formerly developed Adaptively Shifted Integration (ASI) technique
for the linear Timoshenko beam element, which computes highly accurate elasto-plastic solutions even with the minimum number of elements per member. The ASI-Gauss technique gains
still higher accuracy especially in elastic range, by placing the numerical integration pointsof the two consecutive elements forming an elastically deformed member in such a way thatstresses and strains are evaluated at the Gaussian integration points of the two-element member. Moreover, the technique can be used to express member fracture, by shifting the numericalintegration point to an appropriate position and by releasing the resultant forces in the element
simultaneously. The numerical code can be used to analyze dynamic behaviors of framed structures, initiating from elastic range to a total collapse. Simple numerical tests are conducted toverify the proposed code. Moreover, a seismic collapse analysis considering member fracture
and contact is performed on a large-scale framed structure, and practical results are obtainedin a reasonable short calculation time.

1

D. Isobe and N. Katahira

1

INTRODUCTION

In the conventional design of a building, only static analysis in the horizontal and uniaxial
directions is commonly carried out, in order to minimize calculation costs. This approach canensure the structural strength of the building, if there is sufficient strength to support the load in
the vertical direction. For a similar reason, the mass system model replaces the building layerin dynamic analysis, and the complicated dynamic behavior of the structure at the member level
becomes difficult to be examined sufficiently. Therefore, the development of a more precise andmore efficient dynamic analysis code is strongly desired.
Recently, significant advances in the field of computers have been eliminating the calculation
cost restrictions, and various dynamic analysis codes are being developed. Among those codes,
there are some codes applicable to dynamic collapse problems which contain strong nonlinearities and discontinuities, such as the Distinct Element Method (DEM) [1] or the Discontinuous
Deformation Analysis (DDA) [2]. These codes have been applied to demolition analyses andseismic collapse analyses [3-8]. However, the above-mentioned discrete numerical methods arecomputationally intensive and need detailed modeling. Therefore, they are suitable only fordetailed analyses of either two-dimensional or small three-dimensional models.
On the other hand, the Finite Element Method (FEM), which is based on continuum material
equations, has been successfully applied to a wide range of engineering problems including
structural analyses of large-scale structures. However, the FEM is generally limited to analysesof relatively small displacements and it needs complicated modifications to simulate fracture
occurring in structural members or joints. The main objective of this study is to develop afinite element code, which can simulate dynamic behaviors with such strong nonlinearities anddiscontinuities and is efficiently applicable to seismic collapse analyses of framed structures.
Toi and Isobe developed the Adaptively Shifted Integration (ASI) technique [9,10] for the
linear Timoshenko beam and the Bernoulli-Euler beam elements, which can be easily implemented into the existing finite element codes. In this technique, the numerical integration pointis shifted immediately after the occurrence of a fully plastic section in the element so that aplastic hinge is formed exactly at that section. As a result, the ASI technique gives more precise elasto-plastic solutions than the conventional schemes, and has become able to simulate
dynamic behaviors with strong nonlinearities by using only a small number of elements for amember. Structurally discontinuous problems have also become easily handled, by shifting thenumerical integration point of the linear Timoshenko beam element to an appropriate position,and by releasing the resultant forces in the element simultaneously [11,12]. However, when
the number of elements per member is very small, it still lacks accuracy in the elastic range
compared to the converged solution, due to the low-order displacement function of the linear
Timoshenko beam element.
To improve the accuracy in elastic range, Isobe and Lynn modified the ASI technique into an
ASI-Gauss technique [13]. The accuracy in elastic range is improved by placing the numericalintegration points of the two consecutive elements forming an elastically deformed member
in such a way that stresses and strains are evaluated at the Gaussian integration points of thetwo-element member, where the accuracy of bending deformation is mathematically guaranteed
for two-point integration. The rest of the shiftings are done according to the flow of the ASI
technique.
In this study, the ASI-Gauss technique is implemented into the finite element code in orderto develop a more precise and less-calculation-time-consuming seismic collapse analytical tool.
The purpose of this study is to verify the validity of the ASI-Gauss technique in seismic collapse2

D. Isobe and N. Katahiraanalysis and to construct a highly efficient structural design tool for framed structures. Simple
numerical tests are carried out to show the validity of the code, and an analysis involving member fracture and contact is performed on a large-scale framed structure to show the practicabilityof the scheme.
2

NUMERICAL METHOD

General concept of the ASI-Gauss technique in comparison with the ASI technique is explained in this section. Moreover, algorithms considering member fracture and elemental contact, and incremental kinematic equation for excitation at fixed points, are described.
2.1

ASI-Gauss techniqueFigure 1: Linear Timoshenko beam element and its physical equivalent.

Figure 1 shows a linear Timoshenko beam element and its physical equivalence to the rigid
bodies-spring model (RBSM). As shown in the figure, the relationship between the locations of
the numerical integration point and the stress evaluation point where a plastic hinge is actuallyformed is expressed as [14]r = −s(1)
In the above equation, s is the location of the numerical integration point and r the location
where stresses and strains are actually evaluated. We refer to r as the stress evaluation point
later in this paper. s is a nondimensional quantity, which takes a value between -1 and 1.
In both the ASI and ASI-Gauss techniques, the numerical integration point is shifted adaptively when a fully plastic section is formed within an element to form a plastic hinge exactly
at that section. When the plastic hinge is determined to be unloaded, the corresponding numerical integration point is shifted back to its normal position. Here, the normal position means
the location where the numerical integration point is placed when the element acts elastically.
By doing so, the plastic behavior of the element is simulated appropriately, and the converged
solution is achieved with only a small number of elements per member. However, in the ASI
technique, the numerical integration point is placed at the midpoint of the linear Timoshenko
beam element, which is considered to be optimal for one-point integration, when the entire region of the element behaves elastically. When the number of elements per member is very small,solutions in the elastic range are not accurate since one-point integration is used to evaluate thelow-order displacement function of the beam element.
The main difference between the ASI and ASI-Gauss techniques lies in the normal positionof the numerical integration point. In the ASI-Gauss technique, two consecutive elements form3

D. Isobe and N. Katahira

ing a member are considered as a subset, and the numerical integration points of an elastically
deformed member are placed such that the stress evaluation points coincide with the Gaussian integration points of the member. This means that stresses and strains are evaluated at theGaussian integration points of elastically deformed members. Gaussian integration points are
optimal for two-point integration and the accuracy of bending deformation is mathematicallyguaranteed [15]. In this way, the ASI-Gauss technique takes advantage of two-point integration
while using one-point integration in actual calculations.

Figure 2: Locations of numerical integration and stress evaluation points in elastic range.

Figure 2 shows the locations of the numerical integrations points of elastically deformed
elements in the ASI and ASI-Gauss techniques. The elemental stiffness matrix, generalizedstrain and resultant force increment vectors for the ASI technique are given by Eqs. (2) andthose for the ASI-Gauss technique by Eqs. (3).
[K] = L[B(0)]T [D(0)][B(0)](2a){∆ε(0)} = [B(0)]{∆u}(2b){∆σ(0)} = [D(0)]{∆ε(0)}(2c)[K] = L[B(sg )]T [D(rg )][B(sg )](3a){∆ε(rg )} = [B(sg )]{∆u}(3b){∆σ(rg )} = [D(rg )]{∆ε(rg )}

(3c)

2
2
sg = 1 − √ , rg = −1 + √
3
3

(3d)

Here,{∆ε},{∆σ} and {∆u} are the generalized strain increment vector, generalized stress (resultant force) increment vector and nodal displacement increment vector, respectively. [B] is
the generalized strain-nodal displacement matrix, [D] the stress-strain matrix and L is the lengthof the element.
The plastic potential used in this study is expressed by
µ

Mx
f=
Mx0

¶2

Ã

My
+
My0

!2

µ

N
+
N0

¶2

4

µ

Mz
+
Mz0

¶2

− 1 = fy − 1 = 0(4)

D. Isobe and N. Katahira

Here, fy is the yield function, and Mx , My , N and Mz are the bending moments around thex- and y- axes, axial force and torsional moment, respectively. Those with the subscript 0
are values that result in a fully plastic section in an element when they act on a cross section
independently. The effect of shear force is neglected in the yield function.
2.2

Member fracture and contact algorithm

A plastic hinge is likely to occur before it develops to a member fracture, and the plastic
hinge is expressed by shifting the numerical integration point to the opposite end of the fullyplastic section. Accordingly, the numerical integration point of the adjacent element forming
the same member is shifted back to its midpoint where it is appropriate for one-point integration.
Figure 3 shows the locations of numerical integration points for each stage in the ASI-Gauss
technique.

Figure 3: Locations of numerical integration points in each stage.

L’

L

r

θ

D

L'
− 1 = ε Bx or ε By
L

Figure 4: Definition of bending strain.

Member fracture is determined, in this study, by means of bending strain and axial strain
occurred in the elements. By assuming that a pin-supported element is deformed as shown in
Fig. 4, the initial length L and the maximum deformed length L0 of the element are expressed
as
θ
(5a)
L = 2r sin
25

D. Isobe and N. Katahira

D
)θ
(5b)
2
where D is the member width, r and θ are the radius and angle of the bending deformation,
respectively. Bending strain εB is defined using Eqs. (5) as
L0 = (r +

εB =

Dθ
θ
L0
−1=
+
−1
L
2L 2 sin(θ/2)(6)

Here, two angles of bending deformation around the x- and y- axes can be written asθx = |κx |L(7a)θy = |κy |L(7b)where κx and κx are the curvatures around the x- and y- axes, respectively. Substituting Eqs.
(7) into Eq. (6) yields to the following equations that calculate bending strains around the xand y- axes.
D|κx |LεBx = |κx | +
−1
(8a)
2
2 sin(|κx |L/2)
εBy =|κy |LD|κy | +
−1
2
2 sin(|κy |L/2)(8b)

Member fracture is determined when the sum of bending strains εBx , εBy and axial strain εzsatisfies the following condition.
εz0 ≤ εBx + εBy + εz

(9)Here, εz0 is the ultimate tensile strain of the member material.
Contact determination is done by examining the geometrical locations of the elements andonce two elements are determined to be in contact, they are bound with a total of four gap
elements between the nodes.
2.3

Incremental kinematic equation

The incremental kinematic equation for a structure under excitation at fixed points, which isused in this paper, is as follows:[M1 ]{∆ü} + [M2 ]{∆üb } + [K1 ]{∆u} + [K2 ]{∆ub } = 0(10)

The subscript ’1’ indicates the coupled terms between free nodal points, ’2’ indicates the coupled terms between free nodal and fixed nodal points, and ’b’ indicates the components at fixed
nodal points. [M ] is the mass matrix, and vectors {∆ü} and {∆u} are the nodal acceleration
increment and the nodal displacement increment, respectively.
Under the assumption that the displacements at free nodal points are estimated by adding
quasi-static displacement increments {∆us } and dynamic displacement increments {∆ud }, thedisplacements at free nodal points are given as{∆u} = {∆us } + {∆ud }(11){∆us } is evaluated, by neglecting inertia force, as follows:{∆us } = −[K1 ]−1 [K2 ]{∆ub }6(12)

D. Isobe and N. Katahira

Substituting Eqs. (11) and (12) into Eq. (10), the following equation is obtained:[M1 ]{∆üd } + [K1 ]{∆ud } = ([M1 ][K1 ]−1 [K2 ] − [M2 ]){∆üb }(13)

In this scheme, equivalent forces are calculated by substituting nodal acceleration increments atfixed points into the right side of the above equation. Damping matrices are not considered inthis study.
3

NUMERICAL ESTIMATIONS

To verify the numerical code using the ASI-Gauss technique, simple elasto-plastic response
analyses under seismic excitation are performed. In the analyses, we use the following threeschemes to verify the accuracy of the ASI-Gauss technique: (a) a conventional finite element
scheme in which the numerical integration point of each element is fixed at its midpoint; (b) the
ASI technique; and (c) the ASI-Gauss technique. Figure 5 shows the analyzed space frame withits dimensions, physical properties and material constants. An EW component of JMA Kobe
seismic wave is multiplied by 1.5 for the amplitude as shown in Fig. 6, and is subjected at the
fixed points according to Eqs. (10) to (13). An implicit solution scheme using Newmark’s βmethod (β = 4/9 and δ = 5/6 is used to consider numerical damping), updated Lagrangian
formulation, and a consistent mass matrix are adopted in the analyses.

Figure 6: Input seismic wave (EW, JMA-Kobe).

Figure 5: Analyzed model.

Figure 7: Comparison of the numerical accuracies between three schemes.

Figure 7 shows the results of the analyses. Since plastic hinges are expressed only at the
midpoints of the elements, the conventional scheme shows a very slow convergence and it
still does not converge even when the 64-element modeling is used. The ASI technique gives
comparatively better results than the conventional scheme. However, the two-element modeling
yet does not give converged solution, since the elements lack with accuracy in the elastic range.
7

D. Isobe and N. Katahira

On the other hand, the ASI-Gauss technique shows a very fast convergence and a converged
solution is obtained with the minimum subdivision, since the stress evaluation points in theelements are adaptively controlled in both elastic and plastic ranges. The three-dimensional
excitation gave the same tended results.
4

SEISMIC COLLAPSE ANALYSIS OF A LARGE-SCALE FRAMED STRUCTURE

Figure 8: Analyzed model.

A seismic collapse analysis is carried out on a 22-story steel framed structure as shown in
Fig. 8. Member list of the structure is shown in Table 1. We assumed a structure with lower
stiffness and longer natural period to see the effect of a large seismic wave in such condition.
The natural periods of the structure are 13 s to the X-direction, and 20 s to the Y-direction. EW
and NS components of JMA Kobe seismic wave (Fig. 6 and Fig. 9), also multiplied by 1.5
for the amplitude, are subjected orthographically (EW to X-direction, NS to Y-direction) at the
fixed points on the ground floor. Dead load for each floor is assumed as 200 kg/m2 , and a value
8

D. Isobe and N. Katahiraof 0.22 is used for εz0 in Eq. (9). The time increment is set to 5 ms and the calculation is done
for a total of 8000 steps. The calculation takes approximately 30 CPU minutes with a personal
computer (Core2 DUO E6300, 1.86 GHz CPU×2, 3 GB RAM).
MemberColumn (1F-11F)
Column (11F-22F)
BeamSpec. (material)□-250*250*9 (SN490)□-200*200*6 (SN490)
H-200*200*8*12 (SN490)

Table 1: Member list.

Figure 9: Input seismic wave (NS, JMA-Kobe).

The numerical result is shown in Fig. 10. A total of six colors are used to demonstrate thetendency of the elements to yielding. The numbers in the color palette indicate the values of fy
in Eq. (4). Red color (fy = 1.0) means a plastic hinge is formed within the element. Fractured
elements are excluded from the figure. Although the maximum amplitude of the input seismic
wave is observed between 5 s to 12 s, only the ground level columns are yielded during theperiod. However, a so-called incremental collapse continues to proceed in the columns due tosmall excitations after the period, and finally, the structure begins to collapse totally after 26s when excitations are decreased to an appearingly small amount. In this case, a difference
between the natural period of the structure and the seismic wave led to a stress concentration
on ground level columns, which reduced the total vertical resistant force of the structure itself.
The collapse process can be confirmed in detail by observing the story-shear and drift angle
relations of one of the columns on the ground floor (Fig. 11) and the drift angle time histories(Fig. 12). By observing the drift angle values, it is clear that the EW component (X-direction)of the seismic wave played a main role in the total collapse behavior of the structure. It may also
be said that the member fracture and contact algorithms have succeeded to express a realistic
phenomenon of the total collapse.
5

CONCLUSIONSIn this paper, a nonlinear finite element code using the ASI-Gauss technique was developed,
in order to analyze seismic collapse problems including structural discontinuities. The fractureof a section was modeled by shifting the numerical integration point with simultaneous releaseof the resultant forces. The proposed code was improved by considering the contact between
members in order to obtain results that agree more closely with the actual behavior. The code
gives accurate and practically realistic results in a short calculation time, which may be valid
9

D. Isobe and N. Katahira

Figure 10: Seismic collapse analysis of a large-scale framed structure.

10

D. Isobe and N. Katahira

Figure 11: Story-shear and drift angle relations on the ground floor.

Figure 12: Time histories of drift angles.

when used in a numerical estimation of the seismic design of large-scale framed structures.
However, the damping matrices, as well as more detailed member properties, should also be
considered in the future studies to improve the accuracy.
6

ACKNOWLEDGEMENT

This research is done partially as a trustee project of National Research Institute for Earth
Science and Disaster Prevention (NIED), Japan.
REFERENCES[1] P.A. Cundall, A Computer Model for simulating Progressive, Large-scale Movement inBlocky Rock System, Proceedings of the International Symposium on Rock Mechanics,
II-8, 129-136, 1971.
[2] G.H. Shi and R.E. Goodman, Discontinuous Deformation Analysis, Proceedings of 25th
U.S. Symposium on Rock Mechanics, 269-277, 1984.
11

D. Isobe and N. Katahira[3] K. Meguro and M. Hakuno, Simulation of Collapse Process of Structures due to Earthquake, Proceedings of Symposium on Computational Methods in Structural Engineeringand Related Fields, 15, 325-330, 1991. [in Japanese][4] M. Itoh, N. Yoshida, N. Utagawa and I. Kondo, Simulation of Blast Demolition of Reinforced Concrete Buildings, Proceedings of the Third World Congress on Computational
Mechanics, 1152-1153, Chiba, Japan, August 1-5, 1994.
[5] M.Y. Ma, P. Barbeau and D. Penumadu, Evaluation of active thrust on retaining walls
using DDA, Journal of Computing in Civil Engineering, 1, 820-827, 1995.
[6] I. Kondo, N. Utagawa, M. Ito and N. Yoshida, Numerical Method to Simulate Collapse
Behavior in Blasting Demolition of Space Framed Structure, Proceedings of the 13th Symposium on Computer Technology of Information, Systems and Applications, 49-54, 1990.
[in Japanese][7] E. Yarimer, Demolition by Controlled Explosion as a Dynamical Process, Structures under
Shock and Impact, 411-416, 1989.
[8] N. Tosaka, Y. Kasai and T. Honma, Computer Simulation for Felling Patterns of Building,
Demolition Methods and Practice, 395-403, 1988.
[9] Y. Toi and D. Isobe, Adaptively Shifted Integration Technique for Finite Element Collapse
Analysis of Framed Structures, International Journal for Numerical Methods in Engineering, 36, 2323-2339, 1993.
[10] Y. Toi and D. Isobe, Finite Element Analysis of Quasi-static and Dynamic Collapse Behaviors of Framed Structures by the Adaptively Shifted Integration Technique, Computers
and Structures, 58(5), 947-955, 1996.
[11] D. Isobe and Y. Toi, Analysis of Structurally Discontinuous Reinforced Concrete Building
Frames Using the ASI Technique, Computers and Structures, 76(4), 471-481, 2000.
[12] D. Isobe and M. Tsuda, Seismic Collapse Analysis of Reinforced Concrete Framed Structures Using the Finite Element Method, Earthquake Engineering and Structural Dynamics, 32(13), 2027-2046, 2003.
[13] D. Isobe and K.M. Lynn, A Finite Element Code for Structural Collapse Analyses ofFramed Structures Under Impact Loads, Proceedings of the 4th European Congress
on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2004),
Jyväskylä, Finland, July 24-28, 2004.
[14] Y. Toi, Shifted Integration Technique in One-Dimensional Plastic Collapse Analysis Using Linear and Cubic Finite Elements, International Journal for Numerical Methods inEngineering, 31, 1537-1552, 1991.
[15] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in
FORTRAN: The Art of Scientific Computing, 2nd Edition. Cambridge University Press,
1992.

12