European Congress on Computational Methods in Applied Sciences and Engineering
ECCOMAS 2004
P. Neittaanmäki, T. Rossi, K. Majava, and O. Pironneau (eds.)
R. Owen and M. Mikkola (assoc. eds.)
Jyväskylä, 24—28 July 2004

A FINITE ELEMENT CODE FOR STRUCTURAL COLLAPSEANALYSES OF FRAMED STRUCTURES UNDER IMPACT LOADS
Daigoro Isobe* and Kyaw Myo Lynn†*

Institute of Engineering Mechanics and Systems, University of Tsukuba
1-1-1 Tennodai Tsukuba-shi, Ibaraki 305-8573, Japane-mail: isobe@kz.tsukuba.ac.jp, web page: http://www.kz.tsukuba.ac.jp/~isobe/†

Graduate school, University of Tsukuba
1-1-1 Tennodai Tsukuba-shi, Ibaraki 305-8573, Japane-mail: kyaw@kz.tsukuba.ac.jp

Key words: ASI-Gauss Technique, ASI Technique, Structural Collapse Analysis, Framed
Structures, Impact Loads, Finite Element Method.
Abstract. In this paper, we describe a new finite element code, which can be efficiently
applied to structural collapse analyses of framed structures under impact loads. The code isdeveloped by using the ASI-Gauss technique, 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 inelastic range, by placing the numerical integration points of the two consecutive elements
forming an elastically deformed member in such a way that stresses and strains are evaluatedat the Gaussian integration points of the two-element member. Moreover, the technique canbe used to express member fracture, by shifting the numerical integration point to an
appropriate position and by releasing the resultant forces in the element simultaneously. An
impact phenomenon between structures is expressed by means of contacts between theelements involved, and the elemental contact algorithm is verified from the point ofconservation of energy. Impact analyses considering member fracture with different sets ofparameters are performed on a high-rise framed structure collided by a small aircraft. From
the results obtained, we can observe propagation of impact loads and shock waves in theframed structure. Also, the results indicate that the larger the pre-impact kinetic energy, thelarger the impact damage. In the cases of the same pre-impact kinetic energy, slightdifferences are observed when the parameters of the aircraft are varied. Moreover, soon after
impact, tensile stresses are observed in the columns that were compressed by dead loads
before impact.

1

Daigoro Isobe and Kyaw Myo Lynn.

1 INTRODUCTION
The progressive collapse of World Trade Center towers in New York, USA revealed thestructural vulnerability of tall steel buildings to the impact of a fast moving object. Technicalinvestigations are necessary to clarify the technical conditions that contributed to theprogressive collapse and to improve the ways buildings are designed, constructed, maintained
and used. Generally, numerical simulation is considered to be one of the means to investigate
such problems.
Advancements in computational technology have helped the development of various
numerical analysis methods and some of them have been applied to analyze dynamic
behaviors with strong nonlinearities and discontinuities. The Distinct Element Method (DEM)[1] follows the discontinuum approach, that is, the analytical model is formed by a number ofdistinct interacting general-shaped elements to simulate the perfectly discrete behavior. The
Discontinuous Deformation Analysis (DDA) [2] models a discontinuous material as a system
of individually deformable blocks that move independently without interpenetration. DDA is
known to be able to cope with large deformations and strong discontinuities. The Applied
Element Method (AEM) [3], which can predict crack initiation and propagation in thematerial and can also follow the total behavior from zero loading to complete collapse, hasalso been developed and applied to detailed nonlinear analysis of reinforced concrete
structures. However, the above-mentioned discrete numerical methods are computationallyintensive and need detailed modeling. Therefore, they are suitable only for detailed analysesof 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 analyses of relatively small displacements and it needs complicated modifications to apply
to strongly nonlinear and discontinuous problems. The main objective of this study is todevelop a finite element code, which can simulate dynamic behaviors with strong
nonlinearities and discontinuities and is efficiently applicable to impact collapse analyses offramed structures.
Toi and Isobe developed the Adaptively Shifted Integration (ASI) technique [4,5] for the
linear Timoshenko beam and the Euler-Bernoulli beam elements, which can be easily
implemented into the existing finite element codes. In this technique, the numericalintegration point is shifted immediately after the occurrence of a fully plastic section in theelement so that a plastic hinge is formed exactly at that section. As a result, the ASI techniquegives 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 numberof elements for a member. Structurally discontinuous problems have also become easily
handled, by shifting the numerical integration point of the linear Timoshenko beam element toan appropriate position, and by releasing the resultant forces in the element simultaneously[6,7]. 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-degree displacement
function of the linear Timoshenko beam element.

2

Daigoro Isobe and Kyaw Myo Lynn.

