Modal Analysis of a 3D Part

Eigenfrequencies and mode shapes with 3D finite elements

Run AllReset

Every elastic structure has a set of natural frequencies at which it likes to vibrate, each with a characteristic mode shape. Finding them is modal analysis, and it is the foundation of structural dynamics — from tuning a guitar body to keeping a bridge away from resonance.

Mathematically, modal analysis is the generalized eigenproblem

$K\,\phi = \lambda\,M\,\phi,$

where $K$ is the stiffness matrix, $M$ the mass matrix, the eigenvalues are the squared angular frequencies $\lambda = \omega^2$ , and the eigenvectors $\phi$ are the mode shapes. This whole pipeline — build a 3D mesh, assemble $K$ and $M$ , solve the eigenproblem, and visualize the modes — runs entirely in your browser.

Problem: Find the lowest vibration modes of a steel cantilever bracket — a $1.0 \times 0.1 \times 0.1\,\text{m}$ block, clamped on one end and free everywhere else.

This extends the linear elasticity tutorial from 1D to a full 3D continuum: displacements are vectors, the strain is a tensor, and each element contributes a $24 \times 24$ block. The assembly pattern — element matrices, sparse triplets, boundary conditions — is exactly the same.

Step 1: Material and the Elasticity Matrix

For an isotropic material, stress and strain are linked by Hooke's law $\sigma = D\,\varepsilon$ . In 3D, using Voigt notation with the six strain components $(\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{xy}, \gamma_{yz}, \gamma_{xz})$ , the constitutive matrix $D$ is built from the Lamé parameters $\lambda$ and $\mu$ .

Code [1]Run

# Steel
E = 210e9        # Young's modulus [Pa]
nu = 0.3         # Poisson's ratio
rho = 7800.0     # density [kg/m^3]

# Lamé parameters
lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E / (2.0 * (1.0 + nu))

# 6x6 isotropic elasticity matrix (Voigt notation)
D = zeros(6, 6)
D[1,1] = lam + 2*mu
D[1,2] = lam
D[1,3] = lam
D[2,1] = lam
D[2,2] = lam + 2*mu
D[2,3] = lam
D[3,1] = lam
D[3,2] = lam
D[3,3] = lam + 2*mu
D[4,4] = mu
D[5,5] = mu
D[6,6] = mu

println("Shear modulus mu = ${mu} Pa")

Step 2: The Mesh

We discretize the bracket with a structured grid of trilinear hexahedral elements (VTK cell type 12 — the eight-node brick). With nx, ny, nz elements per direction there are $(n_x{+}1)(n_y{+}1)(n_z{+}1)$ nodes, and every node carries three displacement degrees of freedom $(u_x, u_y, u_z)$ .

Nodes are numbered so that node $(i_x, i_y, i_z)$ gets global index $g = 1 + i_x + (n_x{+}1)\,i_y + (n_x{+}1)(n_y{+}1)\,i_z$ . We store each node's reference coordinates for later warping.

Code [2]Run

# Bracket dimensions [m]
Lx = 1.0
Ly = 0.1
Lz = 0.1

# Elements per direction
nx = 12
ny = 2
nz = 2

hx = Lx / nx
hy = Ly / ny
hz = Lz / nz

nnode = (nx+1)*(ny+1)*(nz+1)
ndof = 3*nnode
nel = nx*ny*nz

# Reference node coordinates
X0 = zeros(nnode)
Y0 = zeros(nnode)
Z0 = zeros(nnode)
for iz in 0:nz {
  for iy in 0:ny {
    for ix in 0:nx {
      gg = 1 + ix + (nx+1)*iy + (nx+1)*(ny+1)*iz
      X0[gg] = ix*hx
      Y0[gg] = iy*hy
      Z0[gg] = iz*hz
    }
  }
}

println("=== Mesh ===")
println("Elements: ${nel}")
println("Nodes: ${nnode}")
println("Degrees of freedom: ${ndof}")

Step 3: Element Stiffness and Mass

Each hexahedron is integrated with $2\times2\times2$ Gauss quadrature. At every Gauss point we evaluate the eight trilinear shape functions $N_a$ and their spatial derivatives, assemble the $6 \times 24$ strain–displacement matrix $B$ , and accumulate

