METHOPS IN APPLIED MECHANICS AND ENGINEERING 20 (1979) 105-131

pUBLISHlNG COMPANY

ON TfrE GEOMETRICAL STIFFNESS OF A BEAM IN SPACE-

A CONSISTENT V.W. APPROACH

J.H. ARGYRIS, 0. HILPERT, G.A. MALEJANNAKIS, D.W. SCHARPF*

Insbtur fiir St@ik und Dynamik der Lufr- und Raumfahrtkonstruktionen, Unioersity of StuUgarr.
Germany

Manuscript received 26 September 1978

The application of the standard virtual work expressions to the large displacement-small strain
domain merely requires the replacement of the standard linear strain-displacement relations by the
quadratic ones. In fact. the geometrical stiffness matrix of an arbitrary finite element can be derived
timed+tely fropl the virtual work of the second order terms in the strains.
A correct gqqetrical stiffness is, however, only obtained if the prerequisites of the energy theorems
are strictly observed. This proves to be a rather difficult task as soon as the elements contain rotational
freedom. In @is case, the solution demands a comprehensive understanding of the true nature of the
strains, stresses and nodal displacements, rotations as well as forces, moments,
qfter an extensive study of the relevant entities the above principle is successfully applied to the
derivation of the geometrical stiffness of a beam element in space. The consistency of the present
approech is degloastrated by the full agreement with the prior results of the authors based on the
natural mde technique [2.3].
Spme numerical examples demonstrate the practical importance of the present development for the
geometrically ponlinear analysis of three-dimensional frame structures.

I. In$ruduction

in this paper we tackle once more the problem of rotational freedoms in finite elements
subject to large displacements- but small strains [2], [3]. In particular, we discuss the questions
which arise when applying the standard energy theorems to finite rotations. In what follows we
present a method for the derivation of the geometrical stiffness matrices of arbitrary elements
which is based 01) the standard virtual work expression of the linear theory in conjunction with
nonlinear strain-~&placement relations. To this purpose we require only the well-known
quadratic strain terms given in Love’s classical text book [S]. The virtual work of these
qvadratic strains directly yields an expression for the geometrical stiffness of the element.
Although the principle of the method is extremely simple, its application to elements
subject to finite rotations in space proves rather intricate. A correct geometrical stiffness of an
element specifically designed for general three-dimensional structures is only obtained if the
prerequisites for the validity of the virtual work principle are strictly observed. This entails in a
large displacement analysis a consideration not only of the initial position of the structure, but
also of tb$ displaced one.

lResearcij Consultant, Stuttgart.

106 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

In particular, the above-mentioned quadratic strain terms actually represent strains arising
in a displaced position of the system. Therefore it follows that the corresponding stresses,
entering into the respective virtual work expression, must also be referred to the displaced
position. For finite elements with rotational freedoms the latter statement leads to a depen-
dency of the nodal moments on the nodal rotations, since these moments must be understood
as resultants of the stresses which accompany the structure as it deforms.
This brief discussion demonstrates that the correct interpretation of the strains, stresses,
nodal displacements and loads is a most difficult task; it forms the main subject of the present
paper. For this reason we restrict the formulation of the theory to the simplest finite element
possible, which is evidently given by a straight cylindrical beam with a doubly symmetrical and
solid cross-section. This restriction appears acceptable for the purpose of this paper. A
subsequent publication will deal with thin-walled closed and open (cylindrical) beams. We also
note that most beams are anyhow not cylindrical and require an analysis as folded plate
systems, for which e.g. the TRUMP element [4] may be applied.
In the following section we first develop the basic relations required for the calculation of
the geometrical element stiffness. As mentioned above, these formulae immediately derive
from the virtual work of the quadratic strain terms. The respective quadratic strain-displace-
ment relations are applied in the third section to the beam element. The resulting quadratic
forms for the components of strain are critically analysed in section 4 and lead to two
important conclusions. Firstly, we observe that the quadratic terms of the strains arise in a
beam as a consequence of the assignment of the displacements to a cross-section subject to an
antecedent rotation. Secondly, it can be stated that the quadratic terms do not depend on a
precise definition of the rotations.
Section 5 presents the straightforward stress distribution in our simple beam element. These
stresses and, in particular, the resultant moments are carefully studied in the next section 6. It
is found thereby, that the moments resulting from the stresses on a cross-section are modified as
the cross-section is rotated. Moreover, it becomes evident that the bending and torsion
moments behave differently under rotations of the cross-section. Using the nomenclature of
Ziegler [7], we confirm that the bending moments are quasitangential but that the torsional
moment is semitangential.
Clearly, different definitions of nodal moments in finite elements are unacceptable in a large
displacement analysis of three-dimensional systems. This means that the geometrical stiffness
of a finite element with rotational degrees of freedom cannot be based exclusively on nodal
moments which correspond to the element stresses, but must observe the extended require-
ments of the large displacement theory.
The main difficulty in the actual evaluation of the geometrical stiffness is hence a correction
to the element nodal moments which is inserted into the fundamental virtual work expression.
This procedure leads to a symmetrical matrix which must be added to the matrix derived
directly from the virtual work of the quadratic strain terms. The numerical examples of the
final section extracted and developed further from [3] clearly demonstrate that this additive
correction to the geometrical stiffness is of great practical importance and completely indis-
pensable for beam structures in space.
Finally, we should explain briefly the relations between the present paper and some recent
publications. We mention first that similar methods for the derivation of the geometrical
stiffness matrices are described in [9] and [ll]. However, these publications do not consider
the principal problem of the present paper, viz. rotations about different axes in space. The
 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach 107

