1.2 Buckling Analysis — Part I

Second-order axial strain

Consider a straight, prismatic two-node truss element. Although the displacements are small, the terms of second degree in the nodal displacements are retained.

Initial and displaced configurations of a two-node 2D truss element — Figure 1.Initial and displaced configurations in the global frame.

With $\Delta u=u_2-u_1$ , $\Delta v=v_2-v_1$ , and initial element angle $\theta$ , the exact engineering strain is:

\varepsilon=\sqrt{\left(\cos\theta+\frac{\Delta u}{L}\right)^2+\left(\sin\theta+\frac{\Delta v}{L}\right)^2}-1

This expression can be written as:

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

a=2\cos\theta\frac{\Delta u}{L}+2\sin\theta\frac{\Delta v}{L}+\frac{(\Delta u)^2}{L^2}+\frac{(\Delta v)^2}{L^2}

The binomial expansion [1] is:

\sqrt{1+a}=1+\frac{1}{2}a-\frac{1}{8}a^2+\frac{1}{16}a^3-\frac{5}{128}a^4+\frac{7}{256}a^5-\cdots

For $|a|\ll1$ , the first three terms are required because $a$ itself already contains terms of second degree:

\sqrt{1+a}\approx1+\frac{1}{2}a-\frac{1}{8}a^2

After discarding terms of degree three and higher in $u$ and $v$ ,

\varepsilon\approx\frac{\Delta u}{L}\cos\theta+\frac{\Delta v}{L}\sin\theta+\frac{(\Delta u)^2}{2L^2}\sin^2\theta+\frac{(\Delta v)^2}{2L^2}\cos^2\theta-\frac{\Delta u\,\Delta v}{2L^2}\sin2\theta

Compact matrix form

Introduce the element displacement vector:

\boldsymbol{u}_{el}=\begin{pmatrix}u_1&v_1&u_2&v_2\end{pmatrix}^{T}

The second-order strain can then be written as:

\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}

\varepsilon_0=\boldsymbol{B}\boldsymbol{u}_{el}

\boldsymbol{B}=\frac{1}{L}\begin{pmatrix}-\cos\theta&-\sin\theta&\cos\theta&\sin\theta\end{pmatrix}

Using the column vector:

\boldsymbol{q}=\begin{pmatrix}-\sin\theta&\cos\theta&\sin\theta&-\cos\theta\end{pmatrix}^{T}

the symmetric matrix in the quadratic term is:

\boldsymbol{C}=\frac{1}{L^2}\boldsymbol{q}\boldsymbol{q}^{T}

Element deformation energy

For a linear elastic truss element,

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

Substitution of the second-order strain would generate terms up to degree four. Neglecting terms of degree four gives:

