A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering

Tome XVI [2018] | Fascicule 4 [November]

1.
Mariano T. BOCOVO, 1.Ernesto C. HOUEHANOU, 1.Gérard GBAGUIDI AÏSSE, 1.Gérard DEGAN

STIFFNESS MATRIX METHOD FOR NONLINEAR ANALYSIS OF

PLANE FRAMES
Laboratory of Applied Energetics and Mechanics (LEMA), University of Abomey-Calavi, BENIN
1.

Abstract: In this paper, geometric nonlinear analysis of plane frames was performed by the stiffness matrix method
using stability functions. At first, the argument of the stability functions (υ) was set as 0.01. The stiffness matrix of the
frame, [K], has been assembled, as well as the nodal load vector of the frame, {F}. The boundary conditions (support
restraint and windbracing restraint) were introduced for the reduction of this matrix and the nodal load vector. At
this stage, the determinant of the reduced stiffness matrix, [Kr], and the reduced nodal displacement vector, {Dr}, are
calculated. The argument of the stability functions is incremented by 0.01 and the operations are repeated until the
determinant of the reduced stiffness matrix, │Kr│, changes sign. The argument υ of the iteration preceding the sign
change of │Kr│ and corresponding to the positive value of the determinant is taken and refined by a process
described in the paper. The buckling loads of the frame members are determined at this stage. Windbracing and
the increase of supports stiffness increase the value of the critical load (less sensitivity to the phenomena of elastic
instability) and have been identified as factors of stability.
Keywords: geometric nonlinear analysis, stiffness matrix method, stability functions, determinant, buckling load

1. INTRODUCTION
The geometric nonlinear analysis of the plane frames will be made by the stiffness matrix method using stability
functions. Geometric nonlinear analysis of frames is a problem of elastic instability. Although the structure is in
the elastic range, the effects of large displacements (large deformations) modify the geometry of the structure;
which renders inapplicable the theory of linear elasticity. Examples of these phenomena of instability are
buckling, local buckling and lateral buckling [1].
Stability analysis in solid mechanics began with Euler's solution of buckling of an elastic column (Euler, 1744).
Most basic linear elastic problems of structural stability were solved by the end of the 19th century, although
further solutions have been appearing as new structural types were being introduced. The twentieth century
has witnessed a great expansion of the stability theory into nonlinear behavior, caused either by large
deflections or by nonlinearity of the constitutive law of the material [2].
Several authors have studied stiffness matrix method and nonlinear analysis of structures. Some of these studies
include: Sezer (1995) did some works on nonlinear analysis of reinforced concrete plates by finite element
method [3]; Gummandi and Palazotto (1997) performed nonlinear analysis of beams and arches by finite
element method [4]; Bhatt (1999) worked on the elastic stability of structural systems by stiffness matrix method
[5]; Ozer (2006) used various methods for nonlinear analysis of structural systems [6]. Orumu (2013) has
proposed mathematical formulae for estimating the critical loads of rigid sway frames with pinned or fixed base
restraints [7].
The aim of the study is to find the critical loads of the frame members. The determinant of the reduced stiffness
matrix of the frame, as well as the argument of the stability functions are used to determine these critical loads.
A computer program is used to carry out the various operations.
2. GEOMETRY, LOADING AND SUPPORT CONDITIONS
This paper is devoted to calculating the critical load of rectangular plane frames made of straight members. The
posts of the frame, which are vertical, rise continuously from the foundation to the top of the structure. Similarly,
the beams (horizontal) are continuous from the left to the right of the structure. Thus, there is no missing bar in
an interior panel. The posts of the first storey can be of different lengths. The bars are connected to the nodes,
which are considered non-deformable, by rigid assemblies. The feet of the posts can be pinned or fixed in a
rigid manner to the foundation.
The external forces considered in the present study are such that, before the loss of stability, the members
undergo only axial compression (or traction). By way of example, Figure 1 shows such a system of forces. Here,
the loads are only concentrated loads applied to the nodes, acting in the axis of the posts.
A frame of a building is, in practice, subject to loads due to self-weight and overloads. The transformation of
dead weights and overloads acting on the horizontal members into nodal forces is done by a simple process.
Each beam is considered between two consecutive nodes as a beam on simple supports. The equivalent
concentrated loads acting at the nodes are equal and of opposite sign to the reactions of this simple beam.
When there are horizontal forces acting at the nodes, the axial forces in the members are found by a preliminary
calculation.

