FEA TipsFinite Element Theory & MATLAB Programs
HomeThree-Dimensional Truss Structures2.3

Section 2.3 of Chapter 2: Three-Dimensional Truss Structures

Elastic and Geometric Stiffness Matrices

A direct comparison of two- and three-dimensional truss elements, showing the common structure of their elastic and geometric stiffness matrices and the complete second-order strain expansion from which they follow.

2D and 3D trussesSmall displacementSecond-order strainTheory

01

Two-dimensional truss element

We begin with the plane truss element and write its elastic and geometric stiffness matrices explicitly.

For a plane element, let c=cos⁡θc=\cos\theta and s=sin⁡θs=\sin\theta. Then

d={cs},p={−c−scs}\boldsymbol{d}=\begin{Bmatrix}c\\s\end{Bmatrix},\qquad \boldsymbol{p}=\begin{Bmatrix}-c\\-s\\c\\s\end{Bmatrix}

The elastic stiffness matrix is:

kel=EAL[c2cs−c2−cscss2−cs−s2−c2−csc2cs−cs−s2css2]\boldsymbol{k}_{el}=\frac{EA}{L}\begin{bmatrix}c^2&cs&-c^2&-cs\\cs&s^2&-cs&-s^2\\-c^2&-cs&c^2&cs\\-cs&-s^2&cs&s^2\end{bmatrix}

The transverse projector removes the component parallel to the element axis. For any vector v\boldsymbol{v}, the axial component is ddTv\boldsymbol{d}\boldsymbol{d}^{T}\boldsymbol{v} and the transverse component is Pv\boldsymbol{P}\boldsymbol{v}, where

P=I2−ddT=[s2−cs−csc2]\boldsymbol{P}=\boldsymbol{I}_{2}-\boldsymbol{d}\boldsymbol{d}^{T}=\begin{bmatrix}s^2&-cs\\-cs&c^2\end{bmatrix}

In two dimensions this projection leaves the direction perpendicular to the bar; in three dimensions it leaves the plane perpendicular to the bar. This is why the projector appears in the geometric stiffness matrix.

and the geometric stiffness matrix becomes:

kG,el=NL[s2−cs−s2cs−csc2cs−c2−s2css2−cscs−c2−csc2]\boldsymbol{k}_{G,el}=\frac{N}{L}\begin{bmatrix}s^2&-cs&-s^2&cs\\-cs&c^2&cs&-c^2\\-s^2&cs&s^2&-cs\\cs&-c^2&-cs&c^2\end{bmatrix}

This is the same matrix used in Sections 1.2 and 1.3.

02

Three-dimensional truss element

Let ll, mm, and nn be the direction cosines of the spatial element axis:

d={lmn},l2+m2+n2=1\boldsymbol{d}=\begin{Bmatrix}l\\m\\n\end{Bmatrix},\qquad l^2+m^2+n^2=1

The elastic stiffness matrix is:

kel=EAL[l2lmln−l2−lm−lnlmm2mn−lm−m2−mnlnmnn2−ln−mn−n2−l2−lm−lnl2lmln−lm−m2−mnlmm2mn−ln−mn−n2lnmnn2]\boldsymbol{k}_{el}=\frac{EA}{L}\begin{bmatrix}l^2&lm&ln&-l^2&-lm&-ln\\lm&m^2&mn&-lm&-m^2&-mn\\ln&mn&n^2&-ln&-mn&-n^2\\-l^2&-lm&-ln&l^2&lm&ln\\-lm&-m^2&-mn&lm&m^2&mn\\-ln&-mn&-n^2&ln&mn&n^2\end{bmatrix}

The geometric stiffness matrix is:

kG,el=NL[1−l2−lm−lnl2−1lmln−lm1−m2−mnlmm2−1mn−ln−mn1−n2lnmnn2−1l2−1lmln1−l2−lm−lnlmm2−1mn−lm1−m2−mnlnmnn2−1−ln−mn1−n2]\boldsymbol{k}_{G,el}=\frac{N}{L}\begin{bmatrix}1-l^2&-lm&-ln&l^2-1&lm&ln\\-lm&1-m^2&-mn&lm&m^2-1&mn\\-ln&-mn&1-n^2&ln&mn&n^2-1\\l^2-1&lm&ln&1-l^2&-lm&-ln\\lm&m^2-1&mn&-lm&1-m^2&-mn\\ln&mn&n^2-1&-ln&-mn&1-n^2\end{bmatrix}

This is the matrix used for the spatial buckling analysis in Section 2.2.

03

Common structure and matrix decomposition

The preceding two- and three-dimensional results have the same algebraic structure. Let r=2r=2 for a plane truss and r=3r=3 for a spatial truss. Denote by Ir\boldsymbol{I}_{r} the identity matrix of order rr, and define

d∈Rr,dTd=1,p={−dd}\boldsymbol{d}\in\mathbb{R}^{r},\qquad \boldsymbol{d}^{T}\boldsymbol{d}=1,\qquad \boldsymbol{p}=\begin{Bmatrix}-\boldsymbol{d}\\\boldsymbol{d}\end{Bmatrix}
B=1LpT,kel=EALBTB=EALppT\boldsymbol{B}=\frac{1}{L}\boldsymbol{p}^{T},\qquad \boldsymbol{k}_{el}=EAL\boldsymbol{B}^{T}\boldsymbol{B}=\frac{EA}{L}\boldsymbol{p}\boldsymbol{p}^{T}

