Journal of Natural Disaster Science, Volume 24, Number 1, 2002, pp25-34

Applied Element Method Used for Large Displacement Structural Analysis

1)

Kimiro MEGURO1) and Hatem Sayed TAGEL-DIN2)
Associate Professor, International Center for Urban-Safety Engineering (ICUS/INCEDE),
Institute of Industrial Science, The University of Tokyo
(4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan)
2)
Assistant Professor, Structural Engineering Department, Faculty of Engineering,
Cairo University, Giza, Egypt(Received for 19. Jun., 2002)

ABSTRACTA new extension of the Applied Element Method (AEM) for structural analysis is introduced. A brief overview
of the method’s formulation is presented. Then, modifications needed to analyze the behavior of structures subjected to large displacements under static loading are introduced. As no geometric stiffness matrix is needed, theformulation is simple, general and applicable to any type of structural configuration or material. A series ofexamples that verify the applicability of the proposed technique are presented. The AEM is shown to be an efficient tool for structural analysis in both the small and large displacement ranges.

1. INTRODUCTION
Studies of structure collapse behavior are very important inorder to reduce the number of casualties that occur during earthquakes. From the numerical point of view, structures undergo
small and large displacements before collapse. In the small displacement range, the variation of the structure geometry duringloading can be neglected. The buckling of columns, however, amain reason of structural failure, can only be detected by large displacement theory.
The Finite Element Method (FEM) has been considered themain tool for studying structural buckling behavior. By its use, thebuckling mode and load, as well as post-buckling behavior, can befollowed (Waszczyszyn et al. 1994; Szabo et al. 1986).
Unfortunately, the FEM assumes that a material is a continuum,
and special techniques must be adopted to consider the separationof structural members. In most cases, the fracture plane is arbitrary and unknown before analysis. Modeling of the separation
using joint elements is reliable only when the crack location can bepredicted. This condition is common to numerical methods thatconsider the structure as a continuum.
A recent method that deals with structure failure analysis is
the Modified or Extended Distinct Element Method, MDEM or
EDEM (Meguro and Hakuno, 1989, 1994). It can follow highly
nonlinear geometric changes in the structure during failure. Its
main disadvantage when compared to the FEM, however, is that it
is less accurate in the small displacement range and requires relatively longer CPU time. The applicability of the EDEM or othermethods that adopt rigid elements, such as the Rigid Body and
Spring Model, RBSM (Kawai, 1980; Kikuchi et al, 1992), has notbeen verified for buckling analysis. Discontinuous Deformation

Analysis (Amadei et al, 1996; Sitar et al., 1997) has been used forthe large deformation of rock blocks, but it too is limited to analysis of elements that are separated from the beginning of analysis.
It is not verified in cases in which deformation is small.
Our proposed method, the Applied Element Method (AEM),
is based on division of the structural members into virtual elements
connected through springs. Each spring entirely represents thestresses, strains, deformations, and failure of a certain portion of
the structure. The main advantage of this method is that it can follow structural behavior from the initial loading stages until complete collapse with reliable accuracy in reasonable CPU time. The
applicability and accuracy of the AEM in various fields have been
discussed elsewhere (Tagel-Din and Meguro, 1998, 1999, 2000;
Meguro and Tagel-Din, 1997,1998).
The current areas of AEM application are given in Table 1.
Numerical results were compared with theoretical and experimental results whenever possible for each field of application. TagelDin and Meguro (1998) and Meguro and Tagel-Din (1997, 1998)
verified the accuracy of the AEM in small deformation ranges.
Effects of Poisson’s ratio, normally neglected in methods based on
rigid body elements, also are successfully handled (Tagel-Din and
Meguro, 1998). The AEM can track complicated nonlinear behavior such as crack initiation, propagation, and opening and closure,
as well as estimate failure loads (Meguro and Tagel-Din, 1997,1998).
To follow structural behavior until complete collapse, accuracy of the method in the large displacement range must be verified
quantitatively. We here present the AEM basic formulation anddiscuss the adjustments needed for analysis of structures subjectedto large displacements. Unlike other methods, in the AEM there isno need to determine the geometrical stiffness matrix. This makes

