Aircraft Impact Analyses of the World Trade Center Towers
D. Isobe and Z. Sasaki
University of Tsukuba, Japan

Abstract
Although some official statements have been released by the Federal Emergency Management
Agency (FEMA) and the National Institute of Standards and Technology (NIST), the 9/11
tragedy of the New York World Trade Center (WTC) towers still remains as an unresolved
mystery since the progressive collapse sequence of the towers was unnaturally rapid. We
sought the true cause of the collapse by performing full-model, dynamic collapse analyses ofthe WTC towers due to aircraft collision. The numerical code used in the analyses is based on
the ASI-Gauss technique, a modified version of the formerly developed adaptively shifted
integration (ASI) technique, which computes highly accurate elastoplastic solutions even withthe minimum number of elements per member. The ASI-Gauss technique gains still a higher
accuracy particularly in the elastic range, by placing the numerical integration points of the twoconsecutive elements forming an elastically deformed member in such a way that stresses andstrains are evaluated at the Gaussian integration points of the two-element member. Moreover,
the technique can be used to express member fracture, by shifting the numerical integration
point to an appropriate position and by releasing the resultant forces in the element
simultaneously. The numerical results show instant disconnections between the core columns
and the main double trusses during the impact, and a spring-back phenomenon of the corecolumns due to a very rapid unloading is observed. The spring-back phenomenon is fatalenough to produce a gigantic tensile force on the columns, which might have triggered the total
fracture of some connections, eventually leading to the total collapse of the towers.
Keywords: World Trade Center towers, 9/11, dynamic unloading, spring-back, ASI-Gauss
technique, finite element analysis______________________________________
Daigoro Isobe
University of Tsukuba
1-1-1 Tennodai Tsukuba-shi,
Ibaraki 305-8573
Japan
Email: isobe@kz.tsukuba.ac.jp
Tel: +81-29-853-6191

1.0 Introduction
The 9/11 terrorist attack on the New York World Trade Center (WTC) towers caused an
unprecedented tragedy in the history of architecture. The twin towers WTC 1 and 2 stood inflames caused by jet fuel until finally, both collapsed totally to the ground with thousands ofpeople trapped in the buildings. Both towers collapsed at an unnaturally high speed, which wasobserved to be nearly equal to that of free fall. Official statements have already been released by
the Federal Emergency Management Agency (FEMA) in 2002 [1], and also by the National
Institute of Standards and Technology (NIST) in 2005 [2], regarding the incident. FEMA [1]concluded that the heat of burning jet fuel induced additional stresses into the damaged structural
frames while simultaneously softening and weakening these frames, and this additional loading
and the resulting damage were sufficient to induce the collapse of both structures. Many detailed
numerical analyses were carried out in [2], and the report was concluded by stating “the WTC
towers likely would not have collapsed under the combined effects of aircraft impact damage andthe extensive, multifloor fires, if the thermal insulation had not been widely dislodged or had been
only minimally dislodged by aircraft impact.” However, did the fire really spread to such a wide
range in the buildings to cause the total collapse? Moreover, was the heat high enough to reduce
the strengths of structural members? What really caused the free-fall total collapse? Many
questions still remain unresolved. Also, it is to be noted that although many detailed analyses were
performed in [2], not even a single result can be found showing the dynamic behaviors of thetowers during aircraft impact. It is clear that the towers experienced an extreme dynamic load thatno other high-rise building has ever experienced in history. Therefore, we performed some aircraft
impact analyses of a full-model WTC tower, to determine the possibility of the impact itself acting
as a fatal cause of the collapse.
The numerical code used in the analyses is based on the ASI-Gauss technique [3], which wasdeveloped and modified by Isobe et al., from the adaptively shifted integration (ASI) technique [4]for the linear Timoshenko beam element. It enables the computation of highly accurate
elastoplastic solutions even with the minimum number of elements per member. In this paper, wefirst describe an outline of the technique and the modeling of the WTC tower, followed by somenumerical results on WTC 2, for which some observed data are available and easy to verify.

2.0 Numerical Code
2.1 ASI-Gauss Technique
Toi and Isobe developed the ASI technique [4-6] for the linear Timoshenko beam and BernoulliEuler beam elements, which can easily be incorporated into existing finite element codes. In thistechnique, 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 elastoplastic solutions than the conventional schemes, andenables us to simulate dynamic behaviors with strong nonlinearities using only a small number of
elements for a member.

