Structural Collapse Analysis of Framed Structures under Impact Loads
Using ASI-Gauss Finite Element Method
Kyaw Myo Lynn1* and Daigoro ISOBE2
1*

Graduate School, University of Tsukuba, 1-1-1 Tennodai, Tsukuba-shi, Ibaraki 305-8573, Japan(kyaw@kz.tsukuba.ac.jp)
2
Institute of Engineering Mechanics and Systems, University of Tsukuba, 1-1-1 Tennodai, Tsukuba-shi, Ibaraki
305-8573, Japan (isobe@kz.tsukuba.ac.jp)
Abstract
The main objective of this study is to devise a technique, which, when implemented into finite element
codes, is efficiently applicable to impact collapse analyses of framed structures. In this study, the formerly
developed Adaptively Shifted Integration (ASI) technique for the linear Timoshenko beam element is modified
into the ASI-Gauss technique 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 evaluated at the Gaussian
integration points of the two-element member. On comparison with the ASI technique, the ASI-Gauss techniqueproves its higher accuracy and efficiency in elastic range. Instead of applying impact loads in the form of nodal
forces, this study simulates the impact phenomenon by means of contacts between the elements involved and theelemental contact algorithm is verified from the point of conservation of energy. Impact analyses considering
member fracture with different sets of parameters are performed using a high-rise framed structure and a small
aircraft. From the results obtained, we can observe propagation phenomena of impact loads and shock waves.
Different parameters of the aircraft produce proper difference in impact damage. Results also indicate that the
influence of mass of the aircraft has a stronger influence on impact damage than its velocity. Moreover, soonafter impact, tensile stresses are observed in the columns that were compressed by dead loads before impact.
Keywords: ASI-Gauss technique; ASI technique; Impact; Elemental contact; Member fracture; Finite element method

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 the progressive
collapse and to improve the ways buildings are designed, constructed, maintained and used. And
numerical simulation is considered to be one of the means to investigate such problems.
Advancements in computing technology have helped the development of various numericalanalysis methods and some of them have been applied to analyze dynamic behaviors with strong
nonlinearities and discontinuities. The Distinct Element Method (DEM) [1] follows thediscontinuum approach, that is, the analytical model is formed by a number of distinct 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 tocope with large deformations and strong discontinuities. The Applied Element Method (AEM),
which can predict crack initiation and propagation in the material [3] and can also follow the total
behavior from zero loading to complete collapse, has also been developed and applied to detailed
nonlinear analysis of reinforced concrete structures. 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, has been successfully applied to a wide range of engineering problems including
structural analyses of large-scale structures. However, FEM is generally limited to analyses ofrelatively small displacements and it needs complicated modifications to apply to nonlinear anddiscontinuous problems. The main objective of this study is to devise a technique, which, when
implemented into finite element codes, can simulate dynamic behaviors with strong nonlinearities
and discontinuities and is efficiently applicable to impact collapse analyses of framed structures.
Toi and Isobe developed the Adaptively Shifted Integration (ASI) technique [4][5][6] for the
linear Timoshenko beam element, which can be easily implemented into the existing finite element
codes. In this technique, the numerical integration point is shifted immediately after the occurrenceof a fully plastic section in the element so that a plastic hinge is formed exactly at that section. As aresult, the ASI technique gives a more precise solution than the conventional finite element method,and has become able to simulate dynamic behaviors with strong nonlinearities by a small numberof elements. However, when the number of elements per member is very small, it still lacks
accuracy in the elastic range compared to the converged solution. If the model to be analyzed is alarge structure, it is desirable, from the point of calculation costs, to use a small number of elements
per member. So it is very important to have a numerical scheme that ensures accuracy even when
the number of elements per member is small. In this study, the ASI technique is modified into the
ASI-Gauss technique 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 evaluated at theGaussian integration points of the two-element member. And the accuracy of the ASI-Gauss
technique is verified by comparing with the ASI technique and the conventional FEM.
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. This study tries to simulate the impact phenomenon by means of contacts
between the elements involved. The elemental contact algorithm [7][8] is verified from the point ofconservation of energy. Furthermore, dynamic collapse analyses considering member fracture are
performed to simulate 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 theimpact damage.
2. Methods
Incremental finite element formulation for linear Timoshenko beam element, which is usedthroughout this study, is summarized in this section. Then, the existing ASI technique, the newly
proposed ASI-Gauss technique and algorithms considering member fracture and elemental contact
are explained.
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)

