Section 2.1 of Chapter 2: Three-Dimensional Truss Structures

Large Displacements of 3D Truss Structures

A total-potential-energy formulation for straight, prismatic three-dimensional truss elements, including exact displacement kinematics, a consistent tangent matrix, Newton–Raphson solution, and a comparison with ANSYS.

Linear elastic materialSmall strainLarge displacementMATLAB tested in R2017b and R2026a

01

Overview and assumptions

Consider a spatial structure assembled from straight, prismatic, two-node truss elements. Each element carries axial force only.

MaterialLinear elastic
StrainSmall
DisplacementPotentially large

Notation

u\boldsymbol{u}
Nodal-displacement vector in the global reference frame xyzxyz
F\boldsymbol{F}
Nodal-force vector in the global reference frame xyzxyz
U, V, ΠU,\ V,\ \Pi
Strain energy, potential of the applied forces, and total potential energy
σ, ε\sigma,\ \varepsilon
Axial stress and engineering strain
E, A, LE,\ A,\ L
Young's modulus, cross-sectional area, and initial element length
I3\boldsymbol{I}_3
Three-by-three identity matrix

02

Three-dimensional element kinematics

For the initial nodal coordinates(x1,y1,z1)(x_1,y_1,z_1) and(x2,y2,z2)(x_2,y_2,z_2), the initial element length is:

L=(x2x1)2+(y2y1)2+(z2z1)2L=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2}

Let(ui,vi,wi)(u_i,v_i,w_i) be the displacement components of nodeii. The deformed length is:

Lf=(x2x1+u2u1)2+(y2y1+v2v1)2+(z2z1+w2w1)2L_f=\sqrt{(x_2-x_1+u_2-u_1)^2+(y_2-y_1+v_2-v_1)^2+(z_2-z_1+w_2-w_1)^2}

The engineering strain is therefore:

ε=LfLL\varepsilon=\frac{L_f-L}{L}

The six element degrees of freedom are collected in the column vector:

uel=(u1v1w1u2v2w2)T\boldsymbol{u}_{el}=\begin{pmatrix}u_1&v_1&w_1&u_2&v_2&w_2\end{pmatrix}^{T}

03

Total potential energy and equilibrium

For a linear elastic material,σ=Eε\sigma=E\varepsilon. The total potential energy of the assembled structure is:

Π=U+V=12elVelεσdVuTF=12elEALε2uTF\Pi=U+V=\frac{1}{2}\sum_{el}\int_{\mathcal{V}_{el}}\varepsilon\sigma\,d\mathcal{V}-\boldsymbol{u}^{T}\boldsymbol{F}=\frac{1}{2}\sum_{el}EAL\varepsilon^2-\boldsymbol{u}^{T}\boldsymbol{F}

Here Vel=AL\mathcal{V}_{el}=AL, and the summation includes every truss element. The principle of minimum total potential energy requires

Ψ=Πu=UuF=0\boldsymbol{\Psi}=\frac{\partial\Pi}{\partial\boldsymbol{u}}=\frac{\partial U}{\partial\boldsymbol{u}}-\boldsymbol{F}=\boldsymbol{0}

Substitution of the element energy gives:

Ψ=elEALεεuF=0\boldsymbol{\Psi}=\sum_{el}EAL\varepsilon\frac{\partial\varepsilon}{\partial\boldsymbol{u}}-\boldsymbol{F}=\boldsymbol{0}

Because the strain depends nonlinearly on the nodal displacements, these equilibrium equations are nonlinear.

04

Element internal force and tangent stiffness

The strain energy of one truss element is:

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

Differentiation with respect to the six element degrees of freedom gives the internal nodal-force vector and the consistent tangent stiffness matrix:

fel=Ueluel=EALεεuel\boldsymbol{f}_{el}=\frac{\partial U_{el}}{\partial\boldsymbol{u}_{el}}=EAL\varepsilon\frac{\partial\varepsilon}{\partial\boldsymbol{u}_{el}}
kel=feluel=EAL[ε2εuel2+εuel(εuel)T]\boldsymbol{k}_{el}=\frac{\partial\boldsymbol{f}_{el}}{\partial\boldsymbol{u}_{el}}=EAL\left[\varepsilon\frac{\partial^2\varepsilon}{\partial\boldsymbol{u}_{el}^2}+\frac{\partial\varepsilon}{\partial\boldsymbol{u}_{el}}\left(\frac{\partial\varepsilon}{\partial\boldsymbol{u}_{el}}\right)^T\right]

