thThe 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China

DEVELOPMENT OF MACRO-MODEL SEISMIC COLLAPSE SIMULATORFOR FRAMED STRUCTURES USING ASI-GAUSS TECHNIQUE
123

N. Katahira , D. Isobe , T. Ine and K. Kajiwara

4

1

Graduate Student, University of Tsukuba, Japan
Associate Professor, Dept. of Engineering Mechanics and Energy, University of Tsukuba, Japan
3
Invited Research Fellow, National Research Institute for Earth Science and Disaster Prevention, Japan
4
Senior Researcher, National Research Institute for Earth Science and Disaster Prevention, Japan
Email: e0311273@edu.esys.tsukuba.ac.jp, isobe@kz.tsukuba.ac.jp, ine@bosai.go.jp, kaji@bosai.go.jp
2

ABSTRACT :The main objective of this study is to develop a macro-model finite element code, which can simulate dynamic
behaviors with strong nonlinearities and discontinuities and is efficiently applicable to seismic collapse analysesof framed structures. The applied technique is called the ASI-Gauss technique for the linear Timoshenko beam
element, which computes highly accurate elasto-plastic solutions even with the minimum number of elements
per member. It gains high accuracy especially in elastic range, by placing the numerical integration points of thetwo consecutive elements forming an elastically deformed member in such a way that stresses and strains areevaluated at the Gaussian integration points of the two-element member. Moreover, the technique can be used toexpress member fracture, by shifting the numerical integration 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 in this paper to verify the proposed code. Moreover, a three-dimensional seismic collapse analysis
considering member fracture and contact is performed on a large-scale framed structure, and practical results are
obtained in a reasonable short calculation time.
KEYWORDS: Seismic collapse simulator, Framed structures, ASI-Gauss technique, Linear Timoshenko beam
element, Non-linear FEM
1. INTRODUCTION
In the conventional design of a building, only static analysis in the horizontal and uniaxial directions iscommonly carried out, in order to minimize calculation costs. This approach can ensure the structural strengthof 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 layer in dynamic analysis, and the complicated dynamic behaviorof the structure at the member level becomes difficult to be examined sufficiently. Therefore, the developmentof a more precise and more 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 applicableto 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 and seismic collapse analyses [3-6]. However, the above-mentioned discrete
numerical methods are computationally intensive and need detailed modeling. Therefore, they are suitable only
for detailed 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, hasbeen successfully applied to a wide range of engineering problems including structural analyses of large-scale
structures. However, the FEM is generally limited to analyses of 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 a finite element seismic collapse simulator, which can simulate dynamic behaviors withsuch strong nonlinearities and discontinuities and is efficiently applicable to seismic collapse problems offramed structures.

th

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
Toi and Isobe developed the Adaptively Shifted Integration (ASI) technique [7,8] 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 point is shifted immediately after the occurrence of a fully
plastic section in the element so that a plastic 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 tosimulate dynamic behaviors with strong nonlinearities by using only a small number of elements for a member.
Structurally discontinuous problems have also become easily handled, by shifting the numerical integration
point of the linear Timoshenko beam element to an appropriate position, and by releasing the resultant forces inthe element simultaneously [9,10]. However, when the number of elements per member is very small, it stilllacks 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 [11]. The accuracy in elastic range is improved by placing the numerical integration points of the twoconsecutive elements forming an elastically deformed member in such a way that stresses and strains areevaluated at the Gaussian integration points of the two-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 a finite element code in order to develop a more
precise and less-calculation-time-consuming seismic collapse simulator. The purpose of this study is to verify
the validity of the ASI-Gauss technique in seismic collapse problems and to construct a highly efficient
structural design tool for framed structures under three-dimensional excitation. Simple numerical tests arecarried out to show the validity of the simulator, and an analysis involving member fracture and contact isperformed on a large-scale framed structure to compare the effects of three-dimensional excitation against
two-dimensional excitation.
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 equation of motion
for excitation at fixed points, are described.

Figure 1 Linear Timoshenko beam element and its physical equivalent2.1. ASI-Gauss TechniqueFigure 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 andthe stress evaluation point where a plastic hinge is actually formed is expressed as [12]r = Äs(2.1)