KEY WORDS: Applied Element Method, buckling, post-buckling behavior, large displacement, computer simulation

26

K. MEGURO AND H. S. TAGEL-DIN

Table 1.

Fig. 1

Areas of Application of AEM

Modeling of the structure in the AEM

the method general and applicable to different loading and structure types.
Its accuracy in the large displacement range is verified bycomparing the numerical results with theoretical ones for the buckling and post-buckling stages; agreement is very good. The aim of
the proposed method is to simulate the behavior of a structure fromthe stage in which no load is applied until complete collapse.
Results presented here show the applicability of the method tolarge deformation analysis. Note that this research is one of thefirst attempts to utilize rigid body elements with three degrees of
freedom for the study of structure buckling behavior.
2.

APPLIED ELEMENT METHOD FORMULATIONIn the AEM, the structure is divided in small elements asshown in Fig. 1. Two elements are assumed to be connected at
discrete points along their edges by a pair of normal and shear
springs. Spring stiffness is determined by
GdT
Kn = EdT
a and Ks = a(1)where d is the distance between springs, T the element thickness, a
the length of the representative area, E the material Young’s modulus, and G the material shear modulus. Equation (1) implies thateach spring represents the stiffness of a volume with the dimensions d, T, and a. When reinforcement is present, rebar stiffness is
added to the material stiffness found with Equation (1). In thetwo-dimensional model, each element has three degrees of freedom (DOF), which represent the element’s rigid body motion.
Although elements move as rigid bodies (element shape does notchange), internal deformations are represented by spring deformations (the element assembly is deformable). Tagel-Din and
Meguro (1998) introduced the Poisson’s ratio effect to the AEM.
Stiffness matrix components are determined by applying aunitary displacement to a DOF while keeping the remaining DOF’s
fixed. The forces needed to generate this configuration are thestiffness matrix components, which are equal to the summation of
the contributions of the springs surrounding the element. The contribution of the contact spring shown in Fig. 2 to DOFs u1, u2, andu3 isSin2(θ+α)Kn+Cos2(θ+α)Ks−KnSin(θ+α)Cos(θ+α)+KsSin(θ+α)Cos(θ+α)
Cos(θ+α)KsLSin(α)−Sin(θ+α)KnLCos(α)−KnSin(θ+α)Cos(θ+α) Cos(θ+α)KsLSin(α)+KsSin(θ+α)Cos(θ+α) −Sin(θ+α)KnLCos(α)
Sin2(θ+α)Ks
Cos(θ+α)KnLCos(α)+Cos2(θ+α)Kn+Sin(θ+α)KsLSin(α)
Cos(θ+α)KnLCos(α)
L2Cos2(α)Kn+Sin(θ+α)KsLSin(α)+L2Sin2(α)Ks(2)

where all the terms are illustrated in Fig. 2. Equation (2) showsAPPLIED ELEMENT METHOD USED FOR LARGE DISPLACEMENT STRUCTURAL ANALYSIS

27

Fig. 2 Element shape, contact point, and degrees of freedom

one-quarter of the element stiffness matrix. The global stiffness
matrix, K, is found by summing up the contributions of all thesprings in the system.
In spite of its relatively simple formulation, the AEM has beenused successfully for the nonlinear analysis of RC structures subjected to small deformations (Meguro and Tagel-Din, 1998). It
cannot, however, be used for large deformation analysis unlessgeometrical changes in the structure’s shape are considered.
3. LARGE DISPLACEMENT ANALYSIS WITH AEM
3.1 Numerical Procedure
Whereas in the FEM the effects of large displacements are
considered by adopting a geometrical stiffness matrix(Waszczyszyn et al., 1994), in the AEM there is no need for such amatrix. This makes the method more general and applicable to any
type of loading or structure. One limitation of the AEM is that thedirection of the applied forces is assumed to be constant. Becauseof this, loading conditions (Waszczyszyn et al., 1994) in which theforce direction changes, as present when a member buckles, cannot
be analyzed by the AEM.
To adapt the formulation in the previous section for static
large deformation analysis, the following modification is introduced:
ΔU = Δf + Rm + RG
KΔ

