119

PROC. OF JSCE,
No. 312, August 1981

ON DERIVATION OF TIMOSHENKO BEAM

STIFFNESS EQUATION

By Tetsuo IWAKUMA*, Masahiro AI** and Funiio NISHINO***

While energy principles are used, only strain

1. INTRODUCTION energy and energy due to inertia force are
evaluated and used to derive the stiffness and
Since shear deformation is generally small in mass matrices. Since the work done by boundary
comparison with bending deformation, it is forces is not included in energy expressions in
usually neglected in static analysis of beams. The these references, energy principles do not reveal
effect of shear deformation increases with increase components of boundary forces corresponding to
of the ratio of cross-sectional dimensions to beam the assumed interpolation functions, nor boundary
length. Shear deformation plays an important conditions. Instead, they are derived from physical
role in such dynamic problems as beam vibration considerations, which have led to errors in bound-
of higher modes and response of a beam to impact ary conditions and in the assembly of a global
forces. Because of this, the influence of shear matrix from the local matrices4.
deformation has been extensively studied in It has been reported?, for problems governed
relation to dynamic analysis of beams. Beams by ordinary differential equations, that the finite
showing shear deformation are usually called element method results in the exact stiffness
Timoshenko beams. Timoshenko beam theory is matrix by employing the homogeneous solution
based on the assumption that the plane normal as an interpolation function. In spite of the fact
to the beam axis before deformation is not normal that the homogeneous solution of a Timoshenko
to the axis after deformation but that it remains beam equation is a polynomial, the finite element
a plane. With a further increase in the ratio of method using polynomials with a sufficient num-
the cross-sectional dimensions to length, the ber of degrees of freedom as interpolation func-
assumption that the plane section remains a tions does not necessarily result in an exact
plane no longer holds and, as a result, the so-called stiffness equation, but results in a variety of
shear lag phenomenon starts to play a significant stiffness matricesis. This is due to the energy
role. Timoshenko beam theory is applicable only expression for the Timoshenko beam, for which
for beams in which shear lag is insignificant. This the selection of the polynomial identical to the
implies that Timoshenko beam theory considers homogeneous solution of the governing equation
shear deformation, but that it should be small in is not obvious. The exact stiffness equation for
quantity. a prismatic Timoshenko beam can be easily
A number of finite element analyses have been obtained by integrating the governing differential
reported for vibration of Timoshenko beamsls> equations8'9, and thus it is not important to
based on energy principles. Although several kinds derive the stiffness matrix by the finite element
of polynomials and the homogeneous solution for technique. However, the energy expression for
static problems are used as interpolation functions the Timoshenko beam is a good example in
resulting in different stiffness and mass matrices, selection of an interpolation function for the
there has been no discussion of the accuracy of finite element method, since the resulting stiffness
the stiffness matrices themselves. Errors, how- matrix can be compared with the exact stiffness
ever, have been discussed for natural frequencies matrix. This paper discusses this selection of
resulting from approximations in both matrices. interpolation functions. The energy expression
for the Timoshenko beam derived from basic
* Member of J SCE, M. Eng., Research Assistant, assumptions on kinematic field and energy princi-
Dept. of Civil Eng., Northwestern Univ. ples for a three-dimensional elastic continuum is
** Member of JSCE, Dr. Eng., Lecturer, Dept. also presented, and is used to clarify the confu-
of Civil Eng., Hosei Univ. sion present in the literature in boundary condi-
*** Member of JSCE, Ph. D., Professor, Dept. of tions for the finite element formulation of the
Civil Eng., Univ. of Tokyo Timoshenko beam.
 120 T. IWAKUMA, M. AI and F. NISHINO

a<z<b as

2. GOVERNING EQUATIONS Tgdz-6uTxdz-[13uT T]z -a- [BUTT]z=b

The in-plane behavior of straight prismatic =0 ...... (6)

