Computers & Swuctures Vol. 34. No. 2. pp 239-250. 1990 0045.7949190 13.00 + 0.

00
Printed in Great Britain. Q 1990 Pergamon Press plc

TIMOSHENKO BEAM THEORY WITH EXTENSION

EFFECT AND ITS STIFFNESS EQUATION FOR
FINITE ROTATION
TETSUO IWAKUMA
Department of Civil Engineering, Tohoku University, Aoba, Sendai 980, Japan

(Received 7 February 1989)

Abstract-The nonlinear effect of axial extension is taken into account in the elastic Timoshenko beam
theories and corresponding stiffness equations for finite rotation. Elastic buckling problems of a column
and a beam on an elastic foundation are solved to examine the effect. This effect enhances the critical load
of a column, and plays an important role in the buckling loads of higher modes. In the formulation of
theories, two constitutive models for the shear resultant force are introduced, leading to two different
buckling formulas, one of which corresponds to Engesser’s formula and the other to the so-called modified
one. The stiffness equations for these two models with finite rotation are derived by superposition of the
rigid finite rotation and substantial deformation of a finite element. In order to show the consistency
between the theories and stiffness equations, the FEM solutions with finite rotation are compared with
direct solutions of the differential equations by numerical integration. A few more illustrative examples
are solved to ensure the potentiality of the FEM solutions with finite rotation.

1. INTRODUCTION The effect of shear deformation becomes more

significant in proportion to the square of the
Within the framework of the Timoshenko beam thickness parameter [8]. However, it is known that
theory, there exist two distinct constitutive models for the effect of extension (shortening) of the axis
the shear force. In the buckling problems, these two becomes relatively larger for shorter bars. Hence, the
models lead to two formulas, Engesser’s formula and difference between the constitutive models mentioned
the so-called modified one. Ziegler [l] has shown that above becomes more evident. Ziegler [l] and
Engesser’s approach holds for bars, and that the Iwakuma and Kuranishi [3] have pointed out that the
modified approach applies for helical springs. The shortening effect on the buckling load becomes
difference between these models has also been important for short columns. As is clear from govern-
discussed in terms of stress resultants, in conjunction ing equations in nondimensional form, this effect is
with the strain energy, by Reissner [2]. Furthermore also proportional to the square of the thickness
Iwakuma and Kuranishi [3] have shown that parameter. Therefore both effects must be taken into
Engesser’s formula is obtained in the case of bars, account simultaneously for a more rigorous analysis
when the constitutive laws are specified in terms of of bars.
the physical components of the 2nd Piola-Kirchhoff As for the Timoshenko beam element, several
stress tensor and the corresponding strain measures, stiffness matrices have been derived in many manners,
using a similar approach by Reissner [4]. but some based on the FEM approach do not reduce
Such a rigorous approach for bars indicates the to the conventional one of the Bernoulli-Euler beam
contribution of the normal stress component to the element when the shear deformation is neglected.
shear force on a cross section. This leads to a Iwakuma et al. [9] have reviewed such stiffness
nonlinear constitutive equation for the shear force matrices, and have shown that the appropriate one
coupled with the extension of the axis. The difference can be obtained by choosing displacement functions
between the two models essentially stems from the with at least five degrees of freedom for bending in
treatment of this nonlinear term. one element. The corresponding geometric stiffness
Recently Cheng and Pantelides [S] have examined matrix can also be obtained by the same method [lo].
the differences between these constitutive models in A similar approach was taken by Mourelatos and
the dynamic problem of beam<olumns on an elastic Parsons [6] to obtain a Timoshenko beam-column
medium. In the stability analysis of a beam-column stiffness matrix on an elastic foundation.
on an elastic foundation, the Timoshenko beam However, there seems to exist no report on the
theory must be employed, because a buckling load stiffness matrices of a Timoshenko beam-column
of the higher mode may become smaller than that of with the extension effects. EzechilS [l l] derived such
the first one. Mourelatos and Parsons [6] formulated a matrix of thick arcs but for small displacements by
a finite element of such a beam-column, and directly solving the governing differential equation.
Yokoyama [7] examined its parametric instability. Dvorkin et al. [12] employed an incremental
239
240 T. IWAKUMA

approach to solve curved Timoshenko beams with Reference state

finite rotation, but neglected the nonlinear terms of 4 X
extension. II
The objectives of the present work are to summa- 9,(X,Y 1
rize the Timoshenko beam theories in conjunction
~W)
with the extension effect and the difference between
constitutive models, to derive a corresponding -tf-
beam-column stiffness matrix by the FEM approach, v(X)1I
Current configuration

L*
and to deduce a stiffness equation for finite rotations.