Figure 1. Linear Timoshenko beam element and its physical equivalent.
Figure 1 shows a linear Timoshenko beam element and its physical equivalence to the rigid
bodies-spring model (RBSM). As shown in the figure, the relationship between the locations of
the numerical integration point and the stress evaluation point where a plastic hinge is actuallyformed is expressed as [7]r = Äs .

(1)

In the above equation, s is the location of the numerical integration point and r is the location
where stresses and strains are actually evaluated. s is a nondimensional quantity, which takes avalue 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 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. By
doing so, the plastic behavior of the element is simulated appropriately, and the converged
solution is obtained 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 theelement behaves elastically. When the number of elements per member is very small, solutions in
the elastic range are not accurate since one-point integration is used to evaluate the low-orderdisplacement function of the beam element.
The main difference between the ASI and ASI-Gauss techniques lies in the normal position of thenumerical integration point. In the ASI-Gauss technique, two consecutive elements forming amember are considered as a subset, and the numerical integration points of an elastically deformed
member are placed such that the stress evaluation points coincide with the Gaussian integration
points of the member. This means that stresses and strains are evaluated at the Gaussian
integration points of elastically deformed members. Gaussian integration points are optimal fortwo-point integration and the accuracy of bending deformation is mathematically guaranteed [8].
In this way, the ASI-Gauss technique takes advantage of two-point integration while using onepoint 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 elementsin the ASI and ASI-Gauss techniques. The elemental stiffness matrix, generalized strain andresultant force increment vectors for the ASI technique are given by Eqs. (2a)-(2c) and those for
the ASI-Gauss technique by Eqs. (3a)-(3c).
[K] = L[B(0)]T [D(0)][B(0)]fÅ"(0)g = [B(0)]fÅug
fÅõ(0)g = [D(0)]fÅ"(0)g[K] = L[B(sg )]T [D(rg )][B(sg )]fÅ"(rg )g = [B(sg )]fÅugfÅõ(rg )g = [D(rg )]fÅ"(rg )g

(2a)
(2b)
(2c)
(3a)
(3b)
(3c)

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. sg and rg
are the locations of the numerical integration and stress evaluation points for the ASI-Gauss
technique in the elastic range. They take values as shown in Fig. 2. [B] is the generalized strainnodal displacement matrix, [D] is 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

+

í

Mz
Mz0

ì2

Ä 1 = fy Ä 1 = 0 .

(4)Here, fy is the yield function, and Mx , My, N and Mz are the bending moments around the x- and yaxes, axial force and torsional moment, respectively. Those with the subscript 0 are values that
result in a fully plastic section in an element when they act on a cross section independently. The

effect of shear force is neglected in the yield function.
2.2 Member Fracture and Contact Algorithm
Structurally discontinuous problems have also become easily handled in the ASI and ASI-Gauss
techniques, 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 [3,6].
Elemental contact is considered in the code by introducing gap elements between two elements
determined to be in contact by geometrical relations. More details on the code can be found in [36].

North face

East face

West faceCore beamSouth faceCore column

Main double
truss(a) Upper view

Figure 3. Analyzed model.

(b) South view(c) East view

PerimetercolumnDiagonal brace
memberIntermediate
support angle

Transverse truss

SpandrelFigure 4. Cross section of the WTC towers.

3.0 Numerical Model and Analytical Conditions
3.1 Numerical Model
A full model of the WTC tower is constructed to investigate the overall effect of aircraft impact on
the building. Details on the structural members and construction data are extracted and adopted
from the reports [1,2]. The tower and the B767-200ER aircraft are both modeled with linear
Timoshenko beam elements, of which all members are subdivided into two elements. Figure 3
shows a schematic view of the WTC 2 tower model and the aircraft model. The tower model
contains 604,780 elements, 435,117 nodes and 2,608,686 degrees of freedom, while the aircraftmodel contains 4,322 elements, 2,970 nodes and 17,820 degrees of freedom. Figure 4 shows across-sectional view of the towers. A floor is roughly divided into two sections, external and core
structures, each with different roles from the structural point of view. The external structure ismainly designed to sustain wind loads, while the main double trusses connected to the core
structure support all the floor loads. The whole tower is subjected to a designed dead load of 2,890