beams of which the cross-sectional deformation
where
is small and can be neglected is considered in this
QT=[N V MJ, ET=Le y K]
paper. A rectangular Cartesian coordinate system
(xi, x2, x3) with a reference frame of unit vectors uT=LW v A], XT=Lpz Pr ill]
(i1, i2, i3) along the coordinate axes is defined TT =LN V M]
with the x3-axis parallel to the beam axis.
...... (7ae)
From this selection of the coordinate and the
and the components of these vectors are defined
problem statement, loading and displacement by
vectors having components only on the x2-x3
plane are considered. The assumptions for the e=w', y=A+v', K=A' ....... (g.a-c)
kinematics of the Timoshenko beam can be
N=o'zzdA, V=o'ndA, M=cry dA
expressed in strain tensor terms asio
.......(g-ac)
e23(xi, x2, x3)=e32(xi, x2, x3)- y(z)
Pz=pzdA, P=' pdA, n= pzydA
eii(xi, x2, x3)=e22(xi, x2, x3)=ei2(xi, x2, x3)=0

...... (l-a-e) ...... (10a-c)

where y(z) denotes a function of z, and represents N=0 tzdA, V= tdA, M=tzydA
a shear strain constant over a cross section at z.
The notation (x, y, z) is used in lieu of x2, together ...... (11a-c)
with (x1, x2, x3) whenever convenient. The A denotes cross-sectional area.
strain-displacement relations are Integrating Eq. (6) by parts yields equilibrium
equations and boundary conditions. The result-
ing equilibrium equations in (a, b) can be expressed
as
where, ui= ii component of displacement vector
u, and ( ), i= derivative with respect to x 2. N-FPz=O, V'+Pv=O, M'-V+m=0

The displacement components satisfying Eq. (2) ...... (12a-c)

and the strain conditions of Eq. (1) can be written and the boundary conditions at z=a and b can
as be expressed as
91(x1, x2, x3)= 0, u2(xl, x2, x3)= v(z) w=w or nzN=N
93(x1, x2, x3)=W(z)+yA(z), A(z)= y(z)-v'(z) v=v or nzV=V ...... (13a-c)
...... (3a-d) or nzM=M
where ( )'=d( )/dz, and v and w = u2 and u3, where w, v, and A are the prescribed displace-
respectively, at x2= 0 and x3=z. A(z), defined ment components at the boundaries, and nz is
in Eq. (3.d), denotes the change of angle of the the direction cosine between iz and the outward
beam axis due to bending. Substitution of Eq. normal vector of the cross sections of the ends.
(3) into Eq. (2) yields -1, z=a
n2= ...... (14)
e33(x1,X2,x3)=e2z(y,z)=w'+yA' ........ (4) 1, z=b

If the ii component of body force vector is When a concentrated force Tc'=LNc Vc Mc

denoted by pi, and the surface force vector acts at z=c between a and b, addition to Eq. (6)
working on both ends by ti, the virtual work of the virtual work done by this force results in
equation can be expressed as the same equilibrium equation as Eq. (12) for
the domains, a to c, and c to b, and the identical
o-ijeijd V-13uid V-i uid S=0 boundary conditions at z= a and b as Eq. (13).
In addition, the virtual work equation gives
...... (5) the internal boundary conditions at z= c as
where V is the body and S its surface of con- we=wc+ and Nc-Nc+=Nc
tinuum, oij are the components of the stress
vc=vc+ and Vc-Vc+=Vc
tensor, and' denotes virtual kinematic quantities.
Substituting Eqs. (1) to (4) into Eq. (5) gives A=Ac+ and 1Itc -Mc+=Mc
the virtual work equation for a beam lying on ...... (15a-vc)
 On Derivation of Timoskenko Beam Stiffness Equation 121

