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

Section 2.2 of Chapter 2: Three-Dimensional Truss Structures

Buckling Analysis of 3D Truss Structures

A second-order strain formulation for straight, prismatic three-dimensional truss elements, leading to the elastic and geometric stiffness matrices, a generalized eigenproblem, and four verified buckling modes.

Linear elastic materialSmall displacementSmall strainMATLAB tested in R2017b and R2026a

01

Second-order axial strain

The analysis assumes small displacements, small strains, and a linear elastic material. Terms of second degree in the nodal displacements are nevertheless retained.

For the two-node spatial truss introduced in Section 2.1, the exact engineering strain is:

ε=(x2−x1+u2−u1L)2+(y2−y1+v2−v1L)2+(z2−z1+w2−w1L)2−1\varepsilon=\sqrt{\left(\frac{x_2-x_1+u_2-u_1}{L}\right)^2+\left(\frac{y_2-y_1+v_2-v_1}{L}\right)^2+\left(\frac{z_2-z_1+w_2-w_1}{L}\right)^2}-1

Let ll, mm, and nn be the direction cosines of the initial element axis. The same expression can be written as:

ε=(l+u2−u1L)2+(m+v2−v1L)2+(n+w2−w1L)2−1\varepsilon=\sqrt{\left(l+\frac{u_2-u_1}{L}\right)^2+\left(m+\frac{v_2-v_1}{L}\right)^2+\left(n+\frac{w_2-w_1}{L}\right)^2}-1

Using l2+m2+n2=1l^2+m^2+n^2=1, define

a=2lu2−u1L+2mv2−v1L+2nw2−w1L+(u2−u1)2L2+(v2−v1)2L2+(w2−w1)2L2a=2l\frac{u_2-u_1}{L}+2m\frac{v_2-v_1}{L}+2n\frac{w_2-w_1}{L}+\frac{(u_2-u_1)^2}{L^2}+\frac{(v_2-v_1)^2}{L^2}+\frac{(w_2-w_1)^2}{L^2}

Then

ε=1+a−1\varepsilon=\sqrt{1+a}-1

Although the displacements are small, the quadratic terms in uu, vv, and ww must be kept. The binomial expansion [1] is limited to its first three terms:

1+a≈1+12a−18a2\sqrt{1+a}\approx1+\frac{1}{2}a-\frac{1}{8}a^2

02

Compact matrix form

Neglecting terms of degree three and higher in all three displacement components gives:

ε≈Buel+12uelTCuel=ε0+12uelTCuel\varepsilon\approx\boldsymbol{B}\boldsymbol{u}_{el}+\frac{1}{2}\boldsymbol{u}_{el}^{T}\boldsymbol{C}\boldsymbol{u}_{el}=\varepsilon_0+\frac{1}{2}\boldsymbol{u}_{el}^{T}\boldsymbol{C}\boldsymbol{u}_{el}
uel=(u1v1w1u2v2w2)T\boldsymbol{u}_{el}=\begin{pmatrix}u_1&v_1&w_1&u_2&v_2&w_2\end{pmatrix}^{T}

where

ε0=Buel\varepsilon_0=\boldsymbol{B}\boldsymbol{u}_{el}
B=1L(−l−m−nlmn)\boldsymbol{B}=\frac{1}{L}\begin{pmatrix}-l&-m&-n&l&m&n\end{pmatrix}
C=1L2[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{C}=\frac{1}{L^2}\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}

If d=(lmn)T\boldsymbol{d}=\begin{pmatrix}l&m&n\end{pmatrix}^{T} and P=I3−ddT\boldsymbol{P}=\boldsymbol{I}_3-\boldsymbol{d}\boldsymbol{d}^{T}, the same matrix has the compact form:

C=1L2[P−P−PP]\boldsymbol{C}=\frac{1}{L^2}\begin{bmatrix}\boldsymbol{P}&-\boldsymbol{P}\\-\boldsymbol{P}&\boldsymbol{P}\end{bmatrix}

03

Element deformation energy

The deformation energy of one linear elastic truss element is:

Uel=12EALε2U_{el}=\frac{1}{2}EAL\varepsilon^2

Substitution of the second-order strain generates terms up to degree four. After neglecting the fourth-degree terms,

Uel≈12EALε02+12EALε0uelTCuelU_{el}\approx\frac{1}{2}EAL\varepsilon_0^2+\frac{1}{2}EAL\varepsilon_0\boldsymbol{u}_{el}^{T}\boldsymbol{C}\boldsymbol{u}_{el}

Since the displacements are small, the signed axial force from the preliminary linear analysis is:

N=EAε0N=EA\varepsilon_0

Therefore,

Uel=12uelT(kel+kG,el)uelU_{el}=\frac{1}{2}\boldsymbol{u}_{el}^{T}\left(\boldsymbol{k}_{el}+\boldsymbol{k}_{G,el}\right)\boldsymbol{u}_{el}

04

Elastic and geometric stiffness matrices

Define the column vector:

p=(−l−m−nlmn)T\boldsymbol{p}=\begin{pmatrix}-l&-m&-n&l&m&n\end{pmatrix}^{T}

The elastic stiffness matrix is:

kel=EALBTB=EAL[l2lmln−l2−lm−lnlmm2mn−lm−m2−mnlnmnn2−ln−mn−n2−l2−lm−lnl2lmln−lm−m2−mnlmm2mn−ln−mn−n2lnmnn2]\boldsymbol{k}_{el}=EAL\boldsymbol{B}^{T}\boldsymbol{B}=\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=NLC=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}=NL\boldsymbol{C}=\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}

