FEA TipsFinite Element Theory & MATLAB Programs
HomeTwo-Dimensional Truss Structures1.3

Section 1.3 of Chapter 1: Two-Dimensional Truss Structures

Buckling Analysis — Part II

A second derivation of the elastic and geometric stiffness matrices, performed first in the element's local reference frame and then transformed to global coordinates.

Small displacementSmall strainLinear elastic materialEquations verified in MATLAB R2026a and R2017b

01

Local reference frame

Consider a two-node truss element under the hypotheses of small displacements, small strains, and linear elastic material behaviour.

The element is initially inclined at an angleθ\theta in the global reference framexyxy. Its local axis xˉ\bar{x} is aligned with the initial element.

Inclined truss element shown in the global frame with its local axes
Figure 1.Global and local reference frames.
Initial and displaced configurations of the truss element in its local frame
Figure 2.Element displacements measured in the local frame.

02

Second-order axial strain

As established in Section 1.1, the exact engineering strain in global coordinates is:

ε=(x2−x1+u2−u1L)2+(y2−y1+v2−v1L)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}-1

Equivalently,

ε=(cos⁡θ+u2−u1L)2+(sin⁡θ+v2−v1L)2−1\varepsilon=\sqrt{\left(\cos\theta+\frac{u_2-u_1}{L}\right)^2+\left(\sin\theta+\frac{v_2-v_1}{L}\right)^2}-1

In the local frame,xˉ2−xˉ1=L\bar{x}_2-\bar{x}_1=Landyˉ2−yˉ1=0\bar{y}_2-\bar{y}_1=0. Thus the initial element angle is zero in that frame, and the strain becomes:

ε=(1+uˉ2−uˉ1L)2+(vˉ2−vˉ1L)2−1\varepsilon=\sqrt{\left(1+\frac{\bar{u}_2-\bar{u}_1}{L}\right)^2+\left(\frac{\bar{v}_2-\bar{v}_1}{L}\right)^2}-1
ε=1+2uˉ2−uˉ1L+(uˉ2−uˉ1)2L2+(vˉ2−vˉ1)2L2−1\varepsilon=\sqrt{1+2\frac{\bar{u}_2-\bar{u}_1}{L}+\frac{(\bar{u}_2-\bar{u}_1)^2}{L^2}+\frac{(\bar{v}_2-\bar{v}_1)^2}{L^2}}-1

The barred displacements are measured in the local reference frame. Introduce

a=2uˉ2−uˉ1L+(uˉ2−uˉ1)2L2+(vˉ2−vˉ1)2L2a=2\frac{\bar{u}_2-\bar{u}_1}{L}+\frac{(\bar{u}_2-\bar{u}_1)^2}{L^2}+\frac{(\bar{v}_2-\bar{v}_1)^2}{L^2}

For ∣a∣≪1|a|\ll1, use the binomial expansion [1]

1+a=1+12a−18a2+116a3−5128a4+7256a5−⋯\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

Retaining the first three terms gives:

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

After neglecting terms of degree three and higher in the nodal displacements, the second-order strain is:

ε≈uˉ2−uˉ1L+(vˉ2−vˉ1)22L2\varepsilon\approx\frac{\bar{u}_2-\bar{u}_1}{L}+\frac{(\bar{v}_2-\bar{v}_1)^2}{2L^2}

03

Compact matrix form

Define the local element displacement vector:

uˉel=(uˉ1vˉ1uˉ2vˉ2)T\bar{\boldsymbol{u}}_{el}=\begin{pmatrix}\bar{u}_1&\bar{v}_1&\bar{u}_2&\bar{v}_2\end{pmatrix}^{T}

The second-order strain can be written as:

ε≈Bˉuˉel+12uˉelTCˉuˉel=ε0+12uˉelTCˉuˉel\varepsilon\approx\bar{\boldsymbol{B}}\bar{\boldsymbol{u}}_{el}+\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{C}}\bar{\boldsymbol{u}}_{el}=\varepsilon_0+\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{C}}\bar{\boldsymbol{u}}_{el}

where

Bˉ=1L(−1010)\bar{\boldsymbol{B}}=\frac{1}{L}\begin{pmatrix}-1&0&1&0\end{pmatrix}
Cˉ=1L2[0000010−100000−101]=1L2{010−1}{010−1}T\bar{\boldsymbol{C}}=\frac{1}{L^2}\begin{bmatrix}0&0&0&0\\0&1&0&-1\\0&0&0&0\\0&-1&0&1\end{bmatrix}=\frac{1}{L^2}\begin{Bmatrix}0\\1\\0\\-1\end{Bmatrix}\begin{Bmatrix}0\\1\\0\\-1\end{Bmatrix}^{T}

04

Element deformation energy

For one linear elastic truss element,

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

Substituting the second-order strain and neglecting the fourth-degree term gives:

Uel≈12EALε02+12EALε0uˉelTCˉuˉelU_{el}\approx\frac{1}{2}EAL\varepsilon_0^2+\frac{1}{2}EAL\varepsilon_0\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{C}}\bar{\boldsymbol{u}}_{el}
Uel=12uˉelTkˉeluˉel+12uˉelTkˉG,eluˉelU_{el}=\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{k}}_{el}\bar{\boldsymbol{u}}_{el}+\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{k}}_{G,el}\bar{\boldsymbol{u}}_{el}

or