th

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
In the above equation, s is the location of the numerical integration point and r the location where stresses andstrains are actually evaluated. We refer to r as the stress evaluation point later in this paper. s and r are
nondimensional quantities, which take values 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 theelement acts elastically. By doing so, the plastic behavior of the element is simulated appropriately, and theconverged 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 behaveselastically. When the number of elements per member is very small, solutions in the elastic range are notaccurate since one-point integration is used to evaluate the low-order displacement function of the beam
element.
The main difference between the ASI and ASI-Gauss techniques lies in the normal position of the numericalintegration point. In the ASI-Gauss technique, two consecutive elements forming a member are considered as asubset, and the numerical integration points of an elastically deformed member are placed such that the stressevaluation points coincide with the Gaussian integration points of the member. This means that stresses andstrains are evaluated at the Gaussian integration points of elastically deformed members. Gaussian integration
points are optimal for two-point integration and the accuracy of bending deformation is mathematicallyguaranteed [13]. 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, generalized strain and resultant force increment
vectors for the ASI technique are given by Eqn. 2.2 and those for the ASI-Gauss technique by Eqn. 2.3.
[K] = L[B(0)]T [D(0)][B(0)](2.2a)

fÅ"(0)g = [B(0)]fÅug

(2.2b)

fÅõ(0)g = [D(0)]fÅ"(0)g(2.2c)[K] = L[B(sg )]T [D(rg )][B(sg )](2.3a)

fÅ"(rg )g = [B(sg )]fÅug

(2.3b)

fÅõ(rg )g = [D(rg )]fÅ"(rg )g(2.3c)

22
sg = 1 Ä p ; rg = Ä1 + p
33(2.3d)

thThe 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
Here, fÅ"g , fÅõg and fÅug 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 length of the element.
The plastic potential used in this study is expressed by
f=

í

Mx
Mx0

ì2

+

í

My
My0

ì2

+

í

N
N0

ì2

Ä 1 = fy Ä 1 = 0(2.4)

Here, fy is the yield function, and Mx , My and N are the bending moments around the x- and y- axes andaxial force, 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 torsion and shear force is neglected in the yieldfunction.
2.2. Member Fracture and Contact Algorithm

Figure 3 Locations of numerical integration points in each stage
A plastic hinge is likely to occur before it develops to a member fracture, and the plastic hinge is expressed byshifting the numerical integration point to the opposite end of the fully-plastic section. Accordingly, thenumerical 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 foreach stage in the ASI-Gauss technique.
Member fracture is determined, in this study, by means of curvatures, shear strains and axial tensile strain
occurred in the elements as shown in the following equation.
å
å
å
å
å
å
å
å
å îx å
å
å
å
å
å
å
å
å Ä 1 ï 0 or å îy å Ä 1 ï 0 or å çxz å Ä 1 ï 0 or å çyz å Ä 1 ï 0 or èz Ä 1 ï 0
åîx0 å
åîy0 å
åçxz0 ååçyz0 åèz0

(2.5)

where κx, κy, γxz, γyz, εz, κx0, κy0, γxz0, γyz0 and εz0 are the curvatures around the x- and y-axes, the shear strains
for the x- and y-axes, the axial tensile strain and the critical values for these strains, respectively. We fixed thecritical values using the information from actual experimental data [14]. Contact determination is done byexamining the geometrical locations of the elements and once two elements are determined to be in contact, they
are bound with a total of four gap elements between the nodes [10,11].
2.3. Incremental Equation of Motion
The dynamic equilibrium equation at time step t=t can be formulated as[M ]f°
ugt = fEgt Ä fF gt

(2.6)

th

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
ugt , fEgt and fF g t are the mass matrix, acceleration vector, nodal external force vector andwhere [M], f°internal force vector at time step t=t, respectively.
The following equation is substituted into Eqn. 2.6 at t=t+1 in the implicit code:fF gt+1 = fF gt + [K ]fÅug

(2.7)

Then the following incremental stiffness equation is evaluated:[M ]f u°g t+1 + [K ]fÅug = fE g t+1 Ä fF g t(2.8)

where [K] is a stiffness matrix at time step t=t. By neglecting residual forces, an implicit code is obtained byevaluating the following incremental equation of motion:[M ]fÅ u°g + [K]fÅug = 0(2.9)

Consequently, the incremental equation of motion for a structure under excitation at fixed points, which is usedin this paper, yields to:[M1 ]fÅ u°g + [M2 ]fÅ u°b g + [K1 ]fÅug + [K2 ]fÅub g = 0(2.10)The subscript '1' indicates the coupled terms between free nodal points, '2' indicates the coupled terms between
ug andfree nodal and fixed nodal points, and 'b' indicates the components at fixed nodal points. Vectors fÅ°fÅug 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 fÅus g and dynamic displacement increments fÅud g , the displacements at free nodal
points are given asfÅug = fÅus g + fÅud g(2.11)

fÅus g is evaluated, by neglecting inertia force, as follows:fÅus g = Ä[K1 ]Ä 1 [K2 ]fÅub g(2.12)