(3)

where K is the nonlinear stiffness matrix, ΔU the incremental displacement vector, Δf the incremental load vector, Rm the residual
force vector due to cracking or incompatibility between springstrains and stresses, and RG the residual force vector due to geometrical changes in the structure during loading.
Application of the AEM:
1. Assume that Rm and RG are null and solve Equation (3) to getΔU.
2. Modify the structural geometry according to the calculated
incremental displacements.
3. Modify the direction of the spring force vectors according to the

new element configuration. The geometrical changes generate
incompatibility between the applied forces and internal stresses.
4. Verify whether cracking occurred and calculate Rm . In elastic
analysis Rm is zero.
5. Calculate the element force vector, Fm , by summing the forces ofthe springs around each element.
6. Calculate the geometrical residuals around each element with
Equation (4).
RG = f - Fm

(4)where f is the applied force vector. Equation (4) implies thatgeometrical residuals account for the incompatibility between
the external applied and internal forces due to modification of
the structure’s geometry.
7. Small deformations are assumed during each increment.
8. Calculate the stiffness matrix for the structure with the new configuration considering stiffness changes due to cracking or
yielding.
9. Repeat the entire process.
Residuals calculated in the previous increment can be incorporated in the solution of Equation (3), which reduces the calculation time.
Although the technique presented is simple, the numericalresults show high accuracy. The following limitations were noted:
1. Complete symmetry of the structure and loading must be avoided in buckling analysis. The symmetry can be broken by slightchanges in the material parameters of a part of the structure.
2. Small displacement theory is assumed during each increment.
The load or displacement increments therefore should be small.
3. In many applications, the apparent structure stiffness decreases
after buckling; hence, the applied load also should decrease. In
load control analysis with a positive load increment, the cumulative difference between the applied and internal loads appears
after buckling. Applying additional load increments after buckling results in the divergence of geometrical residuals, produc-

28

K. MEGURO AND H. S. TAGEL-DIN

ing large geometrical changes within a few increments.
Because of this, in some cases it is impossible to obtain a solution. Displacement control analysis is suggested to overcome
this problem, but that method is limited to cases in which there
are few prescribed displacements. For a more general solution,
the proposed technique can be modified to adopt the energy or
arc length method (Kleiber, 1989).
3.2 Simulation of a simply supported rubber beam
By the small displacement theory
This example is introduced to show why the AEM cannot deal
with large displacement analysis unless the modifications listed in
section 3.1 are executed. The beam shown in Fig. 3, span 12.0 m,
and square cross section side equal to 1.0 m, is simply supported(pinned support at one side, roller at the other) and subjected to aconcentrated load at its center. The material is elastic. Young’s
modulus is 210 MPa. The following are seen from the deformed
shape shown in Fig. 3:
1. Beam deformation is acceptable in the small displacement
range.
2. The volume of the beam increases when analysis continues up tovery large deformations.
3. The roller does not move even for large deformations. This is
not realistic.
By the large displacement theory
The same beam under the same loading was analyzed by large
displacement formulation to check the accuracy of the proposed

Fig. 3 Deformed shape of a simple supported beam (Small displacement analysis)Fig. 4 Deformed shape of a simple supported beam (Large displacement analysis)

APPLIED ELEMENT METHOD USED FOR LARGE DISPLACEMENT STRUCTURAL ANALYSIS

adjustments. In that case, the realistic results shown in Fig. 4 wereobtained. Load displacement relationships before and after modification are shown in Fig. 5. DY and DY\ respectively represent displacements at the beam center when small and large displacement
analyses are performed. Similarly DX and DX\ represent roller
displacement in those cases. The following are seen from theresults:
1. In the small displacement range, DY and DY\ are similar.
2. DX\ increases as the applied load increases, whereas DX is29