Uel=12uˉelT(kˉel+kˉG,el)uˉel=12uˉelTkˉtotal,eluˉelU_{el}=\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\left(\bar{\boldsymbol{k}}_{el}+\bar{\boldsymbol{k}}_{G,el}\right)\bar{\boldsymbol{u}}_{el}=\frac{1}{2}\bar{\boldsymbol{u}}_{el}^{T}\bar{\boldsymbol{k}}_{\mathrm{total},el}\bar{\boldsymbol{u}}_{el}

Because the displacements are small, the axial force is evaluated from the first-order strain

N=EAε0N=EA\varepsilon_0

05

Stiffness matrices in the local frame

The elastic stiffness matrix is:

kˉel=EALBˉTBˉ=EAL{−1010}{−1010}T\bar{\boldsymbol{k}}_{el}=EAL\bar{\boldsymbol{B}}^{T}\bar{\boldsymbol{B}}=\frac{EA}{L}\begin{Bmatrix}-1\\0\\1\\0\end{Bmatrix}\begin{Bmatrix}-1\\0\\1\\0\end{Bmatrix}^{T}
kˉel=EAL[10−100000−10100000]\bar{\boldsymbol{k}}_{el}=\frac{EA}{L}\begin{bmatrix}1&0&-1&0\\0&0&0&0\\-1&0&1&0\\0&0&0&0\end{bmatrix}

The geometric stiffness matrix is:

kˉG,el=NLCˉ=NL{010−1}{010−1}T\bar{\boldsymbol{k}}_{G,el}=NL\bar{\boldsymbol{C}}=\frac{N}{L}\begin{Bmatrix}0\\1\\0\\-1\end{Bmatrix}\begin{Bmatrix}0\\1\\0\\-1\end{Bmatrix}^{T}
kˉG,el=NL[0000010−100000−101]\bar{\boldsymbol{k}}_{G,el}=\frac{N}{L}\begin{bmatrix}0&0&0&0\\0&1&0&-1\\0&0&0&0\\0&-1&0&1\end{bmatrix}

This energy-based separation into elastic and geometric contributions is the standard stress-stiffening construction [2] [3]

06

Transformation to the global frame

For an element inclined at angleθ\theta, the local and global displacement vectors satisfy

uˉel=Ruel\bar{\boldsymbol{u}}_{el}=\boldsymbol{R}\boldsymbol{u}_{el}
R=[cs00−sc0000cs00−sc]c=cos⁡θs=sin⁡θ\boldsymbol{R}=\begin{bmatrix}c&s&0&0\\-s&c&0&0\\0&0&c&s\\0&0&-s&c\end{bmatrix}\qquad c=\cos\theta\qquad s=\sin\theta

The total element stiffness matrix is therefore:

ktotal,el=RTkˉtotal,elR\boldsymbol{k}_{\mathrm{total},el}=\boldsymbol{R}^{T}\bar{\boldsymbol{k}}_{\mathrm{total},el}\boldsymbol{R}
kel+kG,el=RT(kˉel+kˉG,el)R\boldsymbol{k}_{el}+\boldsymbol{k}_{G,el}=\boldsymbol{R}^{T}\left(\bar{\boldsymbol{k}}_{el}+\bar{\boldsymbol{k}}_{G,el}\right)\boldsymbol{R}

In the global frame, 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}
kel=EAL{−cos⁡θ−sin⁡θcos⁡θsin⁡θ}{−cos⁡θ−sin⁡θcos⁡θsin⁡θ}T\boldsymbol{k}_{el}=\frac{EA}{L}\begin{Bmatrix}-\cos\theta\\-\sin\theta\\\cos\theta\\\sin\theta\end{Bmatrix}\begin{Bmatrix}-\cos\theta\\-\sin\theta\\\cos\theta\\\sin\theta\end{Bmatrix}^{T}

The geometric stiffness matrix is:

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}
kG,el=NL{−sin⁡θcos⁡θsin⁡θ−cos⁡θ}{−sin⁡θcos⁡θsin⁡θ−cos⁡θ}T\boldsymbol{k}_{G,el}=\frac{N}{L}\begin{Bmatrix}-\sin\theta\\\cos\theta\\\sin\theta\\-\cos\theta\end{Bmatrix}\begin{Bmatrix}-\sin\theta\\\cos\theta\\\sin\theta\\-\cos\theta\end{Bmatrix}^{T}

These are the same elastic and geometric stiffness matrices obtained directly in global coordinates in Section 1.2.

07

Verification

Independent matrix check

Both identitieskel=RTkˉelR\boldsymbol{k}_{el}=\boldsymbol{R}^{T}\bar{\boldsymbol{k}}_{el}\boldsymbol{R}andkG,el=RTkˉG,elR\boldsymbol{k}_{G,el}=\boldsymbol{R}^{T}\bar{\boldsymbol{k}}_{G,el}\boldsymbol{R}were checked numerically in MATLAB R2026a and MATLAB R2017b. The transformed matrices agree with their explicit global forms to machine precision.

No additional program package

This section provides an alternative derivation of the matrices used in Section 1.2. It introduces no new MATLAB program, so there is no separate download package.

08

References

  1. Wikipedia contributors, “Binomial theorem”
  2. Ansys, Inc., Ansys Mechanical APDL Theory Reference, Release 2026 R1, Section 3.4, “Stress Stiffening,” 2026.
  3. 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.