INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL.

14,961-986 (1979)
LARGE DISPLACEMENT ANALYSIS OF
THREE-DIMENSIONAL BEAM STRUCTURES
KLAUS-JURGEN BATHEt AND SAID BOLOURCHIt
Department of Mechanical Engineering, Massachusetts institute of Technology, Massachusetts, USA.
SUMMARY
An updated Lagrangian and a total Lagrangian formulation of a three-dimensional beam element are
presented for large displacement and large rotation analysis. It is shown that the two formulations yield
identical element stiffness matrices and nodal point force vectors, and that the updated Lagrangian
formulation is computationally more effective. This formulation has been implemented and the results of
some sample analyses are given.
INTRODUCTION
The possibility of practical static and dynamic nonlinear analysis of structures has during recent
years progressed substantially, due to the effective use of digital computers operating on finite
element representations of the structures. To enable general nonlinear analysis the develop-
ment of versatile geometric and material nonlinear finite elements is in much need, and among
these elements the use of an effective three-dimensional beam element is very important.
Since the first applications of computers to nonlinear analysis of structures, various nonlinear
beam elements have been presented.*-" The large number of publications on nonlinear analysis
of beam structures is, at least partially, due to the fact that various kinematic nonlinear
formulations can be employed, and that at this time it is not clear which formulation is most
effective. The difficulty of obtaining effective solutions is particularly pronounced in the analysis
of three-dimensional beam structures. Namely, considering a beam element it is noted that a
general three-dimensional nonlinear beam formulation is not a simple extension of a two-
dimensional formulation, because in three-dimensional analysis large rotations have to be
accounted for that are not vector quantities.
In the development of a geometrically nonlinear beam element, basically an updated
Lagrangian or a total Lagrangian formulation can be employed. 12,13 These formulations must be
implemented using appropriate displacement interpolation functions. Considering the choice of
these functions it is recognized that for a beam of constant cross-section in small displacement
analysis the Hermitian functions should be employed to interpolate the transverse bending
displacements, and linear interpolation must be used to interpolate the torsional and longi-
tudinal displacements. Therefore, in the search for a beam element that can undergo large
rotations (with small strains), it is natural to employ the same functions but referred to the beam
convected co-ordinate axes. In this way the usual beam kinematic assumptions are used referred
to the current beam geometry.
t Associate Professor of Mechanical Engineering.
$ Research Assistant.
0029-598 1/79/07 14-096 1$0 1 .OO
@ 1979 by John Wiley 8c Sons, Ltd.
Received 17 April 1978
961
962 K. J. BATHE AND S. BOLOURCHI
Considering the formulation of a large displacement beam element, once specific beam
assumptions have been made and the interpolation functions have been selected, basically the
same element stiffness matrices and nodal point force vectors should be calculated using any one
formulation. Therefore, the response predicted using different formulations must be the same, if
the same number of beam elements are employed to model a structure. Indeed, the choice for a
total Lagrangian or an updated Lagrangian formulation should be decided only by the relative
numerical effectiveness of the formulations. However, considering large displacement beam
formulations using the Hermitian interpolations to describe bending deformations and linear
interpolations to specify axial and torsional displacements, a moving co-ordinate formulation
appears quite natural. Namely, in a total Lagrangian formulation for large rotation analysis, the
fact that the different displacement components are interpolated using different order poly-
nomials establishes an interpolation directionality that requires special attention.
In another approach to formulate a beam element that includes large rotation effects, the
transverse displacements, axial displacements and the rotations are interpolated independently.
If the same interpolation functions are employed for all these kinematic variables, the problem
of interpolation directionality under large rotations does not arise in the total Lagrangian
f or m~l at i on. ' ~ However, to obtain the same accuracy as with the beam elements based on
Hermitian functions, in this element formulation about twice as many degrees of freedom are
needed. It can be concluded that, for straight beams, it is more efficient to employ the
conventional beam interpolation functions, but to formulate more general and curved beam
elements the independent interpolation of displacements and rotations is effective.
The objective in this paper is to present two consistent large rotation nonlinear three-
dimensional beam formulations: an updated Lagrangian (U.L.) and a total Lagrangian (T.L.)
formulation. The formulations are derived from the continuum mechanics based Lagrangian
incremental equilibrium equations." The beam elements are assumed to be straight, and the
conventional beam displacement functions are employed to express the displacements of the
elements in convected co-ordinates. In the paper the two formulations are evaluated, and it is
shown that the governing incremental equilibrium equations of the beam elements are identical
but that the updated Lagrangian-based element is computationally more effective. This element
is a very efficient three-dimensional nonlinear beam element. The element has been implemen-
ted for use in elastic, elastic-plastic, static and dynamic analysis and in the paper a few
demonstrative sample solutions are presented.
INCREMENTAL T.L. AND U.L. CONTINUUM MECHANICS FORMULATIONS
The beam element formulations are based on the general incremental T.L. and U.L. continuum
mechanics equations," which are briefly summarized below.
Consider the motion of a body in a fixed Cartesian co-ordinate system, as shown in Figure 1.
Assume that the solutions measured in the co-ordinate system corresponding to all time points
0, At, 2At, . . . , t are known. It is required to solve for the unknown static and kinematic variables
in the configuration at time t + At. In static analysis and implicit time integration the equilibrium
of the body at time t + At is expressed and used to solve for the static and kinematic variables
corresponding to time t + At. On the other hand, in explicit time integration, the equilibrium at
time t is employed to solve for the displacements at time t +
Total Lagrangian (T.L.) formulation
In the total Lagrangian formulation all static and kinematic variables are referred to the initial
configuration at time 0. Considering the equilibrium of the body at time t + At, the principle of
THREE-DIMENSIONAL BEAM STRUCTURES
963
CONFIGURATION
AT TIME t f A t
/ CONFIGURATION
AT TI ME t
x i = O X I + + UI
t + At , I = o x I + t t At UI
CON FI G U R AT I0 N
AT TI ME 0 I 1
u I = t + At U. - t U.
Figure 1. Motion of body in Cartesian co-ordinate system
virtual displacements gives
where
components
is the total external virtual work expression due to the surface tractions with
and body forces with components A&fk,
In equations (1) and (2), SUk is a (virtual) variation in the current displacement components
uk, &r+Ai ~i i is a (virtual) variation in the Cartesian components of the Green-Lagrange strain
tensor in the configuration at time t + At referred t o the initial configuration, and are the
Cartesian components of the 2nd Piola-Kirchhoff stress tensor in the configuration t + At and
[ +At
964 K. J. BATHE AND S . BOLOURCHI
measured in the configuration at time 0:
(3)
(4)
t +Ar 1 r+At i CAf
O&ij = 5( 0ui.j + Ouj, i + I f A; uk , i f +Ab Uk , j )
0
t +A& = - P 0 i +Af 0
11 t +At r+ArXi , k Tk f f +At Xj , l
P
0 0 t t A t
where t+AtXi,j = d xi/d
t + At .
effects."
tal decompositions are used:
xj and the t +ArTkl are the components of the Cauchy stress tensor at time
In dynamic analysis the body force components in equation (2) include the mass inertia
Since the stresses f+A$ij and strains are unknown, for solution the following incremen-
( 5 1
(6)
t +At &.. = 6s.. +, S. .
11 g1 11
t +Ar
OEij = 6Ei j + OEij
where ASij and
the configuration at time t. It follows from equation (6) that St +A; Ei j =
components can be separated into linear and nonlinear parts
are the known 2nd Piola-Kirchhoff stresses and Green-Lagrange strains in
The strain increment
OEij = Oe i j + 0Vij (7)
where
(8)
(9)
1
oeij = d(0ui.j +OUj , i ) + ( 6 uk . i 0 u k , j + 6 uk . j 0 Uk , i ) l
and
1
0 Vi i = d 0 u k . i 0 u k . i )
Finally, the constitutive relations with tensor components ,,CijrS can be used to relate incremental
2nd Piola-Kirchhoff stresses to incremental Green-Lagrange strains
&.. 11 = &.. ilrs OErs (10)
Using equations (5)-(lo), equation (1) can now be transformed to
Equation (1 1) is nonlinear in the incremental displacements ui, and can be linearized by using
the approximations OSij = OCijrs Oers and = SOeii. We thus obtain
which is a linear equation in the incremental displacements.
integration, equation (1) is used corresponding to time t.
Equation (12) is employed in static analysis or implicit time integration. In explicit time
Updated Lagrangian (U.L.) formulation
In the U.L. formulation the same incremental stress and strain decompositions as in the T.L.
formulation are employed, but all variables are referred to the configuration at time t, i.e. the last
known configuration. Thus, corresponding to equation (12), the linearized equilibrium equation
THREE-DIMENSIONAL BEAM STRUCTURES 965
is in the U.L. formulation
fCl,rs ters &e,, du + f ~ , , &q2, dv = +B - 1, T,, &err du
I,, I,
(13)
where the f ~ t l are the Cartesian components of the Cauchy stress tensor at time t ; tell and fqrf are
the Cartesian components of the linear and nonlinear strain increments, respectively, and the
tCLf,s are the components of the tangent constitutive tensor relating small strain increments to the
corresponding stress increments.
U.L. AND T.L. FORMULATIONS OF BEAM ELEMENT
The general three-dimensional straight beam element is formulated based on the continuum
mechanics theory summarized above. The element has two nodes with 6 degrees-of-freedom
per node, and can transmit an axial force, two shear forces, two bending moments and a torque.
Figure 2 shows a typical beam element.
oxp. + x 2
I
4
0
I , 2 = B E A M ELEMENT END NODES
3 = A UX I L I A RY NODE
) O t
X I X I
Figure 2. Schematic view of the three-dimensional beam element local co-ordinate axes
The element is assumed to be straight and of constant cross-section. It is assumed that plane
sections of the beam element remain plane during deformation, but not necessarily perpendic-
ular to the neutral axis, i.e. a constant shear is allowed. The element can undergo large
deflections and rotations, but small strains are assumed. Thus, the cross-sectional area and the
length of the beam element do not change during deformation.
The principal moment of inertia axes of the beam element define the local co-ordinate system
r, s, t, as shown in Figure 2. The two end nodes of the element, 1 and 2, plus a third auxiliary
966 K.J. BATHE AND S. BOLOURCHI
node, 3, are used to define these axes, where it should be noted that in the computations the r-s
plane is defined by nodes 1 , 2 and 3.
Incremental equilibrium equations
In equations (12) and (13) the incremental equilibrium equations of a body in motion are
given corresponding to the global co-ordinate system xi , 7 = 0 or t. Considering a typical beam
element it is more effective to first evaluate the finite element matrices corresponding to the local
principal axes Ti?i of the element (see Figure 3), and then transform the resulting matrices to
correspond to the global Cartesian co-ordinate axes prior to the element assemblage process.I6
AREA
A
AREA / CONFI GURATI ON
AT TI ME t
A
0-
1, XI
t +
s,
At L ENGT H L
ON FI GURATI ON
T TI ME t + A t
. ,
A
CONFI GURATI ON
_____) AT T I ME 0
o x , , t x , !+A+ x ,
0 t t t o i
X3 X 3 x 3
Figure 3. Motion of the three-dimensional beam element and its local co-ordinate axes shown in global co-ordinate
system.
THREE-DIMENSIONAL BEAM STRUCTURES 967
The finite element matrices corresponding to the axes T2i are simply obtained by measuring all
static and kinematic quantities in this co-ordinate system. These new quantities are denoted by a
bar placed over them. Thus, using equations (12) and (13) we obtain for a single beam element,
using the U.L. formulation and considering only static analysis
( :KL + :KNL)u = i+AiR - :F
(14)
and using the T.L. formulationand considering only static analysis
(bKL + ~KNL) u = ttAtR - bF
(15)
where AKL, :KL are linear strain incremental stiffness matrices; AKNL, :KNL are nonlinear strain
(geometric or initial stress) incremental stiffness matrices; R is the vector of externally
applied element nodal loads at time t + At ; AF, :F are vectors of nodal point forces equivalent to
the element stresses at time t ; and u is the vector of incremental nodal displacements.
In dynamic analysis using implicit time integration the inertia forces corresponding to time
t + At are added to the left-hand sides of equations (14) and (15), whereas in dynamic analysis
using explicit time integration the stiffness effect is not included, the inertia forces corresponding
to time t are added to the left-hand sides of equations (14) and (15), and the applied external
loads correspond to time t.
The element matrices in equations (14) and (15) are evaluated using the displacement
interpolation functions of the beam element. Table I summarizes these calculations. The
following notation is used in Table I with all quantities referred to the co-ordinate systems T%,
T=Oor t :
ABL, :BL = linear strain-displacement transformation matrices
:BNL, :BNL = nonlinear strain-displacement transformation matrices
i +Ai
oC, ,C = incremental stress-strain material property matrices
1- i r
r, r = matrix and vector of Cauchy stresses
AS,
= matrix and vector of 2nd Piola-Kirchhoff stresses.
Table I. Finite element matrices
968 K.J. BATHE AND S . BOLOURCHI
It should be noted that the elements of the stress matrices and vectors in the T.L. and the U.L.
formulations are equal, because small strain conditions are assumed.
Interpolation functions for incremental displacements
To describe the motion of the beam elements the incremental displacement field within the
elements as a function of the incremental nodal point displacement components is required,
k = l
where the ,hl are the interpolation functions corresponding to the local axes 'f,, and the rGk are
the nodal point displacement increments measured in the local axes at time t (see Figure 3).
The interpolation functions in equation (16) are constructed assuming cubic bending dis-
placement variations and a linear variation in the axial and torsional displacements. In order to
include shear effects, constant shear deformations can be included. Using the usual beam
incremental nodal displacements (these are shown for time 0 in Figure 2) and leaving the shear
deformations as independent variables, we obtain the incremental displacement interpolation
functions given in Table 11. In this table the variable N in equation (16) is equal to 12 if no shear
deformations are included; otherwise N is equal to 14. If shear deformations are included the
element stiffness matrices and nodal point force vectors are of order 14, and are reduced to order
12 by static condensation prior to the element assemblage process.
Strain-displacement transformation matrices in the U.L. formulation
The kinematic assumptions used in defining the interpolation functions of Table I1 hold for
small strains, small rigid body incremental rotations in each solution step, but any size
translational displacements. These assumptions are appropriate for the updated Lagrangian
formulation of beams, because the kinematic variables are linearized about the last-known body
position.
Using the interpolation functions in Table 11, the strain-displacement matrices of the U.L.
formulation can directly be evaluated. Table III(B) summarizes the calculation of the matrices
(Br. and :BNL that are required to evaluate the tangent stiffness matrix and nodal point force
vector of an element corresponding to co-ordinate axes 'f, ( i = 1, 2, 3) . The element matrices
have to be transformed to the global co-ordinate system prior to their assemblage into a system
of beam elements.
t -
Strain-displacement transformation matrices in the T.L. formulation
In the discussion of the total Lagrangian formulation and the comparative study of the total
and updated Lagrangian formulations we do not include, for clarity, the effect of shear
deformations. Referring to the definitions in Figures 2 and 3, the displacement increments
within the element at time t measured in the local axes at time 0 are related to the nodal point
displacement increments of the element in its local axes using
12
k = l
OJj = c "hl " Uk
where the " Uk are the element nodal point displacement increments at time t, but measured in
the "ii (i = 1, 2, 3) co-ordinate system. The functions oh: ( i = 1, 2, 3) are the interpolation
THREE-DIMENSIONAL BEAM STRUCTURES
Table 11. Beam interpolation functions
969
We define:
Incremental displacement interpolation matrix
r
L
0 *z t -*zs -
( 1- i ) s -*5L 0 0
,
0 --t
L
r
46
0 $6
where
L = length of the beam element
,hi = vector of interpolation functions in 'f, direction
r, s, t = beam convected co-ordinate axes (see Figure 3)
The incremental displacement vector is
r i i T =[,GI rU2 . . . r ~ 1 2 j ,1113 , ~ ' 4 1
where rU l 3 = p 1 ,d l 4 = p2 (shear deformations)
functions corresponding to the convected axes r, s, t and measured in the co-ordinate system O f i
( j = 1,2,3).
The interpolation functions are obtained using
970 K. J. BATHE AND S. BOLOURCHI
Table 111. Matrices used in beam analysis
A. TOTAL LAGRANGIAN FORMULATION
1. Incremental strains
- - -
2. Linear strain-displacement transformation matrix
Using
where
oG = vector of incremental nodal displacements measured in O f i ( i = 1,2,3) co-ordinate system
u = vector of incremental nodal displacements in global co-ordinate system
O R = transformation matrix
THREE-DIMENSIONAL BEAM STRUCTURES 97 1
Table 111-continued
where
4. 2nd Piola-Kirchhoff stress matrix and vector
symmetric] ; $= [ : 1] ASl2
0 i s 1 3
where 1 3 is a 3 x 3 identity matrix.
2. Linear strain-dispfacement transformation matrix
Using
$5 = :BL
where
and
972 K. J. BATHE AND S. BOLOURCHI
- -
1 -
71 1
f -
0 71 I symmetric
0 0 7 1 1
g-= f - 7 1 2 0
0 0 0 0
0 713 0 0 0 0 0
r -
0 0
0 0 0 0
t -
I -
713
f -
-
Table III-continued
where
ru = vector of incremental nodal displacements measured in the !fi ( i = 1, 2, 3) co-ordinate system
u = vector of incremental nodal displacements in the global co-ordinate system
'R = transformation matrix between the local co-ordinate system at time t and the global co-ordinate
system
rh;,1 h i . 1 . . .
(rh;.2+th:,1) (thi.z+rh;,l) .. .
(thl.3 +th:,r) (rhi.3 +hi , , ) . . .
a,h',
ar f j
N = 12 if shear effects are neglected
N = 14 if shear effects are included
fh ;.j = -
where
3. Nonlinear strain-displacement transformation matrix
:BNL =
where 'Rim is the element (i, rn) of the matrix 'R, which transforms displacements measured in
the co-ordinate system '.ti (i = 1, 2, 3) to displacements measured in the system ().ti (i = 1, 2, 3) as
defined in equation (22).
A typical derivative required in the calculation of the strain-displacement transformation
matrix is (a,ui/ao2j) = C k = , (dohL/aofj)oak. To evaluate these derivatives it should be noted that
the axes 02j ( j = 1, 2, 3) correspond to the convected co-ordinates axes r, s, t at time 0 (i.e.
0 - n -
12
x l - r ; x 2 = s ; ' 2 3 Et ) . Using Eq. (18) we have
Therefore, double transformations are needed for the strain calculations in the T.L. formula-
tion. In comparison, the U.L. formulation does not require the above transformation.
THREE-DIMENSIONAL BEAM STRUCTURES 973
In order to evaluate the strain increments it is also necessary to calculate the derivatives of the
total displacements. The kinematics of the rigid body rotations of the beam give
i = 1, 2, 3
bUi,i = 'R.. - 6..
1 1 1 I J ( j = 1, 2, 3 )
where the 'Rij are the direction cosines of the ' Xi axes with respect to the "Xi axes, as defined in
equation ( 23) , and Sii is the Kronecker delta.
Transformation between current and original beam co-ordinate axes
In the U.L. and T.L. formulations a transformation matrix 'R that relates displacements
measured in the current configuration to displacements measured in the original configuration is
needed.
The 'R transformation matrix is evaluated using Euler angles which define the rotations
of the beam. These angles are shown in Figures 4(a) and 4(b). To arrive at the information given
in this illustration it is required first to evaluate the relative translational displacements of nodes
' Q = ROTATION OF COORDINATE AXES ABOUT O x 2 AXI S
' p = ROTATION OF COORDINATE AXES ABOUT 7 AXI S
-
(OR, . O R 2 ,Ox,) TO ( I J 1, '% , 7 )
( GI Xp , 7 ) TO (7, ? , T )
0-
PLANE PI I S PERPENDICULAR TO PLANE P 2
Figure 4(a). Rotation of beam element co-ordinate axes in large displacement analysis (first stepj
974 K. J. BATHE AND S. BOLOURCHI
0.
\
' y = ROTATION OF COORDINATE AXES ABOUT 7 AXI S
(7, T, 7) TO ( ' 21, ' 2 2 , 'i?,)
PLANE P 3 IS PERPENDI CULAR TO 7 A X I S
Figure 4(b). Rotation of beam element co-ordinate axes in large displacement analysis (final step)
1 and 2 measured in the beam original co-ordinate system. Denoting for clarity nodes 1 and 2 as
nodes I and J, these relative displacements are evaluated as
where the ' u are the element nodal point displacements measured in the global co-ordinate
system, and the ORij are components of the matrix OR that transforms the global nodal point
displacements to the element local axes at time 0.
The components of the matrix 'R are then constructed from the direction cosines of the axes
',fi ( i = 1, 2, 3) with respect to the axes O X j ( j = 1, 2, 3) . We have
1 -
R =
where 'R is a matrix of order 3 X 3,
THREE-DIMENSIONAL BEAM STRUCTURES 975
In equation (23) 'Rd is the transformation matrix due to the relative translational displacements
of nodes J and I, and 'R" is the transformation matrix that takes into account the axial rotation
of the beam,
The components of the matrix 'Rd are the direction cosines of the axes ?' , g, ;with respect to '4
(i = 1, 2, 3) . These $omponents are
cos ('a) cos (' p) sin (' p) sin ('a) cos
-sin (' a) 0 cos (' a)
(24)
('a) sin (' p) cos ('p> -sin ('a) sin ('PI
where the angle 'a represents the rotation about the (negative) O f 2 axis
O L +&;I
cos(' a)= -
IJi
51= {(OL + OUJI ) + ( OUJ I ) }
and O L is the original length of the beam,
'-1 2 ' - 3 2 1f2
The angle ' p represents the rotation about the positive t" direction;
:us,
sin (' p) = 7
L
where
The components of the 'Ra matrix, which are the direction cosines between ' . C (i = 1, 2, 3) and
the ?, g, t" axes are computed using
(29)
O I
0
@" = 0 cos ( ' y) sin ( ' y)
' - " o -sin ( ' y) cos ('7)
where ' y is the rigid body rotation of the beam about the ?-axis in the configuration at time t. This
angle is calculated using
y =${'fi4+',1o}
= (<)i i 4 + o ~ ' o ) + '@?2 ( ( ) i i 5 + o ~ l ') + + oii'Z)) (30)
and then
(3 1)
t y = '--Aty +
Substituting the relations in equations (23)-(31) into equation (22) we obtain the transformation
matrix between the beam local axes at times t and 0.
Calculation of beam element stresses
In the development of the incremental U.L. and T.L. equilibrium equations corresponding to
time t, we assumed that the stress components corresponding to the configuration at time t are
known (see Table I). The solution of the incremental equations (14) and (15) will then yield
nodal point displacement increments, from which the corresponding stress increments must be
calculated. These stress increments are added to the stress components at time t to obtain the
stress components corresponding to time t + At.
976 K. J . BATHE AND S. BOLOURCHI
To evaluate the stress increments accurately it is realized that the tangent approximation in the
strain-displacement relation for the normal strain, as employed in the :BL and bBL matrices,
does not yield an increment in normal strain if the element deflects transversely without
bending. However, for large displacement analysis, the corresponding extension of the element
is taken into account in the incremental strains using in the U.L. formulation
where the iBL,, are the components of the linear strain-displacement matrix given in Table III(B)
and Si j is the Kronecker delta; also
2 1 1 = ('L -"L)/OL
Using the T.L. formulation the corresponding calculations are
where the constraints
should be imposed to evaluate the appropriate normal strains.
With the incremental strains known, the corresponding stress increments can be calculated as
~i sua1. l ~ In general large displacement and elastic-plastic analysis, the stiffness matrices and
nodal point force vectors must be evaluated using numerical integration. Also, to improve the
solution accuracy it may be necessary to employ equilibrium iterations in the incremental
solution. The iterative equations are directly obtained from equations (14) and (15) in the usual
manner.
12
COMPARISON OF T.L. AND U.L. FORMULATIONS
In the 'T.L. formulation the reference co-ordinate system used is given by the element principal
axes of inertia in the configuration at time 0, O X i (i = 1, 2, 3) . Therefore, the complete stiffness
matrix (including the linear and nonlinear strain stiffness matrices), the nodal point force vector
and the local displacement increments are referred to this co-ordinate system and must be
transformed to the global co-ordinate system "xi (i = 1, 2, 3) :
where OR is the transformation matrix that expresses the nodal point displacements measured in
the beam local co-ordinate system OXi (i = 1, 2, 3) in terms of the global nodal point displace-
ments.
The reference co-ordinate system used in the U.L. formulation is defined by the principal
axes of the beam element in the position at time t, i.e. Ifi (i = 1, 2, 3) . Therefore, the local
stiffness matrix and the nodal point force vector are referred to this co-ordinate system. These
THREE-DIMENSIONAL BEAM STRUCTURES 977
matrices are transformed to the global co-ordinate system using
where R is the transformation matrix relating the co-ordinate systems Xi and Xi (i = 1, 2, 3), as
defined in equation (22).
The principal difference between the U.L. and the T.L. formulations is that in the T.L.
formulation the transformation on the interpolation functions in equation (18) is carried out to
refer the displacement interpolations to the original configuration, and the ABLl matrix is
included in the calculations. The transformations on the interpolation functions and the use of
the ABL1 matrix in the T.L. formulation together are equivalent to the additional transformation
matrix R that is employed in equation (36) in the U.L. formulation. Indeed, as shown in more
detail in the Appendix, using these formulations the same element stiffness matrices and nodal
point force vectors are obtained.
Although the same final element stiffness matrices and nodal point force vectors are
generated in the two formulations, it is noted that using numerical integration the trans-
formation on the interpolation functions in equation (18) and the evaluation of the ABL1 matrix is
carried out at each integration point. Therefore, the U.L. formulation is computationally more
effective.
SAMPLE SOLUTIONS
The updated Lagrangian-based beam element was implemented in the computer program
ADINA17 and a number of sample analyses were carried out. We report here the results of some
of the analyses. In these analyses the beam linear strain stiffness matrices :KL were evaluated in
closed form, and the nonlinear strain stiffness matrices :KNL and force vectors :F (see Table I)
were evaluated using Newton-Cotes integrationI6. Also, in all analyses beam shear deformations
were neglected.
Large deflection analysis of a shallow arch
The clamped circular arch with a single static load at the apex was analysed for buckling using
the beam element, as shown in Figure 5. The material of the arch was assumed to be isotropic
linear elastic. One half of the arch was idealized using 6, 12 and 18 equal beam elements. The
same arch was also analysed using eight six-node isoparametric elements with 2 x 2 Gauss
integration.
This arch was also analysed by Mallet and Berke, who used four equilibrium-based
elements. Dupuis et al. analysed the same arch using curved beam elements, and used this
example to demonstrate the convergence of their Lagrangian and updated formulations,
Figure 5 shows the predicted load-deflection curve of the arch. It is observed that in this
analysis the use of the beam elements is quite effective.
Large deflection and rotation analyses of a cantilever beam
The objective in this analysis was to investigate the performance of the beam element in large
displacement and rotation problems. Two problems were analysed. First, a large deflection and
978 K. J. BATHE AND S. BOLOURCHI
-
n
c
n
a
0
0
J
R = I33 114 in
h = 3/16 in
b = I 0 i n (WI DTH)
L = 34 O i n
H = I 0 9 i n
B = 7 3397O
A = 0188 in2
I = 000055hn4
E = 10x l O61b/ i n2
I u = 0 2
FOR HALF OF ARCH
- 18 BEAM ELEMENTS
0 12 BEAM ELEMENTS
0 12 BEAM ELEMENTS
A 6 BEAM ELEMENTS
NO EQUILIBRIUM ITERATION
NO EQUILIBRIUM ITERATION
WI TH EQUILIBRIUM ITERATION
NO EQUILIBRIUM ITERATION
ISOPARAMETRIC ELEMENTS
NO EQUILIBRIUM ITERATION
x EI GHT 6NODE
VERTI CAL DI SPLACEMENT AT APEX w, [in.J
Figure 5 . Large deflection analysis of shallow arch under concentrated load.
moderate rotation analysis of a clamped cantilever with a concentrated end load was carried out
as shown in Figure 6. The second problem was the large displacement and large rotation analysis
of a cantilever beam subjected to a concentrated end moment (Figure 7).
In the cantilever analysis subjected to the concentrated tip load, the objective was to
demonstrate the effects of the aspect ratio of an element on its performance in the geometric
nonlinear range. Figure 6 shows the response predicted by ADINA using four different models
and an analytical solution. It is noted that the cantilever models using beam elements and
two-dimensional isoparametric elements (2 X 2 Gauss integration), with an aspect ratio of 2,
predict responses quite close to the analytical solution, and that the performance of the
THREE-DIMENSIONAL BEAM STRUCTURES 979
I SO PAR AM ETRl C ELEMENT BEAM ELEMENT MODEL
MODEL
BEAM ELEMENTS, +=4 8 + = 4 0
--- ISOPARAMETRIC ELEMENTS , $- = 2
- _ _ _ _ - ISOPARAMETRIC ELEMENTS , f = 100
LARGE DISPLACEMENT
AN A LY T I C A L S 0 L U T I 0 N
0.8 -
-I
3
' 0. 7
2 0. 6 -
a
(z 0. 5
2
0 0.4 -
t
0
w 0.3 -
-
I-
-
0 1.0 2 .o 3.0 4.0
0.8 -
-I
3
' 0. 7
2 0. 6 -
a
(z 0. 5
2
0 0.4 -
t
0
w 0.3 -
-
I-
-
LARGE DISPLACEMENT
AN A LY T I C A L S 0 L U T I 0 N
PL*
LOAD PARAMETER k = -
1
Figure 6. Large deflection analysis of a cantilever subjected to a concentrated load.
beam element does not change with its aspect ratio. However, as is well known, the predicted
response using two-dimensional isoparametric elements deviates from the analytical solution
with increasing element aspect ratios.
Figure 7 shows the results obtained in the analysis of the cantilever subjected to an end
moment. The cantiIever was modelled using 5 and 20 beam elements. The figure shows that
the predicted response compares well with the analytical solution up to 90 degrees rotation.'* It
is also seen that as the number of elements increases the numerically predicted response
improves. This increase in accuracy is due to the fact that the geometry of the deformed
cantilever is defined more accurately with a larger number of elements.
980
1.0
-slk
K. J. BATHE AND S. BOLOURCHI
-
L = 100 IN.
I = 0.01042 I N 4
A = 0.5 IN.*
M= CONCENTRATED
END MOMENT
E = 1.2 x 104PSI
0 .o
- -- ANALYTICAL SOLUTION
ADl NA , 90 STEPS, NO EQUI LI BRI UM I TERATI ON
- EOBEAM ELEMENTS
0 0 A 5 BEAM ELEMENTS
/
0.0 0.1 0. 2 0.3 0. 4
M L
MOMENT PARAMETER 7 = __
2 T E I
Figure 7. Moment-deflection curve.
Large displacement three-dimensional analysis of a 45 -degree bend
The large displacement response of a cantilever 45-degree bend subjected to a concentrated
end load, as shown in Figure 8, was calculated. The bend has an average radius of 100in,
cross-sectional area 1 in2 and lies in the X-Y plane. The concentrated tip load is applied into the
2-direction.
The bend was idealized using 8 equal straight beam elements and 16 sixteen-node three-
dimensional solid elements. For the beam elements the Newton-Cotes formula of order
3 X 3 X 3 was used and, for the isoparametric elements, Gauss integration of order 2 x 2 x 2 was
employed.16 The material was assumed to be linear elastic.
THREE-DIMENSIONAL BEAM STRUCTURES
t
98 1
NODE LAY-OUT FOR ONE OF THE
SI XTEEN EQUAL SOLID ELEMENTS
b Y
16 EQUAL SI XTEEN NODE SOLI D EL EMENTS
Z
R = 100"
8 EQUAL STRAI GHT BEAM EL EMENT S
-.
L
k 1 ' ' 4
B E A M
CROS S - SECTI ON
Figure 8. Finite element modelling of a 45-degree circular bend.
Figure 9 shows the tip deflection predicted by ADINA using the two finite element models. To
the accuracy that can be shown in the illustration, the same response is predicted using the beam
element idealization and the isoparametric element discretization. The deflected shapes of the
bend at various load levels are shown in Figure 10.
CONCLUSIONS
To develop capabilities for large displacement and large rotation analysis of beam structures, an
updated Lagrangian and a total Lagrangian formulation of a geometric nonlinear beam element
982 K.J. BATHE AND S. BOLOURCHI
t
FIX ED
, Y
X
T,
It
R=l OO. O IN
u= 0.
E - I 0 p s i
7
BEAM CROSS -SECTION
/
J
LL
u
04-
n
* LARGE DISPLACEMENT RESPONSE
_I 60 EQUAL LOAD STEPS
a
z 03- 8 BEAM ELEMENTS
I6 THREE-DIMENSIONAL ELEMENTS
z
0 0 10 2 0 3.0 4.0 5 0 6 0 7 0
P R 2
LOAD PARAMETER k __
E l
Figure 9. Three-dimensional large deflection analysis of a 45-degree circular bend.
have been presented. The incremental displacement fields within the straight two-noded beam
element are defined using the usual beam displacement functions. It has been shown that the two
formulations yield identical element stiffness matrices and nodal point force vectors, and that
the updated Lagrangian formulation is computationally more effective. This formulation can be
used efficiently for the general nonlinear analysis of beam structures.
ACKNOWLEDGEMENT
The work reported in this paper has been supported financially by the ADINA users group. We
would like to acknowledge gratefully this support.
THREE-DIMENSIONAL BEAM STRUCTURES 983
X J
Figure 10. Deformed configurations of a 45-degree circular bend.
APPENDIX: DETAILED COMPARISON OF BEAM T.L. AND U.L.
FORMULATIONS
In the text we showed that the U.L. formulation is more effective than the T.L. formulation for
beam analysis. The objective in this appendix is to compare the T.L. and U.L. formulations
presented in detail. Assume that the beam is deformed to the configuration at time t, as shown in
Figure 3. It is shown in this appendix that all element matrices are identical in both formulations.
Linear strain stiffness matrices
equation (20). Thus we have corresponding to Table III(A)
Consider first the T.L. formulation. The initial displacement effect is taken into account using
i = 1, 2, 3
j = 1, 2, 3 0
-
I.. = R. . - 8..
11 I I1 11 (37)
The Pij are the direction cosines of the & axes with respect to the Xi axes defined in equation
(23), and 8ij is the Kronecker delta. Using equation (37) the 6BL1 matrix defined in Table III(A)
is
1
(PI, - 1) ohfi + g2i oh:i + g31 oh!i
~ B L I = (811 - 1) oh,\ + 8 1 2 oh,*i + 821 ohfz + ( g22 - 1) oh;i + @310h:2 + 8 3 2 oh:1
i (g11 - 1) oh,; + 8 1 3 oh:i + 821 oh;3 + R 2 3 ohf1 + Rg i oh:3 + (g33 - 1) ohTi
(38)
984 K. J. BATHE AND S. BOLOURCHI
where (not considering shear deformations)
oh:= [ohi,j oh;,+ . . oh;z.j] (39)
Adding the JBLI matrix of equation (38) to the matrix JBLo defined in Table III(A) yields the
linear strain-displacement matrix,
IF::] (40)
t -
8 3 1 o o o o o
1 - 8 3 2 1 - ell 1 - 8 2 1 f - 8 3 1 0 0
g33 0 0 0 feii K21 K31 o H. 3
where we define
Th,i
TH,i=[;-] ( T =O, t )
The derivatives of the interpolation functions defined in equation (18) are
i = 1, 2, 3
j = l , 2 , 3 0
3 12
m 1 -
0hi.j = 1 c Ri m 1hn.j R n k
m= l n = l
where the incremental interpolation functions & are defined in Table 111. Equations (42) may
be rewritten in matrix form
where R and R are defined in equations (22) and (23), respectively. Substituting equations (43)
into equation (40) and simplifying gives,
(44)
In the U.L. formulation the geometric linear strain-displacement matrix is given in Table
III(B)
:BL = [ (:hi?h:2)] (45)
(1h;i + 1h.i)
Comparing equations (44) and (45) yields
4BL = :BLfR (46)
Substituting the above relation into Table I to evaluate the linear strain stiffness matrices in both
formulations, and comparing, we obtain
(47)
f - T i - t -
; KL= R tKL R
Therefore the two formulations lead to identical linear-strain stiffness matrices corresponding
to the global co-ordinate system.
THREE-DIMENSIONAL BEAM STRUCTURES 985
Nodal point force vectors
The components of the Cauchy stress tensor referred to the Xi axes are numerically equal to
the components of the 2nd Piola-Kirchhoff stress tensor referred to the OXi axes, i.e. the stress
vectors 4 and [as defined in Table 111) are equal. It follows from equation (46), Table I and the
above fact that the nodal point force vectors corresponding to the global axes are equal in both
formulations.
Nonlinear strain stiffness matrices
The nonlinear strain-displacement matrix in the T.L. formulation is defined in Table III(A):
Substituting equation (43) into (48) gives
The geometric nonlinear stiffness matrix is evaluated as defined in Table I,
0 0
where the 2nd Piola-Kirchhoff stress matrix is given in Table III(A). We also have
0 0 R 0 0
and the rh:2 and h:3 are null vectors. Thus equation (50) can be written as
6KNL = RT( lv;BcL $ :BNL du }a
where the matrix :BNL is defined in Table III(B), and
1 -
0s=
(49)
986 K. J. BATHE AND S. BOLOURCHI
The geometric nonlinear stiffness matrix based on the U.L. formulation is evaluated by using
the matrices of Table 111,
JV
Since the stress vectors and 68 are numerically identical for the beam element we have
bKNL = ,RT : K N L z R ( 55)
Therefore the two formulations lead to identical nonlinear strain stiffness matrices correspond-
ing to the global co-ordinate system.
REFERENCES
1. J. H. Argyris and P. C. Dunne, A simple theory of geometrical stiffness with application to beam and shell
problems, 2nd Znt. Symp. Computing Meth. Appl . Sci. & Engng, Versailles, France (1975); ZSD-Report No. 183,
Univ. of Stuggart (1975).
2. Z. P. Bazant and M. E. Nimeiri, Large-deflection spatial buckling of thin-walled beams and frames, ASCE, J. Eng.
Mech. Diu., 10247-10281 (1973).
3. C. Oran and A. Kassimaili, Large deflection of framed structures under static and dynamic loads, Comp. & Struct.
4. R. W. H. Wu and E. A. Witmer, Nonlinear transient responses of structures by the spatial finite element method,
5. T. Y. Yang, Matrix displacement solution to elastica problems of beams and frames, Znt. J. Sol. Struct. 9,828-842
6. C. Orean, Tangent stiffness in plane frames, J. Struct. Din, ASCE, 99 (ST6), 973-985 (1973).
7. T. Belytschko, L. Schwer and M. J. Klein, Large displacement transient analysis of space frames, Znt. J. Num. meth.
8. S. Ranganath and R. J. Clifton, Finite deflection dynamics of elastic beams, Znt. J. Sol. Struct. 10,557-568 (1974).
9. J. T. Holden, On the finite deflections of thin beams, Znt. J. Sol. Struct. 8, 1051-1055 (1972).
10. R. H. Mallet and L. Berke, Automated method for the large deflection and instability analysis of three-dimensional
truss and frame assemblies, AFFDL-TR-66-102, (1966).
11. G. A. Dupuis, H. D. Hibbitt, S. F. McNamara and P. V. Marcal, Nonlinear material and geometric behavior of shell
structures, Comp. & Struct. 1, 223-239 (1971).
12. K. J. Bathe, E. Ramm and E. L. Wilson, Finite element formulations for large deformation dynamic analysis, Znt. J.
num. Merh. Engng, 9,353-386 (1975).
13. K. J. Bathe, An assessment of current finite element analysis of nonlinear problems in solid mechanics, Symp.
Num. Solutions of Partial Difi Eqns, Univ. Maryland (1975).
14. K. J. Bathe, Finite element formulation, modeling and solution of nonlinear dynamic problems, in Numerical
Methods for Partial Differential Equations (Ed. S . V. Parter), Academic Press, 1979.
15. K. J. Bathe, Static and dynamic geometric and material nonlinear analysis using ADINA, Report 82448-2,
Acoustics and Vibration Lab., Mechanical Engineering Dept., MIT (1977).
16. K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, N.J.,
1976.
17. K. J. Bathe, ADINA - A finite element program for automatic dynamic incremental nonlinear analysis, Report
82448-1, Acoustics and Vibration Lab., Mechanical Engineering Dept., MIT (1977).
18. E. Ramm, A plate/shell element for large deflections and rotations, in Formulations and Computational
Algorithms in Finite Element Analysis, (Eds. K. J. Bathe, J. T. Oden and W. Wunderlich), MIT Press, 1977.
6, 536-547 (1976).
AI AA J. 11(8), 1110-1117 (1973).
(1973).
Engng, 11,65-84 (1977).