almost zero. Fig. 5 shows that DX\ is negative at the beginning.
This is because of the rotation of the beam’s cross section about
the neutral axis.
3. During loading, the shape of the beam changes from a straightline to an arch. This shows that beam stiffness increases whengeometrical changes are considered.
3.3 Simulation of buckling and post-buckling behavior
Fixed base column
The first case study of buckling and post-buckling behavior isFig. 5 Load versus displacement of a simple supported beam analyzed by considering or not geometrical residuals

Fig. 6 Post buckling behavior of an elastic column

30

K. MEGURO AND H. S. TAGEL-DIN

that of a fixed base elastic column subjected to axial load. The
load direction is assumed to be constant during the analysis.
Column height is 12.0 m, and the side of the square cross section is
1.0 m. Young’s modulus is 840 MPa. Analysis is performed for
300 elements, and vertical displacement is applied to the top of thecolumn at a constant rate. To break the system’s symmetry, thestiffness of the left edge element at mid-height is increased by 1%with respect to the other elements. Figure. 6 shows the deformed
shape of the column before and after buckling. Obviously, thehighly nonlinear geometrical changes have been successfully followed.
Figure. 7 shows the horizontal and vertical displacements at
the loading point obtained with the AEM, with and without considering geometrical residuals, and the theoretical load-displacement
relationships (Timoshenko and Gere, 1961). In the theoreticalanalysis, effects of axial and shear deformations are neglected
whereas in AEM analysis these effects, although being relatively
small, are taken into account.
Figures. 6 and 7 indicate the following:

1. Load displacement relationship obtained when geometricalresiduals are considered is close to the theoretical values, evenunder very large displacements. This shows that the proposed
method is accurate and numerically stable.
2. The buckling load calculated using only geometrical modifications, not geometrical residuals, is approximately 470 KN, larger than the theoretical value, 78 KN. This shows that considering only geometrical modifications is not sufficient to obtainaccurate results.
3. The calculated load displacement relation is tangential to thehorizontal line at the buckling load. This is in good agreementwith theoretical results.
4. A slight increase in the load after buckling results in very large
displacements. This indicates that applying the load control
technique after buckling produces very large deformations in
only a few increments.
5. When vertical displacement is about 9 m, horizontal displacement starts to decrease.
6. The column shape changes after buckling, increasing the speci-

Fig. 7 Load displacement relationship of an elastic column subjected to vertical loading

Fig. 8 Load versus vertical normal stress at point AAPPLIED ELEMENT METHOD USED FOR LARGE DISPLACEMENT STRUCTURAL ANALYSIS

men’s stiffness.
Figure. 8 shows the relationship between the load and verticalnormal stress at point A under the load application point. Before
buckling, stress is mainly uniform compression and increases linearly. When the buckling load is reached, compression stresses
decrease until zero is reached, stage (8) in Fig. 6, when the load
direction is parallel to the column edge. After this, tension stresses
develop and increase.
Internal stress changes at the intermediate sections are shownin Fig. 9. Before buckling, the stresses are mainly uniform compression and only axial deformation occurs. After reaching thebuckling load, although the applied load is constant (P=78 KN),

Fig. 9Fig. 1031

buckling bending moments appear and cause large deformation.
The strong point of the AEM is that it accurately follows the
behavior of any point in the structure even if large deformations
occur.
Snap through buckling of a two-member truss
The second case study is that of simulation of the buckling
behavior of a two-member truss (Szabo et al., 1986). Figure. 10shows the truss layout and loading point location. The material
Young’s modulus is 210 MPa. Due to structural symmetry, only
half the truss is analyzed using 0.1m-side square elements. The
load was applied at a constant rate.

Variation in the internal stress distribution during buckling

