"""BEAM188 — 3D 2-node linear beam.
femorph-solver implements the slender-beam (Euler–Bernoulli) form of BEAM188:
axial, torsional, and two independent bending planes combine into a
12×12 local stiffness matrix rotated into global coordinates. This is
exactly the limit BEAM188 converges to when the Timoshenko shear
parameter Φ → 0, i.e. for slender sections where bending dominates over
shear deformation. A future revision can introduce Φ from an explicit
shear area or section integration.
Local DOF ordering at each node:
(u, v, w, θx, θy, θz)
so global assembly fills rows/cols 0..5 (node i) and 6..11 (node j) in
MAPDL's canonical ``(UX, UY, UZ, ROTX, ROTY, ROTZ)`` order.
Real constants
--------------
real[0]: AREA — cross-sectional area A
real[1]: IZZ — bending inertia about local z (resists y-deflection)
real[2]: IYY — bending inertia about local y (resists z-deflection)
real[3]: J — Saint-Venant torsion constant
Material
--------
EX : Young's modulus
PRXY : Poisson's ratio (→ G = E/(2(1+ν)) for torsion)
DENS : mass density (required for M_e)
Orientation
-----------
The local x-axis points from node I → node J. Without an explicit
orientation node, the local y-axis is chosen so that local z stays
"up" (world +Z) whenever the beam is not parallel to the global Z-axis,
falling back to world +Y for vertical beams. For axisymmetric sections
(``IYY == IZZ``) the orientation is immaterial.
Stiffness (Theory Ref § 14.188 slender limit)
---------------------------------------------
K_local = [ K_axial · · · ]
[ · K_bend_xy · · ]
[ · · K_bend_xz · ]
[ · · · K_torsion ]
with the familiar Hermite-cubic 4×4 bending blocks. A block-diagonal
direction-cosine matrix ``Λ`` rotates it: ``K_global = Λᵀ K_local Λ``.
Mass (consistent, Hermite cubic)
--------------------------------
Axial / torsional: ``ρ A L / 6 · [[2,1],[1,2]]`` (resp. ``ρ I_p L / 6``
with ``I_p ≈ I_y + I_z``). Bending uses the standard 4×4 block
``ρ A L / 420 · (…)`` in both planes.
Lumped mass puts ``ρ A L / 2`` on each translational DOF and zeros the
rotational ones — a simple diagonal approximation useful for explicit
dynamics but less accurate for modal work than the consistent form.
References
----------
* **Euler--Bernoulli beam** (slender-beam limit, shear neglected):
Timoshenko, S.P., *Vibration Problems in Engineering*, 4th ed.,
Wiley, 1974, §5.3. Beam kinematics and derivation of the
12-DOF stiffness / consistent-mass blocks: Cook, Malkus, Plesha,
Witt, *Concepts and Applications of Finite Element Analysis*,
4th ed., Wiley, 2002, §2.4 (axial), §2.5 (bending Hermite
cubics), §2.6 (torsion).
* Hermite cubic shape functions for transverse displacement and
slope continuity at the nodes: Zienkiewicz, O.C. and Taylor, R.L.
*The Finite Element Method: Its Basis and Fundamentals*, 7th
ed., Butterworth-Heinemann, 2013, §2.5.1 eqs. (2.26)–(2.27).
* Consistent mass matrix from Hermite cubics (coefficients
``ρAL/420`` on the 4×4 bending block and ``ρAL/6`` on the
axial / torsional block): Cook et al. Table 16.3-1.
* Timoshenko-beam form (non-zero shear, Φ = 12EI/(κAGL²) — future
revision) and comparison with the Euler-Bernoulli limit:
Bathe, K.J., *Finite Element Procedures*, 2nd ed., Prentice Hall,
2014, §5.4.
MAPDL compatibility — specification source
------------------------------------------
* Ansys, Inc., *Ansys Mechanical APDL Element Reference*,
Release 2022R2, section "BEAM188 — 3-D 2-Node Beam".
Short factual summary (paraphrased): 2-node straight beam; 6
DOFs per node (3 translations + 3 rotations); cross-section
properties supplied via ``REAL`` constants (``AREA``, ``IZZ``,
``IYY``, ``J``); optional local-z orientation node.
femorph-solver reproduces the slender-beam (Euler-Bernoulli)
limit; the transverse-shear (Timoshenko) form is roadmap. The
Ansys Element Reference is the compatibility spec only — theory
is Cook / Timoshenko / Bathe cited above.
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
def _beam_frame(coords: np.ndarray) -> tuple[float, np.ndarray]:
"""Return ``(L, R)`` — beam length and 3×3 direction-cosine matrix."""
d = coords[1] - coords[0]
L = float(np.linalg.norm(d))
if L == 0.0:
raise ValueError("BEAM188: coincident nodes, element length is zero")
ex = d / L
world_z = np.array([0.0, 0.0, 1.0])
if abs(np.dot(ex, world_z)) > 0.99:
ref = np.array([0.0, 1.0, 0.0])
else:
ref = world_z
ey = np.cross(ref, ex)
ey /= np.linalg.norm(ey)
ez = np.cross(ex, ey)
R = np.stack([ex, ey, ez])
return L, R
def _transform_block(R: np.ndarray) -> np.ndarray:
"""Return the 12×12 block-diagonal transformation matrix ``Λ``."""
T = np.zeros((12, 12), dtype=np.float64)
for i in range(4):
T[3 * i : 3 * i + 3, 3 * i : 3 * i + 3] = R
return T
def _k_local(L: float, E: float, G: float, A: float, Iy: float, Iz: float, J: float) -> np.ndarray:
k = np.zeros((12, 12), dtype=np.float64)
# Axial: u_i (0), u_j (6)
ka = E * A / L
k[0, 0] = k[6, 6] = ka
k[0, 6] = k[6, 0] = -ka
# Torsion: θx_i (3), θx_j (9)
kt = G * J / L
k[3, 3] = k[9, 9] = kt
k[3, 9] = k[9, 3] = -kt
# Bending in local x-y plane (v, θz) — uses I_z
# DOFs: v_i (1), θz_i (5), v_j (7), θz_j (11)
a = E * Iz / L**3
b = E * Iz / L**2
c = E * Iz / L
dofs_z = [1, 5, 7, 11]
kb = np.array(
[
[12 * a, 6 * b, -12 * a, 6 * b],
[6 * b, 4 * c, -6 * b, 2 * c],
[-12 * a, -6 * b, 12 * a, -6 * b],
[6 * b, 2 * c, -6 * b, 4 * c],
]
)
for ii, I in enumerate(dofs_z):
for jj, J_ in enumerate(dofs_z):
k[I, J_] += kb[ii, jj]
# Bending in local x-z plane (w, θy) — uses I_y, sign flips from dw/dx = -θy
# DOFs: w_i (2), θy_i (4), w_j (8), θy_j (10)
a = E * Iy / L**3
b = E * Iy / L**2
c = E * Iy / L
dofs_y = [2, 4, 8, 10]
kb = np.array(
[
[12 * a, -6 * b, -12 * a, -6 * b],
[-6 * b, 4 * c, 6 * b, 2 * c],
[-12 * a, 6 * b, 12 * a, 6 * b],
[-6 * b, 2 * c, 6 * b, 4 * c],
]
)
for ii, I in enumerate(dofs_y):
for jj, J_ in enumerate(dofs_y):
k[I, J_] += kb[ii, jj]
return k
def _m_local_consistent(L: float, rho: float, A: float, Iy: float, Iz: float) -> np.ndarray:
m = np.zeros((12, 12), dtype=np.float64)
mass = rho * A * L
# Axial
ka = mass / 6.0
m[0, 0] = m[6, 6] = 2 * ka
m[0, 6] = m[6, 0] = ka
# Torsional rotary inertia: I_p ≈ I_y + I_z per unit length
kt = rho * (Iy + Iz) * L / 6.0
m[3, 3] = m[9, 9] = 2 * kt
m[3, 9] = m[9, 3] = kt
# Bending x-y (v, θz)
base = mass / 420.0
dofs_z = [1, 5, 7, 11]
mb = base * np.array(
[
[156, 22 * L, 54, -13 * L],
[22 * L, 4 * L**2, 13 * L, -3 * L**2],
[54, 13 * L, 156, -22 * L],
[-13 * L, -3 * L**2, -22 * L, 4 * L**2],
]
)
for ii, I in enumerate(dofs_z):
for jj, J_ in enumerate(dofs_z):
m[I, J_] += mb[ii, jj]
# Bending x-z (w, θy) — sign flips on rotation coupling
dofs_y = [2, 4, 8, 10]
mb = base * np.array(
[
[156, -22 * L, 54, 13 * L],
[-22 * L, 4 * L**2, -13 * L, -3 * L**2],
[54, -13 * L, 156, 22 * L],
[13 * L, -3 * L**2, 22 * L, 4 * L**2],
]
)
for ii, I in enumerate(dofs_y):
for jj, J_ in enumerate(dofs_y):
m[I, J_] += mb[ii, jj]
return m
def _m_local_lumped(L: float, rho: float, A: float) -> np.ndarray:
m = np.zeros((12, 12), dtype=np.float64)
half = 0.5 * rho * A * L
for i in (0, 1, 2, 6, 7, 8):
m[i, i] = half
return m
[docs]
@register
class BEAM188(ElementBase):
name = "BEAM188"
n_nodes = 2
n_dof_per_node = 6
dof_labels = ("UX", "UY", "UZ", "ROTX", "ROTY", "ROTZ")
vtk_cell_type = 3 # VTK_LINE
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(
f"BEAM188 requires REAL[0..3] = AREA, IZZ, IYY, J; got {r.size} values"
)
E = float(material["EX"])
nu = float(material.get("PRXY", 0.3))
G = E / (2.0 * (1.0 + nu))
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
J = float(r[3])
k_loc = _k_local(L, E, G, A, Iy, Iz, J)
T = _transform_block(R)
return T.T @ k_loc @ T
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
L, R = _beam_frame(np.asarray(coords, dtype=np.float64))
r = np.asarray(real, dtype=np.float64)
if r.size < 4:
raise ValueError(
f"BEAM188 requires REAL[0..3] = AREA, IZZ, IYY, J; got {r.size} values"
)
rho = float(material["DENS"])
A = float(r[0])
Iz = float(r[1])
Iy = float(r[2])
if lumped:
m_loc = _m_local_lumped(L, rho, A)
else:
m_loc = _m_local_consistent(L, rho, A, Iy, Iz)
T = _transform_block(R)
return T.T @ m_loc @ T