In this study, the ASI technique using the linear Timoshenko beam element is modified
into the ASI-Gauss technique to improve its accuracy especially in elastic range. It is handledby placing the numerical integration points of the two consecutive elements forming anelastically deformed member in such a way that stresses and strains are evaluated at theGaussian integration points of the two-element member, where the accuracy of bending
deformation is mathematically guaranteed for two-point integration. The numericalintegration points are adaptively shifted according to the change in material property (elastoplasticity) of the member. The validity of the ASI-Gauss technique is verified, in this paper,by comparing with the ASI technique and the conventional scheme. The impact phenomenon
between structures is expressed by means of contacts between the elements involved, and theelemental contact algorithm is verified from the point of conservation of energy. Generally, it
is difficult to determine the loads resulted in the structures due to an impact, and applying
impact loads to an analytical model in the form of nodal forces may not well simulate theimpact phenomenon. But by using the elemental contact algorithm, the phenomenon can besuitably expressed even when the parameters for the colliding structure are varied. Dynamic
collapse analyses considering member fracture and elemental contact are performed tosimulate the impact of a small aircraft into a ten-storied steel framed structure, andobservations are made to evaluate the influence of mass, velocity and size of the aircraft on
the impact damage.
2 METHODS
In this section, the ASI technique using the linear Timoshenko beam element issummarized. It is followed by the explanation on the newly proposed ASI-Gauss technique,and the algorithms considering member fracture and elemental contact.
2.1 Adaptively Shifted Integration (ASI) technique
Elasto-plastic analyses of framed structures using the linear Timoshenko beam element are
based on the following incremental form of the principle of virtual work:

δＶ－δＷ＝∫

L/2−L / 2{δ∆ε }T {∆σ }dz − {δ∆u}T {∆f } = 0(1)where V and W are the internal work and the external work, {∆ε}, {∆σ}, {∆u}, and {∆f} are
the generalized strain increment vector, the generalized stress (resultant force) increment
wz

L/2
L/2

1

2
θy

θx

θz

u
x

v
y
Figure 1: Definitions of local coordinates and displacements.

3

Daigoro Isobe and Kyaw Myo Lynn.

vector, the nodal displacement increment vector and the external force increment vector,
respectively. The symbols δ and ∆ denote variation and increment. The local coordinates,
nodal numbers and displacement components (u, v, w, θx, θy, θz) are shown in Fig.1.
The relation between the generalized strain increment and nodal displacement increment
vectors are given by the following equation:{∆ε (r )} = [ B( s )]{∆u}(2)

where [B(s)] is the generalized strain-nodal displacement matrix. s is the location of thenumerical integration point and r the location where stresses and strains are actually evaluated.
We refer r as the stress evaluation point later in this paper. s is a non-dimensional quantity,
which is defined as z/(L/2) where L is the length of the element. s takes a value between –1
and 1. It should be noted here that the linear Timoshenko beam element uses one-point
integration to avoid shear locking.
When the element behaves elastically, the relation between resultant force increment
vector and generalized strain increment vector is given by the following equation:{∆σ (r )} = [ De (r )]{∆ε (r )}(3)where [De] is the stress-strain matrix for elastic deformation. Substitution of Eqs. (2) and (3)into Eq. (1) leads to the following form of elemental stiffness matrix:[ K ] = L[ B( s )]T [ De (r )][ B( s )](4)

In the ASI technique, numerical integration points are shifted immediately after theformation of a fully plastic section in the element in order to form a plastic hinge exactly atthat section. When a plastic hinge is judged to be unloaded, the corresponding numericalintegration point is shifted back to its normal position. Here, the normal position means thelocation where the numerical integration point is placed when the element acts elastically.
Figure 2 shows a linear Timoshenko beam element and its physical equivalence to a rigidbody-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 follows:
r1 = − s1(5)

When the entire region of an element behaves elastically, the numerical integration point isplaced at the midpoint of the element (s1=0) which gives r1=0 from Eq. (5) in this case. The
elemental stiffness matrix for the element, the generalized strain increment vector and theresultant force increment are calculated as follows:
-1
θ1
u1

s1

0
s

Linear Timoshenko beam element

1

θ1

θ2

r1
r

u1

u2

u2
RBSM

Figure 2: Linear Timoshenko beam element and its physical equivalent.

4

θ2

Daigoro Isobe and Kyaw Myo Lynn.

[ K ] = L[ B(0)]T [ De (0)][ B(0)](6a){∆ε (0)} = [ B(0)]{∆u}{∆σ (0)} = [ De (0)]{∆ε (0)}

(6b)
(6c)

The resultant forces calculated from Eq. (6c) are physically those at the midpoint of theelement.
According to the elementary beam theory, the relations between the bending moments andshear forces are given as follows:
Vx = −