The order of the two gradient factors is important: their outer product is a 6×66\times6 matrix.

05

Exact strain derivatives

Introduce the deformed coordinate differences

qx=x2x1+u2u1,qy=y2y1+v2v1,qz=z2z1+w2w1q_x=x_2-x_1+u_2-u_1,\qquad q_y=y_2-y_1+v_2-v_1,\qquad q_z=z_2-z_1+w_2-w_1

and the compact quantities

a={qxqyqzqxqyqz},G=[I3I3I3I3]\boldsymbol{a}=\begin{Bmatrix}-q_x\\-q_y\\-q_z\\q_x\\q_y\\q_z\end{Bmatrix},\qquad \boldsymbol{G}=\begin{bmatrix}\boldsymbol{I}_3&-\boldsymbol{I}_3\\-\boldsymbol{I}_3&\boldsymbol{I}_3\end{bmatrix}

The complete six-component strain gradient and the exact Hessian are:

εuel=1LLfa\frac{\partial\varepsilon}{\partial\boldsymbol{u}_{el}}=\frac{1}{LL_f}\boldsymbol{a}
2εuel2=1LLfG1LLf3aaT\frac{\partial^2\varepsilon}{\partial\boldsymbol{u}_{el}^2}=\frac{1}{LL_f}\boldsymbol{G}-\frac{1}{LL_f^3}\boldsymbol{a}\boldsymbol{a}^{T}

When axial deformation is very small, the approximation LfLL_f\approx L may be used. This gives:

εuel1L2a\frac{\partial\varepsilon}{\partial\boldsymbol{u}_{el}}\approx\frac{1}{L^2}\boldsymbol{a}
Exact derivative used in the program

Replacing LfL_f by LL does not by itself reduce the exact Hessian to(1/L2)G(1/L^2)\boldsymbol{G}. That shorter expression also omits the outer-product term. The supplied deriv.m uses the exact gradient and exact Hessian shown above.

06

Newton–Raphson solution and MATLAB implementation

At each iteration, the program assembles the structural tangent matrix and residual vector, applies the displacement constraints, and solves

Δu(i)=(Kt(i))1Ψ(i)\Delta\boldsymbol{u}^{(i)}=-\left(\boldsymbol{K}_t^{(i)}\right)^{-1}\boldsymbol{\Psi}^{(i)}
u(i+1)=u(i)+Δu(i)\boldsymbol{u}^{(i+1)}=\boldsymbol{u}^{(i)}+\Delta\boldsymbol{u}^{(i)}
deriv.mExact derivatives
% Exact engineering strain and derivatives
ep = Lf/L - 1;
q  = [dx+du; dy+dv; dz+dw];
a  = [-q; q];
G  = [eye(3) -eye(3); -eye(3) eye(3)];

ep1 = a/(L*Lf);
ep2 = G/(L*Lf) - (a*a')/(L*Lf^3);
stiff.mElement contribution
fel = EA*L*ep*ep1;
kel = EA*L*(ep*ep2 + ep1*ep1');

K(ip,ip) = K(ip,ip) + kel;
F(ip)    = F(ip) + fel;
MATLAB compatibility

Tested in MATLAB R2017b and MATLAB R2026a. Both versions produced the same displacement history and converged at all five load steps. 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 example in gen.m contains 18 repeated spatial modules, 76 nodes, 234 truss elements, and 228 structural degrees of freedom. The complete load is applied in five equal steps.

At load step 4, the three displacement components of node 76 are compared below. The values are expressed in millimetres.

Solutionuuvvww
MATLAB1760.92508.7−2639.6
ANSYS1775.32520.0−2665.3
Relative difference0.81%0.45%0.96%

The agreement is close, with every reported component differing by less than 1%. The supplied ANSYS macro uses BEAM188 with a very small bending inertia to approximate an axial member; therefore, it is not mathematically identical to the pure truss model. ANSYS did not converge for the fifth load level, whereas the MATLAB program converged in seven iterations.

Deformed configurations of the spatial truss during the five load steps, with node 76 identified
Figure 1. Successive deformed configurations of the 3D truss structure. Node 76 is identified at the top.
Scope

Within the stated assumptions, the formulation uses exact displacement kinematics and a consistent tangent. The number of load steps controls the numerical continuation path, not the converged equilibrium state. Local buckling of individual truss members is not included.

08

Program files

The MATLAB archive contains seven commented source files and a short README. The ANSYS macro is supplied separately so that the comparison model remains transparent.