In Eq.(1), δV and δW are the internal work and the external work respectively. {∆ε}, {∆σ},{∆u}, and {∆f} are the generalized strain increment vector, the generalized stress (resultant force)increment 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)Here, [B(s)] is the generalized strain-nodal displacement matrix. s is the location of the numerical

w z

L/2
L/2

θz

2
u

θy

1

θx
y

x

v

Fig.1. Definitions of local coordinates and displacements

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 asz/(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)Here, [De] is the stress-strain matrix for elastic deformation. [De] will change to [Dp], when plastic
deformation is determined by the plastic flow theory using the following form of yield function:f y (σ (r )) − 1 = 0(4)

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 )](5)

In the ASI technique, numerical integration points are shifted immediately after the formationof a fully plastic section in the element in order to form a plastic hinge exactly at that section.
When a plastic hinge is judged to be unloaded, the corresponding numerical integration point isshifted back to its normal position. Here, the normal position means the location where thenumerical integration point is placed when the element acts elastically.
Figure 2 shows a linear Timoshenko beam element and its physical equivalence to arigid-body-spring-model (RBSM). As shown in Fig.2, the relationship between the locations of thenumerical integration point and the stress evaluation point where a plastic hinge is actually formed
is expressed as follows:

r1 = − s1(6a)s1 = −r1(6b)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.(6) in this case. The
elemental stiffness matrix for the element is given by the following equation:

-1
θ1
u1

s1

0
s

Linear Timoshenko beam element
Numerical integration point

1
θ2
u2

θ1

r1
r

θ2

u1

u2
RBSM

Rotational and shear spring connecting rigid bars
(Plastic hinge including the effect of shear force)Fig.2. Linear Timoshenko beam element and its physical equivalent[ K ] = L[ B (0)]T [ De (0)][ B (0)](7a)

The generalized strain increment vector {∆ε(0)} and the resultant force increment {∆σ(0)} are
calculated as follows:{∆ε (0)} = [ B (0)]{∆u}{∆σ (0)} = [ De (0)]{∆ε (0)}

(7b)
(7c)

The resultant forces calculated from Eq.(7c) are physically those at the midpoint of theelement. However, the resultant forces except for bending moments have constant values in theelement. According to the elementary beam theory, the relations between the bending moments andshear forces are given as below.

Vx = −

dM y

dz
dM x
Vy = −
dz

(8a)
(8b)

Here, Mx and My are the bending moments around x- and y axes while Vx and Vy are the shear
forces along x- and y- axes. The bending moments in an 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 )

(9a)
(9b)

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 right endof an element (r1= 1), the numerical integration point is shifted to the left end (s1= −1) and viceversa. 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 )](10a){∆ε (r1 )} = [ B( s1 )]{∆u}{∆σ (r1 )} = [ D p (r1 )]{∆ε (r1 )}

(10b)
(10c)

Here, [Dp] is the stress-strain matrix for plastic deformation.
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 isshifted 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 numerical integration points.
The relationship between the locations of the numerical integration point and the stressevaluation point is expressed as follows:

rg = − s g(11a)s g = − rg

(11b)

When the entire region of an element behaves elastically, the numerical integration point is locatedat (s1= sg). As shown in Fig.3, the normal position of the numerical integration point in theMember
-1

-1/√3

1/√3

0

1

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

rg 0 sg

sg 0 rg-1 -1

Element 1

1

Stress evaluation point

Element 2

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

ASI-Gauss technique is sg= 1-(2/√3) and the stress evaluation point is rg= −sg = −1+(2/√3).
When two consecutive elements forming a member are considered as in Fig.3, the two stressevaluation points coincide with Gaussian integration points of the member. This means that stresses
and strains are evaluated at Gaussian integration points when the element behaves elastically.
Gaussian integration points are known to be optimal for two-point integration and the accuracy ofbending deformation is mathematically guaranteed [9]. 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 bending moments in an elastically deformed element are determined by the following
approximations:L(s + sg )2
L∆M y (s) = ∆M y (sg ) − ∆V x (sg ) (s + sg )2
∆M x (s) = ∆M x (sg ) − ∆V y (s g )