$k_e = \int_{\Omega_e} B^\top D\, B \; dV, \qquad m_e = \int_{\Omega_e} \rho\, N^\top N \; dV.$

Because every element is an axis-aligned box of the same size, the Jacobian is diagonal and constant ( $\partial\xi/\partial x = 2/h_x$ , etc.), so the derivative mapping is a simple scaling and $\det J = \tfrac{h_x}{2}\tfrac{h_y}{2}\tfrac{h_z}{2}$ . All elements are identical, so we build one reference $k_e$ once and reuse it.

For the mass we use a lumped (diagonal) mass matrix: each node receives $\rho\,V_e/8$ of the element mass in each direction. Lumping preserves the total mass, keeps $M$ diagonal and positive-definite, and is standard practice for modal analysis.

Code [3]Run

# Reference hexahedron nodes in natural coordinates (VTK ordering)
xi  = [-1.0, 1.0, 1.0,-1.0,-1.0, 1.0, 1.0,-1.0]
eta = [-1.0,-1.0, 1.0, 1.0,-1.0,-1.0, 1.0, 1.0]
zta = [-1.0,-1.0,-1.0,-1.0, 1.0, 1.0, 1.0, 1.0]

# 2-point Gauss rule (weights are 1)
g = 1.0 / sqrt(3.0)
gp = [-g, g]
detJ = (hx/2) * (hy/2) * (hz/2)

# Element stiffness by 2x2x2 quadrature
Ke = zeros(24, 24)
for iq in 1:2 {
  for jq in 1:2 {
    for kq in 1:2 {
      s = gp[iq]
      t = gp[jq]
      u = gp[kq]

      # Shape-function derivatives w.r.t. x, y, z at this Gauss point
      dNdx = zeros(8)
      dNdy = zeros(8)
      dNdz = zeros(8)
      for a in 1:8 {
        dNdx[a] = 0.125*xi[a]*(1+eta[a]*t)*(1+zta[a]*u) * (2.0/hx)
        dNdy[a] = 0.125*(1+xi[a]*s)*eta[a]*(1+zta[a]*u) * (2.0/hy)
        dNdz[a] = 0.125*(1+xi[a]*s)*(1+eta[a]*t)*zta[a] * (2.0/hz)
      }

      # Strain-displacement matrix B (6 x 24)
      B = zeros(6, 24)
      for a in 1:8 {
        c1 = 3*(a-1)+1   # ux column
        c2 = c1+1        # uy column
        c3 = c1+2        # uz column
        B[1,c1] = dNdx[a]
        B[2,c2] = dNdy[a]
        B[3,c3] = dNdz[a]
        B[4,c1] = dNdy[a]
        B[4,c2] = dNdx[a]
        B[5,c2] = dNdz[a]
        B[5,c3] = dNdy[a]
        B[6,c1] = dNdz[a]
        B[6,c3] = dNdx[a]
      }

      Ke = Ke + (transpose(B) * D * B) * detJ
    }
  }
}

# Lumped element mass: rho * V_e / 8 per node, per direction
mlump = rho * (hx*hy*hz) / 8.0
Me = zeros(24, 24)
for a in 1:8 {
  for d in 1:3 {
    Me[3*(a-1)+d, 3*(a-1)+d] = mlump
  }
}

println("Element stiffness Ke[1,1] = ${Ke[1,1]}")
println("Nodal lumped mass = ${mlump} kg")

Step 4: Global Assembly

We loop over elements, map each element's 8 local nodes to their global indices (in VTK ordering, so the connectivity can be reused for visualization), and scatter the $24\times24$ blocks of $k_e$ and $m_e$ into global sparse triplets. As always, sparse sums duplicate $(i,j)$ entries — exactly what assembly needs, since shared nodes accumulate contributions from every adjoining element.

Code [4]Run

# Local-node offsets within an element (VTK hexahedron ordering)
ox = [0,1,1,0,0,1,1,0]
oy = [0,0,1,1,0,0,1,1]
oz = [0,0,0,0,1,1,1,1]

# 24 x 24 = 576 triplets per element
ntrip = nel * 576
Ki = zeros(ntrip)
Kj = zeros(ntrip)
Kv = zeros(ntrip)
Mv = zeros(ntrip)

# Element -> global node connectivity (kept for the viewer)
conn = zeros(nel, 8)