221 | F a s c i c u l e 4
A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

a)

b) c) d)
Figure 1: Categories of structures studied
3. ASSUMPTIONS
This study is based on the following assumptions:
— The materials that make up the bars are supposed to be perfectly elastic.
— The nodes are supposed to be rigid.
— Forces are expected to maintain their original direction and initial point of application during buckling.
— In addition, the external forces are such that, before the loss of stability, the bars undergo only an axial
force (compression or traction).
— Deformations (in the displaced state) are considered small.
— The case of buckling studied is that of bending buckling in the plane of the frame. It is assumed that the
spatial buckling accompanied by twisting and local buckling of the walls are prevented.
In addition, the elastic length variations
due to axial forces are neglected. All the
nodes of a crossbar thus have the same
transverse displacement.
4. STABILITY STUDY
Nonlinear analysis of structures can be
formulated using Figure 2 and Equation 1.
The stiffness matrix ([ki] in local coordinate
system) and that of rotation
transformation, [Ti], of a member << i >>
are modified during the increase of the
axial load. However, in this study, [Ti] is
assumed unaffected.
Figure 2: Displaced position and degrees of freedom (d.o.f.) of a member

{ fi } = [ ki ].{di } (1)
where:
 Pi   ui 
V  v 
 i   i
 Mi   θi 
={ fi } =  and {di }  
 Pj  u j 
V  v 
 j   j
M j  θ j 
   
{fi} is the nodal load vector of the member; Pi, Vi and Mi respectively represent axial force, shear force and
bending moment at node i.
{di} is the vector of nodal deformations of the member; ui, vi and θi respectively represent axial displacement ,
lateral displacement and rotation at node i.
The stability functions [8] of the members (equation 2) cannot be calculated because the axial forces acting on
them are unknown. An iterative method is used to overcome this difficulty. One starts by calculating these
functions using as load argument υ = 0.01.

222 | F a s c i c u l e 4
 A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

 AE AE 
 L 0 0 − 0 0 
L
 
 0 12 EI 6 EI 12 EI 6 EI
.η2 (υ ) .ϕ4 (υ ) 0 − .η 2 (υ ) .ϕ4 (υ ) 
 L3 L2 L3 L2 
 
 0 6 EI 4 EI 6 EI 2 EI
.ϕ 4 (υ ) .ϕ2 (υ ) 0 − .ϕ 4 (υ ) .ϕ3 (υ ) 
 L2 L L2 L 
[ ki ] =  AE AE

− 0 0 0 0 
 L L 
 
 0 12 EI 6 EI 12 EI 6 EI
− .η2 (υ ) − .ϕ4 (υ ) 0 .η2 (υ ) − .ϕ4 (υ ) 
 L3 L2 L3 L 2 
 
 0 6 EI 2 EI 6 EI 4 EI
.ϕ4 (υ ) .ϕ3 (υ ) 0 − .ϕ4 (υ ) .ϕ2 (υ ) 
 L2 L L2 L 
