Stiffness and mass matrices

From a continuous element to nodal matrices

A real member has a continuous displacement field: its value can be found at any point along the length. A finite element replaces this infinite set of unknowns with a few nodal degrees of freedom q_e. Displacements between nodes are interpolated using predefined shape functions:

u(x) = N(x) q_e

The choice of N(x) determines the displacement distribution inside the element. The matrices follow by substituting this approximation into strain energy and kinetic energy expressions:

K_e = ∫_V B^TDB dV,    M_e = ∫_V ρN^TN dV,    ε = Bq_e

For a two-node type 112 member, axial displacement is interpolated linearly. With ξ = x/L:

u(x) = (1 − ξ)u′₁ + ξu′₂

The axial strain ε = du/dx is therefore constant within the element. Here the operator B contains first derivatives of the shape functions. For a type 113 frame, axial displacement is also linear, but transverse deflection uses cubic Hermite functions:

v(x) = H₁(x)v′₁ + H₂(x)θ₁ + H₃(x)v′₂ + H₄(x)θ₂

In the bending part, the operator B contains second derivatives of the functions H. Thus the curvature κ = d²v/dx² and, for constant EI, the bending moment vary linearly along the length.

Two points therefore do not reproduce an arbitrary deformed shape: they define the selected approximation. Subdividing a structure into more elements adds nodes and allows a more accurate description of displacement and strain distributions.

Plane truss member, type 112

In local axes, x′, y′ the element resists only changes in the distance between its nodes. The implementation retains both translational components at each node, so the local vector has four components.

k = EA/L

K_local = k · [ 1  0 -1  0
                0  0  0  0
               -1  0  1  0
                0  0  0  0 ]

Consistent mass distributes translational inertia between the two nodes:

M_local = ρAL/6 · [ 2  0  1  0
                    0  2  0  1
                    1  0  2  0
                    0  1  0  2 ]

The coefficients 2 and 1 come from the integral ∫N^TN dx. Off-diagonal terms are not artificial coupling: the motion of every interior point depends on the displacements of both nodes. This mass matrix is therefore consistent with the chosen shape functions.

For a free member, the rank of K is 1: only one linear combination of displacements produces axial strain. The other three independent modes—two translations and one infinitesimal rigid-body rotation—cause no axial elongation in the linear formulation. Their elastic strain energy is therefore zero.

Plane Euler–Bernoulli frame, type 113

A frame element combines axial stiffness and bending. The Euler–Bernoulli assumption neglects transverse shear deformation, so rotation is related to the derivative of deflection.

Integrating derivatives of these shape functions gives the characteristic bending coefficients 12, 6, 4 and 2 in the stiffness matrix:

a = EA/L,   b = 12EI/L³,   c = 6EI/L²,   d = 4EI/L,   e = 2EI/L

K_local = [ a   0   0  -a   0   0
            0   b   c   0  -b   c
            0   c   d   0  -c   e
           -a   0   0   a   0   0
            0  -b  -c   0   b  -c
            0   c   e   0  -c   d ]

In createStructuralElement.m the standard consistent mass matrix is used:

M_local = ρAL/420 ·
[ 140    0       0     70    0        0
    0   156     22L     0   54      -13L
    0   22L     4L²     0   13L     -3L²
   70    0       0    140    0        0
    0    54     13L     0  156      -22L
    0  -13L    -3L²     0  -22L      4L² ]

The same cubic deflection interpolation is used in kinetic energy, giving coefficients 156, 54, 22L, 13L and the other terms of the consistent mass matrix.

A free frame element has three rigid-body modes: two translations and one rotation. Its 6×6 stiffness matrix therefore has rank 3.

Units and coefficient meanings

Matrices with rotational DOFs do not have one common physical unit for all entries. For example, force–displacement coupling has units N/m, while moment–rotation coupling has units N·m. Check each generalized-force/generalized-displacement pair individually.

Five checks before assembly

1. Size: 4×4 for type 112, 6×6 for type 113.
2. Finiteness: no coefficient is NaN or Inf.
3. Symmetry: K = K^T and M = M^T.
4. Rank: 1 for a free truss member, 3 for a free frame.
5. Mass: under unit translation, the quadratic form r^TMr equals ρAL.

A free element stiffness matrix must be singular. This is not an error: the element can move without deforming. After assembly and applying supports, the global matrix should become nonsingular; otherwise the model still has an unrestrained motion or an internal mechanism.

Laboratory 01. Stiffness and mass of an inclined element

The aim is to relate matrix coefficients to geometry and explain why the axial stiffness of an inclined member couples X and Y displacements.

Input data and notation

Nodes: (0, 0) and (3, 4) m. Truss element: A = 0.01 m², E = 200 GPa, ρ = 7850 kg/m³. Frame: A = 0.02 m², E = 210 GPa, ρ = 7850 kg/m³, I = 8·10⁻⁶ m⁴. No supports or loads are specified: the matrices of a single free element are studied.

Before running: hand calculations

1. Find L, c = Δx/L, s = Δy/L and the axial member stiffness EA/L.
2. Calculate three coefficients of its global matrix: K₁₁ = (EA/L)c², K₁₂ = (EA/L)cs, K₂₂ = (EA/L)s².
3. Find the physical mass ρAL of each model.

Run in Octave

Run from the project root. The supplied function prints the results; you do not need to edit its code for this exercise.

addpath(pwd);
setup;
addpath(fullfile(pwd, "reference", "examples"));
r = example_01_element_matrices();

Working with the results

1. Compare the length and direction cosines with the beginning of the output. In Truss local stiffness find the axial stiffness, and in Truss global stiffness find the three calculated coefficients. Explain the nonzero K₁₂.
2. Using the labels Truss DOFs and Frame DOFs label the matrix rows. In a frame element, the third and sixth components are rotations.
3. Examine rows 2 and 4 of the local truss stiffness and the diagonal of its mass matrix: why does zero transverse stiffness not imply zero mass?
4. Run the commands below to find the mass when both nodes undergo equal X translation. Compare with ρAL. Use the formulas to predict how mass and axial stiffness change if the length doubles while the direction stays the same.

tx = [1; 0; 1; 0];
fx = [1; 0; 0; 1; 0; 0];
disp(tx' * r.truss.mass * tx);
disp(fx' * r.frame.mass * fx);

What to include in the report

Geometry calculations, a table of the three coefficients and two masses, and an explanation of X–Y coupling and the difference between transverse stiffness and mass.