dM y

dz
dM x
Vy = −
dz

(7a)
(7b)

Here, Mx and My are the bending moments around x- and y-axes while Vx and Vy are theshear forces along x- and y- axes. Thus the distributions of the bending moment increments inan elastically deformed element are determined by the following approximations:
Ls
2
Ls∆ M y ( s ) = ∆ M y (0 ) − ∆ V x (0 )2

∆ M x ( s ) = ∆ M x (0) − ∆ V y (0 )

(8a)
(8b)

The above equations show that the bending moments are subject to a linear change in anelement and are likely to take the maximum value on either ends (s=±1). As other resultant
forces have constant values in the element, a fully-plastic state can be determined with anyield function by comparing the calculated distributions from Eqs. (8).
After the formation of a fully plastic section, the numerical integration point is shifted
immediately to a new point (s1= −r1). For instance, if a fully plastic section occurs at the rightend of an element (r1=1), the numerical integration point is shifted to the left end (s1= −1) andvice versa. In this case, the elemental stiffness matrix, the generalized strain and resultant
force increment vectors are given as follows:[ K ] = L[ B( s1 )]T [ D p (r1 )][ B( s1 )](9a){∆ε (r1 )} = [ B( s1 )]{∆u}{∆σ (r1 )} = [ D p (r1 )]{∆ε (r1 )}

(9b)
(9c)

Here, [Dp] is the stress-strain matrix for plastic deformation. It is determined by the plastic
flow theory using the following form of yield function:f y (σ (r )) − 1 = 0

5

(10)

Daigoro Isobe and Kyaw Myo Lynn.

Member
-1

1

-1/√3

rg 0 sg
Element 1

0

1/√3

-1 -1

sg 0 rg1

sg = 1 - (2/√3)
rg = -1 + (2/√3)
Numerical integration point

1

Stress evaluation point

Element 2

Figure 3: Locations of numerical integration and stress evaluation points in ASI-Gauss technique.

2.2 ASI-Gauss techniqueIn this paper, we propose the ASI-Gauss technique for elasto-plastic collapse analyses offramed structures using the linear Timoshenko beam element. As in the ASI technique, thenumerical integration point is shifted immediately after the formation of a plastic hinge and it
is shifted back to its normal position when the element is unloaded. The main difference
between the ASI and the ASI-Gauss techniques lies in their normal positions of the numericalintegration points.
The relationship between the locations of the numerical integration point and the stressevaluation point is expressed as follows:
rg = − s g(11)When the entire region of an element behaves elastically, the numerical integration point islocated at (s1=sg). As shown in Fig.3, the normal position of the numerical integration point in
the ASI-Gauss technique is sg=1−(2/√3) and the stress evaluation point is rg= −sg = −1+(2/√3).
The two stress evaluation points coincide with Gaussian integration points of the member.
This means that stresses and strains are evaluated at Gaussian integration points when theelement behaves elastically. Gaussian integration points are known to be optimal for twopoint integration and the accuracy of bending deformation is mathematically guaranteed [8].
In this case, the elemental stiffness matrix, the generalized strain and resultant force increment
vectors are given as follows:[ K ] = L[ B( s g )]T [ De (rg )][ B( s g )](12a){∆ε (rg )} = [ B( s g )]{∆u}(12b){∆σ (rg )} = [ De (rg )]{∆ε (rg )}(12c)

The distributions of the bending moment increments in an elastically deformed element are
determined by the following approximations:∆M x (s) = ∆M x (sg ) − ∆V y (s g )6L(s + sg )2(13a)Daigoro Isobe and Kyaw Myo Lynn.

∆M y (s) = ∆M y (sg ) − ∆V x (sg )L(s + sg )2(13b)

Immediately after the occurrence of a fully plastic section, the numerical integration pointis shifted to a new point (s1= −r1) as in the ASI technique. For instance, if a fully plastic
section occurs at the right end of an element (r1=1), the numerical integration point is shiftedto the left end (s1= −1) and vice versa. The elemental stiffness matrix, the generalized strain
vector and resultant force increment vector are given by Eqs. (9a), (9b) and (9c). It is to be
noted here, that in this case, the numerical integration point of the adjacent element forming
the same member is shifted back to the midpoint of the element where its accuracy isguaranteed against one-point integration. If the plastic hinge is judged to be unloaded, both
numerical integration points in the two-element member are, again, shifted to the locations
according to Fig.3. In other words, the locations of numerical integration points are adaptively
shifted according to the change in material property (elasto-plasticity) of the member.
The yield function assumed in this study is expressed by the following equation:
2

2

2

2

 M   My   N   Mz 
 + + 
 = 1f y ≡  x  +   M x0   M y 0   N 0   M z0 