The following steps should be implemented for the nonlinear analysis of a frame:
Step 1: Idealize the structure and establish global axes,
Step 2: Number the nodes (A, B, C, D, ...) and the degrees of freedom (d.o.f.) of the structure (0 for inactive d.o.f.,
and 1, 2, 3, 4, 5, 6, 7, …, for the others),
Step 3: Number the members and assign an arrow to each member so that ends i and j are defined,
Step 4: Enter the geometric characteristics (area, inertia, length) of each member and its orientation angle, θ, as
well as the mechanical properties of the materials (YOUNG’s modulus, COULOMB’s modulus),
υ 2 tan υ
ϕ1 (υ ) =
3(tan υ − υ )
υ (tan υ − υ )
ϕ 2 (υ ) =
υ υ
8 tan υ (tan − )
2 2 (2)
υ (υ − sin υ )
ϕ3 (υ ) =
υ υ
4sin υ (tan − )
2 2
υ
ϕ 4 (υ ) = ϕ1 ( )
2
υ3
η1 (υ ) =
3(tan υ − υ )
υ
η 2 (υ ) = η1 ( )
2
Step 5: Initialize the argument υ to 0.01 and the iterations counter << Count. >> to 1,
Step 6: Enter the expressions of the stability functions and form the stiffness matrix, [ki], and the nodal load
vector, {fi}, of each member in the local coordinate system; for uncompressed members,
ϕ=
(υ ) η=
(υ ) 1 (3)
Step 7: Form the rotation transformation matrix, [Ti], of each member and calculate its stiffness matrix, [Ki], and
its nodal load vector, {Fi}, in the global coordinate system,
 cos θ sin θ 0 0 0 0
 − sin θ cos θ 0 0 0 0 

 0 0 1 0 0 0 (4)
[Ti ] =  0 0 0 cos θ sin θ 0


 0 0 0 − sin θ cos θ 0
 
 0 0 0 0 0 1 

[ Ki ] = [Ti ]T .[ ki ].[Ti ] (5)

Step 8: Assemble the stiffness matrix, [K], and the nodal load vector, {F}, of the structure,
Step 9: Form the reduced stiffness matrix, [Kr], and the reduced nodal load vector, {Fr}, of the structure by
ignoring the rows and columns corresponding to inactive d.o.f. (restraint of supports and restraint of
windbracing) and of zero order number << 0 >>,
Step 10: Calculate the reduced nodal deformation vector, {Dr}, of the structure:
{Dr } = [ K r ]−1 .{Fr } (6)
Step 11: Calculate the determinant of the reduced stiffness matrix, │Kr│,
» If │Kr│ sign changes, go to step 12
» If not, increment υ and Count., and return to step 6:
υ= υ + 0.01 and Count. =
Count. + 1 (7)

223 | F a s c i c u l e 4
A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

Step 12: Continue the calculations with the penultimate value of the argument υ, with a small increment (0.001)
to improve the accuracy in determining the critical argument, υcr, and therefore, the critical load, Pcr; stop
operations as soon as │Kr│ changes sign, and take υ for which │Kr│> 0, as critical argument, υcr.
5. FLOWCHART OF THE COMPUTER PROGRAM
The algorithm required for the nonlinear analysis of plane frames (by stiffness matrix method) is represented in
the form of a flowchart (Figure 3) showing the important tasks to be performed:

Figure 3: Flowchart of nonlinear analysis of frames

