.
8
The TL
Bar Element
8–1
Chapter 8: THE TL BAR ELEMENT 8–2
TABLE OF CONTENTS
Page
§8.1. Introduction 8–3
§8.2. The Two-Dimensional Bar Element 8–3
§8.2.1. Element Kinematics . . . . . . . . . . . . . . . 8–3
§8.2.2. Strain Measure . . . . . . . . . . . . . . . . 8–4
§8.2.3. Stress Measure . . . . . . . . . . . . . . . . . 8–7
§8.2.4. Total Potential Energy and Residual Force Equations . . . . 8–8
§8.2.5. The Tangent Stiffness Matrix . . . . . . . . . . . . 8–8
§8.3. FEM Coding Using Mathematica 8–10
§8.3.1. The Example Structure . . . . . . . . . . . . . . 8–10
§8.3.2. Forming the Internal Force . . . . . . . . . . . . . 8–11
§8.3.3. The Equilibrium Equations . . . . . . . . . . . . 8–14
§8.3.4. Plotting the Equilibrium Paths . . . . . . . . . . . . 8–15
§8.3.5. Having Some Fun: Snapshots and Animation . . . . . . 8–15
§8.3.6. Forming the Tangent Stiffness Matrix . . . . . . . . . 8–18
§8.3.7. Critical Point Study . . . . . . . . . . . . . . . 8–18
§8. Exercises . . . . . . . . . . . . . . . . . . . . . . 8–25
8–2
8–3 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT
§8.1. Introduction
In the present Chapter the key concepts of nonlinear continuum mechanics reviewed in Chapter 7
are applied to the development of governing equations of bar (truss) elements based on the Total
Lagrangian (TL) kinematic description. For brevity these will be referred to as “TL bar elements.”
There are two ways to construct TL elements:
1. The Standard Formulation (SF)
2. The Core Congruential Formulation (CCF).
The first method is easier to describe and will be presented in this Chapter through examples. The
second one is more flexible and powerful but it is more difficult to teach because it proceeds in
stages. As such it will be relegated to Chapters 10-11.
§8.2. The Two-Dimensional Bar Element
The element developed in this Chapter is a prismatic bar element that can be used to model pin-
jointed plane truss structures of the type sketched in Figure 8.1. These structures undergo large
displacements and rotations but their strains are assumed to remain small so that the material
behavior stays in the linear elastic range. These assumptions allows us to consider only geometric
nonlinear effects.
A two-node bar element appropriate to model members of such truss structures is shown in Figure
8.2. The element moves in the (X, Y ) plane. In the reference (base) configuration the element has
cross section area A0 (constant along the element) and length L0. In the current configuration the
cross section area and length become A and L , respectively. The material has an elastic modulus
E that links the axial-stress and axial-strain measures defined below.
Because thisChapter deals primarilywith the formulation of an individual element, the identification
superscript (e) will be omitted to reduce clutter until assemblies are considered.
The element has four node displacements and associated node forces. These quantities are collected
in the vectors
u =
uX1
uY1
uX2
uY2
, f =
fX1
fY1
fX2
fY2
, (8.1)
The loads acting on the nodes will be assumed to be conservative.
§8.2.1. Element Kinematics
In accordance with bar theory, to describe the element motion it is sufficient to consider a generic
point of coordinates X located on the longitudinal axis of the reference configuration C0. That point
maps to point x in the current configuration C. The bar remains straight in any configuration. These
coordinates can be parametrically interpolated from the end nodes as
X(ξ) =
[
X (ξ)
Y (ξ)
]
=
[ 1
2 (1 − ξ)X1 + 12 (1 + ξ)X2
1
2 (1 − ξ)Y1 + 12 (1 + ξ)Y2
]
x(ξ) =
[
x(ξ)
y(ξ)
]
=
[ 1
2 (1 − ξ)x1 + 12 (1 + ξ)x2
1
2 (1 − ξ)y1 + 12 (1 + ξ)y2
] (8.2)
8–3
Chapter 8: THE TL BAR ELEMENT 8–4
Motion
Reference configuration�
(same as base)
Current configuration
Figure 8.1. A plane truss structure undergoing large displacements
while its material stays in the linear elastic range.
Here ξ is the usual isoparametric coordinate that varies from −1 at node 1 to +1 at node 2. The
displacement field is obtained by subtracting the foregoing position vectors:
u(ξ) = x(ξ) − X(ξ) =
[
uX (ξ)
uY (ξ)
]
=
[ 1
2 (1 − ξ)uX1 + 12 (1 + ξ)uX2
1
2 (1 − ξ)uY1 + 12 (1 + ξ)uY2
]
, (8.3)
which expressed in matrix form is
u(ξ) =
[
uX (ξ)
uY (ξ)
]
=
[ 1
2 (1 − ξ) 0 12 (1 + ξ) 0
0 12 (1 − ξ) 0 12 (1 + ξ)
]
uX1
uY1
uX2
uY2
= N(ξ) u.
(8.4)
The element kinematic defined by these equations is depicted in Figure 8.3.
§8.2.2. Strain Measure
As discussed in Chapter 7, in the Total Lagrangian (TL) description the Green-Lagrange (GL)
strains and the second Piola-Kirchhoff (PK2) stresses are frequently used as conjugate measures in
the formulation of the internal energy. The only GL strain that appears in the energy expression is
8–4
8–5 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT
X, x
Y, y
X1
u
u
u
u
Y1
X2
Y2
f X1
fY1
f X2
fY2
1(X1, Y1)
2(X2, Y2)
1(x1, y1)
2(x2, y2)
Reference configuration
Area A , length L
C0
PK2 stress ,s0 GL strain e0 = 0
Current configuration C
PK2 stress ,s GL strain e
0 0
Area A, length L
�
Figure 8.2. The geometrically nonlinear, two-node, two-dimensional bar element
in Total Lagrangian description. This element may be used to
model members of a plane truss such as that shown in Figure 8.1.
the axial strain e1 ≡ e which is most expediciously defined using the length change as
e = L
2 − L20
2L20
, (8.5)
rather than through displacement gradients. Because of the linear displacement assumptions (8.4)
the strain e is constant over the element.1
This expression can bemaneuvered into a matrix function of the node displacements as follows. Let
X21 = X2−X1, Y21 = Y2−Y1, aX = X21/L0, aY = Y21/L0, uX21 = uX2−uX1, uY21 = uY2−uY1,
umX = (uX2 − uX1)/2 and umY = (uY2 − uY1)/2. Some of these quantities can be geometrically
interpreted as illustrated in Figure 8.4. Then
1 This is in fact the only use of the displacement interpolation (8.4) in the following derivations.
8–5
Chapter 8: THE TL BAR ELEMENT 8–6
X
X
, x
Y, y 1(X1, Y1)
2(X2, Y2)
1(x1, y1)
2(x2, y2)
X(ξ) =
{ X (ξ)
Y (ξ)
}
=
{ 1
2 (1 − ξ)X1 + 12 (1 + ξ)X2
1 (1 − ξ)Y + 1 (1 + ξ)Y
}
2 1 2 2
x(ξ) =
{
x(ξ)
y(ξ)
}
=
{ 1
2 (1 − ξ)x1 + 12 (1 + ξ)x2
1
2 (1 − ξ)y1 + 12 (1 + ξ)y2
}
u(ξ) = x(ξ)−X(ξ) =
{
uX (ξ)
uY (ξ)
}
=
{ 1
2 (1 − ξ)vX1 + 12 (1 + ξ)vX2
1
2 (1 − ξ)vY1 + 12 (1 + ξ)vY2
}
Bar longitudinal axis
C0
C
_
Figure 8.3. The definition of displacement field for the
two-dimensional TL bar element.
L2 = (X21 + uX21)2 + (Y21 + uY21)2,
e = L
2 − L20
2L20
= 1
L0
(aXuX21 + aY uY21) +
1
2L20
(
u2X21 + u2Y21
)
= 1
L0
[−aX −aY aX aY ]u +
1
L20
[−umX −umY umX umY ]u
= (Bl + Bn(u))u.
(8.6)
Observe that the GL strain e has been separated into two parts: e = el + en , where
1. el = Blu, where Bl is constant, depends linearly in the node displacements u. This is called
the linear part of the strain.
2. en = Bnu, in which Bn is a function of the node displacements, depends quadratically on
the displacements. This is called the nonlinear part of the strain, because Bn(u) vanishes if
u → 0. Note that both matrices are constant over the element.
8–6
8–7 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT
X, x
Y, y 1(X1, Y1)
2(X2, Y2)
C0
C
X21
Y21
L0
L
ψ
ψ0
aX =
X21
L0
= cosψ0
aY =
Y21
L
= sinψ0
0
ay =
y21
L0
= L
L0
sinψ1(x1 = X1+ u uX1, y1 = Y1+ Y1)
u u2(x2 = X2+ X2, y2 = Y2+ Y2)
ux21 = X21 + X21
uy21 = Y21 + Y21
ax =
x21
L0
= L
L0
cosψ
X
_
Figure 8.4. Geometric interpretation of quantities used in the study of element kinematics.
The variation of e induced by variations δ u in the node displacements is
δe = Bl δu + δ(Bnu) = B δu, (8.7)
The matrix B that links δe to δu is
B = ∂(B� + Bnu)
∂u
= 1
L0
[−ax −ay ax ay ] (8.8)
in which
ax = aX +
uX21
L0
= x2 − x1
L0
= x21
L0
, ay = aY +
uY21
L0
= y2 − y1
L0
= y21
L0
. (8.9)
For a geometric interpretation of ax and ay see Figure 8.4.
8–7
Chapter 8: THE TL BAR ELEMENT 8–8
§8.2.3. Stress Measure
The stress measure conjugate to GL strains is the second Piola-Kirchhoff (PK2) stress tensor. The
only component that appears in the internal energy is the axial stress s, which is related to e through
the constitutive equation
s = s0 + Ee, (8.10)
where s0 is the axial stress in the reference configuration, and E is the elastic modulus.
The axial force based on this stress is
N = A0s (8.11)
Note that this is not the true axial force in the current configuration C , which would be
Ntrue = Aσ (8.12)
where σ is the true or Cauchy stress in C and A is the actual area.
§8.2.4. Total Potential Energy and Residual Force Equations
In what follows it is assumed that the element is subjected only to node forces f that are conservative
and proportional, so that f = λq. The Total Potential Energy (TPE) of the element in the current
configuration is
� = U − P =
∫
V0
(s0e + 12 Ee2) dV0 − fT u =
∫
L0
A0(s0 e + 12 Ee2) d X¯ − λqT u. (8.13)
where X¯ is directed along the bar longitudinal axis in C0, as shown in Figure 8.4. This energy
expression is separable because the internal energy depends only on u through e and not on λ.
The finite element residual equations are obtained by differentiating this equation or, equivalently,
making (8.13) stationary with respect to virtual displacement variations δ u:
δ� = δU − δP = ( p − f)T δu = 0. (8.14)
The first variation of U is given by
δU = p T δu =
∫
L0
A0(s0 δe + Eeδe) d X¯ =
∫
L0
A0 s δe d X¯ =
∫
L0
N0B δu d X¯ = NL0B δu
(8.15)
where N = A0s is the axial force in the current configuration measured per area of the reference
configuration.2 It follows that the internal force vector is
p = N0L0BT = N
−ax
−ay
ax
ay
. (8.16)
This equation admits of a simple geometric interpretation; see Figure 8.5. The relation between
N0 = A0s and the true axial force N = Aσ can be worked out from inspection of this diagram.
The load potential variation δP simply generates f = λq as can be expected.
2 This is the PK2 axial force; cf (8.11).
8–8
8–9 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT
1
2
N0
fX1 = −N0ax
fY1 = −N0ay
fX2 = N0ax
fY2 = N0ay
Figure 8.5. Geometrical interpretation of the internal force vector.
The axial force N0 = A0s would be positive as shown.
§8.2.5. The Tangent Stiffness Matrix
Because the residual equations are separable the tangent stiffness matrix is obtained simply by
differentiating the internal force with respect to the node displacements u:
K = ∂p
∂u
= ∂(NL0B
T )
∂u
= A0L0BT ∂s
∂u
+ A0L0s ∂B
T
∂u
= KM + KG . (8.17)
The above expression shows that K splits naturally into two parts: KM and KG , which are called
the material stiffness matrix and geometric stiffness matrix, respectively, in the FEM literature.
To get KM note that
∂s
∂u
= ∂(s0 + Ee)
∂u
= E ∂e
∂u
= EB. (8.18)
Consequently
KM = E A0L0BT B. (8.19)
Inserting the expression (8.8) for B yields
KM = E A0L0
a2x axay −a2x −axay
axay a
2
y −axay −a2y
−a2x −axay a2x axay
−axay −a2y axay a2y
(8.20)
This component of K looks formally similar to the stiffness matrix of a linear bar element,3 except
that B now depends on u. The dependence of KM on the material properties (here the elastic
modulus E) explains the name “material stiffness” given in the FEM literature.
3 To which it reduces if u = 0. In that case ax and ay become the sine and cosine of the angle ψ0 shown in Figure 8.4.
8–9
Chapter 8: THE TL BAR ELEMENT 8–10
The other component can be obtained by differentiating B with respect to the node displacements,
the result being a constant 4 × 4 matrix:
∂BT
∂u
= 1
L20
1 0 −1 0
0 1 0 −1
−1 0 1 0
0 −1 0 1
. (8.21)
Inserting this into (8.17) one gets
KG = NL0
1 0 −1 0
0 1 0 −1
−1 0 1 0
0 −1 0 1
. (8.22)
This component ofK depends only on the stress state in the current configuration, because N = A0s.
No material properties appear. Thus the name “geometric stiffness” applied to KG .4
Remark 8.1. Assuming that E , A0 and L0 are nonzero, the rank of KM is obviosly one because B is a 1 × 4
matrix. On the other hand the rank of the numerical matrix in (8.21) is 2 (because its eigenvalues are 2, 2, 0, 0).
Consequently KG has rank 2 if s is nonzero and 0 otherwise. Combining these results it can be shown that the
rank of K = KM +KG is 1 if the configuration is stressed and 2 otherwise. In other words, the rank deficiency
is 3 and 2, respectively. The implications of this property in the analysis of stability are studied later.
Remark 8.2. The addition of KG increases the bar stiffness if the current configuration is in tension (s > 0),
but it reduces it if the current configuration is in compression (s < 0). This is in accord with physical intuition.
The main effect of this stiffness is on the rotational rigid-body motions of the bar about the Z axis.
§8.3. FEM Coding Using Mathematica
This section illustrates the use of Mathematica as a quick way to do some simple problems with the TL bar
element. The use of a computer symbolic system as prototyping language has several advantages:
(i) The code is very compact, because of higher lever array notation. Because of its compactness, it can be
easily translated into other high-level languages such as, for example, Matlab or Maple.
(ii) Debugging is often quicker than with programming languages because of the interpretative nature of the
language front end. This is of course platform dependent; the implementation of Mathematica on a
Macintosh or PC is friendler than on a Unix workstation.
(iii) Graphic output is part of the language. On the Macintosh, graphics can be transported into documents
such as this one via PostScript and Adobe Illustrator.
(iv) The language can manipulate both algebraic (symbolic) and numerical expressions. The later may be
done exactly or in floating-point arithmetic. This flexibility is convenient for many situations, as the
examples below illustrate.
4 In the pre-1970 FEM literature, the name “initial stress stiffness” was used for KG by some authors.
8–10
8–11 §8.3 FEM CODING USING MATHEMATICA
��
��
��
��
α
S
H
X, x
Y, y
(2)(1)
1
2
3
uX
uY
fY = λ
E, A0E, A0
Here λ > 0 means crown�
load is applied upwards
Figure 8.6. An arch as a two-bar FE model. This model has two degrees of
freedom: uX and uY , which are the displacements of node 2. This
structure will be used as example for use of Mathematica in §8.3.
§8.3.1. The Example Structure
The two-spring arch studied in Exercise 6.3 will be used as a structure that illustrates the use of Mathematica
in conjunction with finite elements. The structure is shown in Figure 8.6. It is modeled with two TL bar
elements labeled (1) and (2), which join nodes 1-2 and 2-3, respectively. The modulus of elasticity and the
bar areas are unity. The arch geometry is defined by the span S and height H . The reference configuration is
assumed stress free.
Because nodes 1 and 3 are pinned to the ground, the only degrees of freedom are the horizontal and vertical
displacements of node 2. To avoid clutter these will be denoted in the sequel as uX and uY rather than uX2
and uY2, respectively. The vertical motion uY is positive upwards. The node-2 components of the applied
and internal force vector will be also denoted simply by ( fX , fY ) and (pX ,pY ), respectively. Displacement and
force components for nodes 1 and 3 will be explicitly removed from the governing equations.
The arch is loaded by a vertical force fY = λ applied at node 2, positive upwards. The applied horizontal
force fX will be always zero. Recall that in Exercise 6.3 only symmetric motions of the arch under a vertical
load were considered. That constraint effectively reduces the model to one degree of freedom: uY , and allows
only limit point behavior. When two degrees of freedoms are considered, however, a richer set of possibilities
emerges; notably lateral bifurcation under a vertical load.
§8.3.2. Forming the Internal Force
The Mathematicamodule FormIntForce2DTwoNodeBar listed in Figure 8.7 forms the internal force vector
of an individual 2D bar finite element handled by the TL description. The arguments of this module are
XY1 Coordinates {X1, Y1} of end node 1 in C0.
XY1 Coordinates {X2, Y2} of end node 2 in C0.
uX1 Displacements {uX1, uY1} of end node 1.
8–11
Chapter 8: THE TL BAR ELEMENT 8–12
FormIntForce2DTwoNodeBar[XY1_,XY2_,uXY1_,uXY2_,Em_,A0_,s0_]:=
Module[{X1,Y1,X2,Y2,X21,Y21,uX1,uY1,uX2,uY2,uX21,uY21,L0,L,
e,s,ax,ay,ifv},
{X1,Y1}=XY1; {X2,Y2}=XY2; X21=X2-X1; Y21=Y2-Y1;
{uX1,uY1}=uXY1; {uX2,uY2}=uXY2; uX21=uX2-uX1; uY21=uY2-uY1;
L0=Sqrt[X21^2+Y21^2]; L=Sqrt[(X21+uX21)^2+(Y21+uY21)^2];
e=(L-L0)*(L+L0)/(2*L0^2); s=s0+Em*e;
ax=(X21+uX21)/L0; ay=(Y21+uY21)/L0;
ifv=A0*s*{{-ax}, {-ay}, {ax}, {ay}};
Return[Simplify[ifv]]
];
pe=FormIntForce2DTwoNodeBar[{2,3},{5,7},{1,0},{1,0}, 20,12, 0];
Print["Internal elem force pe=",pe];
pe=FormIntForce2DTwoNodeBar[{2,3},{5,7},{1,0},{1,0}, 20,12, s0];
Print["Internal elem force pe=",pe];
Internal elem force pe={{0}, {0}, {0}, {0}}
-36 s0 -48 s0 36 s0 48 s0
Internal elem force pe={{------}, {------}, {-----}, {-----}}
5 5 5 5
Figure 8.7. A Mathematica module that computes the internal force vector p of an individual
2-node 2D bar element. Output from two test statements is shown.
uX2 Displacements {uX2, uY2} of end node 2.
Em Elastic modulus E .
A0 Cross sectional area A0.
s0 PK2 stress in C0
Two test statements, one with purely numeric arguments and the other one with one symbolic argument (s0
for the bar prestress), are shown in the figure along with their output.
Figure 8.8 shows two modules, AssembleMasterIntForceOfTwoBarArch and
MergeElemIntoMasterIntForce that, together with FormIntForce2DTwoNodeBar, form the internal force
vector p for the two-bar example structure. The arguments of the former are the element internal force vector
eif, the element freedom table eftab that maps element-level freedom to master-level freedoms and the
master internal force vector mif before the element is merged. The function returns the updated internal force
vector. The arguments of the second module are:
span Span S
height Height H .
Em Elastic modulus E .
A0 Cross sectional area A0.
uX,uY Crown displacements uX , uY .
8–12
8–13 §8.3 FEM CODING USING MATHEMATICA
MergeElemIntoMasterIntForce[pe_,eftab_,p_]:=
Module[{i,ii,neldof,fmaster}, fmaster=p;
neldof=Dimensions[eftab][[1]];
For[i=1, i<=neldof, i++, ii=eftab[[i]];
If [ii>0,fmaster[[ii,1]]+=pe[[i,1]]]
]; Return[fmaster]
];
AssembleMasterIntForceOfTwoBarArch[span_,height_,Em_,A0_,uX_,uY_]:=
Module[{f1,f2,fmaster},
fmaster=Table[0,{2},{1}];
f1=FormIntForce2DTwoNodeBar[{-span/2,0},{0,height},
{0,0},{uX,uY},Em,A0,0];
fmaster=MergeElemIntoMasterIntForce[f1,{0,0,1,2},fmaster];
f2=FormIntForce2DTwoNodeBar[{0,height},{span/2,0},
{uX,uY},{0,0},Em,A0,0];
fmaster=MergeElemIntoMasterIntForce[f2,{1,2,0,0},fmaster];
Return[Simplify[fmaster]]
];
ClearAll[Em,A0,S,H,uX,uY,f];
p=AssembleMasterIntForceOfTwoBarArch[2,2.5,10,.75,-0.4,0.25];
Print["Master Int Force p= ",p];
p=AssembleMasterIntForceOfTwoBarArch[S,H,Em,A0,uX,uY];
Print["Master Int Force p= ",p];
Master Int Force p= {{-0.5336499821957073}, {1.555758744891104}}
Master Int Force p= {{(4*A0*Em*uX*
(S^2 + 2*uX^2 + 4*H*uY + 2*uY^2))/(4*H^2 + S^2)^(3/2)},
{(8*A0*Em*(H + uY)*(uX^2 + 2*H*uY + uY^2))/
(4*H^2 + S^2)^(3/2)}}
Figure 8.8. Two Mathematica modules that assemble the internal force vector p
for the two-bar arch structure of Figure 8.6. Two test statements, one
numeric and the other symbolic, are shown along with output.
Two test statements are also shown in the figure. The first one is purely numeric and sets the following
properties for the structure:
S = 2, H = 2.5, E = 10, A0 = 0.75, uX = 0.4, uY = −0.25 (8.23)
and the module AssembleMasterIntForceOfTwoBarArch returns
p =
[ pX
pY
]
=
[−0.53365
1.55576
]
(8.24)
The second one specifies symbolic values for all of the arguments: S for span S, H for the height H , etc. The
result (dis