(13a)
(13b)

Immediately after the occurrence of a fully plastic section, the numerical integration point isshifted 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 shifted to 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.(10a), (10b) and (10c).
The yield function assumed in this study is expressed by the following equation:
2

2

2

2

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

(14)Here, Mx, My, N and Mz are the bending moments around x- and y- axes, the axial force and thetorsional 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 of shear forces isneglected in the yield function.
The dynamic collapse analysis is made possible by adding inertial force terms to thequasi-static computational procedure. That is, the following equation for the complete system is
integrated by an appropriate time integration scheme:..

[ M ]{u}t + ∆t + [ K ]{∆u} = {F }t + ∆t − {R}t

(15)

Here, [M] is the total mass 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 and {∆u} the nodal

displacement increment vector between t and t+∆t. In this study, we use consistent mass as anelemental mass matrix which is assumed to be unaffected by the shifting of numerical integration
points. Newmark’s β method is used as a numerical integration scheme for Eq.(15). And conjugate
gradient (CG) method is used as a solver to reduce the need of memory resources.
2.3 Algorithm considering member fracture
In this study, member fracture is determined by examining axial and bending strains in twoelements forming a member. Member fracture is expressed by shifting the numerical integration
point to the opposite end of the fractured section and by releasing the resultant forces exerting on
the section simultaneously as shown in Fig.4. This causes a virtual plastic hinge to form at thefractured section. For example, if fracture occurs at the right end of the element (r=1), thenumerical integration point is shifted to the left end (s=-1) and the virtual plastic hinge is formed at
the right end.
Fracture is judged to occur in an element when at least one the conditions in the following
equation are satisfied for the element.



κ 
 κx 
 − 1 ≥ 0,  y  − 1 ≥ 0,  ε z  − 1 ≥ 0

ε 
κ 
κ 
 fz 
 fy 
 fx 

(16)

In the above equation, κx, κy, εz, κfx, κfy and εfz are the bending strains around x- and y- axes, theaxial strain and the critical values for each strain, respectively. The values κfx=κfy=0.0004 andεfz=0.3 are used in the analyses.
In the Finite Element Method, the mass of an element is distributed at its nodes. So, whenmember fracture occurs, the mass matrix [M] needs to be recalculated according to the newconnecting conditions of the elements. Figure 5 shows the distribution of mass before and after theoccurrence of fracture.
Releasing resultant forces exerting on the fractured section may cause longitudinal waves in
the model. Therefore, a gradual release of resultant forces or introducing artificial viscosity is
considered necessary to reduce this effect. However, an instant release technique is used in thisstudy to preserve the simplicity of the algorithm.
It should also be noted that new virtual node points for the fractured sections are added at theMember fractureMember fracture(a) Linear Timoshenko beam element(b) RBSM

Numerical integration point

Rotational and shear spring connecting rigid bars
(Plastic hinge including the effect of shear force)

Fig.4. ASI-Gauss technique dealing with member fracture

Element 1
m/2

Element 2

Element 1
m/2

m
Before fracture

m/2

Element 2
m/2

m/2

After fracture