6. NUMERICAL EXAMPLES
— Example 1
Let's start with a digital application for a basic construction. Let us study the stability of the simple symmetrical
frame and symmetrically loaded shown in Figure 1a for the numerical values:
Lp = 12.00 m, Lb = 10.00 m, Ip = Ib = 18260 cm4, E = 21000 kN/cm2, A = 331 cm2, P = 1000 kN.
The following table presents the computation of the determinant │Kr│ and the nodal displacements {Dr} of the
frame, for increasing values of υ:
Table 1: Determinant and nodal displacements of the frame 1
Load P Nodal displacements {Dr} of the frame
No. υ │Kr│
(kN) d.o.f. Node B Node C
u (mm) 1.34102E-17 1.34068E-17
1 0.01 2.66E-02 6.65824E+22 v (mm) -4.59718E-05 -4.59718E-05
θ (rad) -5.98879E-22 -6.00533E-22
u (mm) 5.36951E-17 5.36951E-17
2 0.02 1.07E-01 6.65792E+22 v (mm) -0.000183887 -0.000183887
θ (rad) -2.39882E-21 -2.39882E-21
u (mm) . .
. . . . v (mm) . .
θ (rad) . .
u (mm) 0 0
282 2.773 2.05E+03 1.37024E+19 v (mm) -3.535015094 -3.535015094
θ (rad) 0 0
u (mm) 1.81899E-12 4.54747E-13
283 2.774 2.05E+03 -2.20118E+19 v (mm) -3.53756515 -3.53756515
θ (rad) -2.77556E-17 -2.77556E-17
The reduced stiffness matrix of the frame, [Kr], is a function of the load because of the geometrical non-linearity;
this is used as a stability criterion:
K r > 0 ⇒ stable equilibrium
K r < 0 ⇒ instable equilibrium
(8)
K r = 0 ⇒ neutral equilibrium and Pcr

224 | F a s c i c u l e 4
 A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

The following table presents the calculation of some characteristics of the frame:
Table 2: Critical parameters of the frame 1
Element Critical value υcr Pcr (kN) λcr ρ μ
Post 2.773 2.05E+03 2.047657493 0.77911218 1.132921981
The value υcr is close to that given by the table of Bleich's book [9]; (υcr = 2.775).
We have:
2
υ .E.I
Pcr = cr : critical load
L2
2
υcr .E.I Pcr
λcr
= = : critical load parameter (9)
P.L2 P
Pcr π 2 .E.I
=ρ = where PE : Euler ' s load
PE L2
π
µ= : effective length coefficient
υcr
— Example 2
Now let us consider the same construction as before, but for which the transverse displacement of the nodes
is prevented (Figure 1b). The following table presents the computation of the determinant │Kr│ and the nodal
displacements {Dr} of the frame, for increasing values of υ:
Table 3: Determinant and nodal displacements of the frame 2
Load P Nodal displacements {Dr} of the frame
No. υ │Kr│
(kN) d.o.f. Node B Node C
v (mm) -4.59718E-05 -4.59718E-05
1 0.01 2.66E-02 2.45727E+20
θ (rad) 8.27181E-25 -4.1359E-25
v (mm) -0.000183887 -0.000183887
2 0.02 1.07E-01 2.45724E+20
θ (rad) 1.65436E-24 -1.65436E-24
. . . . . . .
v (mm) -11.92446452 -11.92446452
514 5.093 6.91E+03 2.50735E+15
θ (rad) 0 0
v (mm) -11.92914767 -11.92914767
515 5.094 6.91E+03 -8.98958E+16
θ (rad) 8.67362E-19 -4.33681E-19
The following table presents the calculation of some characteristics of the frame:
Table 4: Critical parameters of the frame 2
Element Critical value υcr Pcr (kN) λcr ρ μ
Post 5.093 6.91E+03 6.907246073 2.628134619 0.61684521
The value υcr is close to that given by the table of Bleich's book [9]; (υcr = 5.095).
— Example 3
Now let us consider the same construction as that of Example 1, but for which the supports are pinned (Figure
1c). The following table presents the computation of the determinant │Kr│ and the nodal displacements {Dr}
of the frame, for increasing values of υ:
Table 5: Determinant and nodal displacements of the frame 3
Nodal displacements {Dr} of the frame
Load P
No. υ (kN) │K r│ d.o.f.
u (mm) v (mm) θ (rad)
Node A -3.3087E-24
Node B 3.3881E-20 -4.5972E-05 -1.6544E-24
1 0.01 2.66E-02 2.01414E+36
Node C 6.0986E-20 -4.5972E-05 -2.4815E-24
Node D -6.6174E-24
Node A -2.9779E-23
Node B 2.4395E-19 -1.8389E-04 -3.3087E-24
2 0.02 1.07E-01 2.01375E+36
Node C 1.8974E-19 -1.8389E-04 -9.9262E-24
Node D -2.6470E-23
. . . . . . . .
Node A -1.0608E-13
Node B 9.0540E-10 -8.7549E-01 -2.0027E-14
140 1.38 5.07E+02 2.28503E+33
Node C 9.0540E-10 -8.7549E-01 -2.0022E-14
Node D -1.0610E-13
Node A 7.5403E-11
Node B -6.4349E-07 -8.7676E-01 1.4232E-11
141 1.381 5.08E+02 -3.21836E+30
Node C -6.4349E-07 -8.7676E-01 1.4227E-11
Node D 7.5417E-11