genuine difficulties due to the rotational freedoms are also avoided in [6], since assemblies of
beams with a common straight axis are only studied. On the other hand, this restriction
permits the consideration of general beam cross-sections, including warping, which becomes
most problematic in frameworks.
We should also mention here our prior papers [2] and [3], which are based on the natural
mode technique initiated in [l]. In contrast to [2], the present developments start from
considerations of strains and stresses and not of natural displacements and forces. The natural
mode technique is hence not an essential requirement of the present theory. Any natural
formulation used in the present paper serves only the purpose of abbreviating the presentation
and simplifying the comparison with our previous results.
Naturally, completely different solutions are possible for the problem of the rotational
freedoms. Whilst the above-cited references use rotations which agree at least in the limit of
infinitely small rotations with the standard engineering definition, the rotations can also be
described by Euler angles, which permit a unique definition of the rotated positions of a rigid
body. The reader is referred to [lo] for this method of solution.

2. Principle of the method

Consider a state of finite displacements and small strains in a deformable body idealized by
an assembly of finite elements. The derivation of the elastic as well as the geometrical stiffness
matrix of an arbitrary finite element can be based on the familiar virtual work expression

& ‘ud V = Sp’P. (2.1)

”

Here which we write form

The strains with the nodal displacements p whilst the

stresses with the nodal p and P
must be corresponding in order to yield a meaningful virtual work expression on the
right-hand side.
For finite displacements it is helpful to split the vector E into two components

E =&L+&Q (2.4)

where eL contains the linear terms and cQ the quadratic terms. Classical theory [5] yields for
the strain components, the linear terms

(2.5)
108 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach

and the quadratic terms

