17 — Modal Analysis of a Cantilever Beam¶

Curriculum slot: Tier 5, slot 17.
Prerequisite: 02 — 2D Cantilever Beam, 06 — Sections Catalog.
Strategy used for results: spec.capture(...) — modal stages only live on this path (the classic recorder path can't drive ops.eigen()).

Purpose¶

Tier 5 covers the four major analysis types beyond static-linear: modal, buckling, pushover, time-history. This is the modal entry point.

A linear-elastic cantilever of length $L$ with mass per unit length $\bar{m} = \rho A$ has natural circular frequencies

$$ \omega_n \;=\; (\beta_n L)^{2} \sqrt{\dfrac{E I}{\bar{m}\,L^{4}}}, $$

where $\beta_n L$ satisfies $\cos(\beta L)\cosh(\beta L) + 1 = 0$. First three roots:

Mode	$\beta_n L$	$\omega_n / \omega_1$
1	1.8751	1.000
2	4.6941	6.267
3	7.8548	17.55

Why this notebook is interesting¶

• Modal stages are capture-only. spec.emit_recorders raises immediately if the spec contains a modal record because the classic recorder path can't drive ops.eigen(...). This notebook uses cap.capture_modes(N) — apeGmsh runs the eigenvalue analysis itself and writes one stage per mode.
• Read back via results.modes. Each mode is a stage with kind="mode", exposing mode_index, eigenvalue, frequency_hz, period_s as scoped properties.
• Mode shapes are slabs with T = 1 (one "step" per mode) — same read API, just mode.nodes.get(component="displacement_z").

1. Imports and parameters¶

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import openseespy.opensees as ops

from apeGmsh import apeGmsh, Results
from apeGmsh.results.spec import Recorders

# Geometry / section
L  = 3.0
E  = 2.1e11
nu = 0.3
G  = E / (2.0 * (1.0 + nu))
A, Iy, Iz, J = 1.0e-3, 1.0e-5, 1.0e-5, 2.0e-5
rho   = 7850.0          # steel
m_bar = rho * A          # mass per unit length

# Mesh density — higher than the static cantilever so the higher modes
# get enough segments per half-wavelength.
N_ELEM  = 30
LC      = L / N_ELEM
N_MODES = 3

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import openseespy.opensees as ops from apeGmsh import apeGmsh, Results from apeGmsh.results.spec import Recorders # Geometry / section L = 3.0 E = 2.1e11 nu = 0.3 G = E / (2.0 * (1.0 + nu)) A, Iy, Iz, J = 1.0e-3, 1.0e-5, 1.0e-5, 2.0e-5 rho = 7850.0 # steel m_bar = rho * A # mass per unit length # Mesh density — higher than the static cantilever so the higher modes # get enough segments per half-wavelength. N_ELEM = 30 LC = L / N_ELEM N_MODES = 3

2. Geometry + mesh¶

A line cantilever — same as slot 02 but with a transfinite curve so the element count is exactly N_ELEM.

In [ ]:

g_ctx = apeGmsh(model_name="17_modal", verbose=False)
g = g_ctx.__enter__()

p_base = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC)
p_tip  = g.model.geometry.add_point(L,   0.0, 0.0, lc=LC)
ln     = g.model.geometry.add_line(p_base, p_tip)
g.model.sync()

g.physical.add(0, [p_base], name="base")
g.physical.add(0, [p_tip],  name="tip")
g.physical.add(1, [ln],     name="beam")