2. GOVERNING EQUATIONS
A(X)
2.1. Kinematics u(X) ,”
AM)
Since the formulation of the Timoshenko beam
-7
theory has been given in the literature, e.g. [3], here c$v) q(x,yy(x4yJ
we simply enumerate the necessary governing
equations. Figure 1 shows the kinematics of the
Timoshenko beam, where the shear deformation of a Fig. 1. Kinematics.
cross section is represented by the angle A(X)
between its normal n and the axis. The upper-case
letter X indicates the x-component of the initial
position of the axis. Let u(X) and v(X) denote the
where A is the area of a cross section. While T is the
x- and y-components of the displacement of the axis,
shear stress on a cross section in the G,-direction, u
respectively. The convected unit base vector gi moves
is a physical component of the 2nd Piola-Kirchhoff
to Gi in the deformed configuration, which no longer
stress tensor in the G,-direction. Since a is not acting
forms an orthogonal system. Then the physical
normal to a cross section, the axial force and bending
components of the strain measure can be defined by
moment are defined by its normal component, and
the extension c(X) and the shear deformation y(X) as
the contribution of a to the shear force emerges as the
second term of the integrand in eqn (4c). Let p(X)
6=&cosn-1
and q(X) be the x- and y-components of the dis-
tributed load per unit initial axial length, then the
y =&sinA (1)
equilibrium equations are expressed by
where
(Ncos1 - Vsin1)‘+p=O (5a)
go=(l +u’)2+(u’)2 (2)
(Nsin1 +Vcosh)‘+q=O (5b)
and a prime stands for differentiation with respect to
M’+(l+~)f’-yN=O. (5c)
X. Geometrical consideration leads to the following
kinematic relations:
Correspondingly, the boundary conditions for a
beam with length L can be given by
u’=(l+~)cosl -y sin1 - 1 (3a)
24= given or v(Ncos1 - Vsinl)=f @a)
u’=(l +c)sini +ycos1 (3b)
where k(X) is the rotation of a cross section. v=given or v(N sin). + VcosI)=s (6b)
2.2. Equilibrium and boundary conditions
1=given or VM =c (6~)
A fundamental approach [3] chooses the definition
of the three stress resultants; axial force N(X), bend- where f and s are the x- and y-components of the
ing moment M(X) and shear force V(X) as applied concentrated loads at the end; c is the applied
concentrated moment there, reckoned to be positive
N= acos4dA
when it is clockwise; v is defined by
(44
5A
1 atX=L
v= (7)
M= acos4(-Y)dA (4b) ( -1 atX=O.
sA
2.3. Constitutive relations
V= (7 +a sin4)dA (k) Only the elastic materials are considered, and the
sA strains are assumed to be negligibly small when
 Timoshenko beam theory with extension effect 241

compared with unity. Although the constitutive rela- We call this approximation the ‘small extension ap-
tions have been discussed in terms of the stress proximation’ hereafter. Note that c in eqn (3) of the
resultants in conjunction with the strain energy [2,4], kinematics is not neglected, because otherwise the
it may be straightforward to specify them in terms of theory applies for inextensible beams.
the stress components and their corresponding strain
measures. In the linear elastic case, they may be given
2.5. Nondimensional expressions of field equations
by
For clear comparison, all the field equations, eqns
r~ = EC, T = Gy (‘3) (3), (5) and (9)-(13), can be expressed in terms of the
nondimensional quantities if the beam is uniform in
where E and G are the elastic constants [3]. Although shape and material. Let
substitution of eqn (8) into eqn (4) results in highly
nonlinear expressions, we can linearize these nonlin- z, = (N cos 1 - V sin l)L2/(EZ),
ear constitutive relations [3] using the small-strain
assumption, to obtain z2 = (N sin i + V cos i)L2/(EI),

z, = ML/(EI), zq = u/L, zg = v/L, z6 = 1, (14)

N= EAc (W
y, = z, cos z6 + z2 sin z,, y, = - 2, sin z6 + z2 cos z6,
M= EIl’ (9b) 91 = PLOW), 42 = qL3/(E0.

and
Then the field equations for Model A can be written
I’ = GkAy from eqns (3), (5) (9) and (10) as
(10)

where I is the moment of inertia of a cross section, i, = -4, (15a)

and the shear coefficient [13] k is introduced in eqn i,= -q* (15b)
(IO). This constitutive model is hereafter called the
‘1st order approximation’ or simply ‘Model A’. i3= -{I +P2(1 -a)y,}y, (15c)
Since all the nonlinear terms of strains are neg- i, = (1 + /?‘y,) cos z6 - af12y2sin z6 - 1 (15d)
lected above, the shearing contribution of the normal
stress in eqn (4c) is also discarded in deriving eqn i, = (1 + /I*y, ) sin z6 + x/3’y2cos z6 (15e)
(10). Approximating $J by A in eqn (4a, c), and using i, = 2, (15f)
eqn (1), we can keep a leading term of this contribu-
tion as
where a superposed dot indicates differentiation with
V=(GkA +&), (11)
respect to the nondimensional variable (X/L), and

a = E/(kG), /I = m/L. (16)

in place of eqn (10). Since N is given by eqn (9a), the
second term of the right-hand side of eqn (11) is a
/J’ is called the thickness parameter [8]. Then the
nonlinear term of c and y. The model represented by
equations for Model B are the same as eqn (15) except
eqns (9) and (11) is henceforth called the ‘higher order
eqns (1 %)-(I 5e) which must be replaced by
approximation’ or ‘Model B’.

&= -y2(l +B*Y,)/{~ +aS*y,l(l +B2yl))

2.4. Small extension approximation
i, = (1 + /3*y,) cos zg
The small-strain assumption has been applied in
the constitutive relations where only the most domi- - aP2y2sin z,/{ 1 + aB2y, /( 1 + B2y,)) - I
nant terms are kept in linearization. However, in
i, = (1 + p’y, ) sin z6
particular when relatively slender beams are consid-
ered in elasticity, the axial extension is negligibly +aB*y,cosz,/{l +aP2yll(l +P’Y,)). (17)
small compared with unity. Therefore in such beams,
E can be neglected in the moment equilibrium of eqn
The difference between the constitutive models ap-
(5c), which becomes
pears only in these three equations.
From the definition of eqn (14), the extension is
M’+V-Ny=O. (12)
expressed by c = b*y,. Therefore, when the ‘small
extension approximation’ of the previous section is
Furthermore the constitutive equation of Model B,
applied, eqn (I SC)must be replaced by
eqn (1 l), is also approximated by

V = (GkA + N)y. (13) i3 = -(I - a8*y,)y2 (18)

242 1. IWAK~JMA

for Model A. Similarly eqn (17) of Model B must be sented by eqn (17) has more nonlinear terms than
reduced to Model A, the characteristic equation becomes a cubic
equation for z as
i, = -Yz/(l + H12YI)
z(1 -/32z)’ 7tZ
.?,=(I +/?*y,)cosz,-cr#I*y,sinr,/(l +aj3*y,)- 1 (25)
1 -(l +a)/J*z=-1-’
(19)
and no explicit formula of the buckling load is
i, = (1 + /?*y,) sin zs + @*y, cos z.J( 1 + rQ*y,).
obtained. Equation (25) is exactly the same as equa-
tion (35) of iwakuma and Kuranishi [3], which was
3. BUCKLING ANALYSIS
obtained without any approximation on the constitu-
3.1. Column buckling tive model where only eqn (8) was assumed.
Moreover, when the ‘small extension approxima-
For comparison between theories, the buckling
tion’ is applied, eqns (18) and (19) must be used. Then
analysis of a column is carried out first. Let a column
Model A yields the formula
of length L with one end fixed be subjected to a
compressive force P at the other end. No distributed
load is applied. Then the fundamental solution is (26)
given by the nondimensional quantities as
which is the so-called modified formula. On the other
z, = - PL’/(EI), z* = 0, z3 = 0,
hand, Model B leads to
z, = - PX/(EAL), zs = 0, z, = 0 (201
XL/4
(27)
for all the theories. The buckling analysis can be r = 1 + ?@r*/4
carried out by seeking solutions in the vicinity of the
fundamental solution, eqn (20). To this end, let which is Engesser’s formula. Equations (26) and (27)
zi:=[fundamental solution in eqn (20)] + Azi, can be derived from eqns (24) and (25) respectively by
(i = l-6), and substitute them into the governing setting fi* to zero but keeping a/?‘.
equations. Linearizing with respect to AZ, we eventu- Buckling loads predicted by these formulas are
ally obtain a linear differential equation for bending shown by the solid lines in Fig. 2 for typical values
expressed only by AZ,. For Model A, where the 1st of a and fl; a = 3, 0.05 6 /3 6 0.2. Comparison
order approximation is applied, such a linearized between the 1st or higher order approximations and
equation is obtained from eqn (15) as the two famous formulas indicates that the
shortening effect enhances the buckling loads about
ii& + z(1 - fit{1 - a)z) AZ, = 0 (21) a few percent. On the other hand, the higher order
approximation or Engesser’s formula predicts smaller
with the boundary conditions given by critical loads than the 1st order or the modified one,
because the nonlinear contribution of normal stress
AZ,(O) = 0, A&(O) = 0, in the shear force given by eqn (11) decreases the
apparent shear rigidity in compression. These
A&(L) = 0, differences become much more significant as the
column becomes shorter. In the same figure, the
iii,(L) +z{l -/?*(l -a)z} A&(L) =0 (22) approximate formula with shortening effect by
Ziegler [I] is shown by a dashed line, and gives critical
where loads rather close to the 1st order approximation of
the m-esent work.
z = P,,L’/(EZ) (23) *
3.2. cockling of a simple beam on an elastic foundation
and P,, is the buckling load. The ordinary eigenvalue The Timoshenko beam theory is usually employed
analysis of eqns (21) and (22) leads to the critical to solve a beam on an elastic foundation. Let kJ
loads, the minimum of which is obtained as denote the elastic constant of the Winkler founda-
tion. Then the field equations for such a beam can be
1 - Jl - 7?/32(1 -cc) obtained by addition of ( -kfv) to the left-hand side
r= (24)
282(1 -a) . of eqn (5b), and eqn (15b) must be replaced by

This formula is identical to equation (2.67) of Timo- i, = s zs - q* (28)

shenko and Gere [14] and equation (31) of Iwakuma
and Kuranishi [3]. where
The same manipulation follows when the higher
order approximation is taken. Since Model B repre- s = kfL4/(EI). (29)
 Timoshenko beam theory with extension effect 243

A Km with E=O

Number of Elements
Fig. 2. Buckling formulas of a column and finite element solutions.

The same procedure in the previous section yields a for the n th buckling mode, while Model B again has
linear differential equation for AZ, for each theory no explicit formula but the characteristic equation for
from eqns (15), (17)+ 19) and (28), but the boundary z is given as
conditions must be given by
-s(l -~%)3=(n7r)4(1 -jI2(1 +a)z}
AZ, = 0, AZ, = 0, (30) - (nn)2(1 - /Fz){z(l - 822) - a/?%}. (33)
at both ends. For example, a linearized field equation
Furthermore when the ‘small extension approxima-
of Model A is given by
tion’ is applied, the critical loads are given by

&,+[z{1 -/?2(1 -E)Z} -c$2S]AZs

z = -(n7rY-sB2(l -a)+$
+s(l -j’z){l -/12(1 -cc)z}Az,=O (31) (34)
2cIfi2{(nk)2 + S/I*}

which reduces to eqn (21) when s = 0. R = {@Jr)’ + S/12(l - LX)}’+ 4@32{(nrr)2

When the dynamic effect must be included, the
right-hand side of eqn (5b) is simply replaced by + s/?2}((n71)4
+ s + aS2s(nn)2}
(PA d2v/8r2), where p is mass density. The rotary
inertia is also taken into account by addition of for Model A. and
(pla21/at2) in the right-hand side of eqn (5~). Then
the linearized equations for AZ, of Models A and B (nn)4+ a/9*s(nn)2+ s
under the ‘small extension approximation’ coincide z=(m)2 + ap2(nn)4+ 1% (35)
with those of the second and first approaches by
Cheng and Pantelides [5], respectively. for Model B. If the extension effect on the kinematics,
The eigenvalue analysis of Model A from eqns (30) eqn (3), is also neglected as an approximation for
and (31) leads to the buckling formula as relatively slender beams, the formula eqn (35) be-
comes
(nn)2 + $I(2 - Lx)- JE
(32) (nn)4+ ag2s(nn)2+ s
z = 2/?2(1 - cr){(nn)2 + $2) ’
(36)
’= + ajJ2(n7t)4 ’
(r17r)~
Q = {@>’ + sP2(2 - a,}’ - 4/?2(1 - ‘X){(?nr)2
+ sp2}{(nn)4+ s + ap2s(nn)2} This eqn (36) eventually becomes identical with a

Table 1. Buckling loads, P,,L*/(EI), of a simple beam on an elastic

foundation; a = 13/(5 x 0.85) u 3.06, B = 0.08, s = 2n4
Model

1st order approx. Higher order approximation

A B
Extension finite small finite small neglected
n=l 26.5 25.6 25.7 25.3 28.0
2 33.8 30.9 26.5 26.7 27.2
3 55.1 48.9 31.6 34.4 34.6
4 79.6 69.3 34.3 39.8 39.8
5 105.0 90.6 35.7 43.1 43.1
244 T. IWAKUMA

formula for the Bernoulli-Euler beam [14], when CIis where the quantities in the last term are defined in
set to zero. Fig. 3. For Model A, substitution of the constitutive
For the specific values for each parameter used by eqn (10) into the first term of eqn (39) yields the
Yokoyama [7], the buckling loads calculated from internal virtual work (I.V.W.) term expressed as
eqns (32)-(36) are shown in Table 1. Yokoyama’s
numerical results using the conventional stiffness (I.V.W.) = GkAy 6y + EZI’ 61’
matrices show good agreement with those of eqn (36),
the right-most column of the table, where all the
nonlinear extension effects are neglected. In compari- +& (v’&‘- ydy}]dx. (40)
son with results of Model A, Model B shows very
slow increase of the buckling loads as n increases, On the other hand, substitution of eqn (I 1) into eqn
where the contribution of the normal stress to the (39) leads to
shear force again plays an important role. The higher
the mode is, the more significant this effect becomes. (I.V.W.) = L GkAydy + EIl’&l
S[0

1
N
4. STIFFNESS EQUATION +- 0’6u’ dx, (41)
I+6
4.1. StifJizess matrix
for Model B.
The total Lagrangian formulation with the polar
For simplicity, we here choose y and v to which the
decomposition theorem [ 151 is employed instead of
displacement functions are assumed. Since y itself
the incremental approach. The total large motion of
appears in eqns (40) and (41), the simplest displace-
an element is decomposed into its rigid finite rotation
ment function for y is constant. Then v must be
and the substantial small deformation [16]. Using the
expressed by the third order polynomial as is the case
same approach, Iwakuma et al. [17] deduced a stiff-
for the Bernoulli-Euler beam element, in order to
ness equation of the Timoshenko beam without the
obtain a suitable stiffness matrix [9]. Consequently
extension effect, and showed that the correct solution
one element temporarily has five degrees of freedom
with finite rotation can never be obtained if the
for bending, and the size of the corresponding stiff-
geometric stiffness matrix is neglected.
ness matrix is 5 x 5. Furthermore, following the
Therefore the relatively small displacement theory
procedure given by Iwakuma et al. [9] and Hasegawa
such as the beam<olumn theory is enough to derive
et al. [lo], we condense such a stiffness matrix into a
the stiffness matrix. The linearized equilibrium
4 x 4 matrix by eliminating an extra freedom adopted
equation in the axial direction coincides with that of
in y, and then linearize it with respect to N.
the infinitesimal displacement theory. When the
When the distributed loads are neglected, the final
nonlinear terms due to finite rotation are approxi-
form of the stiffness equation from eqns (40) or (41)
mated, eqns (3b) and (5b) may be rewritten as
becomes with the linear stiffness in the axial direction
r’=(l +t)i fy (37)
f= (K, + K(Nml)d, (m = A, B) (42)
and
where a superscript (m) distinguishes the constitutive
(Nj. + V)’ + q = 0. (38) models, and
From eqns (38) and (5c), we formally write the virtual f’=[fff;], dT=[d:d;]
work equation of one element of its length L for
bending as f;= @‘@I) s,L2/(EZ) c,L/(EZ)] (43)
df=[u,/L v,/L A,] (j= 1,2)
o= - L,6r-{(Ni + V)‘$-q}
s0
+~~{M’+ V(1 +c)-NA}]dx.

For the present, keeping N and .S constant, we Initial ’

74
integrate the above equation by parts. Then using -L Configuratioi7
eqn (37), and substituting eqn (9b), we obtain the
u1
virtual work equation as
YJ t h’ uz

51

Current V2
Configuration -&y
A2

+s,Sr,+c,6iz]=0 (39) Fig. 3. Definitions of nodal quantities of an element.

 Timoshenko beam theory with extension effect 245

each quantity of which is defined in Fig. 3, where

( )’ indicates a transpose. The stiffness matrices are
expressed as follows:

l/B2
(K,), =
Symm. (4 + 12Y)/A 1

1
0 0
-12/{(1 +c)‘A} 6/{(1 +c)A}
-6/((1 + C) A} (2 - 12Y)/A w

1
l/B’ 0
6~13= 12/((1 !Q2 Al -6/{(l + ~1 A}
Symm. (4 + 12Y)/A

which reduce to the conventional linear matrices

when t = 0; the part which corresponds to the geo-
metric stiffness is given by

1
0 0
Pm, =
I Symm.
6 A\m’,{5; + E) A’} A$“’/( 10 A2)
(1 +&)(2/H + A$“‘}/A’

1
0 0 0
{R$yJ},=
I 0
0 - 6 A\“‘#(
- A:“”/( 110+A*)
6) A’} -(I +~){1/30+A’;“‘}/A~
Aim’/( 10 A2) (45)

1
0 0 0
(WI3 =
L Symm.
6 A\‘“‘/{5( 1 + C) A’} -A$“)/(10 A’)
(1 +6)(2/U + A~““}/A*

A=1+12Y Y = o$P/(l + E)2,

1 + 2OY, (m = A) A(,,,)_ 2Y - 24Y2, (m = A)

A\“’ = 2-
1 + 10 A’i”‘, (m = B), 2Y + 12Y2, (m = B) (46)

A@,)_ 2Y + 48Y*, (m = A) 1 - 720Y2, (m = A)

3 -
Ai”’ =
Ai@, (m = B), 1, (m = B).

Since the difference between the two constitutive

models stems from a nonlinear part, it appears only
in eqn (45). The effect of extension is described by E de formations represented by the stiffness matrices
in the elements of matrices. When the ‘small exten- derived in the previous section. Following the proce-
sion approximation’ is applied, 6 must be set to zero. dure of Iwakuma et al. [17], it can be written as
Note that those matrices of Model B with E = 0
become identical with the conventional ones of the f, = T’K2T(d2 - d, - D)
Timoshenko beam element [lo, 181. (47)
f, = TT K, T(d, - d, - D)
4.2. St@ness equation for finite rotation
As mentioned in the previous section, a stiffness where
equation for finite rotation is deduced by the super-
position of the finite rigid rotation and substantial Ki = (K,), + z{R$“),, j = 2,3, m = A, B (48a)
246 T. IWAKUMA

D’= [(cos I, - 1) sin 1, 0] (48b) sin A, 1 0 0

0 -1 -cos1, 0 1 0 (54e)
0 0 -1 0 0 1
T = (IIin’i, z;l $ (48~)
cos1, 0
-sin,& 0 (54f)
and the axial force, z in eqn (48a), can be approxi-
0 0
mated by
R =dK,=@ti}, ~_% 2y %I
z = [{(U, - u,)/L - (cos I, - 1)) cos 1, ’ dc p’+ a(1 +c) i fc ay’
+ {(u, - u,)/L -sin A,} sin l,]/fi*. (49) j=2,3, m =A,B. (54g)

Therefore the extension in the matrices is also com- In eqn (54g), all the elements in K,‘s are formally
puted from considered as functions of three independent vari-
ables, z, (1 + t) and Y. Furthermore, note that the
E = z/P. (50) second and third terms on the right-hand side of eqn
(54g) vanish when the ‘small extension approxima-
The finite rotation is taken into account explicitly in tion’ is applied.
eqn (47) by D and T. Therefore eqn (47) is highly
nonlinear, even when the ‘small extension approxi-
mation’ is applied and t is set to zero in the elements 5. NUMERICAL EXAMPLES

of the matrices. It will be seen from comparison

In order to show the correctness of matrices and
between theories and FEM solutions that the exten-
the applicability of stiffness equations, eqn (47) is
sion in eqn (3) is eventually secure even when 6 is set
solved by the Newton-Raphson method. Since capa-
to zero in matrices, because of D and T.
bility of the present stiffness equation for finite
In order to solve eqn (47), an iterative technique
rotation has been demonstrated by Iwakuma et
may be employed. When the Newton-Raphson
al. [ 171 within the ‘small extension approximation’,
method is used, a solution of the (n + 1)th iterative
we here solve a few examples with emphasis on the
step can be calculated from
effect of axial extension. Note that realistic values are
not always specified for j in order to exaggerate the
variation of solutions.
At each incremental loading step, solutions are
updated iteratively until the condition,
f, - (T?<,T(d2 - d, - D)}‘“’ 1,I’”+ 1)_ ,I(“)1/Id@‘+‘)/ < 10-9, is satisfied. An arc
’ f, - {T?<,T(d, - d, -D)}‘“’ ’ (“)
length method [19] is employed to trace the nonlinear
behavior of structures stably. The skyline
K, is the tangent stiffness matrix defined by
technique [20] is also adopted, so that the numerical
computations can be carried out on a 16-bit
K = TTK2T(d2-d,-D)
I (52) microcomputer.
TTKJT(d2-d, -D) ’
5.1. Column buckling
which becomes nonsymmetric. Some cumbersome First examined is the buckling of columns solved
manipulation leads to its explicit expression as analytically in Sec. 3.1. The bifurcation points are
obtained from the condition of singularity of the
H2 0 0 s, 0 0 0 tangent stiffness matrix of eqn (53). However, since
K, =
O(H3 + 0 0 s3 0 0 0 the tangent stiffness matrix is highly nonlinear with
respect to the axial extension, we here solve a straight
+ F [-cos 1, - sin 1, g cos 1, sin 1, 0] column incrementally and iteratively by eqn (51)
(53)
0 3 until the tangent stiffness becomes singular. The load
at which the matrix first becomes singular is the
where 0 is a 3 x 1 null vector and smallest buckling load of the column.
Numerical results with 256 elements are shown in
H, = TTKjTC (54a) Fig. 2 by circles and triangles. Closed circles indicate
S, = (Q?(,T + T?(,Q)(d, - d, - D) results using the stiffness matrix of Model A, and all
(54b)
of them lie on the solid line obtained by the formula
P, = T%,T(d, - d, - D), 0’ = 2, 3) (54c) of eqn (24), while open circles by Model B lie on that
g = -{(u, - u,)/L - (cos I, - 1)} sin 1, by eqn (25). A typical converging process of the FEM
solutions into an analytical solution is also depicted
+ {(u2 - v, )/L - sin 1, } cos 1, (54d) in the same figure, when CL= 3, /l = 0.2, and the
 Timoshenko beam theory with extension effect 247

matrix of Model B is used. It is seen that eight each figure indicate the solutions of eqns (15) and
elements are enough for a relative error less than (17)-(19) by the direct integration method. Symbols
0.1%. On the other hand, triangles in the figure show ‘A’ or ‘B’ attached indicate the constitutive models in
results wfien L in the matrices is set to zero. They also theories, i.e. eqns (15) or (17), and the ~~nthetj~l
converge to solutions by the formulas, eqns (26) and symbols imply the corresponding solutions under the
(27) which are the modified and Engesser’s formulas. ‘small extension approximation’, i.e. from eqns (18)
Consequently the stiffness matrices derived in Sec. 4.1 or (19). The method of adjoints of shooting
for the two constitutive modeIs can precisely predict methods [21] is employed, and 128 integration
the buckling loads if enough finite elements are used. intervals are used so that the obtained sofutions have
at least six-significant-dirt accuracy. The tolerance
5.2. Cantilever beam subjected to shear or compression for both the integration scheme and the iteration of
Since a bifurcation analysis in the previous shooting is set to lo-“‘.
example can be carried out without large displace- When the stiffness matrix of Model A is used, the
ment, it is necessary to solve a system which FEM solutions converge to those of eqn (15), while
undergoes large rotation in order to check the the solutions of eqn (17) are consecutively obtained
consistency between the theories and numerical from the stiffness equation with the matrix of Model
methods. To this end, the finite displacement B. The same correspondence can be clearly observed
behavior of a cantilever beam subjected to either when the ‘small extension approximation’ is applied.
shear or compression is analyzed. Typical converging Hence it becomes clear that the solutions of either
processes are demonstrated in Figs 4 and 5 by the eqn (15) or eqn (17) can be numerically obtained by
FEM solutions at the specific level of loading for the stiffness eqn (47), if the proper stiffness matrix is
short beams. In the figures, u and u express the axial chosen from eqns (44) and (45). Therefore the stiff-
and transverse dispIacements at the tip respectively. ness matrices derived in relatively small dispia~ment
A symbol K(“‘(m = A, B) near lines distinguishes the theory can fully describe the substantial deformation
constitutive models used for the stiffness matrix. of a small structural element. The finite rotation in
Dashed lines are the results when L in the matrices is plane problems is completely taken into account in
set to zero. Although the converging process is not the stiffness eqn (47).
always uniform, it is clear that these FEM solutions
5.3. Cantilever beam subjected to moment
are consistently approaching the values indicated by
short horizontal lines. Dvorkin et al. [12] solved a cantilever beam
These horizontal lines on the right-hand side of subjected to moment at the tip by an incremental

(a)

0.82

0.80
I
64 256
&nbe?of Elements timber of Elements

fb)

Nuder of ELements kmber of Elements

Fig. 4. (a and b) Convergence of FEM solutions of a cantilever beam subjected to shear.

CAS 34/2-E
248 T. IWAKUMA

-.
(a)
1.14 v/L
t . ---with E=O u==-p
Q95 g=3.7
-:’ k

l-
1.06

1.02

(b)
4
“..p

16 64
--_ ---_-_

K’PJ

256
Number of Elements

\\
-0
-A

I 4

4 16 64 256
Number of Elements

1
--p

‘\ ---with E=O

0.815
0.746. t

0.742 -(A)
-B

-A
I I
4 16 64 256 4 16 64 256
Number of Elements Number of Elements

Fig. 5. (a and b) Convergence of FEM solutions of a cantilever beam subjected to compression.

technique. The same problem is solved here to 5.4. Out-of-plane behavior of a tower of a cable-
illustrate the efficiency of the present method. Since stayed bridge
the beam is in pure bending, no difference emerges The last example is a simulation of the out-of-plane
from the two constitutive models. In Table 2, the behavior of a tower of a cable-stayed bridge. Such a
number of iterations necessary at each loading step is tower can be modeled by a column pulled toward
compared with those of the reference, though the its base by a cable attached at the top. As shown in
convergence criterion is different. The present method Fig. 6, this model is like a fishing pole pulled by the
shows uniform convergence regardless of the loading fishing line near at hand. The load here is applied at
step, but the rate of convergence seems to be the roller, which is initially located at the base but
comparable. moves vertically downward. If the load is always
The accuracy of the solution is again checked by applied through a fixed point at the base, the solution
comparison with the direct numerical solution by a is the same as that of a simple beam [14].
shooting method in Table 3, where three displace-
Figure 6 shows the relation between the applied
ment components at the tip are shown. Particularly,
load and the horizontal displacement at the top of the
the rotation is obtained quite accurately. The lateral
displacement u shows the worst accuracy, but the
relative error becomes about 0.6% when 16 elements
Table 3. Accuracy of an FEM solution of a cantilever
are used. subjected to moment, M(L) = rrEI/L, at the tip; G(= 2.88,
B = 0.0289
Table 2. Number of iterations required in solving a
Number of elements or intervals
cantilever subjected to moment at the tip; a = 2.88,
8 16 80
p = 0.0289
Present method: Dvorkin et al.: Finite element Shooting
4 or 8 elements 4 2-node elements Displacement method method
Error tolerance IO-’ lo-” lo-’ -4L)lL 0.998389 0.999799 1.00000
No. of iterations W)IL 0.652961 0.640709 0.636620
at each step 3 4 3 or 4 W)lL 3.14159 3.14159 3.14159
 Timoshenko beam theory with extension effect 249

Table 4. Governing equations, buckling formulas, and stiffness matrices of each model
Model
Extensible model Small-extension approximation
A B A B
Kinematics u’=(l+c)cosrl-ysin1-1 v’=(l +c)sinl +ycos1
Force equilibrium (NcosI - VsinI)‘+p =0 (Nsin1+ VcosI)‘+q =0

Constitutive laws
for axial force N=EAe M=EIi’
and bending moment
Moment equilibrium M’+V(l+c)-NI+m=O M’+ V-NY +m =0

Constitutive law 1st order Higher order 1st order Higher order
for shear force approximation approximation approximation approximation

V = GkAy .+A +&), V = GkAy V = (GkA + N)y