(14)Here, Mx, My, N and Mz are the bending moments around x- and y- axes, the axial force andthe torsional moment, respectively. Those with the subscript 0 are the values that cause a fully
plastic section in an element when they act on the cross-section independently. The effect ofshear forces is neglected in the yield function.
The dynamic collapse analysis is made possible by adding inertial force terms to the quasistatic computational procedure. That is, the following equation for the complete system is
integrated by an appropriate time integration scheme:[ M ]{u&&}t + ∆t + [C ]{u&}t + ∆t + [ K ]{∆u} = {F }t + ∆t − {R}t

(15)

Here, [M] is the total mass matrix, [C] the total damping matrix, [K] the total stiffness matrix,{F}t+∆t the nodal external force increment vector at t+∆t, {R}t the nodal internal force
increment vector at t. {u&&}t + ∆t , {u&}t + ∆t and {∆u} are the nodal acceleration, velocity anddisplacement increment vector at t+∆t. In this study, we use consistent mass as an elemental
mass matrix which is assumed to be unaffected by the shifting of numerical integration points.
Newmark’s β method [9] is used as a numerical integration scheme for Eq. (15), andconjugate gradient (CG) method is used as a solver to reduce the need of memory resources.
2.3 Algorithm considering member fracture
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 4 shows the locations of numerical integration points for each stage in the7Daigoro Isobe and Kyaw Myo Lynn.

s=1
s=sg

s=0
Plastic hinge

s=sg

s=0

Node
(a) Elastic stage

Fractures=0

Numerical integration point(b) Plastic stage(c) Member fracture

Figure 4: ASI-Gauss technique dealing with member fracture.

ASI-Gauss technique. A member fracture is determined, in this study, by examining axial andbending strains in two elements forming a member, using the following conditions.
 κx κ  − 1 ≥ 0,  y  − 1 ≥ 0,  ε z  − 1 ≥ 0
κ 
κ 
ε 
 fx 
 fy 
 fz 

(16)

Here, κx, κy, εz, κfx, κfy and εfz are the bending strains around x- and y- axes, the axial strain
and the critical values for each strain, respectively. When fracture is judged to occur, a newnode is created at the fractured section and the numerical integration point of the element isshifted back to its midpoint. Resultant forces exerting at the fractured section are released
instantly at the next step after fracture has occurred. The nodal mass is equally redistributed ateach separated node after the step.
2.4 Algorithm considering elemental contact
Generally, it is difficult to determine the loads resulted in the structures due to an impact.
Applying impact loads to the model in the form of nodal forces may not well simulate theimpact phenomenon. Thus, the impact phenomenon is expressed, in this paper, by means ofcontact between the elements involved. Elemental contact algorithm is also useful to simulate
fractured elements colliding with other elements, and to consider the variation of parametersof the colliding structure.
Contact determination is done by examining the following two factors: (i) the distance
between the approaching element and the another one, (ii) the condition under which all four
nodal points lie on the same plane. Figure 5 shows the geometric relation in the contact
determination. Using the coordinates of the nodal points of the colliding element;
A1(xc1,yc1,zc1), A2(xc2,yc2,zc2), and another element Ei; Bi1(xi1,yi1,zi1), Bi2(xi2,yi2,zi2), thecondition under which all four nodal points lie on the same plane can be expressed as follows:8

Daigoro Isobe and Kyaw Myo Lynn.

A1

Bi1

Colliding element Ec

Element Ei

A2
Bi2
Figure 5: Geometric relation in the determination of elemental contact.

f ( x, y, z ) ≡ {( yi1 − y c 2 )( z i 2 − z c 2 ) − ( yi 2 − y c 2 )( z i1 − z c 2 )}( xc1 − xc 2 ) +{( xi 2 − xc 2 )( z i1 − z c 2 ) − ( xi1 − xc 2 )( z i 2 − z c 2 )}( y c1 − y c 2 ) +(17){( xi1 − xc 2 )( yi 2 − y c 2 ) − ( xi 2 − xc 2 )( yi1 − yc 2 )}( z c1 − z c 2 ) = 0
When two elements lie on the same plane and exist within a specific distance from eachother as calculated by the following equation, they are assumed to be in contact:
A1 Bi1 + A1 Bi 2 + A2 Bi1 + A2 Bi 2 ≤ C l (Lc + Li )(18)

Here, A1 Bi1 , A1 Bi 2 , A2 Bi1 , A2 Bi 2 are the distances between the four nodal points and Lc, Li thelengths of colliding element Ec and the element Ei, respectively. Cl is the coefficientconcerning contact length and a value of 1.8 is used throughout this paper.
For the elements that do not lie on the same plane but exist within a certain distancespecified by Eq. (18), another condition is given as expressed by the following equation,
which assumes that the elements nearly form a plane:f ( x, y , z ) ≤ C f(19)

