"""Fixed-free axial-rod natural frequency.
Reference geometry
------------------
Slender uniform rod of length :math:`L`, cross-sectional area
:math:`A`, material density :math:`\\rho`, Young's modulus
:math:`E`. Fixed at :math:`x = 0` (all three translational DOFs
pinned), free at :math:`x = L`.
Closed-form quantities
----------------------
The fixed-free rod satisfies the wave equation
.. math::
\\rho A \\, \\ddot{u} = E A \\, u'',
with boundary conditions :math:`u(0) = 0` and
:math:`E A\\,u'(L) = 0`. Separation of variables yields
mode shapes :math:`u_n(x) = \\sin\\!\\bigl((2n-1) \\pi x /
(2L)\\bigr)` and natural frequencies
.. math::
f_n = \\frac{2n - 1}{4 L} \\sqrt{\\frac{E}{\\rho}},
\\qquad n = 1, 2, 3, \\ldots
The fundamental (n = 1) is :math:`f_1 = \\sqrt{E/\\rho} / (4L)`.
References
----------
* Rao, S. S. *Mechanical Vibrations*, 6th ed., Pearson, 2017,
§8.2 — longitudinal vibration of a bar.
* Meirovitch, L. *Fundamentals of Vibrations*, Waveland, 2010,
§6.6 — axial vibration of a rod with mixed BC.
Cross-references (public verification manuals — fair-use
citation of problem IDs only; no vendor content vendored):
* **Abaqus AVM 1.6.x** bar-natural-frequency family.
Implementation notes
--------------------
Modelled with a string of ``TRUSS2`` elements in line along
the global x-axis. The truss kernel carries three
translational DOFs per node but only axial strain contributes
to the element stiffness; transverse DOFs are suppressed via
``fix(..., dof=["UY", "UZ"])`` on every interior node so the
only surviving mode family is the axial one we want to
verify. Without this suppression the pin-transverse rigid
spring is zero and the solver returns near-zero transverse
frequencies first.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any
import numpy as np
import pyvista as pv
import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.validation._benchmark import (
BenchmarkProblem,
PublishedValue,
)
[docs]
@dataclass
class AxialRodNaturalFrequency(BenchmarkProblem):
"""Fixed-free axial-rod fundamental natural frequency."""
name: str = "axial_rod_nf"
description: str = (
"Fixed-free slender rod, fundamental axial natural "
"frequency f_1 = (1/4L) sqrt(E/rho) (Rao 2017 §8.2)."
)
analysis: str = "modal"
default_n_modes: int = 1
#: Rod length [m].
L: float = 1.0
#: Cross-sectional area [m²] (any positive value works — the
#: axial natural frequency is independent of ``A``).
area: float = 1.0e-4
#: Young's modulus [Pa].
E: float = 2.0e11
#: Poisson's ratio (unused by the truss kernel but required
#: by the material stamp).
nu: float = 0.3
#: Density [kg/m³].
rho: float = 7850.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
f1 = np.sqrt(self.E / self.rho) / (4.0 * self.L)
return (
PublishedValue(
name="f1_axial",
value=f1,
unit="Hz",
source="Rao 2017 §8.2",
formula="f1 = sqrt(E/rho) / (4 L)",
# Linear-Lagrange truss elements converge quickly
# for axial modes — 1 % on a 20-element mesh is
# the usual floor.
tolerance=0.02,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
n_elem = int(mesh_params.get("n_elem", 40))
# Build a line of n_elem 2-node truss elements along +x.
xs = np.linspace(0.0, self.L, n_elem + 1)
pts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))
# VTK_LINE = 3; each cell carries two node indices.
cells = np.empty((n_elem, 3), dtype=np.int64)
cells[:, 0] = 2 # two nodes per cell
cells[:, 1] = np.arange(n_elem)
cells[:, 2] = np.arange(1, n_elem + 1)
grid = pv.UnstructuredGrid(
cells.ravel(),
np.full(n_elem, 3, dtype=np.int64), # VTK_LINE
pts,
)
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.TRUSS2,
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
real=(self.area,),
)
# Fix the left end in all three translations.
m.fix(nodes=1, dof="ALL")
# Suppress transverse translations (UY, UZ) on every other
# node so the axial family is the only mode family that
# survives; truss elements have zero transverse stiffness
# and would otherwise admit rigid-body transverse modes at
# zero frequency ahead of the axial fundamental.
for i in range(2, n_elem + 2):
m.fix(nodes=i, dof="UY", value=0.0)
m.fix(nodes=i, dof="UZ", value=0.0)
m._bench_n_elem = n_elem # type: ignore[attr-defined]
return m