BEAM188 — cantilever tip deflection and first mode

Slender steel cantilever modelled with a line of BEAM188 elements. Two validations:

1. Static tip deflection matches P L³ / (3 E I) (Euler–Bernoulli).

2. First natural frequency matches (β₁ L)² √(E I / (ρ A L⁴)) with β₁ L = 1.87510407.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem data

Steel, square 50 × 50 mm section, 1 m span, 1 kN tip load.

E = 2.1e11  # Pa
NU = 0.3
RHO = 7850.0  # kg/m³
b = 0.050  # m (square side)
A = b * b
I = b**4 / 12.0  # about either principal axis (square)
J = 2.0 * I  # thin-square torsion approximation
L = 1.0  # m
P = 1.0e3  # N (transverse tip load, +y)

N_ELEM = 10  # discretisation along the span

Build the model

BEAM188 has 6 DOFs per node; d(node, "ALL") on the clamped end fixes all six. Hermite-cubic shape functions recover Euler–Bernoulli exactly for prismatic beams, so 10 elements is already machine-precise on the static answer.

m = femorph_solver.Model()
m.et(1, "BEAM188")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
m.r(1, A, I, I, J)

for i in range(N_ELEM + 1):
    m.n(i + 1, i * L / N_ELEM, 0.0, 0.0)
for i in range(N_ELEM):
    m.e(i + 1, i + 2)

m.d(1, "ALL")  # fully clamp node 1
m.f(N_ELEM + 1, "FY", P)  # tip load in +y

/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:46: DeprecationWarning: Model.et(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).et(et_id, name)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material)` for new code.
  m.et(1, "BEAM188")
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:47: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("EX", 1, E)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:48: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("PRXY", 1, NU)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:49: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("DENS", 1, RHO)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:50: DeprecationWarning: Model.r(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).r(real_id, *values)` for line-by-line APDL deck porting, or the native `Model.assign(element, material, real=[...])` for new code.
  m.r(1, A, I, I, J)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:53: DeprecationWarning: Model.n(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).n(num, x, y, z)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.n(i + 1, i * L / N_ELEM, 0.0, 0.0)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:55: DeprecationWarning: Model.e(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).e(*node_nums)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.e(i + 1, i + 2)
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:57: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(1, "ALL")  # fully clamp node 1
/home/runner/_work/solver/solver/examples/elements/beam188/example_beam188.py:58: DeprecationWarning: Model.f(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).f(node, label, value)` for line-by-line APDL deck porting, or the native `Model.apply_force(node, fx=..., fy=..., fz=...)` for new code.
  m.f(N_ELEM + 1, "FY", P)  # tip load in +y

Static solve + analytical comparison

static = m.solve()
dof = m.dof_map()
tip_uy = np.where((dof[:, 0] == N_ELEM + 1) & (dof[:, 1] == 1))[0][0]
u_tip = static.displacement[tip_uy]
u_expected = P * L**3 / (3.0 * E * I)

print(f"BEAM188 tip UY       = {u_tip:.6e} m")
print(f"Analytical PL³/(3EI) = {u_expected:.6e} m")
assert np.isclose(u_tip, u_expected, rtol=1e-8)

BEAM188 tip UY       = 3.047619e-03 m
Analytical PL³/(3EI) = 3.047619e-03 m

Modal solve + analytical comparison

The transverse DOFs in two planes give two (degenerate) first bending modes. Because I_y = I_z here, femorph-solver’s eigensolver will return them both at the same frequency.

modal = m.modal_solve(n_modes=4)
BETA_L = 1.87510407  # dimensionless cantilever eigenvalue (first mode)
omega_expected = BETA_L**2 * np.sqrt(E * I / (RHO * A * L**4))
omega1 = float(np.sqrt(modal.omega_sq[0]))
omega2 = float(np.sqrt(modal.omega_sq[1]))

print(f"Expected ω₁  = {omega_expected:.4f} rad/s")
print(f"Computed ω₁  = {omega1:.4f} rad/s")
print(f"Computed ω₂  = {omega2:.4f} rad/s (bending in the other plane)")
assert np.isclose(omega1, omega_expected, rtol=5e-3)

Expected ω₁  = 262.4853 rad/s
Computed ω₁  = 262.4855 rad/s
Computed ω₂  = 262.4855 rad/s (bending in the other plane)

Plot the deflected shape and the first mode

Scatter the static displacement onto the grid for visualisation, then overlay the first mode shape coloured by transverse amplitude.

grid = m.grid.copy()
disp = np.zeros((grid.n_points, 3), dtype=np.float64)
mode_disp = np.zeros_like(disp)
for i, nn in enumerate(grid.point_data["ansys_node_num"]):
    rows = np.where(dof[:, 0] == int(nn))[0]
    for r in rows:
        d_idx = int(dof[r, 1])
        if d_idx < 3:  # translations only for warping
            disp[i, d_idx] = static.displacement[r]
            mode_disp[i, d_idx] = modal.mode_shapes[r, 0]

# Scale mode to unit peak
peak = float(np.max(np.abs(mode_disp))) or 1.0
mode_disp /= peak

grid.point_data["static_disp"] = disp
grid.point_data["mode1_disp"] = mode_disp

plotter = pv.Plotter(shape=(1, 2), off_screen=True)
plotter.subplot(0, 0)
plotter.add_text("Static: PL³/(3EI)", font_size=10)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=2)
plotter.add_mesh(
    grid.warp_by_vector("static_disp", factor=50.0),
    scalars=np.linalg.norm(disp, axis=1),
    line_width=5,
    scalar_bar_args={"title": "|u| [m]"},
)
plotter.add_axes()

plotter.subplot(0, 1)
plotter.add_text(f"Mode 1: f = {modal.frequency[0]:.1f} Hz", font_size=10)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=2)
plotter.add_mesh(
    grid.warp_by_vector("mode1_disp", factor=0.3),
    scalars=np.linalg.norm(mode_disp, axis=1),
    line_width=5,
    cmap="coolwarm",
    scalar_bar_args={"title": "|φ|"},
)
plotter.add_axes()
plotter.link_views()
plotter.show()