Cf is a coefficient expressing the degree of planarity, which is fixed by considering the
Bi1
Colliding element Ec

A1
A2
Gap elements

Element Ei
Bi2

Figure 6: Four gap elements binding the two elements in contact.

9

Daigoro Isobe and Kyaw Myo Lynn.

sectional properties of constituent members. In this paper, a value of 3.0x10-4 m3 is used for Cf.
Contact determination is skipped for the elements that are not approaching each other. Whentwo elements are determined to be in contact, they are bound with a total of four gap elements
as shown in Fig.6. Gap elements are assumed to have the same geometric and material
properties as the elements in contact. In the case of contact, impact force is delivered throughthese gap elements to the connecting elements. After certain time of contact, the gap elements
are removed from the analytical model. In this paper, the binding time is set to 0.5 ms. It is
assumed that no energy loss takes place during the impact.
3 NUMERICAL EXAMPLESIn this section, the results of the elasto-plastic collapse analyses using a simple space frame,
which is loaded both statically and dynamically, are discussed. In the analyses, we use thefollowing three schemes in order to verify the accuracy of the newly proposed ASI-Gauss
technique: (a) conventional scheme in which the numerical integration point of each element
is fixed at its midpoint; (b) the ASI technique; (c) the ASI-Gauss technique. Moreover, thealgorithm considering elemental contact is verified from the point of conservation of energy
by carrying out an impact analysis of a two-beam model.
3.1 Elasto-plastic collapse analysis of a space frame
Figure 7 shows the analyzed space frame with its dimensions, geometric properties andmaterial constants. As shown in the figure, the frame is fixed at its lower ends and a static
horizontal load is applied to its upper left corner. Members are subdivided into elements
according to the numbers specified in the legends of the figures. As large deformations are notexpected, the total Lagrangian formulation (TLF) is used in the analyses.
Figure 8 shows the results of the elasto-plastic analyses under static loading. As for theconventional scheme, convergence is extremely slow. The ASI technique shows better
convergence, since the numerical integration points are shifted to form plastic hinges atappropriate locations. However, it still shows slow convergence in the elastic range, due to the
low-degree displacement function of the element. On the other hand, the ASI-Gauss techniqueshows that just two-element modeling is needed to obtain the converged solution in the elastoplastic analysis.
L = 2.00 m, A = 2.50 × 10 −3 m 2

L

I x = I y = 5.21 × 10 − 7 m 4
LP,δρ = 7900 kg / m 3 , E = 206 GPaυ = 0.3, σ y = 245 MPa
M x 0 ≡ M y 0 = 7.66 × 10 3 Nm

LZ

N 0 = 6.13 × 10 5 NYX

M z 0 = 1.12 × 10 5 Nm
Figure 7: Analyzed space frame model.

10Daigoro Isobe and Kyaw Myo Lynn.

100

70
P,δ

60

Displacement δ (mm)8050Load P (kN)

2 elements/member
4 elements/member
8 elements/member
16 elements/member
32 elements/member
64 elements/member

P,δ

40
30
2 elements/member
4 elements/member
8 elements/member
16 elements/member
32 elements/member
64 elements/member

20
10
0

60

40

20

0.1

100
200
Displacement δ (mm)

0.2

0.3

100

70
P,δ

2 elements/member
4 elements/member
8 elements/member
16 elements/member

P,δ

60

Displacement δ (mm)8050Load P (kN)

0.5(a) Conventional scheme(a) Conventional scheme

40
30
20

2 elements/member
4 elements/member
8 elements/member
16 elements/member

10
0

60

40

20

0.1

100
200
Displacement δ (mm)

0.2

0.3

0.4

0.5

Time (s)(b) ASI technique(b) ASI technique
100

70

2 elements/member
4 elements/member
8 elements/member

P,δ

P,δ

60

Displacement δ (mm)8050Load P (kN)

0.4Time (s)4030
20
2 elements/member
4 elements/member
8 elements/member

10
0

60

40

20

0.1

100
200
Displacement δ (mm)

0.2

0.3

0.4

0.5

Time (s)(c) ASI-Gauss technique(c) ASI-Gauss techniqueFigure 9: Elasto-plastic response analysis.

Figure 8: Elasto-plastic collapse analysis.

11Daigoro Isobe and Kyaw Myo Lynn.