p = 0
el = 0
for ez in 0:nz-1 {
  for ey in 0:ny-1 {
    for ex in 0:nx-1 {
      el = el+1
      gn = zeros(8)
      for a in 1:8 {
        ix = ex+ox[a]
        iy = ey+oy[a]
        iz = ez+oz[a]
        gn[a] = 1 + ix + (nx+1)*iy + (nx+1)*(ny+1)*iz
        conn[el, a] = gn[a]
      }
      # Scatter Ke, Me into global triplets
      for a in 1:8 {
        for da in 1:3 {
          gi = 3*(gn[a]-1)+da
          la = 3*(a-1)+da
          for b in 1:8 {
            for db in 1:3 {
              gj = 3*(gn[b]-1)+db
              lb = 3*(b-1)+db
              p = p+1
              Ki[p] = gi
              Kj[p] = gj
              Kv[p] = Ke[la,lb]
              Mv[p] = Me[la,lb]
            }
          }
        }
      }
    }
  }
}

K = sparse(Ki, Kj, Kv, ndof, ndof)
M = sparse(Ki, Kj, Mv, ndof, ndof)

println("=== Global system ===")
println("K: ${ndof} x ${ndof}, ${nnz(K)} nonzeros")
println("M: diagonal, total mass = ${rho*Lx*Ly*Lz} kg")

Step 5: Boundary Conditions

The bracket is clamped on the $x=0$ face: every degree of freedom of those nodes is held at zero. As in the elasticity tutorial, we eliminate the fixed DOFs entirely, renumbering the remaining free DOFs into a smaller system. Working from the triplets keeps this simple: drop any entry that touches a fixed DOF and reindex the rest.

Code [5]Run

# Mark fixed DOFs: all nodes on the x = 0 face
isfixed = zeros(ndof)
for iz in 0:nz {
  for iy in 0:ny {
    gnode = 1 + 0 + (nx+1)*iy + (nx+1)*(ny+1)*iz
    isfixed[3*(gnode-1)+1] = 1
    isfixed[3*(gnode-1)+2] = 1
    isfixed[3*(gnode-1)+3] = 1
  }
}

# New (reduced) index for each free DOF
newidx = zeros(ndof)
nfree = 0
for i in 1:ndof {
  if isfixed[i] < 0.5 {
    nfree = nfree+1
    newidx[i] = nfree
  }
}

# Keep only triplets between two free DOFs, reindexed
rKi = zeros(ntrip)
rKj = zeros(ntrip)
rKv = zeros(ntrip)
rMv = zeros(ntrip)
q = 0
for t in 1:p {
  gi = Ki[t]
  gj = Kj[t]
  if isfixed[gi] < 0.5 {
    if isfixed[gj] < 0.5 {
      q = q+1
      rKi[q] = newidx[gi]
      rKj[q] = newidx[gj]
      rKv[q] = Kv[t]
      rMv[q] = Mv[t]
    }
  }
}

# Trim to the used length and build the reduced matrices
a_i = zeros(q)
a_j = zeros(q)
a_kv = zeros(q)
a_mv = zeros(q)
for t in 1:q {
  a_i[t] = rKi[t]
  a_j[t] = rKj[t]
  a_kv[t] = rKv[t]
  a_mv[t] = rMv[t]
}
Kr = sparse(a_i, a_j, a_kv, nfree, nfree)
Mr = sparse(a_i, a_j, a_mv, nfree, nfree)

println("Clamped ${ndof - nfree} DOFs; reduced system is ${nfree} x ${nfree}")

Step 6: Solve the Generalized Eigenproblem

Now the heart of modal analysis. The three-argument form of eigs solves the generalized eigenproblem $K\,\phi = \lambda\,M\,\phi$ and returns the lowest k modes: the eigenvalues $\lambda$ in ascending order and the mode shapes as the columns of $\Phi$ . Each eigenvalue is a squared angular frequency, so the natural frequency in hertz is

$f = \frac{\sqrt{\lambda}}{2\pi}.$

Code [6]Run

nmodes = 6
sol = eigs(Kr, Mr, nmodes)
lambdas = sol[1]   # eigenvalues (omega^2)
Phi = sol[2]       # mode shapes (columns)

twopi = 2 * 3.141592653589793
println("=== Natural frequencies ===")
for m in 1:nmodes {
  f = sqrt(abs(lambdas[m])) / twopi
  println("mode ${m}:  f = ${f} Hz")
}