g.mesh.structured.set_transfinite_curve(ln, N_ELEM + 1)
g.mesh.generation.generate(1)
fem = g.mesh.queries.get_fem_data()
print(f"mesh: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

g_ctx = apeGmsh(model_name="17_modal", verbose=False) g = g_ctx.__enter__() p_base = g.model.geometry.add_point(0.0, 0.0, 0.0, lc=LC) p_tip = g.model.geometry.add_point(L, 0.0, 0.0, lc=LC) ln = g.model.geometry.add_line(p_base, p_tip) g.model.sync() g.physical.add(0, [p_base], name="base") g.physical.add(0, [p_tip], name="tip") g.physical.add(1, [ln], name="beam") g.mesh.structured.set_transfinite_curve(ln, N_ELEM + 1) g.mesh.generation.generate(1) fem = g.mesh.queries.get_fem_data() print(f"mesh: {fem.info.n_nodes} nodes, {fem.info.n_elems} elements")

3. Build the OpenSees model — vanilla openseespy¶

Standard cantilever, plus restraint of every non-bending DOF on every non-base node. Without that, axial / out-of-plane modes interleave with the bending modes and break the ordering vs the analytical Euler-Bernoulli reference.

In [ ]:

ops.wipe()
ops.model("basic", "-ndm", 3, "-ndf", 6)

for nid, xyz in fem.nodes.get():
    ops.node(int(nid), float(xyz[0]), float(xyz[1]), float(xyz[2]))

ops.geomTransf("Linear", 1, 0.0, 1.0, 0.0)
for group in fem.elements.get(target="beam"):
    for eid, conn in zip(group.ids, group.connectivity):
        ops.element(
            "elasticBeamColumn", int(eid),
            int(conn[0]), int(conn[1]),
            A, E, G, J, Iy, Iz, 1,
        )

# fix the base fully
for n in fem.nodes.get(target="base").ids:
    ops.fix(int(n), 1, 1, 1, 1, 1, 1)

# Free DOFs on non-base nodes: u_z (3), r_y (5) — pure in-plane bending
base_set = {int(n) for n in fem.nodes.get(target="base").ids}
for nid, _ in fem.nodes.get():
    if int(nid) in base_set:
        continue
    ops.fix(int(nid), 1, 1, 0, 1, 0, 1)

ops.wipe() ops.model("basic", "-ndm", 3, "-ndf", 6) for nid, xyz in fem.nodes.get(): ops.node(int(nid), float(xyz[0]), float(xyz[1]), float(xyz[2])) ops.geomTransf("Linear", 1, 0.0, 1.0, 0.0) for group in fem.elements.get(target="beam"): for eid, conn in zip(group.ids, group.connectivity): ops.element( "elasticBeamColumn", int(eid), int(conn[0]), int(conn[1]), A, E, G, J, Iy, Iz, 1, ) # fix the base fully for n in fem.nodes.get(target="base").ids: ops.fix(int(n), 1, 1, 1, 1, 1, 1) # Free DOFs on non-base nodes: u_z (3), r_y (5) — pure in-plane bending base_set = {int(n) for n in fem.nodes.get(target="base").ids} for nid, _ in fem.nodes.get(): if int(nid) in base_set: continue ops.fix(int(nid), 1, 1, 0, 1, 0, 1)

4. Lumped-mass assembly¶

Each non-base node gets the consistent-tributary lumped mass $m_{\text{full}} = \bar{m}\,h$; the tip is a half-tribute. Rotational inertia is negligible for a slender beam.

In [ ]:

h = L / N_ELEM
m_full = m_bar * h
m_half = m_full / 2.0

tip_node = int(next(iter(fem.nodes.get(target="tip").ids)))
for nid in fem.nodes.ids:
    nid = int(nid)
    if nid in base_set:
        continue
    m_here = m_half if nid == tip_node else m_full
    ops.mass(nid, m_here, m_here, m_here, 0.0, 0.0, 0.0)

h = L / N_ELEM m_full = m_bar * h m_half = m_full / 2.0 tip_node = int(next(iter(fem.nodes.get(target="tip").ids))) for nid in fem.nodes.ids: nid = int(nid) if nid in base_set: continue m_here = m_half if nid == tip_node else m_full ops.mass(nid, m_here, m_here, m_here, 0.0, 0.0, 0.0)

5. Declare the modal recorder + capture¶

Recorders().modal(n_modes=N) is a modal-only recorder — no components, no selectors. Modal stages can only be executed via spec.capture(...). Inside the context, cap.capture_modes(N) runs ops.eigen(N) itself and writes one stage per mode.

In [ ]:

recorders = Recorders()
recorders.modal(n_modes=N_MODES)

spec = recorders.resolve(fem, ndm=3, ndf=6)
print(spec)

recorders = Recorders() recorders.modal(n_modes=N_MODES) spec = recorders.resolve(fem, ndm=3, ndf=6) print(spec)

In [ ]:

from apeGmsh import workdir
OUT = workdir()
results_path = OUT / "capture.h5"
if results_path.exists():
    results_path.unlink()

with spec.capture(results_path, fem=fem, ndm=3, ndf=6) as cap:
    cap.capture_modes(n_modes=N_MODES)

print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

from apeGmsh import workdir OUT = workdir() results_path = OUT / "capture.h5" if results_path.exists(): results_path.unlink() with spec.capture(results_path, fem=fem, ndm=3, ndf=6) as cap: cap.capture_modes(n_modes=N_MODES) print(f"wrote {results_path} ({results_path.stat().st_size / 1024:.1f} KB)")

6. Read modes via results.modes¶

Each entry is a Results scoped to one mode stage with mode_index, eigenvalue, frequency_hz, period_s as scoped properties — no need to know the underlying HDF5 layout.

In [ ]:

results = Results.from_native(results_path, fem=fem)
print(results.inspect.summary())

results = Results.from_native(results_path, fem=fem) print(results.inspect.summary())

In [ ]:

modes = sorted(results.modes, key=lambda m: m.mode_index or 0)
print(f"{'mode':>4s}  {'ω [rad/s]':>14s}  {'f [Hz]':>14s}  {'T [s]':>14s}")
for m in modes:
    print(f"{m.mode_index:>4d}  {np.sqrt(m.eigenvalue):>14.6e}  "
          f"{m.frequency_hz:>14.6e}  {m.period_s:>14.6e}")

modes = sorted(results.modes, key=lambda m: m.mode_index or 0) print(f"{'mode':>4s} {'ω [rad/s]':>14s} {'f [Hz]':>14s} {'T [s]':>14s}") for m in modes: print(f"{m.mode_index:>4d} {np.sqrt(m.eigenvalue):>14.6e} " f"{m.frequency_hz:>14.6e} {m.period_s:>14.6e}")

7. Verification against analytical $\omega_n$¶

In [ ]:

beta_L = np.array([1.87510407, 4.69409113, 7.85475744])
omega_analytical = (beta_L**2) * np.sqrt(E * Iz / (m_bar * L**4))

print(f"{'mode':>4s}  {'FEM ω':>14s}  {'analyt ω':>14s}  {'err %':>10s}")
for i, m in enumerate(modes):
    w_fem = np.sqrt(m.eigenvalue)
    w_an  = omega_analytical[i]
    err   = abs(w_fem - w_an) / w_an * 100.0
    print(f"{m.mode_index:>4d}  {w_fem:>14.4e}  {w_an:>14.4e}  {err:>9.4f} %")

beta_L = np.array([1.87510407, 4.69409113, 7.85475744]) omega_analytical = (beta_L**2) * np.sqrt(E * Iz / (m_bar * L**4)) print(f"{'mode':>4s} {'FEM ω':>14s} {'analyt ω':>14s} {'err %':>10s}") for i, m in enumerate(modes): w_fem = np.sqrt(m.eigenvalue) w_an = omega_analytical[i] err = abs(w_fem - w_an) / w_an * 100.0 print(f"{m.mode_index:>4d} {w_fem:>14.4e} {w_an:>14.4e} {err:>9.4f} %")

8. Plot the mode shapes¶

Each mode is a stage with T = 1 — pulling displacement_z returns a single-row slab whose shape over $x$ is the mode shape (normalised by OpenSees so the maximum entry has unit absolute value).

In [ ]:

tag_to_idx = {int(t): i for i, t in enumerate(fem.nodes.ids)}

fig, ax = plt.subplots(figsize=(7, 5))
for i, m in enumerate(modes):
    slab = m.nodes.get(component="displacement_z")
    x = np.array([fem.nodes.coords[tag_to_idx[int(n)], 0]
                  for n in slab.node_ids])
    u = slab.values[0]
    order = np.argsort(x)
    sign = np.sign(u[order][-1]) or 1.0   # tip-positive convention
    ax.plot(
        x[order], sign * u[order],
        "o-", ms=3, color=f"C{i}",
        label=f"mode {m.mode_index}  ({m.frequency_hz:.1f} Hz)",
    )

ax.axhline(0.0, color="gray", lw=0.5)
ax.set_xlabel("x  [m]"); ax.set_ylabel("normalised u_z")
ax.set_title("Cantilever — first three bending mode shapes")
ax.legend(); ax.grid(True, alpha=0.3)
fig.tight_layout()

tag_to_idx = {int(t): i for i, t in enumerate(fem.nodes.ids)} fig, ax = plt.subplots(figsize=(7, 5)) for i, m in enumerate(modes): slab = m.nodes.get(component="displacement_z") x = np.array([fem.nodes.coords[tag_to_idx[int(n)], 0] for n in slab.node_ids]) u = slab.values[0] order = np.argsort(x) sign = np.sign(u[order][-1]) or 1.0 # tip-positive convention ax.plot( x[order], sign * u[order], "o-", ms=3, color=f"C{i}", label=f"mode {m.mode_index} ({m.frequency_hz:.1f} Hz)", ) ax.axhline(0.0, color="gray", lw=0.5) ax.set_xlabel("x [m]"); ax.set_ylabel("normalised u_z") ax.set_title("Cantilever — first three bending mode shapes") ax.legend(); ax.grid(True, alpha=0.3) fig.tight_layout()

9. (Optional) Open the post-solve viewer¶

The viewer's stage selector picks per mode — switch through them to see each shape rendered in the deformed-beam view.

In [ ]:

results.viewer(blocking=False)

results.viewer(blocking=False)

What this unlocks¶

• Modal as a first-class apeGmsh strategy. One call to cap.capture_modes(N) writes N stages of kind="mode" — read side surfaces them via results.modes.
• Why modal is capture-only. The classic recorder pipeline can't drive ops.eigen(...), so spec.emit_recorders raises if the spec contains a modal record. Use spec.capture for any analysis with modal stages.
• Stage-scoped properties — mode_index, eigenvalue, frequency_hz, period_s on each scoped Results.
• Lumped-mass tributary pattern — half-mass at the tip, full mass at internal nodes. Generalises to inhomogeneous beams and shells.

Next: notebook 19 — pushover with a fiber section. The first nonlinear example of the curriculum.

In [ ]:

g_ctx.__exit__(None, None, None)

g_ctx.__exit__(None, None, None)