3.2 Elasto-plastic response analysis of a space frame
In the analyses subjected under dynamic loading, the time increment ∆t is set to 0.05 ms
and the same model shown in Fig.7 is used. A step load of 12 kN is applied for a total of 100
steps (5 ms). Damping effect is not considered in the analysis.
Figure 9 shows the results of the elasto-plastic response analyses under dynamic loading.
As in the analyses described in the previous section, the conventional scheme shows very
slow convergence and at least sixteen-element modeling is necessary to obtain the converged
solution in the analysis. The ASI technique gives comparatively better results than theconventional scheme. However, when it comes to the model that uses only two elements per
member, poor convergence is seen and the vibration mode is much different from that of theconverged solution. On the other hand, the ASI-Gauss technique shows very fast convergencein the result and nearly converged solutions are obtained even when the number of elements
per member is two.
From the above analyses, it is confirmed that the ASI-Gauss technique is an efficientmethod for analyses of frame structures both in elastic and plastic range, and even the twoelement modeling is sufficient to acquire the results with high accuracy and reliability.
3.3 Impact analysis of a simple two-beam model
In this section, impact phenomenon of two beams is simulated to verify the creditability ofthe elemental contact algorithm. The analyzed model, its geometrical properties and material
constants are shown in Fig.10. One of the beams is fixed at its both ends while the other one
moves freely. Initial velocity of 50 m/s is applied to the free beam. The beams are modeledwith two elements per each member. The ASI-Gauss technique is applied, and the updated
Lagrangian formulation (ULF) is used since large deformations are expected in the analysis.
Damping effect is not considered in the analysis. The time increment ∆t is set to 10 µs. The
calculation is done for a total of 3000 steps (0.03 s).

L/2

L = 1.00 m , v = 50 m / sv

A = 2.50 × 10 − 3 m 2

L/2

I x = I y = 5.21 × 10 − 7 m 4v
L/2

ρ = 7900 kg / m 3 , E = 206 GPav

υ = 0.3, σ y = 245 MPa

Fixed ends

M x 0 ≡ M y 0 = 7.66 × 10 3 Nm

Z
L/2

Y

X

N 0 = 6.13 × 10 5 N, M z 0 = 1.12 × 10 5 Nm

Figure 10: Analyzed two-beam model.

12Daigoro Isobe and Kyaw Myo Lynn.

KE+PE+SE
KE+PE
SE

Energy (×103J)

40

30
20

10

0
0

5

10
Time (ms)20

15

Figure 11: Kinetic, potential and strain energies in the two-beam model.

Figure 11 shows the sum of kinetic and potential energies, strain energy and the total of all
three for the whole model. Energies for both the colliding and the fixed beams constituted bya total of four elements are calculated by using the following equations:

(

l
1 4
2
2
ρA i v i1 + v i 2∑
2 i =124l
PE = ∑ ρA i g(h i1 + h i 2 )
2
i =1
4
1
SE = ∑ {σ i }{ε i }l i
2 i =1

KE =

)

(20a)
(20b)
(20c)

Here, KE, PE and SE stand for kinetic energy, potential energy and strain energy, respectively.
{σi},{εi}, hi, li represent the resultant force vector, the generalized strain vector, z coordinates
and the length of the element Ei (i=1,2,3,4). Subscripts 1 and 2 denote the nodal points of thecorresponding element.
As seen from the figure, kinetic and potential energies occupy all of the total energy before
the impact takes place. After the impact, both beams deflect and vibrate as seen in Fig.12. As
a consequence, the strain energy is stored up in the beams and the sum of the other two

Z
Y

Z
X

0.000sY

0.010s

Z

Z
X

Y

X

Y

0.018s

0.014sZ
X

Y

0.022sFigure 12: Impact analysis of a simple two-beam model.

13

Z
X

Y

X

0.026sDaigoro Isobe and Kyaw Myo Lynn.

Body 1

Body 2

Unit

L

16

-

m

Lc

-

15

m

h

30

-

m

hc

-

2

m

ha

-

L/4
2Lc/5

L/4
2Lc/5

Lc

L/4
L/4

2Lc/5

L/4

L/4

L/4
Y

L

Lc/5
2Lc/5

3Lc/10

X

h/10
h/10
h/10

hc
Lc

h/10
h/10
h/10

ha

h/10
h/10
h/10

Z

21.50
-2

1.00x10

m
-2

m2

A

2.19x10Ix

6.66x10-4

8.33x10-6m4

Iy

2.24x10-4

8.33x10-6

m4

ρ

7900

2800

kgm-3

E

206

69

Gpa

σy

240

195

Mpa

ν

0.30

0.335

Mxo

7.99x105

4.88x104

Nm

Myo2.69 x105

4.88x104Nm

No

5.25x106

1.95x106N

Mzo

1.12x105

1.12x105

Nm

Body 1: Ten-storied steel framed structure
Body 2: Aircraft (made of duralumin)

X

h/10