U_{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 displacements are small, the axial force obtained from the preliminary linear analysis is:

N=EA\varepsilon_0

Therefore,

U_{el}=\frac{1}{2}\boldsymbol{u}_{el}^{T}\left(\boldsymbol{k}_{el}+\boldsymbol{k}_{G,el}\right)\boldsymbol{u}_{el}

Elastic and geometric stiffness matrices

Define

\boldsymbol{p}=\begin{pmatrix}-\cos\theta&-\sin\theta&\cos\theta&\sin\theta\end{pmatrix}^{T}

The two element matrices are:

\boldsymbol{k}_{el}=EAL\boldsymbol{B}^{T}\boldsymbol{B}=\frac{EA}{L}\boldsymbol{p}\boldsymbol{p}^{T}

\boldsymbol{k}_{G,el}=NL\boldsymbol{C}=\frac{N}{L}\boldsymbol{q}\boldsymbol{q}^{T}

This decomposition is the usual stress-stiffening formulation [2] [3]

The element matrices are assembled in the usual way to obtain the structural elastic stiffness matrix $\boldsymbol{K}$ and geometric stiffness matrix $\boldsymbol{K}_{G}$ .

The buckling eigenproblem

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

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

With the program convention $N>0$ in tension and $N<0$ in compression, the corresponding generalized eigenproblem is:

\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, the smallest positive eigenvalue is the quantity of interest.

Sign convention

The minus sign is required by the signed axial-force convention used in the MATLAB files and produces positive critical load factors for compression.

MATLAB implementation

The calculation has two stages:

truss2d.m solves the small-displacement linear problem, and stress2d.m computes the signed axial forces.
buckling.m assembles $\boldsymbol{K}_{G}$ and solves the generalized eigenproblem.

buckling.mSigned eigenproblem

% Solve for the free degrees of freedom
fixedDofs=2*(cond(:,1)-1)+cond(:,2);
freeDofs=setdiff((1:neq)',fixedDofs);
Kff=K(freeDofs,freeDofs);
KGff=KG(freeDofs,freeDofs);

[Vfree,lambda]=eig(Kff,-KGff,'vector');
valid=isfinite(lambda) & real(lambda)>0;
lambda=sort(real(lambda(valid)));

MATLAB verification

Tested in MATLAB R2017b and MATLAB R2026a. In both versions the package solves only for the unconstrained degrees of freedom, removes nonphysical infinite eigenvalues, and returns the same positive critical load factors. The code uses standard MATLAB syntax and is expected to be compatible with intermediate and newer MATLAB versions.

Example 1 — eigenvalue buckling

The geometry, material data, constraints, and reference loads are specified in gen.m. The first two buckling modes are shown below.

First and second buckling modes of the example truss tower — Figure 2.First and second buckling modes; the undeformed structure is blue and the normalized mode shape is magenta.

Solution	First buckling load	Second buckling load
MATLAB R2026a and R2017b	940,978.7131	6,181,251.434
ANSYS	9.4098 × 10⁵	6.1813 × 10⁶

Example 2 — nonlinear response near buckling

The same structure is analysed with the nonlinear program from Section 1.1. The vertical force increases from zero to $F=940\,980\ \mathrm{N}$ , while a small horizontal force increases from zero to $H$ .

With the residual convention used in the Section 1.1 program, downward and rightward loads are entered with negative signs:

gen.mApplied loads

forze=[11  1  -H
       11  2  -470490
       12  2  -470490];

Nonlinear deformation of the truss tower under vertical and horizontal loads — Figure 3.Representative nonlinear deformation.

Vertical load versus maximum horizontal displacement for four horizontal forces — Figure 4.Vertical load versus maximum horizontal displacement for four values of $H$ .

For $H=2\,000,\ 5\,000,\ 10\,000$ , and $20\,000\ \mathrm{N}$ , the load-displacement curves tend asymptotically to the first buckling load, as expected.

Load sign convention

The negative signs are included so that the data agree with the coordinate directions, the arrows in the figure, and the residual convention implemented in stiff.m.

Why the third binomial term matters

If the expansion is stopped after only two terms,

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

the resulting strain is:

\varepsilon\approx\frac{\Delta u}{L}\cos\theta+\frac{\Delta v}{L}\sin\theta+\frac{(\Delta u)^2}{2L^2}+\frac{(\Delta v)^2}{2L^2}

This is not the complete second-order expansion of the exact engineering strain: the orientation-dependent quadratic terms, including the mixed term $\Delta u\,\Delta v$ , are lost.

The corresponding alternative geometric matrix would be

\boldsymbol{k}_{G,el}^{*}=\frac{N}{L}\begin{bmatrix}1&0&-1&0\\0&1&0&-1\\-1&0&1&0\\0&-1&0&1\end{bmatrix}

This matrix is invariant under a rotation of the global axes, but it is not identical to the consistent geometric stiffness obtained from the complete second-order expansion above. It should therefore be regarded as an alternative formulation, not as an algebraically equivalent form of $\boldsymbol{k}_{G,el}$ .

Program files

MMATLAB buckling package6 commented MATLAB files + README · tested in R2017b and R2026aDownload ZIP AANSYS comparison macro1 Mechanical APDL input fileDownload ZIP

References

Wikipedia contributors, “Binomial theorem”
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.
I. Němec, M. Trcala, I. Ševčík, and H. Štekbauer, “New Formula for Geometric Stiffness Matrix Calculation”, Journal of Applied Mathematics and Physics, 4 (2016), 733–748.