225 | F a s c i c u l e 4
A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

The following table presents the calculation of some characteristics of the frame:
Table 6: Critical parameters of the frame 3
Element Critical value υcr Pcr (kN) λcr ρ μ
Post 1.38 5.07E+02 0.50712585 0.192956062 2.276516416
— Example 4
Now let us consider the same construction as before, but for which the transverse displacement of the nodes
is prevented (Figure 1d). The following table presents the computation of the determinant │Kr│ and the nodal
displacements {Dr} of the frame, for increasing values of υ:
Table 7: Determinant and nodal displacements of the frame 4
Nodal displacements {Dr} of the frame
No. υ Load P │Kr│ d.o.f.
(kN)
u (mm) v (mm) θ (rad)
Node A -1.2408E-24
Node B -4.5972E-05 1.6544E-24
1 0.01 2.66E-02 3.0846E+34
Node C -4.5972E-05 0.0000E+00
Node D -1.0340E-24
Node A 0.0000E+00
Node B -1.8389E-04 -1.6544E-24
2 0.02 1.07E-01 3.08449E+34
Node C -1.8389E-04 -3.3087E-24
Node D -4.1359E-24
. . . . . . . .
Node A 7.4810E-17
Node B -6.1212E+00 -4.8681E-17
375 3.649 3.55E+03 3.569E+30
Node C -6.1212E+00 4.8572E-17
Node D -7.4593E-17
Node A -9.1897E-16
Node B -6.1246E+00 6.0043E-16
376 3.650 3.55E+03 -2.88946E+29
Node C -6.1246E+00 -6.0043E-16
Node D 9.1919E-16
The following table presents the calculation of some characteristics of the frame:
Table 8: Critical parameters of the frame 4
Element Critical value υcr Pcr (kN) λcr ρ μ
Post 3.649 3.55E+03 3.545727066 1.349111926 0.860946192
7. ANALYSIS OF THE RESULTS
The following figure shows the variation of the determinant of reduced stiffness matrix of the structure
according to the upsilon argument for the four structures previously studied.

Figure 4: Determinant-upsilon curves

226 | F a s c i c u l e 4
 A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

