International Symposium on Structures under

Earthquake, Impact, andBlast Loading (IB'08)
October 10-11, 2008, JAPANNUMERICAL SIMULATIONS FOR INVESTIGATION ON TRUE CAUSE OF
THE TOTAL COLLAPSE OF THE WTC TOWERS
Daigoro ISOBE, Dept. of Engineering Mechanics and Energy, University of Tsukuba
1-1-1 Tennodai, Tsukuba-shi, Ibaraki 305-8573, JAPAN
Email: isobe@kz.tsukuba.ac.jp
Keywords: WTC towers, Collapse analysis, FEM, ASI-Gauss technique

SUMMARY
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 towers collapsed totally to the ground, and,in an unnaturally rapid speed. We applied a low-cost, highly accurate collapse analysis code using the ASI-Gauss
technique to several numerical simulations of framed structures, to seek the true cause of the highly-rapid totalcollapse. One of its applications is a full model simulation of the WTC 2 subjected under an aircraft impact.
According to the simulation, “springback” phenomena due to rapid unloading were observed in the core columns
during the impact, which might have caused the destruction of member joints. Another simulation carried out
was that of the fire-induced collapse of a high-rise tower. The results clearly show the effect of the weak member
joints, which were reported to be 20 to 30 % of the member strengths in the WTC towers, and also the effect of
the strength reduction due to elevated temperatures. It is possible that the highly-rapid total collapse of the WTC
towers was caused not only by the buckling and strength reduction of members due to elevated temperatures, butby the original weakness of the member joints and the destruction of those during the collision of the aircraft.

1. 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 in flames caused by jet fuel until
finally, both collapsed totally to the ground withthousands of people trapped in the buildings. Bothtowers collapsed at an unnaturally high speed, which
was observed to be nearly equal to that of free fall.
Official statements have already been released by theFederal Emergency Management Agency (FEMA) in
2002 [ASCE/FEMA, 2002], and also by the National
Institute of Standards and Technology (NIST) in 2005
[NIST NCSTAR1, 2005], regarding the incident.
FEMA concluded in the report that the heat ofburning jet fuel induced additional stresses into thedamaged structural frames while simultaneously
softening and weakening these frames, and thisadditional loading and the resulting damage were
sufficient to induce the collapse of both structures.
Many detailed numerical analyses were carried out in[NIST NCSTAR1, 2005], and the report was
concluded by stating “the WTC towers likely would
not have collapsed under the combined effects ofaircraft impact damage and the extensive, multifloor
fires, if the thermal insulation had not been widelydislodged or had been only minimally dislodged byaircraft impact.” But, did the fire really spread tosuch a wide range in the buildings to cause the totalcollapse? Moreover, was the heat high enough toreduce the strengths of structural members? Whatreally 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 [NIST NCSTAR1, 2005], not even asingle result can be found showing the dynamic
behaviors of the towers during aircraft impact. It isclear that the towers experienced an extreme dynamic
load that no other high-rise building has everexperienced in history. Therefore, we first performed
an aircraft impact analysis of a full-model WTC
tower, to determine the possibility of the impact itself
acting as a fatal cause of the collapse. We then carriedout some fire-induced collapse analyses of a high-rise
tower to verify the influence of fire patterns, member
strengths and member joint strengths against collapse
phenomena. The collapse phase is a very dynamic
phenomenon, but the phase before collapse initiation
is a very static one and the calculation cost may
increase to analyze the whole process. Therefore, thenumerical code must be a highly accurate,less-memory-consuming one, and has to be reliableMember
-1

1

-1/√3

r g 0 sg

0

1/√3

1

-1 -1

sg 0 r g

1

Element 1
Numerical integration point
Stress evaluation pointsg = 1 - (2/√3)
rg = -1 + (2/√3)

Figure 2 Locations of numerical integration andstress evaluation points in ASI-Gauss techniqueFigure 1 Linear Timoshenko beam element
and its physical equivalent
when member fracture and contact are taken into
account.
The numerical code used in the analyses is basedon the ASI-Gauss technique [Lynn and Isobe, 2007],
which was modified from the adaptively shifted
integration (ASI) technique [Toi and Isobe, 1993] forthe linear Timoshenko beam element. It enables thecomputation of highly accurate elastoplastic solutions
even with the minimum number of elements per
member. In this paper, we first describe an outline of
the technique and the modeling of the WTC tower,
followed by an aircraft impact simulation result on
WTC 2, and some fire-induced collapse simulation
results on a high-rise tower.
2. ASI-GAUSS TECHNIQUE
Figure 1 shows a linear Timoshenko beam element
and its physical equivalence to the RBSM. As shownin the figure, the relationship between the locations of
the numerical integration point and the stressevaluation point where a plastic hinge is formed isexpressed [Toi, 1991] asr = ÄsElement 2(1)

In the above equation, s is the location of thenumerical integration point and r is the location
where stresses and strains are evaluated. s and r are
nondimensional quantities that take values between-1 and 1.
In both the ASI and ASI-Gauss techniques, thenumerical integration point is shifted adaptively whena fully plastic section is formed within an element toform a plastic hinge exactly at that section. When theplastic hinge is unloaded, the corresponding
numerical integration point is shifted back to its
normal position. Here, the normal position means thelocation where the numerical integration point isplaced when the element acts elastically. By doing so,
the plastic behavior of the element is simulated

appropriately, and the converged solution is achieved
with only a small number of elements per member.
However, in the ASI technique, the numericalintegration point is placed at the midpoint of thelinear Timoshenko beam element, which isconsidered to be optimal for one-point integration,
when the entire region of the element behaveselastically. When the number of elements per member
is very small, solutions in the elastic range are notaccurate since one-point integration is used toevaluate the low-order displacement function of thebeam element.
The main difference between the ASI and
ASI-Gauss techniques lies in the normal position of
the numerical integration point. In the ASI-Gauss
technique, two consecutive elements forming amember are considered as a subset, as shown in Fig.
2, and the numerical integration points of anelastically deformed member are placed such that thestress evaluation points coincide with the Gaussian
integration points of the member. This means thatstresses 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 [Press et al, 1992]. In
this way, the ASI-Gauss technique takes advantage of
two-point integration while using one-point
integration in actual calculations.
Structurally discontinuous problems also become
easily handled using the ASI and ASI-Gauss
techniques by shifting the numerical integration pointof the linear Timoshenko beam element to an
appropriate position and by releasing the resultant
forces in the element simultaneously. Material strains
are used for the fracture criterion. Elemental contact
is considered in the code by introducing gap elements
between pairs of elements determined to be in contact
by geometrical relations. More details on the code
can be found in [Lynn and Isobe, 2007] and [Isobe
and Sasaki, 2007].

North face

East face

West face
South face(a) Upper view(b) South view(c) East view

Figure 3 Analyzed modelcross-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 ofview. The external structure is mainly designed tosustain wind loads, while the main double trusses
connected to the core structure support the floor loads.
The whole tower is subjected to a designed dead loadof 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 [ASCE/FEMA, 2002]. All
structural members of the B767-200ER aircraft are
assumed to have box-shaped cross sections and thematerial properties of extra super duralumin. The
total mass of the aircraft at the time of impact is142.5 t, which is the sum of the masses of the aircraft(112.5 t) and the jet fuel (30 t). The mass of eachCore beamCore column
Maindouble truss
Spandrel
Perimeter
Diagonal
column brace member

Transverse truss
Intermediate
support angle

Figure 4 Cross section of the WTC towers
3. NUMERICAL INVESTIGATIONS
Several numerical simulations are performed on afull-model WTC tower and a high-rise tower, to seek
the true cause of the highly-rapid collapse of the
WTC towers.
3.1 Aircraft Impact Analysis of the World Trade
Center Tower
We performed an analysis of the aircraft impact witha full-scale model of WTC Tower 2 [Isobe and Sasaki,2007] using the code described in the previous
section. Details on the structural members andconstruction data are extracted and adopted fromprevious reports ([ASCE/FEMA, 2002] and [NISTNCSTAR1, 2005]). The tower and the B767-200ER
aircraft are both modeled using linear Timoshenko
beam elements, of which all members are subdivided
into two elements. Figure 3 shows a schematic viewof 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 theaircraft model contains 4,322 elements, 2,970 nodes
and 17,820 degrees of freedom. Figure 4 shows a(a) 0.12 s(b) 0.28 s (c) 0.56 s after impact
Figure 5 Moment of impact (global view)(a) 0.12 s(c) 0.56 s after impact(b) 0.28 s
Figure 6 Motion of fuselage and engines (upperview)Table 1 Impact timeline of the B767-200ER aircraft
B767-200ER aircraft
Elapsedtime (s)PhaseLeft engine

0.0
A

Right engine

Snout hits outer-perimeter wall
Hits core
structure

0.048
Hits
outer-perimeter
wall

0.072

B

0.076

C

0.128

D

0.168

E

0.504

F

0.512

Snout

(a) Numerical result

Hits
outer-perimeter
wall
Hits core
structureHits
corecolumn No.902
Hits
core
column No.903
Hits east sideouter-perimeter
wall

280

●:Observed data
[ASCE/FEMA, 2002]

240(b) Observed data [NIST NCSTAR1, 2005]

Speed

200

Figure 8 Damage of the WTC 2 south face

160
120
80
40
A B

C

D

EF

0
0

0.1

0.2

0.3
0.4
Elapsed time

0.5

0.6

0.7

Figure 7 Velocity curve of the right engine
engine is assumed to be 19.315 t. The nose of theaircraft is tilted 11.5 degrees to the east and 5 degrees
downward, and its left wing is inclined downward by
35 degrees. It is assumed to collide with the 81stfloor on the south face of WTC 2 with a cruising
speed of 590 mph (262 m/s).
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 [ASCE/FEMA, 2002]. The beam joints of main
double trusses also had simple connections to sustain
static, vertical dead loads, which satisfied the initialstructural design. The structural design of member
joints in Japan, for example, demands the capacity to
be greater than that of the member itself. In thepresent analysis, the vulnerability of member joints
are considered by introducing fracture conditions asfollows [Lynn and Isobe, 2007]:
å
å
å
å
í
ì
å îx åå îy å
èz
å
å
å
å
Ä1 ï0
åîf x å Ä 1 ï 0 or åîf y å Ä 1 ï 0 or
èf z

(2)

where κx, κy, εz, κfx, κfy and εfz are the bending strains

around the x- and y-axes, the axial tensile strain andthe critical values for these three strains, respectively.
In the present analysis, we fixed the critical strain
values by comparing the impact damage and motionof the engines with the observed data. We also
considered the strain rate effect on yield strength
using the Cowper-Symonds equation [Jones, 1989].
An implicit scheme with a consistent mass matrix is
adopted. No damping except for the algorithmic one
is applied. δ is set to 5/6 and β to 4/9 in Newmark’s βmethod [Press et al, 1992]. The time increment ∆t is
set to 0.2 ms. An updated Lagrangian formulation isused since large deformations are expected in theanalysis, and the conjugate gradient (CG) method isused as the solver to reduce the memory requirement.
The analytical code is run on a high-performancecomputer (1.4 GHz Itanium*2, 8GB RAM), and thecalculation takes approximately two months for aphysical time of 0.7 s.
The numerical results are shown in Figs. 5 and 6.
Fractured elements are deleted in the figures. It isobserved that 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 passes through the north-east corner of
the tower. Table 1 shows the impact timeline of thefuselage and both engines. The velocity curve of theright engine is shown in Fig. 7 and a comparison of
the damaged surfaces between the numerical result
and the observed data [NIST NCSTAR1, 2005] isshown in Fig. 8. Both results are in good agreementwith the observed data. Figure 9 shows the transitionof axial forces in typical core columns every ten

AB

C

D

EF

AB

C

D

EF

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

Displacement (mm)

200

100

0
11F
21F
31F
41F
51F
61F
71F
81F
91F
101F

-100

-200

-300
0

0.1

0.2

0.3
0.4
Elapsed time (s)

0.5

0.6

0.7

Figure 10 Axial displacement due to spring-backphenomena (core column No. 1001)

Figure 11 Distribution of fractured main double
trussesstories from the ground floor to the top floor. The
core columns are compressed constantly until thesnout of the aircraft reaches the core structure (Phase
A in the figure), after which a wave due to the impact
and member fracture propagates in the horizontal andvertical directions. In particular, the compression
decreases instantaneously in the fractured corecolumn (No. 1001) at the floors above the impact
point. At the lower levels, compression changes totension immediately after the left engine hits the core
structure (Phase C), and the values of axial forces
continue to vibrate with large amplitude. Note thatthe lower the height, the larger the amplitude
becomes. The core column No.1001 at the 60th floor,
for example, moves vertically for 25 cm in 0.2 s as aresult of this dynamic transition (see Fig. 10).
What actually caused this dynamic unloadinglarge enough to generate such large tensile forces?
We extracted the time history of the distribution of
fractured main double trusses, as shown in Fig. 11, toinvestigate the cause. The propagation of fractured
trusses from the impact point to upper and lower
levels can be observed in the figure. This propagation
speed almost matched the longitudinal wave speed in
the columns. It can be assumed that the shock wave
due to impact propagated through the tower andoccasionally reflected to other columns at the groundor top levels, destroying each member joint as it
moved on.
Our assumed scenario is as follows; the aircraft
impact destroyed several main double trusses, that
supported the floor slabs, and led to a middle-classunloading of core columns. The deformation due tothe unloading caused other trusses to disconnect,
which released further burden from the core columns.
This disconnect-and-release process advanced rapidly
in a chain reaction manner and a spring-backphenomenon might have occurred in many columns
in such a way that a compressed spring is released ina very short moment. It was fatal enough to produce agigantic tensile force on the columns, which might
have triggered the total fracture of some column
joints, eventually leading to the total collapse of the

18m

18m

(a) Elastic modulus
Blockson fire100m(b) Yield strength
Figure 13 Strength reduction curves of steel due
to elevated temperatures

Figure 12 High-rise tower model
towers. It can be assumed that the damage in lower
levels was larger, since the lower the level went, thelarger the tensile force became.
3.2 Fire-Induced Collapse Analyses of High-Rise
TowersSome analyses of fire-induced collapse are carriedout on a 25-story 3-span tower as shown in Fig. 12 todetermine the effect of the strength reduction due toelevated temperatures. We used a model with twocases of column thickness; thick type (□-600×600
×18×18 at the 1st floor) and thin type (□-430×430
×13×13 at the 1st floor), two cases of member joint
strengths; strong type (CM=1.0) and weak type
(CM=0.2), and three cases of fire patterns; 17th ~ 19th
floors with 3 blocks (symmetry), 4 blocks(asymmetry) and 9 blocks (all floor) on fire. The joint
strength ratio CM indicates the ratio of member joint
strength against the member strength itself. The value
is directly used in the yield function as follows, toexplicitly express the influence of member joint
strength:
f =

í

Mx
CM M x0

ì2 í+

MyCM M y0

ì2 í
ì2 í
ì2
N
Mz
+
+
Ä1 = 0
N0
M z0

(3)

The ratio value is assumed on the basis of the factthat the member joint strengths of typical corecolumns were about 20 to 30 % of those of themembers in the WTC towers [ASCE/FEMA, 2002].
Here, the fracture condition of member joints are
considered by examining bending strains, axial
tensile strain and shear strains in a member, as
follows:

å
å
å
å
ì
í
å îx åååèzåå Ä 1 ï 0 or å îy å Ä 1 ï 0 or
Ä1 ï0åîf x ååîf y å
èf z
å
å
å
å
å
å çxz åå
å Ä 1 ï 0 or å çyz å Ä 1 ï 0
or åå
(4)
å
å
çf xz
çf yz åwhere γxｚ, γyｚ, γfxｚ, γfyｚare the shear strains for thex- and y-axes and the critical values for the strains,respectively. The critical strain values used here are
the values actually obtained from some experiments
concerning the high strength joint bolts [Hirashima etal, 2007]. Curves showing the reduced elastic
modulus and yield strength of steel due to elevatedtemperatures [NIST NCSTAR1, 2005], as shown in
Figs. 13(a) and 13(b), are adopted. The elastic
modulus of the columns with no thermal insulation is
reduced to 60 % of the original value, and the yieldstrength to 10 % of the original strength, at 700 ℃,
which is a natural temperature in normal fire. We
assumed the temperature to rise up to 700 ℃ in 7

Figure 14 Assumed time history of temperature
minutes as shown in Fig. 14. Also, we applied a timeincrement control to enable continuous calculationfrom static to dynamic phenomena.
Two of the numerical results are shown in Figs.
15 and 16. Both are performed with the thin-column
type model. Figure 15 shows the result with 3-floors
4-blocks on fire and Fig. 16 with 3-floors 9-blocks on
fire, each performed on strong-member-joint type
(CM=1.0) and weak-member-joint type (CM=0.2). We(a) CM =1.0 (final phase)

can clearly see the difference in collapse modes
between the two models in both figures. The
strong-member-joint models initiate to collapse once
they buckle at the floors that are on fire, but
withstand to collapse any further by redistributing thestresses in the towers. The weak-member-joint
models, on the other hand, fail to stop the progressionof collapse and end in a total destruction. Table 2
shows the sum of all the results. It can be concluded
that the possibility of total collapse becomes high, inspite of thickness of the columns, when the joint
strength ratio is low. Fire patterns do have influenceon the possibility of collapse. Multi-floor andasymmetric fire patterns may increase the possibilityof collapse. And furthermore, the collapse speeds
were much slower than the free-fall or saturated
speed in these cases, which does not explain thehigh-speed collapse of the WTC towers.
4. CONCLUSION
The numerical results shown in this paper might
indicate that the highly-rapid total collapse of the(b) CM =0.2 (progressive collapse sequence)

Figure 15 Fire-induced collapse analysis with 3F-4 blocks fire pattern (thin column model)(a) CM =1.0 (final phase)(b) CM =0.2 (progressive collapse sequence)

Figure 16 Fire-induced collapse analysis with 3F-9 blocks fire pattern (thin column model)

Table 2 Collapse modes occurred in the fire-induced collapse analyses

Columns□-600×600×
18×18→□-400×400×
12×12□-430×430×
13×13→□-230×230×7×7

Joint strength
ratio CM3F 3 blk. fire

3F 4 blk. fire

3F 9 blk. fire

1.0

No collapse

Total collapse（496s）

Partial collapse（421s）

0.2

Total collapse（886s）

Total collapse（463s）

Total collapse（423s）

1.0

No collapse

Partial collapse Partial collapse（433s）（417s）

0.2

Total collapse（516s）

Total collapse（420s）

Total collapse（412s）( ): Collapse initiation time
WTC towers was caused not only by the bucklingand strength reduction of members due to elevatedtemperatures, but by the original weakness of themember joints and the destruction of those during thecollision of the aircraft. Further investigations byperforming full-model fire-induced collapse analysesof the WTC towers are scheduled.
5. REFERENCES
ASCE/FEMA (2002), World Trade Center Building
Performance Study: Data Collection, Preliminary
Observation and Recommendations.
Hirashima, T., Hamada, N., Ozaki, F., Ave, T. andUesugi, H. (2007), “Experimental Study on Shear
Deformation Behavior of High Strength Bolts at
Elevated Temperature”, J. Struct. Constr. Eng., AIJ,
No.621, pp. 175-180, in Japanese.
Isobe, D. and Sasaki, Z. (2007), “Aircraft Impact
Analyses of the World Trade Center Towers”,
CD-ROM Proc. of 1st Int. Workshop on
Performance, Protection, and Strengthening ofStructures
under
ExtremeLoading(PROTECT2007), Whistler, Canada.

Jones, N. (1989), Structural impact, New York:
Cambridge University Press.
Lynn, K.M. and Isobe, D. (2007), “Finite Element
Code for Impact Collapse Problems of Framed
Structures”, Int. J. Numer. Meth. Eng., Vol. 69, No.
12, pp. 2538-2563.
NIST NCSTAR 1 (2005), Federal Building and Fire
Safety Investigation of the World Trade Center
Disaster: Final Report on the Collapse of theWorld Trade Center Towers, National Institute ofStandards and Technology (NIST).
Press, W.H., Teukolsky, S.A., Vetterling, W.T. and
Flannery, B.P. (1992), Numerical Recipes in
FORTRAN: The Art of Scientific Computing, New
York: Cambridge University Press.
Toi, Y. (1991), “Shifted Integration Technique in
One-Dimensional Plastic Collapse Analysis Using
Linear and Cubic Finite Elements”, Int. J. Numer.
Meth. Eng., Vol. 31, pp. 1537-1552.
Toi, Y. and Isobe, D. (1993), “Adaptively Shifted
Integration Technique for Finite Element Collapse
Analysis of Framed Structures”, Int. J. Numer.
Meth. Eng., Vol. 36, pp. 2323-2339.