You should see the modes come in a physically meaningful order: a degenerate pair of first bending modes (the square cross-section bends identically in $y$ and $z$ ), a second bending pair, then a torsion mode, and finally the first axial (stretching) mode.

A quick reality check

Two classical formulas let us sanity-check the numbers:

First bending (Euler–Bernoulli cantilever): $f_1 = \dfrac{(1.875)^2}{2\pi}\sqrt{\dfrac{EI}{\rho A L^4}}$
First axial (clamped–free bar): $f = \dfrac{1}{4L}\sqrt{\dfrac{E}{\rho}}$

Code [7]Run

A = Ly*Lz
I = Ly*Lz*Lz*Lz/12
f_bend = (1.875104^2)/twopi * sqrt(E*I/(rho*A*Lx^4))
f_axial = (1.0/(4*Lx)) * sqrt(E/rho)

println("Euler-Bernoulli first bending: ${f_bend} Hz")
println("Analytical first axial mode:   ${f_axial} Hz")

The axial mode matches the bar formula to better than a percent. The bending modes come out somewhat higher than Euler–Bernoulli: trilinear hexahedra are stiff in bending on a coarse mesh (an effect called shear locking), and the beam theory itself ignores shear deformation. Refining the mesh along the length drives the FEM frequencies down toward — and eventually below — the slender-beam estimate. The physics, though, is already right: mode ordering, the degenerate bending pairs, and the axial frequency are all captured.

Step 7: Prepare the Mesh for Visualization

To draw the part with vtk.femgrid we need a flat node array, a connectivity array ([8, n0, …, n7] per hexahedron, using 0-based indices), and a cell-type array (12 = hexahedron). The connectivity is fixed, so we build it once.

Code [8]Run

# Connectivity + cell types for vtk.femgrid (built once, reused per mode)
elems = zeros(9*nel)
types = zeros(nel)
for e in 1:nel {
  base = 9*(e-1)
  elems[base+1] = 8           # 8 nodes per hexahedron
  for a in 1:8 {
    elems[base+1+a] = conn[e, a] - 1   # 0-based node indices
  }
  types[e] = 12               # VTK_HEXAHEDRON
}
println("Prepared ${nel} hexahedra for rendering")

Step 8: Visualize the Mode Shapes

A mode shape is a relative displacement pattern — its amplitude is arbitrary. To make it visible we warp the geometry: displace every node by its mode displacement, scaled so the largest motion is a fraction of the part length, and color the surface by displacement magnitude.

The eigenvector only covers the free DOFs, so we first scatter it back into the full node list (clamped DOFs stay at zero). Change modeid to explore any of the computed modes.

Mode 1 — First Bending

Code [9]Run

modeid = 1

# Expand the reduced eigenvector to all DOFs (fixed DOFs = 0)
ufull = zeros(ndof)
for i in 1:ndof {
  if isfixed[i] < 0.5 {
    ufull[i] = Phi[newidx[i], modeid]
  }
}

# Nodal displacement magnitude + peak
maxmag = 0.0
mag = zeros(nnode)
for nd in 1:nnode {
  ux = ufull[3*(nd-1)+1]
  uy = ufull[3*(nd-1)+2]
  uz = ufull[3*(nd-1)+3]
  mag[nd] = sqrt(ux*ux + uy*uy + uz*uz)
  if mag[nd] > maxmag {
    maxmag = mag[nd]
  }
}

# Warp the geometry: peak displacement = 20% of the part length
warp = 0.2 * Lx / maxmag
nodes = zeros(3*nnode)
for nd in 1:nnode {
  nodes[3*(nd-1)+1] = X0[nd] + warp*ufull[3*(nd-1)+1]
  nodes[3*(nd-1)+2] = Y0[nd] + warp*ufull[3*(nd-1)+2]
  nodes[3*(nd-1)+3] = Z0[nd] + warp*ufull[3*(nd-1)+3]
}

f1 = sqrt(abs(lambdas[modeid])) / twopi
vtk.colormap("viridis")
vtk.edges("on")
vtk.femgrid(nodes, elems, types, mag)
title("Mode 1 — first bending (${f1} Hz)")

Drag to rotate the bracket: the free end swings up and down while the clamped end stays put — the classic cantilever bending shape, colored by how far each point moves.