MN and 40 % of the allowed capacity load, 300 MN. The ratio of the loads between the core andexternal structures is adjusted to become approximately 6:4 [1].
All structural members of the B767-200ER aircraft are assumed to have box-shaped cross sections
and the material properties of extra super duralumin. The total mass of the aircraft at the time of
impact is 142.5 t, which is the sum of the masses of the aircraft (112.5 t) and the jet fuel (30 t) [1].
The mass of each engine is assumed to be 19.315 t. The nose of the aircraft is tilted 11.5 degrees tothe east and 5 degrees downward, and its left wing is inclined downward by 35 degrees. It is
assumed to collide with the 81st floor, south face of WTC 2 with a cruise speed of 590 mph (262m/s).
3.2 Analytical Conditions
Most of the column joints were connected with 34-mm-thick end plates attached by 4 to 6 bolts.
The joints were weak against flexural and shear forces; the simple moment capacity of the bolt
group was 20 to 30 % of the plastic moment capacity of the column itself [1]. The beam joints ofmain double trusses also had simple connections to sustain static, vertical dead loads, which
satisfied the initial structural design. The structural design of member joints in Japan, for example,
demands the capacity to be greater than that of the member itself. In the present study, thevulnerability of member joints are considered by introducing fracture conditions as follows [3]:
å
å
å
å
ì
í
å îx åååèzåå Ä 1 ï 0 or å îy å Ä 1 ï 0 or
Ä 1 ï 0,åîf x ååîf y åèf z

(5)where κx, κy, εz, κfx, κfy and εfz are the bending strains around the x- and y-axes, the axial tensile
strain and the critical values for these three strains, respectively. We fixed the critical axial tensile
strain values by comparing the impact damage and motion of the engines with the observed data.
We still need to discuss critical bending strains in the future, since member joint characteristics are
currently not well known.
We also considered the strain rate effect on yield strength using the Cowper-Symonds equation [9].
An implicit scheme with a consistent mass matrix is adopted in this study. No damping except for
the algorithmic one is applied. δ is set to 5/6 and β to 4/9 in Newmark’s β method. The timeincrement ∆t is set to 0.2 ms. An updated Lagrangian formulation is used since large deformations
are expected in the analysis, and the conjugate gradient (CG) method is used as the solver toreduce the memory requirement. The analytical code is run on a high-performance computer (1.4
GHz Itanium*2, 8GB RAM), and the calculation takes approximately two months for a physicaltime of 0.7 s.

4.0 Aircraft Impact Analysis of WTC 2
The numerical results are shown in Figs. 5 and 6. Movie files can be downloaded at the website[10]. The left engine reduces its speed rapidly as it directly enters the core structure. The rightengine, on the other hand, glances off the core structure and goes out of the north-east corner ofthe tower. Table 1 shows the impact timeline of the fuselage and both engines. The velocity curveof the right engine is shown in Fig. 7 and a comparison of the damaged surfaces between the
numerical result and the observed data [2] is shown in Fig. 8. Both results are in good agreementwith the observed data.

(a) 0.12 s, (b) 0.28 s and (c) 0.56 s after impact
Figure 5. Aircraft impact analysis of WTC 2 (global view).

(a) 0.12 s, (b) 0.28 s and (c) 0.56 s after impact
Figure 6. Motion of fuselage and engines (upper view).
The dynamic transition of the resultant forces during impact is investigated. Figure 9 shows thetransition of axial forces in two typical core columns, at the site of every ten stories from theground floor to the top floor. Core columns are compressed constantly until the snout reaches thecore structure (Phase A in the figures), after which a wave due to impact and member fracture
propagates in horizontal and vertical directions. In particular, the compression decreases instantlyin the fractured core column (No. 1001) at the floors above the impact point. At the lower levels,
compression changes to tension immediately after the left engine hits the core structure (Phase C),
and continue to vibrate with large amplitude. Note that the lower the level goes, the larger theamplitude becomes. The core column (No.1001) at the 60th floor, for example, is heaved for 25
cm in 0.2 s as a result of this dynamic transition.
Table 1: Impact timeline of the B767-200ER aircraft
B767-200ER aircraft
Phase