Substituting Eqn. 2.11 and Eqn. 2.12 into Eqn. 2.10, the following equation is obtained:[M1 ]fÅ u°d g + [K1 ]fÅud g = ([M1 ][K1 ]Ä 1 [K2 ] Ä [M2 ])fÅ u°b g(2.13)In this scheme, equivalent forces are calculated by substituting nodal acceleration increments at fixed points into
the right side of the above equation. Damping matrices are not considered in this study.
This derivation process above is made under the condition of neglecting residual forces, which consequently
neglects the influence of gravity in collapse phenomena. The numerical analyses in previous work [10] are all
performed by including the residual forces in vertical directions to consider the gravitational effects, under the
restriction that the seismic wave in vertical direction can not be used as an input simultaneously. To enable theconsideration of vertical seismic wave and gravitational effects at the same time, the calculations of theequations of motion under a seismic wave and dead load are separated in this paper; the equation of motion of astructure subjected under dead load is calculated according to Eqn. 2.8 and the obtained residual forces areadded to the external forces calculated by Eqn. 2.13. The residual forces occurred by dead load are only
considered while those occurred by seismic excitation are neglected. By using this approach, a simultaneousconsideration of vertical seismic wave and gravitational effects became possible. All calculations in this paper
are performed using this new approach.

th

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
3. NUMERICAL ESTIMATIONS
To verify the numerical code using the ASI-Gauss technique described in the previous section, elasto-plastic
response analyses of a simple space frame under seismic excitation are performed. In the analyses, we use thefollowing three schemes 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. An EW component of JMA Kobe seismic wave is multiplied by 1.5 for theamplitude and subjected at the fixed points according to Eqn. 2.10 to Eqn. 2.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.
150

100

50

0

-50

0

10

20

Time (s)

30

40

150
2elem
4elem
8elem

100

Displacement (mm)

2elem
4elem
8elem
16elem32elem64elem

Displacement (mm)Displacement (mm)

150

50

0

-50

0

10

20

Time (s)

30

40

2elem
4elem
8elem

100

50

0

-50

0

10

20

Time (s)3040(a) Conventional scheme(b) ASI technique(c) ASI-Gauss techniqueFigure 4 Comparison of the numerical accuracies between three schemes
Figure 4 shows the results of the analyses. Since plastic hinges are expressed only at the midpoints of theelements, the conventional scheme shows a very slow convergence and it finally converges 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 withaccuracy in the elastic range. On the other hand, the ASI-Gauss technique shows a very fast convergence and aconverged solution is obtained with the minimum subdivision, since the stress evaluation points in the elements
are adaptively controlled in both elastic and plastic ranges. The three-dimensional excitation gave the sametended results.
Table 1 Member list of typical columns
Floor
8
7
6
5
4
3
2
1
Material□－450×450×12□－500×500×12□－500×500×16□－550×550×16
BCP325

Table 2 Member list of typical beams

Y

X

Figure 5 Eight-story frame model [15]

Floor
R
8
7
6
5
4
3
2
Material

H－400×200×9×16
H－600×200×12×19

H－600×250×12×19
SN490Bth

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
4. SEISMIC COLLAPSE ANALYSIS OF A LARGE-SCALE FRAMED STRUCTURE

600

Acceleration (gal)JMA-Kobe(EW)2000JMA-Kobe(UD)JMA-Kobe(NS)400

Acceleration (gal)400

Acceleration (gal)

400

800

600

200
0

-200

-200

0

-200

-400

-400

200

-600

-600
0

10

20

Time (s)(a) EW

30

40

-800

-400
0

10

20

Time (s)

30

40

0

10

(b) NS
Figure 6 JMA-Kobe seismic wave (EW, NS, UD components)

20

Time (s)3040(c) UD

A seismic collapse analysis is carried out on a 30 m high, 8-story steel framed structure [15], andthree-dimensional excitation effect in collapse phenomena is verified by comparing with two-dimensional
excitation result. The general view of the model is shown in Fig. 5 and the member lists of the typical columns
and beams are shown in Table 1 and 2. The sectional sizes of the columns at the 3rd floor are deliberately set to
70% of the values in the table, in order to induce member fracture. EW, NS and UD components of JMA Kobe
seismic wave shown in Fig. 6 are multiplied by 2.0 for the amplitude, and are subjected orthographically (EW to
X-direction, NS to Y-direction and UD to Z-direction) at the fixed points on the ground floor. Dead load foreach floor is assumed as 8.0 kN/m2. The time increment is set to 5 ms and the calculation is done for a total of8000 steps. The calculation takes approximately 80 CPU minutes with a personal computer (Core2 DUO E6300,
1.86 GHz CPU Ç 2, 3 GB RAM).

7.20 s30.64 s32.08 sFigure 7(a) Seismic collapse analysis under two-dimensional seismic wave