where superscripts + and - denote the plus and fT=LTT(a) TT(b)J, qT=LuT (a) uT(b)J
minus sides of the point. (23a,-b)
Considering only an elastic material, the follow-
ing constitutive equations can be assumed with
f is the equivalent nodal force vector, and k
the stiffness matrix. The elements of these are
Young's modulus E, and shear modulus G
given in APPENDIX I. These matrices coincide
0'22=Ee22 oyz=2 Gey2 ........ (16. a, b) with the well-known results for the Timoshenko
beam.8 The stiffness matrix has a, form modified
Integrating Eq. (9) with Eq. (16) and making
use of Eqs. (1) and (4) gives the relation between to that for elementary beam theory by the co-
efficient c defined by
generalized stress and strain as
e=Ce ....... (17) q=EI/GkAl2 ........ (24)
where, l is the length of a beam element and
By selecting the z-axis on the centroidal axis,
hence
the matrix C becomes a diagonal matrix as
l=b-a ....... (25)
EA 0 0
C= 0 GkA 0 ...... (18) With increasing GA, the stiffness matrix con-
0 0 El verges to that of elementary beam theory.
where I - moment of inertia of a cross section. 3. ENERGY DERIVATION OF STIFFNESS
k is the so-called shear coefficient, and is a func- EQUATION
tion of cross-sectional shape. The derivation of
Eq. (18) results in k being equal to 1, regardless Substituting Eqs. (8) and (17) into (6), then
of cross sectional shape. Other considerations assuming proper interpolation functions for
give a value of k other than 111). displacement components, the stiffness equation
The governing equations for a beam with can be derived by the finite element technique.
shear deformation are given by Eqs. (8), (12), Selection of interpolation functions and the
(13) and (17). Since the governing equations for resulting stiffness equations are examined and
axial deformation are independent of those for discussed.
bending and shear, and they are not affected by There are two independent generalized dis-
the shear deformation, these equations are omit- placement components, v and A, as in Eq. (7. c).
ted in what follows. Thus, it is natural to assume two interpolation
Substituting Eq. (17) into Eqs. (12b) and (12c) functions for the components. Functions are
gives two simultaneous ordinary differential chosen in general to satisfy geometrical boundary
equations in terms of displacement components. conditions at the ends of an element. It is common
Elimination of A from the equation results in
practice to use polynomials as the interpolation
EI(viv+Pv"/GhA) +Pv+m'0 ......... (19) functions, the orders of which are determined
by the orders of the derivatives appearing in the
Similarly, the boundary conditions are expressed virtual work equation, such that they satisfy
in terms of displacement components by the necessary condition for convergence that
v=v or n(-EI(v"+Pv/GkA)+m)=V internal virtual work should not identically
vanish. Since the highest order derivatives of
v+GkA1 (-EI (v'+py'/GkA)-m}= the displacements appearing in the definition of
generalized strains, y and K, are first order as
or n2 (-EI (v"+ Py/GkA)}= M can be seen in Eqs. (8.b) and (8.c) and since there
...... (20-a, b) are two degrees of freedom for the generalized
displacement, v and A, one at each end, as shown
The general solution of Eq. (19) is expressed as in Eq. (15), the first order polynomial is the
v=Co+Clz+G2z2+G3z3+vp ...... (21) lowest order polynomial interpolation function
that satisfies the necessary condition for conver-
where vp is the particular solution, and Co-C3
gence. Thus expressing the displacement func-
are integration constants. Determining the tion u by
integration constants of this general solution by
use of the boundary conditions of Eq. (20), the u=Nq ....... (26)
stiffness equation of the Timoshenko beam can this lowest order interpolation function can be
be easily obtained, and is written as expressed as
f=kq-fo ....... (22) 1- 0 0
N=
0 1- 0
...... (27)
where
 122 T. IWAKUMA, M. AT and F. NISHINO

where become dependent when y is equal to zero. To

accomplish this, an interpolation function is
assinged for v, though y is not a generalized
Using the interpolation function of Eq. (27), displacement. In addition further higher order
the corresponding stiffness equation is polynomials are selected for v and A. Since no
continuity condition exists for y at the boundary
f=kiq-f1 ...... (29)
and the order of the derivative of y appearing
The elements of kl and f are given in APPENDIX in the expression of internal virtual work is zero
II. the lowest order polynomial satisfying the require-
With the interpolation function of Eq. (27), ments as an interpolation function for y becomes
the generalized strains of Eq. (8) are expressed in a constant, hence
terms of end displacements as
y=4 =const ....... (32)
By selecting interpolation functions for y and
v, the interpolation function for A is automatically
determined by Eq. (8. b). Since two kinematic
continuity conditions are present at each end of
an element as can be understood from Eq. (7.c),
a total of four degrees of freedom are necessary
For arbitrary given constants K and y, there are
for displacement functions. Selecting a third
no end displacements, A(a), A(b), v (a) and v(b),
to satisfy Eq. (30). This implies that K and y order polynomial for v for this reason, and
can not be independent. One example of this is considering Eqs. (32) and (8b) gives
that, when K is not equal to zero, y is also not
equal to zero, and thus even for a constant
moment, the shear strain can not vanish. In
order to improve on this contradiction, higher where

order polynomials have to be considered. By (1-32+23) -l(-22+3) (32-23)

selecting interpolation functions as expressed by 6(-2)/l (1-4+32) 6(2-4)/l
-l(3-2) l(1-)(1 -2)
(l-) -(1-) 2(1-) ...... (34)
(32-2) 6(-g2)
0 (1-a) 0
Using Eq. (33), the stiffness equation is
...... (31)
it is shown that both generalized strain can take
simultaneously arbitrary constant values within
an element. The stiffness matrix k2 derived from The third and fourth terms of Eq. (6), show
Eq. (31) is given in APPENDIX III. Since both there exists no nodal force, f d, corresponding to
Eqs. (27) and (31) are different from the homo- 4, thus
geneous solution of the governing differential
equations, the resulting stiffness matrices kl
and k2 do not agree with the exact stiffness With this, the last row of Eq. (35) can be regarded
matrix k. While the stiffness matrix k reduces as a constraint condition between q and d .
to that of elementary beam theory at the limit Using this constraint, Eq. (35) can be condensed
of increasing GkA, both kl and k2 do not. to
The components of generalized displacement
for this Timoshenko beam, v and A are related
where
by y, as can be seen in Eq. (8.b). For the special
case where y is equal to zero, A is expressed by
v', i.e., A and v can not be independent. For
Timoshenko beam theory based on assumption
Eq. (1.a), i.e., the assumption that a plane sec-
tion remains a plane, y can be regarded as small Eq. (37) agrees with Eq. (22), that is
as has been pointed earlier. With this con-
ka=k, fag=fo ...... (40a, b)
sideration, the selection of interpolation functions
for v and A that are independent even when y For dynamic problems, the inertia force can
is equal to zero may be regarded as inadequate. be treated as distributed force, thus the virtual
Rather, it is better to select functions that work equation can be expressed by replacing
 On Derivation of Timoshenko Beam Stiffness Equation 123

the second term of Eq. (6) by The last two rows of Eq. (49) play the same
roles as the last row of Eq. (35), however, with
inertia term=! [(-pu2)cu2 condensation of matrix sizes, the term cot ap-
+ (- pu3)8u3]dV ........ (41) pears both outside the matrix as a coefficient
and inside the matrices in the components. Since
Substituting Eq. (3) and omitting the term relat- this operation makes numerical computation
ing to axial displacement, Eq. (41) becomes difficult, it is appropriate to superimpose Eq.
inertia term= -m 1 {v(z,t)8v(z,t) (49) in its original form for the entire structure.
It has to be noted that no continuity condition
+ro2A(z,t)8A(z,t)) dz ......... (42) exists for y(a) and y(b), as can be seen in Eq.
(15). Thus the terms corresponding to y(a) and
where p= density, ('')=O( ) /ate, m = mass per
y(b) are not to be superimposed.
unit length, ro,=radius of gyration, and m and
roe are defined by
4. DISCUSSION
m = pA ....... (43)
roe=l/A ...... (44) In a beam problem governed by Eq. (12c)
under the loading condition m=0, bending
Consider a solution of the form moment distribution is uniform when shear force
v(z, t); v(z) ei"t, A(z, t). =A(z) eiwt ...... (45) is zero, and otherwise it is a function of z. With
the employment of the displacement function of
The same notations v and A, used for the func- Eq. (27), however, non-zero shear strain is always
tions of z only, are used for simplicity in Eq. (45),
present when bending strain exists, regardless
and below for functions of z and t. Substituting of its distribution. Similarly, when A(a) is equal
Eq. (45) into Eq. (42) gives virtual inertia work to A(b), shear strain is not necessarily zero in
as spite of the fact that bending strain distribution
is zero inside the element. From these facts, it
inertia term=mco2'(v8v+ro2A8A)dz(46)
...... is obvious that . deformation of an element can
not be properly represented by Eq. (27). Similar
Substituting assumed displacement functions
in Eq. (46) and integrating it results in the so- remarks are applicable to the displacement func-
tions of Eq. (31). It is known that displacement
called consistent mass matrix. Since the accuracy
functions other than the homogeneous solution
of eigen values obtained by a discretized system
of the governing differential equation can not
generally increases with increasing matrix size, represent correct deformation. While the exact
a function which has a larger number of degrees
stiffness matrix is obtained by using, as inter-
of freedom than Eq. (33) is selected by replacing
polation function, the lowest order polynomials
y=A with a linear function with y (a) and y (b) that satisfy the convergence criterion in the case
at both ends. With this function of y, the dis-
of a beam element with no shear deformation,
placement function, in view of Eq. (8b) becomes the use of the lowest order polynomial does not
q lead to the exact stiffness matrix for the Timo-
u=ND y(a) = NDgD ....... (47) shenko beam element.
y(b) Eq. (27) has sufficient degrees of freedom to
where
satisfy the geometrical boundary conditions and
also satisfies necessary conditions of convergence
ND=
(1-32+23) l(-+22-3) (32-23) determined from the order of the derivatives in
(6-62)h (1-4+32) (6g2-6)/l the expression of internal virtual work. Eq. (7)
l(2-3) l(-22+3) l(3-2) indicates that y and K are the two components of
(32-2) (332)-(3-3g2) generalized strains of the Timoshenko beam.
The displacement function of Eq. (27) can not
...... (48) account for the situation in which both y and K
Considering Eq. (46), Eq. (6) after substitution take simultaneously arbitrary constant values,
of Eq. (47) and integration gives the equation of as has been discussed in Eq. (30). Given this
motion of a discretized system fact, there remains uncertainty as to whether
the solution of Eq. (29), which depends on Eq.
kqD-mo2mgD=0=fD ........ (49) (27), does converge to the exact solution of the
problem. To check this, the discrete system of
Eq. (29) is transformed into a consistent con-
The elements of k and m are given in APPENDIX tinuous system which can be compared with
IV. the governing differential equation of the Timo-
 124 T. IWAKTJMA, M. AT and F. NISHINO

shenko beam. {GkA (v'+A)}'+Fy=O

Consider a nodal point O at z=0, where no {EIA'}'-GkA (v'+A) +nn=0 ...... (54a, b)
load is applied, and two neighbouring elements
A and B are of equal length 1. By making use Since Eq. (54) agrees with the governing equa-
of Eq. (29) for elements A and B, the zero nodal tion for the Timoshenko beam, as can be seen
force at O can be expressed in terms of u(-l), substituting Eqs. (17), (8b) and (8.c) into Eqs.
u (O) and u (l). Then, the displacement vectors (12.b) and (12.c), the solution of the stiffness
at nodes z= -1, 0 and 1, respectively, are as equation, Eq. (29), converges to the exact solu-
follows tion. From the nature of Eq. (27) which has a
larger number of degrees of freedom for shear de-
d - GkA /l GkA /2 2GkA /l formation than for bending deformation, the
0 -GkAi2 -EI/l+GkAl/6 0 stiffness equation based on Eq. (27) may result
0 - GkA/l - GkA/2 in a more accurate solution for a beam with
2 {EI/l+GkA1/3 GkA/2 -EI/l+GkAZ/6 larger shear force than bending moment ratio
for the same number of elements. In order to
Iu(-l) I-v PyAzdz+1 VPyB(l-z)dz check this numerically, two cantilever beams of
different lengths are considered. The properties
of the beams are
E=2.1 x 106kg/cm2 (210 GPa)
G=7.0 x 105kg/cm2 (70 GPa)
where subscript A and B are quantities for ele- A =30.0 cm2 , I=250.0 cm4, k=0.8333
ments A and B, respectively. The Taylor expan-
...... (55)
sion of u, P, and m at 0 leads to
with l=40 cm and 100 cm. The shear coefficients
as defined by Eq. (24) are 1.88 x 10-2 and 3.0
x 10-3, respectively, for the beams. The equa-
tions for the beams are solved for boundary
conditions of v=A=0 at the fixed end, and M=0
and V 0 at other end. The displacement at the
loaded end versus the number of elements rela-
tions for the two beams are shown in Fig. 1.
Better accuracy is obtained with the same number
of elements for the beam of shorter span, which
agrees with the above expectation.
As shown by Eq. (21), the homogeneous solu-
tion for the Timoshenko beam is a polynomial,
thus the exact stiffness equation, Eq. (22) an
and so on. Substituting Eq. (51) into Eq. (50)
be easily derived for it. It is also natural to
and integrating terms with distributed forces
results in expect that the exact stiffness equation can

Where the subscript 0 denoting values at O is

omitted for simplification. O( ) denotes Landau's
symbol, and f=O {g(l)} implies

For the limit as l approaches 0 in Eq. (52), the

consistent differential equation to the discrete Fig. 1 Number of elements and end displace-
system expressed by Eq. (29) becomes ment.
 On Derivation of Timoshenko Beam Stiffness Equation 125

also be easily obtained by finite element tech- expressed uniquely by the four conditions. In
nique from energy principles. The use of the order for this to be true, there must be three
displacement function of Eq. (27), though it relations between v and A. Since homogeneous
satisfies all the requirements for the finite element solutions are of interest, if the distributed force
method, does not give Eq. (22), but gives an terms are set equal to zero, there result
approximate equation, Eq. (29). One reason
A'+v"=0, EIA"-GkA (A+v') =0
for this is that the generalized strain y and K
can not simultaneously take arbitrary constant .......(57a, b)
values. This has been improved in the displace- which become the two constraint conditions.
ment function of Eq. (31), however, the use of From the order of the polynomials of the homo-
Eq. (31) also does not result in Eq. (22). In order geneous solutions, one more condition can be
to obtain Eq. (22), a still higher order polynomial obtained by differentiating Eq. (57. a). Given
has to be used for the displacement function. these conditions, all seven coefficients can be
The reason for this may be the presence of the determined for four geometrical boundary con-
Eq. (8b), constraint relation between v and A. ditions. The resulting displacement function
For the case in which y is equal to zero, it may agrees with the homogeneous solutions of Eq.
be inadequate to treat both v and A as completely (54), and use of this function leads directly to
independent functions. In order to derive an Eq. (22)12).
accurate stiffness equation, it may be necessary Timoshenko beam theory plays an important
to select displacement functions for which in- role in dynamic problems. In deriving stiffness
dependence between v and A is lost when y equations for dynamic problems by finite element
equals zero. A displacement function, Eq. (34), technique, attention has to be paid to the selec-
was introduced adding the dummy displacement tion of displacement functions and kinematic
d, as a techniques to obtain this dependence. conditions at nodes. In numerical considerations,
In order to obtain Eq. (22) without introduc- simple examples are solved. In view of the fact
ing a dummy displacement function, it is neces- that, in the analysis of a continuous system as
sary to employ third and second order poly- a discrete system, better accuracy is generally
nomials for v and A, respectively. This can be seen obtained by using a stiffness equation with a
from Eq. (54). Separating v and A, Eq. (54) can larger number of degrees of freedom, numerical
be transformed into two independent equations solutions are obtained using Eq. (49) which has
EIA"+Py+m'=0 the largest degree of freedom of the three stiffness
EI (viv+Py"/GkA)+Py+n2'=0 equations discussed in this paper.
One report4) deals with free vibration of the
...... (56'a, b) Timoshenko beam by the finite element method,
Clearly, the homogeneous solutions for v and A in which the continuity condition is imposed on
are third and second order polynomials. Assum- shear deformation at nodes. This treatment was
ing forms of these functions that satisfy the four not explained, however, it may have been used
geometrical boundary conditions at both ends, because it is continuous when an exact solution
a stiffness equation is obtained that is equivalent of a continuous system is obtained for a prismatic
to Eq. (22) but in a different form. The four Timoshenko beam. There are only four kinds of
geometrical boundary conditions are not suffici- boundary conditions at a node of a Timoshenko
ent to determine the seven coefficients present beam elment as expressed by Eq. (15), in which
in these polynomials, and hence the displace- no continuity condition for y is present. Impos-
ment functions include three dummy displace- ing no continuity condition for y at nodes, Eq.
ments resulting in a stiffness equation consisting (49) is assembled for a whole beam. The solution
of seven simultaneous equations. With static of this system gives y's with different values at
condensation of these three dummy displace- both sides of a node even for a case of free
ments, the stiffness equation reduces to Eq. (22). vibration of a prismatic beam. This discrepancy,
In order to derive Eq. (22) directly, it is neces- compared with the physically exact solution, is
sary to express all seven coefficients of the poly- due to the fact that, contrary to the smooth
nomials for the four existing geometric boundary distribution of inertia force in a continuous
conditions. The fact that the necessary and system, the distributed inertia force of a discrete
sufficient number of boundary conditions to system is treated as an equivalent nodal inertia
solve the governing differential equation for force, as given by the second term on the left
Timoshenko beam are four, and hence the assign- hand side of Eq. (49). By this reason, it is natural
ment of four geometrical boundary conditions, that y is not continuous in a finite element solu-
two at each end, leads to a unique solution, tion based on Eq. (6). By contrast, imposing
implies that all of the seven coefficients can be a continuity condition for y gives a discrete
 126 T. IWAHUMA, M. Al and F. NISHINO

nodal force, corresponding to y in Eq. (47),

shows that g is dependent on y. Thus unnecessary
prescription of y has an unnecessary influence
on q.

5. CONCLUSION

Since the exact stiffness matrix for Timoshenko

beam elements can be obtained easily by solving
the governing differential equation, it is not
significant to derive it by the finite element
technique. In spite of the fact that the homo-
geneous solution of a prismatic Timoshenko
beam problem is given by a polynomial, the
finite element method using polynomials as
displacement functions does not lead to an exact
Fig. 2 Dimension of matrices and natural stiffness matrix unless special attention has been
frequency. paid to the selection of the polynomial functions.
Taking advantage of the known exact stiffness
matrix, discussions on the selection of displace-
system other than that based on Eq. (6). In ment functions have been presented that com-
order to see the numerical difference due to the pare the resulting stiffness matrices with the
difference in conditions of y at nodes, the natural exact matrix. The discussions of this paper may
frequency of a two span continuous beam is be of help in selecting displacement functions for
evaluated. A beam of total length 90 cm is problems for which no exact solutions are known,
supported at a one third point as shown in Fig. 2. and may be applicable directly to plate and shell
The same properties as given in Eq. (55) are problems with shear deformation.
used. Putting external loads equal to zero in The discussions and conclusions of this paper
Eq. (49), the equation for free vibration is, are summarized as follows:
kqD=m(o2mgD ........ (58) (1) In spite of the fact that the homogeneous
solutions for generalized displacement, v and A,
Numerical solutions of o are obtained under are polynomials of third and second order, the
two conditions, one imposing a continuity con- finite element solutions using the lower order
dition for y, as reported in the literature4, and polynomial displacement functions that have
the other without this condition. The accuracy sufficient degrees of freedom to satisfy geometrical
of the natural frequency is plotted in Fig. 2 boundary conditions and do not make strain
for this to matrix size. In view of the facts that energy identically equal to zero, converge to
the finite element method leads to an upper the exact solution of the continuous system.
bound solution for o and that any kinematic For this convergence, the displacement functions
constraint leads to a larger value of w, it is do not necessarily account for the fact that the
concluded from the results shown in Fig. 2 that generalized strains can take arbitrary constants
the accuracy of w is decreased by imposing con- simultaneously.
tinuity conditions for y. (2) As can be seen in the strain displacement
There exists some confusion in the literature relation of Eq. (8.b), the generalized displace-
on the treatment of boundary conditions. For ments, v and A, are dependent when y is equal to
Timoshenko beams, that is, for the continuum zero. Selection of displacement functions which
for which Eq. (1) is assumed as kinematic con- satisfy this dependence in addition to ordinary
straints, Eq. (13) is the boundary conditions requirements in finite element technique leads
compatible with the constraints. Hence, for to an exact stiffness matrix.
example, the boundary conditions for the free (3) While the lowest order polynomial dis-
end of a cantilever are V=0 and M=0, while placement functions yield a 4 x 4 stiffness matrix,
there is a report5> in which y=0 is added as the the function satisfying the dependence of (2)
third condition. In free vibration of a cantilever, yields a 5 x 5 stiffness matrix. In view of the
it is true that y is equal to zero at the free end approximate nature of the former stiffness
since shear force at the free end is equal to zero. matrix, it may be of advantage in numerical
It does not hold true, however, in a discrete computation to use the latter even though the
model, since the equivalent nodal force does number of degrees of freedom is larger than in
not necessarily vanish at the free end. A zero the former.
 On Derivation of Timoskenko Beam Stiffness Equation 127

(4) The most appropriate displacement func- (2) Stiffness Matrix and Equivalent Nodal
tions for a Timoshenko beam element are a third Forces by Eq. (27)
order polynomial for v, and a second order poly-
nomial for A, whose total number of degrees of GhA /l - GkA /2 - GkA /l
freedom is reduced from seven to four by con- EI/l+GkAl/3 GkA/2
sidering the restriction of Eq. (57). This is pos- GkA/l
sible only when the governing differential equa- Sym.
tion is simple and its homogeneous solution is
- GkA/2
a polynomial. In most cases for which the finite
- EI/l + GkA l/6
element method is the only feasible solution
procedure because of the complexity of obtaining GkA/2
analytical solution, the consideration stated in EI/l+GkA l/3
(4) above is impractical but the selection stated
in (3) may be suitable. fOT1 fl f2 f3 f4j

ACKNOWLEDGEMENT

This research reported here was sponsored by

a Grant-in-Aid for Scientific Research of the
Ministry of Education, Science and Culture.

6. APPENDICES
(3) Stiffness Matrix by Eq. (31)
(1) Analytically Exact Matrix of Eq. (22) and GkA/l -GkA/2 -GkA /l
Equivalent Nodal Forces
EI/l + GkA l/4 GkA /2
k2=
GkA/l
Svm.
12 -6l -12 -6l -GkA/2
(4+l2012 61 (2-12c)l2 -EI/l+GkA1/4
12 67
GkA/2
Sym. (4+12q)l2
EI/l +GkA l/4

(4) Stiffness and Mass Matrices by Eq. (48)

12-6-12-6 6 6
4 6 2 -3 -3
12 6 -6 -6
k-El l2
4 -3 -3
3+1/3c 3+1/6
Sym. 3+1/3
...... (A.6)

13/35 -11/210 9/70 13/420

1/105 -13/420 -1/140
13/35 11/210
vn=l2
1/105

Sym.

11/210 -13/420
-1/105 1/140
13/420 -11/210
1/140 -1/105
1/105 -1/140
1/105
 12/ T. IWAKUMA, M. AT and F. NISHINO

6/5 -1/10 -6/5 -1/10 3/5 3/5 291-299, 1975.

2/15 1/10 -1/30 -1/20 -1/20 6) Thomas, D. L.: Comments on "Finite element
6/5 1/10 -3/5 -3/5 model for dynamic analysis of Timoshenko
+ 12fr beam", J. of Sound and Vibration, Vol. 46,
2/15 -1/20 -1/20
pp. 285-290, 1976.
3/10 3/10
7) Pin Tong: Exact solution of certain problems
Sym. 3/10 by finite-element method, AIAA, J. 7, No. 1,
pp. 178.180, 1969.
8) Archer, J. S.: Consistent matrix formulations
for structural analysis using finite-element
techniques, AIAA, J. 3, pp. 1910-1918, 1965.
REFERENCES 9) Przemieniecki, J. S.: Theory of Matrix Structural
Analysis, McGraw-Hill Book Co., New York,
1) Kapur, K. K.: Vibrations of a Timoshenko 1968.
beam, using a finite element approach, J. of the 10) Nishino, F., A. Hasegawa and E. Natori: Thin-
Acoustical Society of America, Vol. 40, walled rectangular beams with shear deform-
pp. 1058.-d063, 1966. ation and cross sectional distortion, Proc. J SCE,
2) Nickel, R. and G. Secor: Convergence of consist- No. 248, pp. 25-40, April 1976 (In Japanese).
ently derived Timoshenko beam finite elements, 11) Cowper, G. R.: The shearr coefficient in Timo-
International J. of Numerical Methods in En- shenko's beam theory, J. of Applied Mechanics,
gineering, Vol. 5, pp. 243-253, 1972. Vol. 28, pp. 335-340, 1966.
3) Davis, R., R. D. Henshell and G. B. Warburton: 12) Oliveira, E. R. A.: Theoretical foundations of
A Timoshenko beam element, J. of Sound and the finite element method, Intenational J. of
Vibration, Vol. 22, pp. 475487, 1972. Solids and Structure, Vol. 4, pp. 929952, 1968.
4) Thomas, D. L., J. M. Wilson and R. R. Wilson: 13) Yamamoto, Y. and V. Yamada: Error in Matrix
Timoshenko beam finite elements, J. of Sound Structural Analysis, Lecture Series II-5-3 on
and Vibration, Vol. 31, pp. 315-330, 1973. Matrix Structural Analysis, Baihukan, 1972
5) Thomas, J. and B. A. H. Abbas: Finite element (In Japanese).
model for dynamic analysis of Timoshenko (Received April 14, 1980)
beam, J. of Sound and Vibration, Vol. 41, pp.