Column buckling
1 -Jl -x*/3*(1 -a) x* z(1 - B’z)’ z=JX@7-1 z= x*/4
formula z=
2fi*(l -a) T= I -(1 +a)j2z 2aB’
z = P,,L*/(EI) (modified)
Local stiffness K, + K’,“1 K, + K$‘: K, + K$ K, + K’,“1
matrix with c = 0 with 6 = 0
Equilibrium Strictly satisfied Approximately satisfied

tower. In the case of slender towers of #I = 10m4 the peaks of the corresponding curves leads to the
without shear deformation, the equilibrium path conclusion that the buckling load is decreased by the
without initial imperfection bifurcates at the buckling pre-buckling displacement of the load application
load of a simple beam, z = x2, and shows an unstable point. Therefore the shortening effect as well as the
post-buckling state, expressed by dash-dotted lines. shear effect reduces the buckling load, unlike a
The buckling loads of a simple beam can be obtained column or a simple beam.
from eqn (33) with s = 0. When tl = 0, the critical
load increases as the beam becomes shorter, owing to 6. CONCLUSION
the shortening effect. On the contrary, the present
example for a shorter tower of /3 = 0.05 and a = 0 The Timoshenko beam theories are reviewed with
shows a decrease of the buckling load, because of the respect to the shear constitutive models and the effect
motion of the load application point due to the of axial extension. The shortening effect plays an
shortening effect. important role in the elastic buckling of a simple
When the shear deformation is taken into account, beam on an elastic foundation, a column, and a tower
the buckling load decreases more, though the post- pulled by a cable. Two kinds of stiffness matrices
buckling behavior is similar. Incidentally, when a = 3 including the geometric stiffness are explicitly derived
and /I = 0.05, eqns (33) and (35) with s = 0 yield the by the FEM technique. A stiffness equation for finite
critical loads of a simple beam as z = 9.379 and 9.189 rotation deduced in plane problem yields the results
respectively. Comparison between these values and which converge to the direct integration solutions of
each corresponding governing equation. Good con-
sistency between theories and numerical methods are
shown by several examples. Comparisons between
models are summarized in Table 4.