Relationships between applied load, member force, and the displacement ratio for vertical displacement applied at constant rate

32

K. MEGURO AND H. S. TAGEL-DIN

Figure. 10 shows the relationship between the applied load,
member force, and displacement ratio. When results are comparedwith theoretical ones, there is almost no difference. It took only
one minute to perform the AEM analysis with a personal computer(CPU Pentium 267 MHz).
The truss passes through the following deformation stages:
1. At first, member length decreases and the compression force
increases. The shortest member length, which corresponds tothe maximum compression stresses, is obtained when truss
members are horizontal.
2. At this instant, the compression force is maximum, and theapplied load is zero in the horizontal direction.
3. Increasing the displacements after the members become horizontal increases the member length; hence, the compression
force is released. The direction of the applied load is reversed.

Fig. 11

Fig. 12

4. When (d/H) (refer to Fig. 10) equals 2.0, the final member
length is the same as the initial length. The member force andapplied load therefore become zero.
5. Increasing the applied displacement increases the tension force
in the members.
Elastic Frame under Different Support Conditions
The third case study is the simulation of the buckling of anelastic frame. The following cases are solved:
1. Sway frame with fixed supports
2. Nonsway frame with fixed supports
3. Sway frame with hinged supports
Two vertical loads are applied to the tops of the columns. The
frame dimensions and location of the loading points are shown in
Fig. 11. The frame has a 0.5-m side square cross section. Young’s

Load displacement relationship for a sway frame with fixed supports subjected to vertical loads

Load displacement relationship for a nonsway frame with fixed supports subjected to vertical loads calculated by the AEM

APPLIED ELEMENT METHOD USED FOR LARGE DISPLACEMENT STRUCTURAL ANALYSIS

Fig. 13

33

Load displacement relationship for a sway frame with hinged supports subjected to vertical loads

modulus is 210 MPa. The load is applied at a constant rate. To
break the system’s symmetry, the stiffness of one of the edge elements is increased by 1% with respect to the other elements. This
analysis cannot be performed under displacement control because,
after buckling, the displacements of the frame corners differ due tochange in axial force in the columns. To carry out the analysis
under load control, it is necessary to use very small load increments after buckling.
Results of the analysis for the first case, 136 elements, areshown in Fig. 11. The buckling load obtained is very close to thetheoretical one (Timoshenko and Gere, 1961). After buckling, displacements markedly increase in a few increments because loading
is applied under load control.
Figure. 12 shows the results for the second case. Due to symmetry, the analysis is performed for only half of the frame. As theload is applied at only one point, analysis can be made under displacement control up to large displacements without loosing stability. The buckling load is higher than in the previous case. The
theoretical buckling load is not available for this case. The
obtained buckling mode, shown in Fig. 12, however, seems realistic. The analysis was stopped before recontact of the elements.
Tagel-Din and Meguro (1999, 2000) discussed the recontact issue
in detail.
Figure. 13 shows the results for the third case. The analysis of
306 elements was performed under load control. The bucklingload obtained is very close to the theoretical one (Timoshenko andGere, 1961) but smaller than that for a frame with fixed supports.
4. CONCLUSIONSWe present a new extension of the AEM for the analysis of
structures subjected to large displacements. Comparison of theresults obtained by the AEM with those obtained by theoreticalformulations, showed the applicability of the proposed techniques.
The advantages of AEM for large displacement analysis are:1. It is relatively simple compared to other numerical techniques.
2. It accurately follows structural behavior even in the range oflarge displacements in which large geometrical changes occur.
The simulated buckling loads, modes, and internal stresses
agree well with those found theoretically.
3. It is general and applicable to any type of structure or material.
4. It is easily extended to follow the large displacement of structures until total collapse, as shown by Tagel-Din and Meguro(1999, 2000).
The following limitations also were identified:
1. Load direction is constant. Subsequent loading conditions(Waszczyszyn et al 1994) and non-conservative loads cannot be
studied by use of the proposed formulation.
2. Although the load can be applied by either load or displacement
control, both have limitations. The load control technique cannot follow post peak behavior, and the displacement control
technique cannot follow cases when the tangent to the load
deformation curve tends to be vertical (Kleiber, 1989). In addition, the displacement control technique cannot be adopted forcases in which a load is applied at many points. In spite of this,the AEM method can be extended to follow other loading methods such as the energy control or arc length control methods(Kleiber, 1989).
REFERENCES
Amadei B., Lin C., and Dwyer J., 1996. Recent extensions to the DDA
method, Proc. of 1 st Int. Forum on Discontinuous Deformation
Analysis (DDA), Berkley, California, Ed. Salami & Banks.
Kawai T., 1980. Some considerations on the finite element method, Int. J.
for Numerical Methods in Engineering, 16, 81-120.
Kikuchi A., Kawai T., Suzuki N., 1992. The rigid bodies-spring models
and their applications to three dimensional crack problems, Comp.
and Struct., 44(1), 469-480.
Kleiber M., 1989. Incremental finite element modeling in non-linear solid