Figure 13: Analyzed ten-storied steel framed structure and the aircraft.

decreases. The main reason for the fluctuations seen after the impact is considered to be theoccurrence of longitudinal and bending waves. However, the total energy remains almost
constant, and the elemental contact algorithm used in this study is considered to be
sufficiently reliable.
4 STRUCTURAL COLLAPSE ANALYSIS OF FRAMED STRUCTURES UNDERIMPACT LOADSIn this section, dynamic collapse analyses are performed using the ASI-Gauss technique tosimulate the impact of a small aircraft into a ten-storied steel framed structure. A total of fourcases are analyzed using different values for three parameters: cruise speed, size and weight

14

Daigoro Isobe and Kyaw Myo Lynn.

of the aircraft. Details of the analyzed model, propagation of impact loads, response of thestructures and the extent of impact damage are discussed in this section.
4.1 Outline of the analyzed model
Figure 13 shows the analyzed ten-storied framed structure and the aircraft. Both structures
are modeled by using two elements per each member. The model contains 1698 elements and7308 degrees of freedom. The gap elements, which bind the elements in contact for a certain
period of time, are assumed to have the same geometrical and material properties as theelements forming the aircraft. The critical values for bending strain and axial strain are set toκfx=κfy=0.04 and εfz=0.3 in these analyses.
The weight of the floor system is considered and nearly twenty times of the own weight ofthe skeletal structure is applied as dead loads. The ten-storied framed structure is fixed at its
lower ends while none of the element is fixed in the aircraft model. It is assumed that the total
weight of the aircraft is uniformly distributed in its elements. Initial velocity is applied to theaircraft with zero degree of inclination from the horizon.
4.2 Parameters set in the analyses
A total of four analyses are performed using different values for cruise speed, size andweight of the aircraft. The combinations of these parameters are shown in Table 1. The
parameters are selected to fix the pre-impact kinetic energy values for Cases I and IV, Cases
II and III to be the same. In all the four analyses, the ASI-Gauss technique with the linear
Timoshenko beam element is used and member fracture and elemental contact are bothconsidered. Damping effect is not considered in all analyses. The updated Lagrangian
formulation (ULF) is used to handle large deformations of the elements. The Conjugate
Gradient (CG) method is used as the solver to reduce the need for memory resources. For
Newmark’s β method [9], 1/6 to 1/4 is the commonly used range for β. However, it is known
that it becomes less precise but more stable with a value of β between 1/4 and 1/2. Sincestrongly nonlinear and discontinuous behaviors of the problem in this study cause instability
in calculations, we used a value of 0.4 for β. The time increment ∆t is set to 10 µs. The
calculation is done for a total of 35000 steps (0.35 s) for each case. The calculation for eachcase takes approximately 20 CPU hours with a personal computer (2.4 GHz CPU, 256 MBRAM).
4.3 Structural collapse analysis of the ten-storied framed structure under impact loads
Figure 14 shows the results of the first analysis (Case I). A total of six colors are used toTable 1: Parameters of the aircraft

Case No.
Case I
Case II
Case III
Case IV

Cruise speed [kmph]
440
311
440
440

Size Lc [m]
15
15
15
12

Weight [kg]
60000
60000
30000
60000

15

Pre-impact kinetic energy [J]448 x106
224 x106
224 x106448 x106

Daigoro Isobe and Kyaw Myo Lynn.

0ms

3ms

50ms

100ms

150ms

200ms

250ms

300ms

350ms

fy

0.00

0.25

0.50

0.75

1.00

fracture

Figure 14: Impact analysis of a ten-storied framed structure collided by a small aircraft (Case 1).

16

Daigoro Isobe and Kyaw Myo Lynn.

(a) Case I
Cruise speed=440kmph, Lc =15m, Weight=60000kg
Pre-impact kinetic energy = 448x106J
Number of fractured members in framed structure=10(b) Case II
Cruise speed=311kmph, Lc =15m, Weight=60000kg
Pre-impact kinetic energy = 224 x106J
Number of fractured members in framed structure =2(c) Case III
Cruise speed=440kmph,Lc =15m, Weight=30000kg
Pre-impact kinetic energy = 224x106J
Number of fractured members in framed structure =5(d) Case IV
Cruise speed=440kmph,Lc =12m, Weight=60000kg
Pre-impact kinetic energy = 448x106J
Number of fractured members in framed structure =16

Figure 15: Impact damage to the ten-storied framed structure (only the upper five floors are shown andfractured elements are excluded).

demonstrate the tendency of the elements to yield or fracture. The numbers in the color palette
indicate the values of fy in Eq. (14). Red color (fy =1.0) means a plastic hinge is formed within
the element and black color is used for fractured elements.
Propagation of impact loads can be confirmed by observing the spread of colors in the
figure. Immediately after the impact, red color propagates into the structures from the point of17