REFERENCES

1. H. Ziegler, Arguments for and against Engesser’s buck-

ling formulas. Ing. Arch. 52, 105-113 (1982).
2. E. Reissner, Some remarks on the problem of column
buckling. Ing. Arch. 52, 115-119 (1982).
3. T. Iwakuma and S. Kuranishi, How much contribution
does the shear deformation have in a beam theory?
Strucr. Engng Eurthq. Engng JSCE 1, 103s-113s (1984).
4. E. Reissner, On one-dimensional finite-strain beam
0 0.25 theory: the plane problem. ZAMP 23, 795-804 (1972).
v/L 0.50
5. F. Y. Cheng and C. P. Pantelides, Dynamic
Fig. 6. Buckling and post-buckling behavior of a column Timoshenko beam- columns on elastic media. J. Sfruct.
pulled vertically by a cable near the base. Engag 114, 15241550 (1988).
250 T. IWAKUMA

6. Z. P. Mourelatos and M. G. Parsons, A finite element 14. S. P. Timoshenko and J. M. Gere, Theory of Elastic
analysis of beams on elastic foundation including shear Stability, 2nd Edn. McGraw-Hill, New York
and axial effects. Comout. Srruct. 27. 323-331 (1987). (1961).
7. T. Yokoyama, Paradetric instabilit; of Timdshenio 15. L. E. Malvern, Introduclion to the Mechanics of a
beams resting on an elastic foundation. Comput. Sfruct. Continuous Medium. Prentice-Hall, Englewood Cliffs,
28, 207-216 (1988). NJ (1969).
8. S. Antman, General solutions for plane extensible 16. G. Wempner, Finite elements, finite rotations and small
elasticae having nonlinear stress-strain laws. Q. uppl. strains of flexible shells. Int. J. Solids Struct. 5, 117-153
Math. 26, 3547 (1968). (1969).
9. T. Iwakuma, M. Ai and F. Nishino, On derivation of 17. T. Iwakuma, A. Hasegawa, F. Nishino and S.
Timoshenko beam stiffness equation. Proc. JSCE 312, Kuranishi, Principle and numerical check of a stiffness
119-128 (1981). equation for plane frames. Strucr. Engng Earthq. Engng
10. A. Hasegawa, T. Iwakuma and S. Kuranishi, A JSCE 4, 73s-83s (1987).
linearized Timoshenko beam theory in finite displace- 18. J. S. Archer, Consistent matrix formulations for struc-
ments. Slruct. Engng Earthq. Engng JSCE 2, 321s-326s tural analysis using finite-element techniques. AZAA Jnl
(1985). 3, 1910-1918 (1965).
Il. J. EzechiBS, Contribution to the calculation of thick arcs 19. T. Hosono, Analysis of elastic buckling problem by arc
with respect to shearing strain and extension of the length method; Part 2. Trans. Arch. Inst. Japan 243,
centroid axis. Comput. Strut/. 29, 645456 (1988). 21-31 (1976) (In Japanese).
12. E. N. Dvorkin, E. Ofiate and J. Oliver, On a non-linear 20. A. E. Elwi and D. W. Murray, Skyline algorithms for
formulation for curved Timoshenko beam elements multilevel substructure analysis. Int. J. Numer. Mefh.
considering large displacement/rotation increments. Int. Engng 21, 465479 (1985).
J. Numer. Merh. Engng 26, 1597-1613 (1988). 21. S. M. Roberts and J. S. Shipman, Two-Point Boundary
13. G. R. Cowper, The shear coefficient in Timoshenko’s Value Problems: Shooting Method.7. Elsevier. New York
.^__
beam theory. J. appl. Mech. 33, 335.-340 (1966). (lY72).