Mode 5 — Torsion

The same block, now selecting the torsion mode. The cross-section twists about the long axis: points near the axis barely move while the corners sweep the farthest, with essentially no motion along the beam.

Code [10]Run

modeid = 5

ufull = zeros(ndof)
for i in 1:ndof {
  if isfixed[i] < 0.5 {
    ufull[i] = Phi[newidx[i], modeid]
  }
}

maxmag = 0.0
mag = zeros(nnode)
for nd in 1:nnode {
  ux = ufull[3*(nd-1)+1]
  uy = ufull[3*(nd-1)+2]
  uz = ufull[3*(nd-1)+3]
  mag[nd] = sqrt(ux*ux + uy*uy + uz*uz)
  if mag[nd] > maxmag {
    maxmag = mag[nd]
  }
}

warp = 0.2 * Lx / maxmag
nodes = zeros(3*nnode)
for nd in 1:nnode {
  nodes[3*(nd-1)+1] = X0[nd] + warp*ufull[3*(nd-1)+1]
  nodes[3*(nd-1)+2] = Y0[nd] + warp*ufull[3*(nd-1)+2]
  nodes[3*(nd-1)+3] = Z0[nd] + warp*ufull[3*(nd-1)+3]
}

f5 = sqrt(abs(lambdas[modeid])) / twopi
vtk.colormap("coolwarm")
vtk.edges("on")
vtk.femgrid(nodes, elems, types, mag)
title("Mode 5 — torsion (${f5} Hz)")

Mode 6 — Axial (Longitudinal) Stretch

The first mode with motion along the beam: the bracket stretches and compresses like a spring, the displacement growing from zero at the clamp to a maximum at the free end. This is the mode whose frequency matched the $c/(4L)$ bar formula.

Code [11]Run

modeid = 6

ufull = zeros(ndof)
for i in 1:ndof {
  if isfixed[i] < 0.5 {
    ufull[i] = Phi[newidx[i], modeid]
  }
}

maxmag = 0.0
mag = zeros(nnode)
for nd in 1:nnode {
  ux = ufull[3*(nd-1)+1]
  uy = ufull[3*(nd-1)+2]
  uz = ufull[3*(nd-1)+3]
  mag[nd] = sqrt(ux*ux + uy*uy + uz*uz)
  if mag[nd] > maxmag {
    maxmag = mag[nd]
  }
}

warp = 0.2 * Lx / maxmag
nodes = zeros(3*nnode)
for nd in 1:nnode {
  nodes[3*(nd-1)+1] = X0[nd] + warp*ufull[3*(nd-1)+1]
  nodes[3*(nd-1)+2] = Y0[nd] + warp*ufull[3*(nd-1)+2]
  nodes[3*(nd-1)+3] = Z0[nd] + warp*ufull[3*(nd-1)+3]
}

f6 = sqrt(abs(lambdas[modeid])) / twopi
vtk.colormap("viridis")
vtk.edges("on")
vtk.femgrid(nodes, elems, types, mag)
title("Mode 6 — axial stretch (${f6} Hz)")

Summary

We ran a complete 3D structural modal analysis in the browser:

Material — the $6\times6$ isotropic elasticity matrix $D$ from $E$ , $\nu$
Mesh — a structured grid of trilinear hexahedral elements, three DOFs per node
Element matrices — $k_e = \int B^\top D B\,dV$ by Gauss quadrature, and a lumped mass $m_e$
Assembly — scatter the $24\times24$ blocks into sparse $K$ and $M$ via triplets
Boundary conditions — eliminate the clamped face by dropping its DOFs
Eigensolve — eigs(K, M, k) for the lowest generalized modes, $f = \sqrt{\lambda}/2\pi$
Visualize — warp the geometry by each mode shape and color by displacement in an interactive VTK.js viewer

The eigs(K, M, k) solver reduces the generalized problem to a standard symmetric one via the $M^{-1/2}$ congruence, so it works for any symmetric $K$ and positive-definite $M$ — not just this example. Swap in a different mesh, material, or set of boundary conditions and the same seven steps give you the modes of any 3D part.

Workbench

Clear

No variables in workbench

Next Steps

Linear Elasticity

The static counterpart: stress, strain, and element assembly

VTK.js Visualization

More on interactive 3D rendering of meshes and fields