Introduce the dimension-independent block matrix:

Qr=[Ir−Ir−IrIr]\boldsymbol{Q}_{r}=\begin{bmatrix}\boldsymbol{I}_{r}&-\boldsymbol{I}_{r}\\-\boldsymbol{I}_{r}&\boldsymbol{I}_{r}\end{bmatrix}

Because p=(−dTdT)T\boldsymbol{p}=\begin{pmatrix}-\boldsymbol{d}^{T}&\boldsymbol{d}^{T}\end{pmatrix}^{T},

kG,el=NL(Qr−ppT)\boldsymbol{k}_{G,el}=\frac{N}{L}\left(\boldsymbol{Q}_{r}-\boldsymbol{p}\boldsymbol{p}^{T}\right)

The dyadic product is related directly to the elastic matrix:

ppT=LEAkel\boldsymbol{p}\boldsymbol{p}^{T}=\frac{L}{EA}\boldsymbol{k}_{el}

Therefore, the same relation can be written as:

kG,el=NLQr−NEAkel\boldsymbol{k}_{G,el}=\frac{N}{L}\boldsymbol{Q}_{r}-\frac{N}{EA}\boldsymbol{k}_{el}
The projection term must be retained

The quantity L/(EA)L/(EA) is dimensional and cannot be compared with unity. Moreover, (L/EA)kel=ppT(L/EA)\boldsymbol{k}_{el}=\boldsymbol{p}\boldsymbol{p}^{T}, whose magnitude does not decrease when EA/LEA/L increases. Thus the second term cannot be neglected on that basis.

The components of the stiffness matrices depend on the orientation of the element. When the element undergoes a rigid rotation, the direction vector and the transverse projector rotate with it. Consequently, the stiffness matrices change their components consistently, while the physical response of the element remains unchanged. This property is called objectivity.

04

Complete second-order strain expansion

Let Δu\Delta\boldsymbol{u} be the relative nodal displacement and define η=1LΔu\boldsymbol{\eta}=\frac{1}{L}\Delta\boldsymbol{u}. The exact engineering strain is:

ε=1+a−1,a=2dTη+ηTη\varepsilon=\sqrt{1+a}-1,\qquad a=2\boldsymbol{d}^{T}\boldsymbol{\eta}+\boldsymbol{\eta}^{T}\boldsymbol{\eta}

A complete expansion through second degree requires the first three terms of the binomial series:

1+a=1+12a−18a2+O(a3)\sqrt{1+a}=1+\frac{1}{2}a-\frac{1}{8}a^2+\mathcal{O}(a^3)

Retaining all terms up to second degree in the nodal displacements gives:

ε≈dTη+12ηT(Ir−ddT)η\varepsilon\approx\boldsymbol{d}^{T}\boldsymbol{\eta}+\frac{1}{2}\boldsymbol{\eta}^{T}\left(\boldsymbol{I}_{r}-\boldsymbol{d}\boldsymbol{d}^{T}\right)\boldsymbol{\eta}

In element nodal coordinates,

ε≈Buel+12uelTCuel\varepsilon\approx\boldsymbol{B}\boldsymbol{u}_{el}+\frac{1}{2}\boldsymbol{u}_{el}^{T}\boldsymbol{C}\boldsymbol{u}_{el}
C=1L2[P−P−PP]\boldsymbol{C}=\frac{1}{L^2}\begin{bmatrix}\boldsymbol{P}&-\boldsymbol{P}\\-\boldsymbol{P}&\boldsymbol{P}\end{bmatrix}
Why the third binomial term matters

If the term −a2/8-a^2/8 is omitted, the axial projection ddT\boldsymbol{d}\boldsymbol{d}^{T} disappears from P\boldsymbol{P}. The resulting quadratic strain and geometric stiffness are incomplete.

05

Energy interpretation

Consider an equilibrium reference state carrying the signed axial force

N=EAε‾N=EA\overline{\varepsilon}

For a small incremental displacement δuel\delta\boldsymbol{u}_{el}, the quadratic part of the incremental element energy is:

ΔUel(2)=12δuelT(kel+kG,el)δuel\Delta U_{el}^{(2)}=\frac{1}{2}\delta\boldsymbol{u}_{el}^{T}\left(\boldsymbol{k}_{el}+\boldsymbol{k}_{G,el}\right)\delta\boldsymbol{u}_{el}
kel=EALBTB,kG,el=NLC\boldsymbol{k}_{el}=EAL\boldsymbol{B}^{T}\boldsymbol{B},\qquad \boldsymbol{k}_{G,el}=NL\boldsymbol{C}

Tension adds transverse stiffness, whereas compression reduces it. Buckling occurs when the assembled tangent stiffness loses positive definiteness and becomes singular.

06

Conclusions

  • The elastic matrix is the axial dyadic product (EA/L)ppT(EA/L)\boldsymbol{p}\boldsymbol{p}^{T}.
  • The geometric matrix contains the transverse projector Ir−ddT\boldsymbol{I}_{r}-\boldsymbol{d}\boldsymbol{d}^{T}.
  • The 2D and 3D expressions differ only through the dimension and components of the direction vector.
  • The complete second-order strain requires the term −a2/8-a^2/8 in the binomial expansion.

No separate program package is required for this section. The matrices are implemented in the buckling programs presented in Sections 1.2 and 2.2.