Elapsedtime (s)
0
0.048
0.072
0.076
0.128
0.168
0.504
0.512

A
B
C
D
E
F

Left engine

Right engine

Snout

Snout hits outer-perimeter wall
Hits core structure
Hits outer-perimeter wall
Hits outer-perimeter wall
Hits core structure
Hits core column No.902
Hits core column No.903
Hits east side outer-perimeter wall

280●:Observed data [1]
240

Speed
(m/s)
Speed

200
160
120
80
40
A B

C

D

EF

0
0

0.1

0.2

0.3
0.4
ElapsedElapsedtime
time(s)

0.5

0.6

0.7

Figure 7. Velocity curve of the right engine.
What actually caused this dynamic unloading large enough to generate such large tensile forces?
We extracted the time history of the distribution of fractured main double trusses, as shown in Fig.
10, to investigate the cause. The propagation of fractured trusses from the impact point to upperand lower levels can be observed in the figure. This propagation speed almost matched thelongitudinal wave speed in the columns. It can be assumed that the shock wave due to impact
propagated through the tower and occasionally reflected to other columns at the ground or top
levels, destroying each member joint as it moved on. Our assumed scenario is as follows; theaircraft impact destroyed several main double trusses, that supported the floor slabs, and led to a
middle-class unloading of core columns. The deformation due to the unloading caused othertrusses to disconnect, which released further burden from the core columns. This disconnect-andrelease process advanced rapidly in a chain reaction manner and a spring-back phenomenon might
have occurred in many columns in such a way that a compressed spring is released in a very shortmoment. It was fatal enough to produce a gigantic tensile force on the columns, which might havetriggered the total fracture of some column joints, eventually leading to the total collapse of thetowers. It can be assumed that the damage in lower levels was larger, since the lower the level
went, the larger the tensile force became.

(a) Numerical result(b) Observed data [2]Figure 8. Damage of the WTC 2 south face.

AB

C

D

EF

AB

C

D

EF

Figure 9. Time histories of axial forces in core columns No.508 and No.1001.

5.0 Conclusion
The assumptions of member joint strength given in this work may not be sufficient for making any
clear statements. However, there is a possibility that the spring-back phenomenon due to rapid
unloading caused by aircraft impact acted as a fatal cause of the total collapse of the WTC towers.

Further numerical investigations performed by explicitly modeling the behaviors of member joints,
and also for the case of WTC 1 are in progress.

Figure 10. Distribution of fractured main double trusses.

6.0 Acknowledgements
This work was partially supported by a Grant-in-Aid for Scientific Research (No.#16206055),
Japan Society for the Promotion of Science. The authors also would like to thank Kyaw Myo Lynn,
a former graduate student of our lab, for contributing greatly to this work.

7.0 References
1. ASCE/FEMA, World Trade Center Building Performance Study: Data Collection, Preliminary
Observation and Recommendations, 2002.
2. NIST NCSTAR 1, Federal Building and Fire Safety Investigation of the World Trade Center
Disaster: Final Report on the Collapse of the World Trade Center Towers, National Institute ofStandards and Technology (NIST), 2005.
3. Kyaw, M.L. and Isobe, D. 2007. “Finite element code for impact collapse problems of framed
structures,” Int. J. Numer. Methods Eng., 69(12): 2538-2563.
4. Toi, Y. and Isobe, D. 1993. “Adaptively Shifted Integration Technique for Finite Element
Collapse Analysis of Framed Structures.” Int. J. Numer. Methods Eng., 36: 2323-2339.
5. 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 andStructures 58(5): 947-955.
6. 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.
7. Toi, Y. 1991. “Shifted Integration Technique in One-Dimensional Plastic Collapse Analysis
Using Linear and Cubic Finite Elements,” Int. J. Numer. Methods Eng., 31: 1537-1552.
8. Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P., Numerical recipes in
FORTRAN: The art of scientific computing, New York: Cambridge University Press, 1992.
9. Jones, N., Structural impact, New York: Cambridge University Press, 1989.
10. http://www.kz.tsukuba.ac.jp/~isobe/wtc-e.html