40.00 s7.20 s30.64 s32.08 s40.00 sFigure 7(b) Seismic collapse analysis under three-dimensional seismic wave

Figures 7(a) and 7(b) show the result under two-dimensional (EW and NS) excitation and the result under
three-dimensional excitation, respectively. Fractured elements are excluded from the figures. Many plastic
sections appear and some beams fracture in the two-dimensional excitation case, however, the structurewithstands the excitation and does not collapse. But in the three-dimensional excitation case, the gravitational
effect and the small vibration after the peak pulse gradually enlarges the deformation on the 3rd floor, where theth

The 14 World Conference on Earthquake Engineering
October 12-17, 2008, Beijing, China
sectional sizes of the columns are deliberately lessened, and finally, failure initiates at approximately 30 s. Thismeans that the UD component of the seismic wave cannot be neglected if the amplitude of the seismic wave
particularly in vertical direction is considerably large.
5. CONCLUSIONSIn this paper, a nonlinear finite element seismic collapse simulator using the ASI-Gauss technique wasdeveloped, in order to analyze seismic collapse problems including structural discontinuities. The proposed
simulator was improved by separating the calculation of the equations of motion for vertical seismic wave anddead load, to enable the consideration of both effects at the same time. The simulator gives accurate andpractically realistic results in a short calculation time, which may be valid when used in a numerical estimationof the seismic design of large-scale framed structures. However, the damping matrices, as well as more detailed
member properties such as local buckling, should also be considered in the future studies to improve thepracticability.
REFERENCES[1] Cundall, P.A. (1971). A Computer Model for simulating Progressive, Large-scale Movement in Blocky
Rock System, Proceedings of the International Symposium on Rock Mechanics II-8, 129-136.
[2] Shi, G.H. and Goodman, R.E. (1984). Discontinuous Deformation Analysis, Proceedings of 25th U.S.
Symposium on Rock Mechanics, 269-277.
[3] Meguro, K. and Hakuno, M. (1991). Simulation of Collapse Process of Structures due to Earthquake,
Proceedings of Symposium on Computational Methods in Structural Engineering and Related Fields, 15,
325-330, in Japanese.
[4] Ma, M.Y., Barbeau, P. and Penumadu, D. (1995). Evaluation of Active Thrust on Retaining Walls using
DDA, Journal of Computing in Civil Engineering, 1, 820-827.
[5] Yarimer, E. (1989). Demolition by Controlled Explosion as a Dynamical Process, Structures under Shockand Impact, 411-416.
[6] Tosaka, N., Kasai Y. and Honma, T. (1988). Computer Simulation for Felling Patterns of Building,
Demolition Methods and Practice, 395-403.
[7] Toi, Y. and Isobe, D. (1993). Adaptively Shifted Integration Technique for Finite Element Collapse Analysis
of Framed Structures, International Journal for Numerical Methods in Engineering, 36, 2323-2339.
[8] Toi, Y. and Isobe, D. (1996). 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.
[9] Isobe, D. and Toi, Y. (2000). Analysis of Structurally Discontinuous Reinforced Concrete Building Frames
using the ASI Technique, Computers and Structures, 76:4, 471-481.
[10] Isobe, D. and Tsuda, M. (2003). Seismic Collapse Analysis of Reinforced Concrete Framed Structures
Using the Finite Element Method, Earthquake Engineering and Structural Dynamics, 32:13, 2027-2046.
[11] Isobe, D. and Lynn, K.M. (2004). A Finite Element Code for Structural Collapse Analyses of Framed
Structures under Impact Loads, Proceedings of the 4th European Congress on Computational Methods in
Applied Sciences and Engineering (ECCOMAS 2004), Jyvaskyla, Finland.
[12] Toi, Y. (1991). Shifted Integration Technique in One-Dimensional Plastic Collapse Analysis using Linear
and Cubic Finite Elements, International Journal for Numerical Methods in Engineering, 31, 1537-1552.
[13] Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (1992). Numerical Recipes in
FORTRAN: The Art of Scientific Computing, 2nd Edition, Cambridge University Press.
[14] Hirashima, T., Hamada, N., Ozaki, F., Ave, T. and Uesugi, H. (2007). Experimental Study on Shear
Deformation Behavior of High Strength Bolts at Elevated Temperature, Journal of Structural and Construction
Engineering, 621, 175-180, in Japanese.
[15] Kaitani, J. et al. (2006). Design Example and Earthquake Response of the Frames with Knee Damper Part
2: Design Examples of 8-Story Office Building, Summaries of Technical Papers of Annual Meeting
Architectural Institute of Japan 2006, C-1, 753-754, in Japanese.