The element matrices are assembled to obtain the structural elastic stiffness matrix K\boldsymbol{K} and geometric stiffness matrix KG\boldsymbol{K}_{G}. This is the stress-stiffening contribution used in linearized eigenvalue buckling [2]

05

The buckling eigenproblem

Let the signed axial forces and KG\boldsymbol{K}_{G} correspond to a reference load. At a load factor λ\lambda, the total tangent stiffness becomes singular:

det⁡(K+λKG)=0\det\left(\boldsymbol{K}+\lambda\boldsymbol{K}_{G}\right)=0

The MATLAB programs use N>0N>0 in tension and N<0N<0 in compression. The corresponding generalized eigenproblem is:

Kϕ=−λKGϕ\boldsymbol{K}\boldsymbol{\phi}=-\lambda\boldsymbol{K}_{G}\boldsymbol{\phi}

Here λ\lambda is the buckling load factor and ϕ\boldsymbol{\phi} is the buckling mode. Usually, only the smallest positive eigenvalue is required.

Sign convention

The minus sign follows directly from the determinant equation and from the signed axial-force convention used in the MATLAB files. It produces positive critical load factors for compression.

06

MATLAB implementation

The calculation is completed in two stages:

  1. truss3d.m solves the small-displacement linear problem, and stress3d.m computes the signed axial forces in the trusses.
  2. buckling.m assembles KG\boldsymbol{K}_{G} and computes the positive buckling load factors and their eigenvectors.
buckling.mGeometric stiffness and signed eigenproblem
% Projection transverse to the element axis
direction = [dx;dy;dz]/L;
P = eye(3) - direction*direction';
kelG = (N(ie)/L)*[P -P;-P P];

% Compression is negative: K*phi = -lambda*KG*phi
[Vfree,D] = eig(Kff,-KGff);
lambdaAll = diag(D);
MATLAB verification

Tested in MATLAB R2017b and MATLAB R2026a. Both versions returned the same four positive load factors. Restrained degrees of freedom are removed before solving, so nonphysical infinite eigenvalues are not included. The code uses standard MATLAB syntax and is expected to be compatible with intermediate and newer MATLAB versions.

07

Numerical example and ANSYS comparison

The data are defined in gen.m. The structure contains 18 repeated spatial modules, 76 nodes, 234 truss elements, and 228 degrees of freedom. A total reference load of 1 N acts downward at the four top nodes, so the numerical load factors also give the critical loads in newtons.

The four values computed with the MATLAB package are compared below with the ANSYS results.

ModeVerified MATLABANSYSDifference from ANSYS
1142,154.8729 N1.41923 × 105 N0.163%
2158,603.4197 N1.5828 × 105 N0.204%
31,150,305.944 N1.1367 × 106 N1.197%
41,268,603.015 N1.2506 × 106 N1.440%

The supplied ANSYS macro uses BEAM188 with very small bending inertias to approximate axial members. It is therefore close to, but not mathematically identical with, the pure truss model used by MATLAB.

First buckling mode of the spatial truss tower
Mode 1.Fcr=1.4215×105 NF_{cr}=1.4215\times10^5\ \mathrm{N}
Second buckling mode of the spatial truss tower
Mode 2.Fcr=1.5860×105 NF_{cr}=1.5860\times10^5\ \mathrm{N}
Third buckling mode of the spatial truss tower
Mode 3.Fcr=1.1503×106 NF_{cr}=1.1503\times10^6\ \mathrm{N}
Fourth buckling mode of the spatial truss tower
Mode 4.Fcr=1.2686×106 NF_{cr}=1.2686\times10^6\ \mathrm{N}

Figure 1. The first four buckling modes. The undeformed structure is blue and the scaled eigenvector is magenta.

08

Why the third binomial term matters

If the expansion is limited to only its first two terms,

1+a≈1+12a\sqrt{1+a}\approx1+\frac{1}{2}a

the resulting strain is:

ε≈lu2−u1L+mv2−v1L+nw2−w1L+(u2−u1)22L2+(v2−v1)22L2+(w2−w1)22L2\varepsilon\approx l\frac{u_2-u_1}{L}+m\frac{v_2-v_1}{L}+n\frac{w_2-w_1}{L}+\frac{(u_2-u_1)^2}{2L^2}+\frac{(v_2-v_1)^2}{2L^2}+\frac{(w_2-w_1)^2}{2L^2}

This expression is not the complete second-order expansion: orientation-dependent quadratic terms are lost. Its associated geometric stiffness matrix would be

kG,el∗=NL[100−1000100−1000100−1−1001000−1001000−1001]\boldsymbol{k}_{G,el}^{*}=\frac{N}{L}\begin{bmatrix}1&0&0&-1&0&0\\0&1&0&0&-1&0\\0&0&1&0&0&-1\\-1&0&0&1&0&0\\0&-1&0&0&1&0\\0&0&-1&0&0&1\end{bmatrix}

This alternative matrix is invariant under rotation of the coordinate axes, but it is not algebraically equivalent to the consistent matrix obtained from the complete second-order expansion.

09

Program files

The MATLAB package contains six commented source files and a README. The ANSYS Mechanical APDL macro is supplied separately.

MMATLAB buckling package6 commented MATLAB files + README · tested in R2017b and R2026aDownload ZIPAANSYS comparison macroMechanical APDL input fileDownload

10

References

  1. Wikipedia contributors, “Binomial theorem”
  2. Ansys, Inc., Ansys Mechanical APDL Theory Reference, Release 2026 R1, sections “Stress Stiffening,” “Eigenvalue and Eigenvector Extraction,” and “Buckling Analysis,” 2026. The supplied comparison macro uses syntax compatible with ANSYS Release 12.0.