Daigoro Isobe and Kyaw Myo Lynn.

impact. The elements near the area of impact turn to black, that is, they start to fracture.
Moreover, it is seen that the colors travel down the columns and move backwards very fast.
This phenomenon is seen mainly in the columns and considered as the propagation of shockwaves. Red-colored elements are temporarily found even in the lowest floor, and this meansthat even the lowest floor is affected, to some extent, by the impact.
Figures 15(a) to 15(d) show the extent of damage to the steel framed structure for all fourcases. In order to show the damaged areas clearly, only the upper five floors are plotted andthe fractured elements are excluded. If we focus on Case I (Fig.15(a)) and II (Fig.15(b)), theimpact damage is seen larger in Case I in which the cruise speed of the aircraft is larger (andso is the pre-impact kinetic energy), and most parts of the aircraft penetrate the framed
structure. Although the pre-impact kinetic energy of Case II and III (Fig.15(c)) are the same,the results show a slight difference between the damage patterns. If we compare the results
between Case I and IV (Fig.15(d)), where the only difference between the two cases are thewidth of the aircraft, we can see that the shorter the width, the deeper the penetration. In CaseIV, most of the inner columns are cut off, however, the total deformation of the framed
structure is observed slightly smaller than in Case I. This result can be interpreted in such a
way that the size of the aircraft has a strong influence on the damage patterns. Moreover, soonafter impact in Case I, tensile stresses are observed in the columns that had been compressedby gravity loads before impact as shown in Fig.16. There is a possibility that these tensilestresses lead to some degree of structural damage especially in the joints of beams andcolumns. Therefore, further detailed analyses are necessary to investigate this phenomenon.
5 CONCLUSION
The formerly developed ASI technique for the linear Timoshenko beam element provides
higher computational efficiency than the conventional finite element scheme, and is known tobe able to cope with dynamic behaviors with strong nonlinearities including such phenomena
as member fracture. However, it still lacked accuracy in elastic range when the number of
elements per member is small. In this study, the ASI-Gauss technique, which deals elastic
members by locating the stress evaluation points of two consecutive elements at Gaussian
points, is proposed. Elasto-plastic analyses under both static and dynamic loading are
performed to verify the proposed technique, and it is confirmed that nearly converged
solutions are obtained even when the number of elements per member is two. It gives resultswith high accuracy at very low calculation cost, which may act efficiently in application to
analyses of large-scale framed structures.
Instead of applying impact loads in the form of nodal forces, impact phenomenon issimulated by means of contacts between the elements involved. This elemental contact
algorithm is also useful to simulate fractured elements colliding with other elements, and to
consider the variation of parameters of the colliding structure. The elemental contact
algorithm is verified from the point of conservation of energy and the result showed its highreliability.
The ASI-Gauss technique with the elemental contact and the member fracture algorithms
are applied to dynamic analyses with strong nonlinearities and discontinuities, where impact

18

Daigoro Isobe and Kyaw Myo Lynn.

Before impact

On impact

25ms after impact

50ms after impact

75ms after impact

100ms after impact

Elements in which tensile stresses are observed
Fractured elements

Figure 16: Tensile stresses observed in the ten-storied framed structure after impact (Case I).

of an aircraft against a ten-storied framed structure is discussed using four different sets ofparameters for the aircraft. From the results, we observed the propagation phenomenon of
impact loads and shock waves. Also, a proper difference in impact damage is obtained bydifferent sets of parameters, and in practically short calculation time. Therefore, the numericalanalysis code used in this study is considered to be an efficient code for impact analyses oflarge-scale framed structures. However, further study is considered necessary to evaluate thestructural damage caused by the tensile stresses of the columns, and damping matrices as wellas the wall and floor slabs should also be considered in the future to improve the accuracy.

19

Daigoro Isobe and Kyaw Myo Lynn.

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.
[3] H. Tagel-Din and K. Meguro. Nonlinear Simulation of RC Structures Using Applied
Element Method. Structural Engineering / Earthquake Engineering, 17(2), 137-148,
2000.
[4] Y. Toi and D. Isobe. Adaptively Shifted Integration Technique for Finite Element
Collapse Analysis of Framed Structures. Int. J. Numer. Methods Eng., 36, 2323-2339,
1993.
[5] 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.
[6] 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.
[7] 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.
[8] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery. Numerical Recipes in
FORTRAN: The Art of Scientific Computing, 2nd ed. Cambridge University Press, 1992.
[9] N.M. Newmark. A Method of Computation for Structural Dynamics. A.S.C.E. J.
Engineering Mechanics, 85, 67-94, 1959.

20