The following table summarizes the critical load values obtained for the four previous examples:
Table 9: Summary of critical load values
Critical loads Pcr (kN)
Ratio
Posts support conditions (P(f)cr/ P(p)cr)
Fixed(f) Pinned(p)
Unbraced structure1 2.05E+03 5.07E+02 4.0434
Braced structure2 6.91E+03 3.55E+03 1.9465
Ratio (P2cr/ P1cr) 3.3707 7.0020
The analysis of the table reveals that windbracing or increase of supports stiffness increases the value of the
critical load (less sensitivity to the phenomena of elastic instability).
The results obtained by Chajes [10] and Timoshenko [11] strengthen our position.
8. CONCLUSIONS
In this paper, we presented a simple and fast method for the nonlinear analysis of plane frames, using stability
functions, and updating the stiffness matrix of the structure at each iteration. The reduced stiffness matrix
determinant, │Kr│, and the stability function argument, υ, were used to check the singularity condition of the
reduced stiffness matrix, [Kr], and the iterations were stopped when │Kr│ changed sign. At this stage, the critical
load and some parameters related to it have been determined.
Through the results obtained for the cases studied, the windbracing and the increase of the stiffness of the
supports were identified as stability factors.
An important event observed during the analysis is that starting from the first iteration, i.e. "υ = 0.01", the
determinant │Kr│ is positive until υcr (Figure 4); which reflects the stability of the structure for loads P < Pcr.
The analysis focused on four frames; the obtained results are similar to those present in the literature.
Nomenclature
A cross-sectional area of a member (cm2)
d.o.f. degree of freedom
{d} vector of nodal deformations
{D} vector of nodal deformations of the structure
E YOUNG’s modulus of elasticity (kN/cm2)
{F} structural nodal load vector
{f} and {F} nodal load vectors
G COULOMB’s modulus of elasticity (kN/cm2)
I moment of inertia of a section (cm4)
[K] structural stiffness matrix
L length of a member (cm)
M bending moment
P axial force
[T] rotation transformation matrix
u axial displacement
v lateral displacement
V shear force
Greek symbols
θ rotation (positive counter-clockwise)
μ effective length coefficient
υ argument of stability functions
φ(υ) and η(υ) stability functions
Subscripts and superscripts
b beam
cr critical
i member index, node index
f fixed
p pinned, post
r reduced
T transpose of a matrix
References
[1] Morel J.: Calcul des structures métalliques selon l’eurocode 3, Paris, Eyrolles, 331 p., 2005.
[2] Bažant Z.P.: Structural stability, International Journal of Solids and Structures, Vol. 37, pp. 55-67, 2000.
[3] Sezer R.: Nonlinear analysis of reinforced concrete plates by finite elements methods, Selcuk University,
Graduate School of Natural and Applied Sciences, PhD, Konya, Turkey, pp. 64-95, 1995.

227 | F a s c i c u l e 4
A NNALS of Faculty Engineering Hunedoara – International Journal of Engineering
Tome XVI [2018] | Fascicule 4 [November]

[4] Gummandi L.N.B., Palazotto F.: Nonlinear analysis of beams and arches undergoing large rotations, J. Struct. Eng.,
123(4): 394-398, 1997.
[5] Bhatt P.: Structures, Longman, London, pp. 406-435, 1999.
[6] Ozer E.: Nonlinear Analysis of Structures Systems, Technical University, Doctorate Lessons Notes, Istanbul, Turkey,
pp. 37-78, 2006.
[7] Orumu S.T.: Interpolative Empirical Model for Stability of Sway Frames, IOSR Journal of Engineering (IOSRJEN),
Vol. 3, Issue 8, ||V6|| pp. 04-09, 2013.
[8] Karnovsky I.A., Lebed O.: Advanced methods of structural analysis, New York, Springer, London, Dordrecht
Heidelberg, 593 p., 2010.
[9] Bleich F.: Buckling strength of metal structures, Mc Graw-Hill Book Company, 1952.
[10] Chajes A.: Principles of Structural Stability theory, Civil Engineering Mechanics series, Prentice-Hall, Inc., Engle-
wood Cliffs, New Jersey, 1974.
[11] Timoshenko S.P., Gere J.M.: Theory of Elastic Stability, 2nd edition, Mc Graw-Hill Book Company, 1961.

ISSN 1584 - 2665 (printed version); ISSN 2601 - 2332 (online); ISSN-L 1584 - 2665
copyright © University POLITEHNICA Timisoara, Faculty of Engineering Hunedoara,
5, Revolutiei, 331128, Hunedoara, ROMANIA
http://annals.fih.upt.ro

228 | F a s c i c u l e 4
Copyright of Annals of the Faculty of Engineering Hunedoara - International Journal of
Engineering is the property of University Politehnica Timisoara, Faculty of Engineering
Hunedoara and its content may not be copied or emailed to multiple sites or posted to a
listserv without the copyright holder's express written permission. However, users may print,
download, or email articles for individual use.