Fig.5. Mass distribution before and after the occurrence of fracturem/2post-processing stage as elements and node points are still treated as a continuous model in thecalculation process. The elements with the virtual node points are then visualized as rigid bars
thereafter.
2.4 Algorithm considering elemental contact
Generally, it is difficult to determine the loads resulted in the structures due to an impact. And
applying impact loads to the model in the form of nodal forces may not well simulate the impact
phenomenon. In this study, we try to simulate the impact phenomenon by means of contacts
between the elements involved. Elemental contact algorithm is also useful to simulate fractured
elements colliding with other elements. If this algorithm is not implemented in the analyses,
unrealistic appearances such as members going through each other could take place.
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 6 shows the coordinates used 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), the condition under which all
four nodal points lie on the same plane can be expressed as follows: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 )( y i 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 each otheras 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 are
A1

Bi1

Colliding element Ec

Element Ei

A2
Bi2
Fig.6. Coordinates used in determination of elemental contact

Bi1
Colliding element Ec

A1
A2
Gap elements

Element Ei
Bi2

Fig.7. Four gap elements binding two elements in contact

the lengths of colliding element Ec and the element Ei respectively. Cl is the coefficient concerning
contact length and a value of 1.8 is used in this study.
For the elements that do not lie on the same plane but exist within a certain distance specified
by Eq.(18), another condition is given as expressed by the following equation, which assumes theelements nearly forming a plane:f ( x, y , z ) ≤ α(19)

In this study, a value of 3.0x10-4 m3 is used for α. Contact determination is skipped for the elementsthat are not approaching each other. When two elements are determined to be in contact, they arebound with a total of four gap elements as shown in Fig.7. Gap elements are assumed to have thesame physical properties and the cross-section as the other elements. After certain time, the gap
elements are removed from the analytical model. In this study, this binding time is set to 0.5ms. It
is assumed that no energy loss takes place during the impact.
3. Numerical examples
In 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 methods in order to verify the accuracy of the newly proposed ASI-Gauss
technique: (a) conventional finite element method in which the numerical integration point of eachelement is fixed at its midpoint; (b) finite element method using the ASI technique; (c) finite
element method using the ASI-Gauss technique. Moreover, the algorithm considering elemental
contact is verified from the point of conservation of energy.
3.1 Elasto-plastic collapse behaviors of framed structure under static loading
Figure 8 shows the analyzed space frame with its dimensions, physical properties and material
constants. As shown in the figure, the frame is fixed at its lower ends and a static horizontal load isapplied 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 not expected in the analyses, thetotal Lagrangian formulation (TLF) is used.
Figures 9 and 10 show the results of the elastic and elasto-plastic analyses under static loading.
In the elastic analysis in which no plastic section is assumed to form within the elements, the ASI
technique shows slow convergence and at least eight elements are needed for a member to obtainthe converged solution. In the elasto-plastic analysis using the ASI technique, the collapse load
obtained from the two-element model shows little difference from that of the converged solution
since the numerical integration points are shifted to form plastic hinges at appropriate locations.
However, it still shows slow convergence in the elastic range. As for the conventional method,
convergence is extremely slow for both analyses. On the other hand, the ASI-Gauss techniqueshows that just two-element modeling is needed to obtain the converged solution both in elastic andelasto-plastic analyses.

L = 2.00 m , A = 2.50 × 10 −3 m 2

L

I x = I y = 5.21 × 10 − 7 m 4
P,δ

L

L

Z

ρ = 7900 kg / m 3 , E = 206 GPaυ = 0.3, σ y = 245 MPa
M x 0 ≡ M y 0 = 7.66 × 10 3 Nm

Y

N 0 = 6.13 × 10 5 NX

M z 0 = 1.12 × 10 5 Nm
Fig.8. Analyzed space frame

3.2 Elastic and Elasto-plastic responses of framed structures under dynamic loading
In the analyses subjected under dynamic loading, the time increment ∆t is set to 0.05ms andthe same model shown in Fig.8 is used. A load of 8kN is applied in the elastic response analysis
70

70
P,δ

60

60

50

50

Load P (kN)

Load P (kN)

P,δ

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

20
10
0

100

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

20
10

200

0

(a) Conventional method(a) Conventional method
70

P,δ

P,δ

60

60

50

50

Load P (kN)

Load P (kN)200

Displacement δ (mm)

70

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

20
10
0

100

40
30
20

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

10

200

0

200

Displacement δ (mm)(b) ASI technique(b) ASI technique
70

P,δ

P,δ

60

50

50

Load P (kN)

60

40
30
20

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

10
0

100

Displacement δ (mm)70

Load P (kN)

100

Displacement δ (mm)

100

200

40
30
20

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

10
0

100

200

Displacement δ (mm)Displacement δ (mm)(c) ASI-Gauss technique(c) ASI-Gauss techniqueFig.9. Load-displacement relations(elastic analysis)Fig.10. Load-displacement relations(elasto-plastic analysis)

and 12kN in the elasto-plastic response analysis. Loading is done for a total of 100 steps (5ms) inboth analyses. β=1/4 is used for Newmark’s β method.
Figures 11 and 12 show the results of the elastic and elasto-plastic response analyses underdynamic loading. As in the analyses described in the previous section, the conventional method
shows very slow convergence and sixteen-element modeling is necessary to obtain the converged
solution in the elasto-plastic response analysis. The ASI technique gives comparatively better
100P,δ

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

P,δ

80

Displacement δ (mm)Displacement δ (mm)60

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

40

20

60

40

20

0.1

0.2

0.3

0.4

0.1

0.5

Time (s)

0.2

0.3

0.4

0.5

Time (s)(a) Conventional method(a) Conventional method
100

2 elements/member
4 elements/member
8 elements/member
16 elements/member
32 elements/memberP,δ

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

P,δ

80

Displacement δ (mm)Displacement δ (mm)

60

40

20

60

40

20

0.1

0.2

0.3

Time (s)(b) ASI technique

0.4

0.1

0.5

0.4

0.5

(b) ASI technique

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

40

20

P,δ

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

80

Displacement δ (mm)Displacement δ (mm)

60

0.3

Time (s)

100
P,δ

0.2

60

40

20

0.1

0.2

0.3

0.4

Time (s)(c) ASI-Gauss techniqueFig.11. Elastic responses of a space frame

0.5

0.1

0.2

0.3

0.4

0.5

Time (s)(c) ASI-Gauss techniqueFig.12. Elasto-plastic responses of a space frameresults than the conventional method. However, when it comes to the model that uses only twoelements per member, poor convergence is seen and the vibration mode is much different from thatof the converged solution.
On the other hand, the ASI-Gauss technique shows very fast convergence in both analyses andnearly 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 efficient method foranalyses of frame structures both in elastic and elasto-plastic range, and even the two-element
modeling is sufficient to acquire the results with high accuracy and reliability.
3.3 Impact analysis using a simple two-beam model
In this study, impact phenomenon of two beams is simulated to verify the creditability of theelemental contact algorithm. The analyzed model, its physical properties and material constants areshown in Fig.13. One of the beams is fixed at its both ends while the other one moves freely. Initial
velocity of 50m/s is applied to the free beam. The beams are modeled with two elements per eachmember. 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 theanalysis. Although using a value between 1/6 and 1/4 for β in Newmark’s β method gives moreaccurate results, 0.4 is set for β in order to stabilize the calculations [6]. The time increment ∆t is
set to 10µs. The calculation is done for a total of 3000 steps (0.03s).
Figure 14 shows the sum of kinetic and potential energies, strain energy and the total of all

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

ZY

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

X

Fig.13. Analyzed two-beam model

Total Energy
Kinetic Energy + Potential Energy
Strain Energy

40
Energy (×103J)

L/2

30
20

10
0
0

5

20
10
15
Time (ms)Fig.14. Kinetic, potential and strain energies in the two-beam modelthree for the whole model. Energies for both the colliding and the fixed beams constituted by a totalof 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 andthe length of the element Ei (i=1,2,3,4). Subscripts 1 and 2 denote the nodal points of thecorresponding element.

0.000s

0.010s

0.014s

0.018s

0.022s

0.026sFig.15. Collision analysis of a simple two-beam model

As seen from the figure, kinetic and potential energies occupy all of the total energy before theimpact takes place. After the impact, both beams deflect and vibrate as seen in Fig.15. As aconsequence, the strain energy is stored up in the beams and the sum of the other two decreases.
However, the total energy remains almost unchanged although fluctuations are seen after the
impact. The main reason for these fluctuations is considered to be the occurrence of longitudinaland bending waves. Nevertheless, the elemental contact algorithm used in this study is considered
to be sufficiently reliable.
4. Structural collapse analysis of framed structures under impact loads
In this section, dynamic collapse analyses are performed using ASI-Gauss technique tosimulate the impact of a small aircraft into a ten-storied steel framed structure. A total of four cases
are analyzed using different values for three parameters: cruise speed, size and weight of theaircraft. Details of the analyzed model, propagation of impact loads, response of the structures andthe extent of impact damage are discussed in this section.
4.1 Outline of the analyzed model
Figure 16 shows the analyzed ten-storied framed structure and the aircraft. Their dimensions,
physical properties and material constants are also shown in the figure. It is assumed that theten-storied framed structure is made of steel and the aircraft duralumin.
The weight of the floor system, which amounts to twenty times the own weight of the whole
structure, is applied to the model in the form of an increase in the density ρ of the beams. The
ten-storied framed structure is fixed at its lower ends while none of the element is fixed in theaircraft model. It is assumed that the total weight of the aircraft is uniformly distributed in its
elements. And this is implemented by increasing the density ρ of all its elements. The aircraft has
zero degree of inclination from the horizon and an initial velocity is applied to it. Both structures
are modeled by using two elements per each member. The total number of elements is 1698, and

L/4
L/4

2Lc/5 Lc

L/4
2Lc/5

L/4

Body 1

Body 2

Unit

L

16

-

m

Lc

-

15

m

h

30

-

m

hc

-

2

m

ha

-

21.50

m

A

2.19x10-22.50x10-3

m2

Ix

6.66x10-4

5.21x10-7m4

Iy

2.24x10-4

5.21x10-7

m4

ρ

7900

2800

kgm-3

E

206

69

Gpa

σy

240

195

Mpa

ν

0.30

0.335

Mxo

7.99x105

6.09x103Nm

Myo2.69 x105

6.09x103Nm

No

5.25x106

4.88x105

N

Mzo

5

5

Nm

2Lc/5

L/4

L/4

L/4

L

Y
Lc/5

2Lc/5

3Lc/10

X

h/10
h/10
h/10

hc
Lc

h/10
h/10
h/10

ha
h/10

Z

h/10

1.12x101.12x10

Body 1: ten-storied steel framed structureX

h/10

Body 2: colliding body (made of duralumin)

h/10

Fig.16. Analyzed ten-storied steel framed structure and colliding bodythe total number of freedom is 7308.
4.2 Parameters set in the analyses
A total of four analyses are performed using different values for cruise speed, size and weightof the aircraft. The combinations of these parameters are shown in Table 1.
Table 1 Parameters of the aircraft

Analyzed cases

Cruise speed [kmph]

Size Lc [m]Weight [kg]Case I

440

15

15000

Case II

320

15

15000

Case III

440

15

12000

Case IV

440

12

15000

In all the four analyses, the ASI-Gauss technique with the linear Timoshenko beam element isused and member fracture and elemental contact are both considered. Damping effect is notconsidered in all analyses. The updated Lagrangian formulation (ULF) is used to be able to copewith large deformations of the elements. The Conjugate Gradient (CG) method is used as the solver

to reduce the need for memory resources. β=0.4 is used for Newmark’s β method and the timeincrement ∆t is set to 10µs. The calculation is done for a total of 25000 steps (0.25s) for each case.
The calculation for each case takes approximately 9 CPU hours with a personal computer (2.4GHz
CPU, 256MB RAM).
4.3 Structural collapse analysis of the ten-storied framed structure under impact loads
Figure 17 shows the results of the first analysis (Case I). A total of six colors are used to

0ms

2.5ms

60ms

90ms

120ms

150ms

180ms

210ms

250ms

fy

0.00

0.25

0.50

0.75

1.00

fracture

Fig.17. Impact analysis of a small aircraft against a ten-storied steel framed structure (Case I)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 theelement 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 of impact. And
the elements near the area of impact turn to black, that is, they start to fracture. Moreover, it is seenthat 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 shock waves. And red-colored
elements are temporarily found even in the lowest floor, and this means that even the lowest floor
is affected, to some extent, by the impact.
Figures 18(a) to 18(d) show the extent of damage to the steel framed structure for all four(a) Case I
Cruise speed=440kmph, Lc =15m, Weight=15000kg
Pre-impact kinetic energy = 112x106J
Number of fractured elements=713(b) Case II
Cruise speed=320kmph, Lc =15m, Weight=15000kg
Pre-impact kinetic energy = 59 x106J
Number of fractured elements=700(c) Case III
Cruise speed=440kmph,Lc =15m, Weight=12000kg
Pre-impact kinetic energy = 90x106J
Number of fractured elements=568(d) Case IV
Cruise speed=440kmph,Lc =12m, Weight=15000kg
Pre-impact kinetic energy = 112x106J
Number of fractured elements=521

fy 0.00

0.25

0.50

0.75

1.00

Fig.18. Impact damage to the ten-storied framed structure(only the upper five floors are shown and fractured elements are excluded)

cases. In order to show the damaged areas clearly, only the upper five floors are plotted and thefractured elements are excluded. If we focus on Case I, II and III, the largest damage is seen inCase I in which the damage extends to the sixth floor and most of the inner columns are cut off.
Pre-impact kinetic energy is also the largest in Case I, and so it can be said that pre-impact kinetic
energy of the aircraft is one of the factors that influence the impact damage. However, in spite of
the smaller pre-impact kinetic energy, the impact damage in Case II is considerably larger that thatin Case III. This indicates that the weight of the aircraft has a stronger influence than the velocity.
If we focus on Case I and IV, it is seen that less inner columns are cut off in Case IV although
the kinetic energy of the aircraft is the same. This result can be interpreted in such a way that thesize of the aircraft has a strong influence on the impact damage. However, it is seen that the smaller
the aircraft, the larger the number of its members that pass through the ten-storied framed structure
without collision. Therefore, further study is needed to clarify the influence of the size of theaircraft. Moreover, soon after impact, tensile stresses are observed in the columns that had been
compressed by gravity loads before impact as shown in Fig.19. There is a possibility that thesetensile stresses lead to some degree of structural damage especially in the joints of beams andcolumns. Therefore, detailed analyses are necessary to investigate this phenomenon.

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
Fig.19. Tensile stresses observed in the ten-storied framed structure after impact (Case I)

5. Conclusion
The formerly developed ASI technique for the linear Timoshenko beam element provides
higher computational efficiency than the conventional finite element method, and is known to beable 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 numberof elements per member is two. It gives results with 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 is simulated
by means of contacts between the elements involved. This elemental contact algorithm is also
useful to simulate fractured elements colliding with other elements. 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 areapplied to dynamic analyses with strong discontinuities, where impact of an aircraft against aten-storied framed structure body is discussed using four different sets of parameters 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 by different sets of parameters, and in reasonably
short calculation time. Therefore, the numerical analysis code used in this study is considered to bean efficient code for impact analyses of large-scale framed structures. However, further study is
considered necessary to clarify the effect of the size of aircraft on the impact damage and toevaluate the structural damage caused by the tensile stresses of the columns.
References[1] Cundall, P.A., A Computer Model for simulating Progressive, Large-scale Movement in Blocky
Rock System, Proceedings of the International Symposium on Rock Mechanics, II-8, pp.129-136,
1971.
[2] Shi, G.H. and Goodman, R.E., Discontinuous Deformation Analysis, Proceedings of 25th U.S.
Symposium on Rock Mechanics, pp.269-277, 1984.
[3] Tagel-Din, H. and Meguro, K., Nonlinear Simulation of RC Structures Using Applied Element
Method, Structural Engineering / Earthquake Engineering, Vol.17, No.2, pp.137-148, 2000.
[4] Toi, Y. and Isobe, D., Adaptively Shifted Integration Technique for Finite Element Collapse
Analysis of Framed Structures, Int. J. Numer. Methods Eng., Vol.36, pp.2323-2339, 1993.
[5] Toi, Y. and Isobe, D., Finite Element Analysis of Quasi-static and Dynamic Collapse Behaviors of
Framed Structures by the Adaptively Shifted Integration Technique, Computers and Structures,
Vol.58, No.5, pp.947-955, 1996.
[6] Isobe, D. and Toi, Y., Analysis of Structurally Discontinuous Reinforced Concrete Building
Frames Using the ASI Technique, Computers and Structures, Vol.76, No.4, pp.471-481, 2000.
[7] Isobe, D. and Morishita, M., Debris Impact Analysis of Large-Scale Space Structures by Using
ASI-FEM, Transactions of the Japan Society of Mechanical Engineers, Part A, Vol.64, No.627,
pp.2726-2733, 1998, in Japanese.
[8] Isobe, D. and Tsuda, M., Development of Seismic Collapse Analysis Code for RC Framed
Structures Using Finite Element Method, Journal of Structural Engineering, Vol.48B, pp.385-394,
2002, in Japanese.
[9] Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T., Numerical Recipes in
FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge University Press, pp. 140-155,
1992.