;[($y+($y+(gy], ;[(g)‘+(g+($)‘],
~[(~>z+(~~+(~)z],~[~~+~~+~~I; (2.6)
1 auau av*+awaw 1 auau av av aw aw
2 [&Z+2jat
_I-

I I
ay az 7 Z azax+ZZ+ZZ
--

I*
These define the strain vectors eL and eQ as functions of the displacement vector

u=b(x,YJ) VkYJ) W,YJN (2.7)

which derives from p via the assumed interpolation scheme. Using hence eq. (2.5) we obtain
the linear matrix relation

EL = ap. (2.8)
On the other hand, the set of eq. (2.6) leads to a quadratic form for each component of strain
EOijin EQ which can be written

EQij = $ p’pijp.
(2.9)

Introducing eqs. (2.8) and (2.9) into eq. (2.2) we find the typical strain component

(2.10)

and the typical virtual strain component

(SEij = aijSp + P'pijSp. (2.11)

Evidently (Yijdenotes the row of the matrix (Y corresponding to the strain sip
Substituting the expressions of eq. (2.11) for the virtual strains SE of eq. (2.1) and using also
the linear elasticity law

u = E& (2.12)

we immediately obtain the finite element force-displacement relation

(2.13)

with the elastic stiffness matrix

 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach 109

(2.14)
V

and the geometrical stiffness matrix

(2.15)

The sum of the two stiffnesses forms the tangent stiffness matrix

kT = kE + kc (2.16)

which is valid for an incremental loading step at the beginning of which the stresses uij are
acting on the structure.
The application of the above theory to finite elements is straightforward when the
displacement vector p does only contain translational degrees of freedom or at the most
rotations about a single fixed axis as in the case of beam elements in a plane. If, however, the
element freedoms include rotations about at least two different axes in space, the virtual work
expression of eq. (2.1) must be examined most critically as demonstrated in the following
sections.
Here it suffices to note that the precise definition of the rotations and moments entering
into the element matrices p and P is an indispensable prerequisite for the correct derivation of
the geometrical stiffness kc. In order to emphasize this statement it may be mentioned that the
unqualified application of eq. (2.15) leads to errors up to 60% even in extremely simple
structures as shown in fig. 7.7 of [2].

3. Displacements and strains of the beam element

The correct application of energy principles to the derivation of geometrical stiffness

matrices can be demonstrated most exhaustively using the example of a straight cylindrical
beam element in space. As shown in fig. 3.1, the beam element is referred to a local Cartesian
system Oxyz in which the x-axis coincides with the axis of the beam whilst the axes Oy and
Or agree with the principal axes of the cross-section. To simplify the presentation we also
assume that the cross-section is solid and doubly symmetrical with respect to the planes Oxy
and Oxz.
The nodal freedoms of the beam illustrated in fig. 3.1 are collected in the element
displacement vector

P =h Ul WlrplXl Ghuzv2 wz(p2x2$2z). (3.1)

The displacement functions u of the beam element and in particular the matrices a and /3
introduced in eqs. (2.8), (2.9) can be expressed in a particularly simple form, if the Cartesian
displacement vector p of eq. (3.1) is first transformed into the rigid-body and natural freedoms
110 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

Fig. 3.1 Beam element in space.

(3.2)

(3.3)
with

112 * . . * .,,2. . . . .

* ,,2 . . . . . ,,2 . . . .

. . ,,* . . . . . 112 . . .

. . . ,,2 . . . . * l/2 * *

. . ,,( . . . . . -,,[ . . .

. _,,[ . . . . . ,,, . . . .
a, : _, . . . . . . . . . . .

. . . . . . . . . ..-I

. 2/l . . . 1 .-2/I. . . 1

. . .._. . . . . . 1 m

. . 2/1 . -1 . . . -2/l . -1 .

. ..-l.....l..
(3.4)

Since the natural mode technique is not an essential requirement of the theory developed
here, the reader is referred to section 4 of [2] for an elaboration of the rigid body and natural
modes, ejro and pN, respectively. Moreover, the subsequent formulae, eqs. (3.5), (3.6) for the
displacements can be considered as a sufficient interpretation of the alternative displacement
parameters poi, &Ji.
Keeping in mind the well-known beams modes, the reader may easily verify the displace-
ments
x
u(x)= pOl+-pN1 1

w=poz+wo6- ($-&N2+ ($-;)pN3 (3.5)

ww=po3-wo5- ($-$N4+(+$NS
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 111

for a point at a position x on the axis of the beam. Similarly, the rotations of the cross-section
at x become

(3.6)

The three-dimensional displacement field u(x, y, z) required in eqs. (2.5) and (2.6) derives
from eqs. (3.5) and (3.6) by application of the engineers’ theory of bending (ETB) and St.
Venant’s theory of torsion. In fact,

~(x,Y,~)=~(X)+XZ-~Y+~(Y,~)~

4x, Y, 4 = 4x)- w(x) (3.7)

w(x, Y, 4 = w(x) + YdX).

Here ~(y, z) is the St. Venant warping function of the cross-section which need not be
specified explicitly for the following developments. Computing next the beam strains due to
the above displacements, we observe that, for a beam only the three strain components E,,, E,,
and t;* forming the reduced vector

E = {E,, tiExy A&,,} (3.8)

are relevant. The corresponding stress vector reads

u = {a,,d&y &rJ. (3.9)

Hence we only require for the linear and nonlinear strains cL and Ed the respective three
matrices cu,,, (uXY,
(Y,,‘and P& &, /.L
When applying eq: (2.5) to the modal displacements of eqs. (3.5) to (3.7) we note that only
the natural freedoms p,,_yield linear strains rzL.Thus, eq. (2.8) reduces to

EL = %@N (3.10)
with

(3.11)
112 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

A further advantage of our transformation eq. (3.3) leading to rigid body and natural
freedoms becomes apparent when calculating the &-matrices. In fact, we observe that the
pure translations pal, po2 and po3 yie_‘d no strains at all. For this reason the first three rows and
columns of the (12 x 12)-matrices fiij referred to 6 contain only zeros. Eq. (2.9) may thus be
rewritten

(3.12)

with
(3.13)

Table 3.1 Matrices for the calculation of quadratic strain terms.

.
_--
. . . . - : ; rr-l_--j
11.. .
L__7

11 .
L__,

I. .
L__,
i>;_ -! \

I (31’Ij;j . . -<(3p-+)
L-_--4.
1 52 -@t2- ;) -5?
rymm. L_____
i (3t;2- 17 1(352- 9
L---i_,
1 v2+ t’
L--F

.
r--
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 113

The actual matrices but, /&, and bXz are reproduced in table 3.1. Clearly, the matrices & of
eq. (3.12) can also be used in eq. (2.15) and lead to a geometrical stiffness

(3.14)

which is referred to j5 of eq. (3.3). The Cartesian geometrical stiffness kG based on p of eq.
(3.1) derives then immediately from the congruent transformation

(3.15)

with the a,-matrix of eq. (3.4).

Note: The matrix fi of table 3.1 does not contain the contribution of (&/&)* which becomes
only relevant for large strains.

4. On the definition of the nodal rotations

Having derived the quadratic strain terms cQ of the beam element in space we must
examine whether these strains depend on the actual definition of the nodal rotations. To this
purpose, let us consider a typical quadratic term of the strains, viz

which corresponds to the entry “--i” in the (1,2)-position of the matrix bXy in table 3.1. Since
p04, ~05 represent rigid body rotations cp, x about the axes OX and Oy (see fig. 3.1) we may
rewrite eq. (4.1) for convenience in the form

1
Exy= -29x. (4.2)

This latter relation expresses that successive rigid body rotations about the different axes Ox, Oy
yield a shear strain.
The reader should note that the generation of a shear strain is a consequence of the
displacement functions of eqs. (3.7) which are Zinear in the rotations cp,x, J/L We observe,
moreover, that this linearity is not only presumed in our theory but is also contained in the
general formulae eq. (2.6) which derive from a linear displacement field.
Returning now to the shear strain of eq. (4.2) we note that as soon as any expression
involves products of rotations about different axes, the result may depend on the sequential
order of the rotations and consequently on the precise definition of the rotations. For an
extensive discussion of this problem the reader is referred to our previous papers [2], [3]. Here
it suffices to show that the strain E,, of eq. (4.2) does not depend fortunately on the sequential
order of the rotations rp and x even if they are assumed to be axial, i.e. about the fixed axes OX and
OY.
114 .T.H. Argyris et al., On the geometrical stifiess of a beam in space--a consistent V. W. approach

The reason for the occurrence of the shear strain can be explained in the following manner.
Let us consider a square element of area in the x-y-plane with the nodes x = +-1 and y = +l;
see fig. 4.1 a,b. Evidently, neither of the displacements due to 9 or x leads to a shear strain, if
these displacements are referred to the initial position of the element in the x-y-plane. If,
however, the deflections w = +x, respectively w = QY are referred to an element with- an
antecedent rotation then the in-plane components u = +spx, respectively u = +(px lead to the
shear deformation of eq. (4.2) as illustrated in fig. 4.1.

a) Rotation x after preceding rp

b) Rotation up after preceding ,y

Fig. 4.1 On the derivation of a typical quadratic strain term.

The essential point in our present development is that for either sequential order of the
rotations the same strain E,, is generated although the final positions are different. We
conclude from this fact, that a precise definition of the nodal rotations of the finite element is
nor a prerequisite of the appiicability of the nonlinear strain-displacement relations in eqs.
(2.6).
This result is certainly not surprising if we remember that eqs. (2.6) derive from linear
displacement functions and that the quadratic terms only arise as a result of the second order
terms in the Taylor expansion of the square root required for the computation of the length of
a line element; see also [5].
 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach 115

The independence of eqs. (2.6) from a precise definition of the rotational freedoms may be
considered to be an advantage of the derivation of the geometrical stiffness via eq. (2.15). On
the other hand this partial simplification may mislead an inexperienced .observer to believe
that the same is true for the nodal moments. As is shown in section 6, a precise definition
of the nodal moments is indispensable and enforces a correction of eq. (2.15) for the
geometrical stiffness.

5. Nodal forces and stresses of the beam element

We first define the element force vector P for the beam in space shown in fig. 3.1. As stated
in conjunction with eq. (2.1), P must correspond to the element displacement vector p of eq.
(3.1) and thus becomes

Here Xi, Yi, Zi describe the forces in the X-, y-, z-directions of node i and similarly Mxi,
Myi, A4,$the moments about these axes. We will not consider here any distributed loads on the
beam. Hence we must demand that the forces and moments in P form a self-equilibrating
system of loads.
Our next task is now to find a three-dimensional stress field a@, y, z) which is statically
equivalent to the nodal loads P. The latter condition simply means that the resultants of the
stresses c on the end cross-sections must be equal to the nodal forces and moments in I?
The above conditions can be satisfied most elegantly if we introduce the natural forces of
the beam (see also [l], [2])

PN={P Nl PN2 PN3 P N4 PN5 PN6 } (5.2)

which derive from P via

P,g = ALP
(5.3)
with

-112 * . . . .
. . . . . .
. . . * . .
. . . . * -10
. . .-l/2-1/2 *
112 - l *
A, = ,;‘I2
. * * . *
. . . . . .
. . . . . .
. . . . l 112
l . * Ii2 -112 l

’ -l/2 l/2 * * l

L
(5.4)
116 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

As in the case of the natural displacements p N we note also here that the application of
natural forces PN is not a necessary prerequisite of the present theory. In place of PN, to PNh
one could use any consistent set of six parameters describing a selfequilibrating system of
stresses of the beam in space. A particular advantage of the use of the PNi’Sis however, that
the results of this paper can be easily compared with those of our prior publications [2], [3].
We next specify the cross-sectional resultants of the beam in terms of the natural loads:
normal force

N=P,.,, (5.5)

shearing forces

s,=-$PpN3; s, =-$P,, (5.6)

torsional moment
T = P,.w, (5.7)
and finally the bending moments

M,=P,,+P,,; M,=-(P&P,,). (5.8)

The reader may easily verify that the resultant forces and moments of eqs. (5.5) to (5.8)
fulfill the internal equilibrium conditions for a beam not subject to distributed loads.
Moreover, in conjunction with eqs. (5.2) to (5.4), we can readily prove that the values of the
resultants at the ends of the beam are equal to the loads P. In accordance with engineers’
theory of bending and St. Venant’s theory of torsion, we can now set up the functions for the
essential stress components of the vector CTin eq. (3.9) viz.

(5.9)

= - 5PN3gxy,y- g Ph3gxy,++ PN&*y.T (5.10)

uxz = uxz,y + uxz,z+ uxr,T

= - $PN3&r,y - $PNSgxz,, + PN6gxz.T. (5.11)

In the latter relations A, I, and 1, stand for the area and the principal moments of inertia of
the cross-section. The functions gij,k = describe the distribution of the shear stresses
gij,k(y, t)

due to S,, S, and T over the cross-section and need not be specified explicitly, since they only
enter in the integrals of eq. (2.15). The integration over a cross-section yields either the stress
resultants S,,, S, or T, respectively, or contributes a zero value due to the assumed double
symmetry of the beam cross-section. The last statement applies in particular to the integrals
containing the derivatives w,~ and o,= of the warping function.
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 117

6. On the definition of the nodal moments

The nodal loads in P of eq. (5.1) must evidently be understood as the cross-sectional
resultants of the stresses specified explicitly in eqs. (5.9) to (5.11). We thus possess, fortunately,
precise information on the method by which the nodal forces and, in particular, the nodal
moments are actually generated. There remains the task of studying the true nature of the
nodal moments and of inquiring whether the definition of the moments implied by the given
stress distribution is appropriate for a finite element formulation in the large displacement
domain.
Let us first consider the bending moments M, and M, which result from the integrals over
the cross-section

MY = zu,, dA ; M, = - yo;, dA. (6.1)

I I
A A

The nature of these resultant moments emerges if we analyse their change under small
rotations of the cross-section. The stresses must be considered as distributed surface tractions
acting on the end cross-sections. This means that the surface forces resulting from the stresses
follow the rotations of the end cross-sections. This assumption is necessary in order to obtain
the correct virtual work with the quadratic strain terms eoij which in their own turn must be
understood to take place in a rotated position of the element of area; see fig. 4.1.
The effect of a small rotation cp is illustrated in fig. 6.1. An element of area dA at a position
(y, z) on which the infinitesimal force u,.. dA acts is moving under the rotation cp by 2, = -cpz
and w = (py. This motion entails the incremental moments

dM,d = -cpyuxx dA; dM,, = qczu,, dA. (6.2)

Integrating the latter expressions over the cross-section we obtain

(6.3)

Fig. 6.1 Change of the bending moments due to rotation cp of cross-section.

118 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W, approach

Since small rotations about the axes Oy and Oz evidently do not change the moments M,,
and M,, eqs. (6.3) describe the-total change in the bending moments.
The reader may easily verify that moments produced by the simple mechanisms of fig. 6.2
are subject to -the same changes as the bending moments of the beam, see eqs. (6.3).
Independently of the actual cross-section we may hence interpret the bending moments as
being produced by a couple of forces acting on a rigid lever which is orthogonal to the axis of
the beam as well as to the axis of the moment. The magnitude of the forces is then equal to
the value of the moment if the lever is of unit length, as assumed in fig. 6.2.

a) Bending moment Mu b) Bending moment M,

Fig. 6.2. Simplified representation of quasitangential bending moments.

Moments generated as shown in fig. 6.2 have been denoted by Ziegler [7] as quasitangential,
see also [2]. Thus, we have demonstrated the quasitangential nature of bending moments
resulting from the normal stress distribution of eq. (5.9).
Let us next analyse the behaviour of the torque

T =
I
A
(yuxz - my) dA

which derives from u,.. and a,, of eqs. (5.10) and (5.11). We first note that for arbitrary
cross-sections both components of shear stress yield the same contribution to the torque T, i.e.

I
A
yu,, dA = -
A
zo-,, dA = ; T. (6.5)

This can easily be proved by partial integration of eq. (6.4) in conjunction with the equilibrium
relation

(6.6)
which is valid for shear stresses due to pure torsion.
As in the case of the bending moments we consider the displacements of an infinitesimal
element of area within a cross-section which is subject to small rotations. In contrast to the
case of the bending moments, a rotation cp does not change the torque T. Fig.. 6.3 illustrates
the effect of a rotation x about the Oy-axis. The displacement u = xz is seen to yield only the
incremental moment
dT,, = xzuzY dA (6.7)
 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach 119

Fig. 6.3 Change of torque due to rotation x of cross-section.

since the infinitesimal force u,.. dA is rotated together with the cross-section. Similarly, a
rotation $ leads to

dT,‘, = +yux,, dA. (6.8)

Observing eq. (6.3, the integration of dT,,, dTzd produces the total incremental moments

T,, = ; $T; Ll = +T. (6.9)

As in the case of the bending moments we can easily construct a simple mechanism which
generates a torque changing under rotations as postulated in eqs. (6.9). The respective
mechanism consists of two rigid levers perpendicular to the beam axis and is illustrated in fig.
6.4. Ziegler [7] has introduced the notation semitangential for such moments as the torsional
ones on the beam.
Having now studied the true nature of the three beam moments, it becomes immediately
apparent that a different behaviour of bending and torsion moments proves unacceptable in a
large displacement finite element analysis. Imagine, for example, two beam elements forming

Fig. 6.4 Simplified representation of semitangential torsional moment.

120 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

a right angle. A torsion of one of the elements leads to bending of the other. If the common
node of the two elements is subject to a rotation the moment in one element obeys eqs. (6.3)
whilst the same moment in the other element obeys eqs. (6.9). A similar effect turns out even
if both elements of the right angle are subject to pure bending, since the assumed levers for
the generation of the bending moments are always perpendicular to the beam axes and hence
different in elements meeting at an angle; see also fig. 2.7 of [2].
There follows from the above remarks that a thoughtless application of eq. (2.15) yields a
geometrical stiffness which must entail serious errors when applied to general beam systems in
space.

7. The geometrical stiffness of the beam element

After the preceding discussion of the strains and stresses in a beam element and the
resulting consequences on the definition of the displacement and force vectors of the finite
element there remains the evaluation of the stiffness matrices. We first note that the elastic
stiffness can be deduced without any difficulty from eqs. (2.14) (3.3) and (3.10). A reproduc-
tion of the well-known result is not called for here.
The derivation of the geometrical stiffness clearly requires as a first step the evaluation of
eq. (2.15) or preferably in our formulation that of eqs. (3.14) and (3.15). We only present here
the result of the straightforward integration of eq. (3.14) leading to iG in table 7.1.
Instead we concentrate our attention to the necessary corrections of eqs. (2.15) and (3.14)
which follow from the considerations of the preceding section 6. Remembering our extensive
discussion in [2] we state that in a consistent large displacement theory only two definitions of
the element nodal rotations and moments are permissible, viz. the axial and the semitangential
one. Since the former of these possibilities leads to inconvenient non-symmetrical element
stiffness matrices, it suffices to consider in this paper the semitangential definition.
In accepting the latter specification we face no problems with respect to the torsional
moment, the semitangential nature of which has been recognized in the preceding section.
On the other hand, we must find an appropriate substitute for the quasitangential bending
moments. In order to achieve this, we consider the general matrix relations describing the
dependence of the quasitangential moments MY and M, on the rotations of the levers by
which they are generated. For the moments illustrated in fig. 6.2 we deduce from eqs. (6.3) the
relations

M; = M’: -I- M,,A;q; M; = M: + M,A:q (7.1)

Here M;, Mi are the modified (quasitangential) moments which are obtained if the levers of
the initial moments

MO,= (0 My 0); MO,={0 0 M,} (7 -2)

are rotated through

Y, = b x $1. (7.3)
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 121

The quasitangential nature of the moments is characterized in eq. (7.1) by the matrices

A:$ a 81; A+; i ;].

(7.4)

Using the A-matrices of eqs. (2.29) in [2] we can also formulate matrix relations cor-
responding to eqs. (7.1) which hold for semitangential moments M& AC, viz.
M;=M;+M&;q; M; = My+ M,A;+o (7.5)
with

n;= [H i J; A:;[! i EJ
(7.6)

Resolving eqs. (7.5) for MO,and Mt, respectively, and introducing the resulting expression in
eqs. (7.1) leads to
M; = M; + M,A~~; Mf = M; + M,Ayq

00; 0-;
0
(7.7)

q=A;-A;=
000; Ay=h;-A;=
-;
00.
where

;001 [ I [ 0 0 0
Thus, we have established in eqs. (7.7) and (7.8) relations which permit the replacement of
(7.8)

the quasitangential moments Mz, Mz by the semitangential ones MC, Ms in conjunction with
an additive correction proportional to the rotation 4p.
The same substitution of the quasitangential moments is clearly also possible in the element
force vector P. We remember that as a consequence of the virtual work expression eq.
(2.1) the loads P must be resultants of the stresses u. Hence in eq. (5.1) the moments MyI, Mzl,
My2 and Mz2 must be understood as quasitangential, but M,, and Mx2 as semitangential
moments.
Introducing a new element load vector P” based on purely semitangential moments we
immediately derive from eqs. (7.7)
P = P” + k,p (7.9)
in which
kM = rO3 M,& + M,,A: 03 M&y + Mz2A:j. (7.10)
We use here the symbol kM since this matrix establishes a relation between the element load
and displacement vectors as does the element stiffness matrix. In this context the reader
should note that the constituent matrices AqSof eqs. (7.8) and kM are symmetrical as it must be
the case for finite element stiffness matrices based on conservative moments.
The element load vector P of eq. (2.1) can now be replaced by that of eq. (7.9) and we
obtain

SE’O d V = Sp’(P” + k@). (7.11)

I
V
122 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

Since neither the linear nor the quadratic strain terms, E= and cc,, respectively, depend on a
precise definition of the nodal rotations, the element displacement vectors p and 8~ remain
unchanged. For this reason we can also insert eq. (7.9) directly into eq. (2.13) which leads to
p”+k,p=&+k&. (7.12)
The matrices kG and kM can now be combined into
ks, = kc - k,,,,. (7.13)
We may understand k& as a geometrical stiffness matrix based on semitangential moments
and rotations. The new geometrical stiffness k & satisfies hence all demands set up in [2] and
can be successfully applied to general systems of beams in space.
There remains only to reproduce in table 7.2 the final geometrical stiffness matrix k& It is
preferable thereby to replace Myl, Mzl, MY2,Mr2 of eq. (7.10) by the expressions
My, = - (Pm + Pm); Mz, = Pm + Pm; My2 = Pm - P,w; Mzz = Pm - Pm (7.14)
which immediately follow from eqs. (5.8). Due to the latter relations the geometrical stiffness
k,$ can be written as a function of the natural forces PN1 to PN6 as in [2]. The reader may easily
verify the agreement between the present geometrical stiffness k& in table 7.2 and the
corresponding matrix of table 5.2 in [2] but for one difference arising from the additional
entries
I
k4.4= kto,lo = -k,,,o = Iy;IL PrJ,=GPN’ (7.15)

in the new k& This difference can be easily explained by the fact that in [2] we restricted

Table 7.1 Geometrical stiffness matrix & of straight beam referred to rigid-body and natural freedoms.
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 123

ourselves to solid cross-sections in which the above terms are evidently ignorable as long as
only small strains are considered. Since the present paper does not account either for the effect
of warping the new theory is seen to apply strictly to solid cross-sections only. The entries
listed in eq. (7.15) can hence be ignored readily, which leads to complete agreement with our
previous analysis in [2]. As stated in section 1, we hope to return to the subject with an extension of
the theory to arbitrary thin-walled closed and open profiles.

Table 7.2 Cartesian geometrical stiffness matrix kb of straight beam in space.

8. Applications

Examples illustrating the theory have been reported previously in [2], [3]; some additional
information is presented here.
In the first example we consider the linear initial buckling and the geometrically nonlinear

Geometrical Data:

Length I = 240 mm

Width b = 30 mm

lhicknorr t = 0.6mm

Section A-A: Material D&n:

Young’s Mod”lUS E= 71240 N/mm2

Poisson’s Ratio Y = 0.31

Loading:

Fig. 8.1 Right-angle frame under end load. Data of the problem.
124 J.H. Argpis et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

Idealization TRUMP BEAM

Number of 86 20
elements

Number of 384 120

unknowns

Fig. 8.2 Right-angle frame under end load. Idealizations.

Lwd P positive

BEAM

Number of p,, [Nl

Y1 elements -
per ride
% %
x
E- ’1 21 0.5531
0.5604 I
1 .llOl
.I754

P- : j ,“::::: ::r?
6 0.5509 1 .0897

10 0.5507 1.0880

TRUMP

Pcr = 1.1628 IN1

Locd P negative
BEAM

Number
elements
of P,,INl 1
per I ide
kG -r %
0.4269 0.7085

0.4227 0.6873

0.4220 0.6833

0.4217 0.6819

0.4215 0.6809

0.4214
1 0.6804

TRUMP

PC, = 0.7005 [ N1

Fig. 8.3 Right-angle frame under end load. Convergence of linear critical load.
 J.H. Argyti et al., On the geome~a~ stifness of a beam in space-a consistent V. W. approach 125

1.1 Linear buckling F& = 1.088 N

._.-_-. -.-_---

ith perturbation load

Reference load :po : 1N Line A -8: Nullification of perturbation lard

Reference perturbation Imd : Ps, i 1O-3 N

0.05 0.1 0.15 0.2 0.25 uoj’

0, 1 , I I I ,
0 10 20 30 40 50 60 u&mm1

Fig. 8.4 Right-angle frame under end load. Beam idealization. Load-displacement diagram of point@.

Zero perturbation load

1.3-

1.2-

linear buckling P,, = I.088 N

with perturbation lwd

Reference iwd :p, =lN Line A -8: Nullification of perturbation load

Reference perturbation load: 5, I IO-~ N

0.05 0.1 0.15 0.20 0.25 - W@/i

1 I I I
I I , , I *
0 10 20 30 40 50 60 -w [mm]
CD

Fig. 8.5 Right-angle frame under end load. Beam idealization. Load-displacement diagram of point 0.
126 J.H. Argyris et al., On the geometrical stifiess of a beam in space-a consistent V. W. approach
Section A -A:

Dimensions and elastic properties

as in Fig. 8.1

Loading cases :
I. Concentrated load Pu along axis of symmetry

2. End moments Mz at points @ and @

of opposite sign

Fig. 8.6 Symmetrical right-angle frame.

Symmetrical modes only considered

mi
Fig. 8.7 Symmetrical right-angle frame. Idealizations (half frame).
 J.H. Argyris et al., On the geometrical stiffnessof a beam in space-a consistent V. W. approach 127

postbuckling behaviour of a right-angled frame. The frame is built-in at one end and subjected
to a transferse load P at the free end (fig. 8.1). The structure is alternatively idealized with
beam and TRUMP elements (fig. 8.2). Fig. 8.3 shows the results of the linear buckling analysis
for positive and negative directions of the load. We observe that agreement between the
TRUMP and beam analyses is only achieved when the correct geometrical stiffness k& of eq.
(7.13) is used.
In the geometrically nonlinear case we restrict the analysis to the discretization with beam
elements and the correct matrix k& Figs. 8.4 and 8.5 show two typical deflection curves, viz.
the u-displacement of the free end (point 0) and the w-displacement of the corner point 0.
In order to overcome the difficulties in the immediate neighbourhood of the critical load PC,,
we apply a perturbation load P, = P/1000 perpendicularly to the plane of the frame. After
reaching the position A of the deflection curve, which lies beyond PO the perturbation load
P, is removed whereupon the structure takes the new equilibrium position at point B. Starting
from the latter position positive and negative load increments Pd are applied. Whilst positive
Lwd Pu positive

BEAM
t pu Number of P,,cr[N]
elements
per side
kc 6
2.4509 4.0474

2.4228 3.9595

2.4327 3.9535

2.4201 3.9507

2.4197 3.9483

TRUMP
2.4194
I 3.9469

P
Y,Cr = 3.9459 [ N]

Lwd Pu negative

BEAM

Number of Py.cr[N]
elements
per side
%
177.698 F 18.4558

29.1519 14.2362

22.1760 12.7725

20.6852 12.2646

19.7744 11.9167

TRUMP
19.3260
I 11.7437
J

P
y, cr = 12.6943 [ Nl

Fig. 8.8 Symmetrical right-angle frame. Concentrated load PY on symmetry line. Convergence of critical load.
128 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach

increments lead to the proper postbuckling domain, negative increments cause an asymptotic
approach to the initial buckling.
Structure, loading and boundary conditions of the right-angled frame of the second example
(fig. 8.6) are symmetrical. The geometrical and material data are identical with those of the
preceding system. For the numerical analysis one half of the frame is idealized with beam and
TRUMP elements, as shown in fig. 8.7.
In the first loading case the system is subjected to a single load P, acting in the axis of
symmetry. First the linear critical loads P,,., are determined for positive and negative directions
of the load. Fig. 8.8 presents some results from which we deduce a satisfactory agreement
between the TRUMP and beam discretization when kk is applied in the latter case. Interes-
tingly enough, even for the relatively high ratio b/t = 50 of the present frame, it is found that
the inclusion of the terms in eq. (7.15) modifies the buckling load of not more than 0.2%.
In the second loading case the frame is loaded by a pair of moments at the ends of the legs.
The results of the linear buckling analysis (fig. 8.9) demonstrate satisfactory agreement
between the beam and TRUMP models when the geometrical stiffness kb is used. The
theoretical result is quoted in [8] as

MC, = % m = 622.2081 [Nmm]

BEAM

I
Number of
I
-r s
elements
per ride M Yl,~%,C.

I
M
zI,cr z2,cr

I 2267.10 343.62 882.62

2 1102.43 322.71 686.08

3 1028.53 318.72 650.88

4 985.38 317.31 638.30

6 953.86 316.31 629.34

TRUMP
IO 937.84
I 315.79 624.77

M = 624.36 [ Nmm 1
z1,cr

M z2,cr = 632.50 [Nmm 1

Analytical

MC, = 622.21 I Nmm 1

Fig. 8.9 Symmetrical right-angle frame. Moment M, at ends of legs. Convergence of critical moment.
 J.H. Argyris et al., On the geometrical stiffness of a beam in space-a consistent V. W. approach 129

and is independent of the angle enclosed by the legs. On the other hand, fig. 8.10 demonstrates
the marked dependence of the critical moment on the angle (Yif the incomplete geometrical
stiffness kc of eq. (2.15) is applied. The maximum deviation (50%) from the correct result
occurs at (Y= 90” and demonstrates impressively the necessity of introducing the corrected
stiffness matrix of section 7.
In the geometrically nonlinear analysis of the frame we use the idealization with 10 beam
elements and substitute a prescribed rotation $ in place of a moment at point 0. In this way, it
proves possible to observe the behaviour of the structure up to a full rotation of $ = 360”.
Point @ on the axis of symmetry undergoes also a full rotation of p = 360”. In order to
overcome the critical behaviour in the initial and final phases of the calculations, we use
geometrical and load imperfections. By imposing a geometrical imperfection in the neigh-
bourhood of point @ only the first singularity can be overcome (for a more detailed discussion
see [3]). We concentrate therefore our attention here to the application of a perturbation load.
In the initial stage of the geometrically nonlinear analysis the singularity at position B in fig.
8.11 is overcome by an incremental application of a perturbation load P, acting in the positive
z-direction at the point @ in conjunction with the prescribed rotation $ at the point 0. After
bypassing the first singularity the perturbation load is nullified and the analysis is continued
until the second singularity at position C in fig. 8.11 is reached. At this stage the value of the
applied moment exceeds the value of the linear MC,. This is a consequence of a higher order
buckling mode which takes place in the neighbourhood of the second singularity (instability).
Immediately prior to the second singularity (position C) a perturbation load P, = O.Ol[N] is

Fig. 8.10 Symmetrical right-angle frame. Dependence of critical moment on angle (Y.
130 J.H. Argyris et al., On the geomehical stifiess of a beam in space-a consistent V. W. approach

Reference load :Mo = IO3 Nmm

Reference perturbation load: Ps, = low2 N

Loading : Prescribed rotation at point @

y
Fig. 8.11 Symmetrical right-angle frame. Load-displacement diagram.

Ai&--%
A
Linear buckling M,_, = 624.77 Nmm
,’ @\ Symmetric
/
\
0
M
4>* x

Reference load : M, = IO3 Nmm

Reference perturbation load : Ps, = lO-2 N

Loading : Prescribed rotation V at point @)

A_
I e
270’ 3600
%

:.il_._._
____Lj Fig. 8.12 Symmetrical right-angle frame. Load-displacement diagram.

applied in the negative z-direction. Due to this additional loading the system moves to the
new equilibrium position D of fig. 8.11. Interestingly enough, the value of the moment A4 at
this instant lies below the linear MC,. Furthermore, the state of deformation corresponds to the
first buckling mode associated with the linear critical moment MC,.
 J.H. Argyris et al,, On the geomehical stifiess of a beam in space-a consistent V. W. approach 131

After reaching the position D the prescribed rotation is increased beyond $ = 360” and at
the same time the perturbation load decreased down to zero. But for small accumulation errors
the system achieves the initial unloaded state A.
To supplement the information the value of the reaction moment M at the point 0 is
plotted in fig. 8.12 against the prescribed rotation @.

References

[l] J.H. Argyris, Continua and discontinua, Matrix methods in structural mechanics, Opening address, Proceed-
ings of the Conference on‘Matrix Methods held at Wright-Patterson Air Force Base, Ohio (26-28 Oct. 1965).
[2] J.H. Argyris, P.C. Dunne, G.A. Malejannakis and D.W. Scharpf, On large displacement-small strain analysis
of structures with rotational degrees of freedom, Comp. Meths. Appl. Mech. Eng. 14 (1978) 401-451; 15 (1978)
99-135.
[3] J.H. Argyris, H. Balmer, J.St. Doltsinis, P.C. Dunne, M. Haase, M. Kleiber, G.A. Malejannakis, H.P.
Mlejnek, M. Miiller and D.W. Scharpf, Finite element method- the natural approach, Opening address,
Conference on Finite Elements in Nonlinear Mechanics (FENOMECH ‘78) Stuttgart (30.8.-2.9.1978) and
Comp. Meths. Appl. Mech. Eng. 17/18 (1979) l-106.
[4] J.H. Argyris, P.C. Dunne, G.A. Malejannakis and E. Schelkle, A simple triangular facet shell element with
applications to linear and nonlinear equilibrium and elastic stability problems, Comp. Meths. Appl. Mech.
Eng. 10 (1977) 371-403, 11 (1977) 97-131.
[5] A.E.H. Love, A treatise on the mathematical theory of elasticity, 4th edition, Dover Publications, N. York.
1944.
[6] R.S. Barsoum, R.H. Gallagher, Finite-element analysis of torsional and torsional-flexural stability problems,
Int. J. Num. Meth. Engng. 2 (1970) 335-352.
[7] H. Ziegler, Principles of structural stability, 2nd edition, Birkhiuser Verlag, Base1 und Stuttgart, 1977.
[8] S.P. Timoshenko, J.M. Gere, Theory of elastic stability, 2nd edition, McGraw-Hill (1961).
[9] J.T. Oden, Calculation of geometric stiffness matrices for complex structures, AIAA-Journal, 4 (1966)
1480-1482.
[lo] J.F. Besseling, Derivatives of deformation parameters for bar elements and their use in buckling and
post-buckling analysis, Comp. Meths. Appl. Mech. Eng. 12 (1977) 97-124.
[ll] G.A. Malejannakis, Lineare und geometrisch nichtlineare Analyse von allgemein versteiften Platten und
Schalen auf der Basis der “Natural Mode Technique”, Dr.-Ing. thesis, University of Stuttgart, 1977.