34

K. MEGURO AND H. S. TAGEL-DIN

mechanics, Ellis Horwood, New York.
Meguro K. and Hakuno M., 1989. Fracture analyses of structures by themodified distinct element method, Struct. Engrg./Earthquake Engrg.,
6(2), 283s-294s.
Meguro K. and Hakuno M., 1994. Application of the extended distinct element method for collapse simulation of a double-deck bridge, Struct.
Engrg./Earthquake Engrg., 10(4), 175s-185s.
Meguro, K. and Tagel-Din, H., 1997. Development of a new fracture analysis method with high accuracy based on discontinuous material modeling, Proc. of 16th Annual Conf. on Natural Disaster Reduction,
Osaka, Japan.
Meguro, K. and Tagel-Din, H., 1998. A new simplified and efficient technique for fracture behavior analysis of concrete structures, Proc. of
3rd Int. Conf. on Fracture Mech. of Concrete and Concrete Struct.
(FRAMCOS-3), Gifu, Japan.
Nicholas S., and Mary M. MacLaughlin, 1997. Kinematics andDiscontinuous Deformation Analysis of Landslide Movement, Invited
Keynote Lecture, 11 Panamerican Symposium on Landslides, Rio deJaneiro, Nov. 10-14th.
Ohno, H., Adachi, H., Nakanishi, M., Shimizu, Y., and Sugiyama, K.,
1976. Study on seismic capacity of reinforced concrete shear wall,
Part 7, Relation between load history and horizontal reinforcement,

Proc. of Annu. Conf. of Arch. Inst. of Japan, 1601-1602. (in Japanese)
Okamura H. and Maekawa K., 1991. Nonlinear analysis and constitutive
models of reinforced concrete, Gihodo Co. Ltd., Tokyo.
Szabo J., Gaspar Z., Tarnai T., 1986. Post-buckling of elastic structures,
Elsevier Science Publishing Co. Inc., Budapest.
Tagel-Din, H. and Meguro, K., 1998. Consideration of Poisson’s ratio
effect in structural analysis using elements with three degrees of freedom, Bull. of Earthquake Resistant Struct. Res. Ctr., University ofTokyo, 31, 41-50.
Tagel-Din, H. and Meguro, K., 1999. Applied element simulation for collapse analysis of structures, Bull. of Earthquake Resistant Struct. Res.
Ctr., University of Tokyo, 32, 113-123.
Tagel-Din, H. and Meguro, K., 2000. Analysis of a small scale RC building
subjected to shaking table tests using the applied element method,
Proc. of 12th World Conf. on Earthquake Engrg., Auckland, New
Zealand.
Timoshenko S. and Gere J., 1961.Theory of elastic stability, McGraw-Hill
Inc.
Waszczyszyn Z., Cichon C., and Radwanska M., 1994. Stability of structures by finite element methods, Elsevier Science Publishing Co. Inc.,
Amsterdam.