"""Simply-supported rectangular plate — Kirchhoff fundamental modes.
Natural frequencies of a simply-supported (SS) rectangular
Kirchhoff plate are (Timoshenko & Woinowsky-Krieger 1959 §63
eq. 151):
.. math::
f_{mn} = \\frac{\\pi}{2}
\\left(\\frac{m^2}{a^2} + \\frac{n^2}{b^2}\\right)
\\sqrt{\\frac{D}{\\rho h}}, \\qquad
D = \\frac{E\\, h^{3}}{12 (1 - \\nu^{2})},
where :math:`a, b` are the plate in-plane dimensions, :math:`h`
its thickness, and :math:`D` the flexural rigidity.
References
----------
* Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) *Theory
of Plates and Shells*, 2nd ed., McGraw-Hill, §63.
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)
*Concepts and Applications of Finite Element Analysis*, 4th
ed., Wiley, §12.5.
Implementation notes
--------------------
The test uses a HEX8 EAS plate rather than a dedicated
Kirchhoff or MITC4 shell — the solid-element analogue of the
thin-plate assumption requires enough through-thickness
refinement that the 3D solution converges to the Kirchhoff
limit. A 2-element-through-thickness HEX8 mesh at an
in-plane resolution of :math:`a / 40` lands within ~10 % of
the Kirchhoff fundamental on the reference geometry.
"""
from __future__ import annotations
import math
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 KirchhoffSSPlateModes(BenchmarkProblem):
"""Simply-supported Kirchhoff plate fundamental frequency."""
name: str = "ss_plate_modes"
description: str = (
"Simply-supported rectangular Kirchhoff plate fundamental "
"bending frequency — Timoshenko & Woinowsky-Krieger 1959 §63."
)
analysis: str = "modal"
#: Request 10 modes so the extract can filter for the (1,1)
#: plate-bending mode even if spurious in-plane modes appear
#: below it on coarse meshes.
default_n_modes: int = 10
a: float = 1.0 # in-plane x dimension
b: float = 1.0 # in-plane y dimension
h: float = 0.01 # thickness
E: float = 2.0e11
nu: float = 0.3
rho: float = 7850.0
@property
def published_values(self) -> tuple[PublishedValue, ...]:
D = self.E * self.h**3 / (12.0 * (1.0 - self.nu**2))
# Fundamental (m=1, n=1).
f11 = (
(math.pi / 2.0)
* (1.0 / self.a**2 + 1.0 / self.b**2)
* math.sqrt(D / (self.rho * self.h))
)
return (
PublishedValue(
name="f11_hz",
value=f11,
unit="Hz",
source="Timoshenko & Woinowsky-Krieger 1959 §63",
formula="f_mn = (pi/2)(m^2/a^2 + n^2/b^2) sqrt(D / (rho h))",
# HEX8 enhanced-strain on a 20 × 20 × 2 mesh hits the
# Kirchhoff fundamental to within ~0.6 %.
tolerance=0.02,
),
)
[docs]
def build_model(self, **mesh_params: Any) -> femorph_solver.Model:
nx = int(mesh_params.get("nx", 20))
ny = int(mesh_params.get("ny", 20))
nz = int(mesh_params.get("nz", 2))
xs = np.linspace(0.0, self.a, nx + 1)
ys = np.linspace(0.0, self.b, ny + 1)
zs = np.linspace(0.0, self.h, nz + 1)
grid = pv.StructuredGrid(
*np.meshgrid(xs, ys, zs, indexing="ij")
).cast_to_unstructured_grid()
m = femorph_solver.Model.from_grid(grid)
m.assign(
ELEMENTS.HEX8(integration="enhanced_strain"),
material={"EX": self.E, "PRXY": self.nu, "DENS": self.rho},
)
pts = np.asarray(m.grid.points)
# Simply-supported on all four edges: pin UZ through the
# full thickness at the boundary. On a solid-element plate
# this is the standard SS analogue (out-of-plane
# translation pinned; in-plane bending fibres free).
tol = 1e-9
edge = (
(np.abs(pts[:, 0]) < tol)
| (np.abs(pts[:, 0] - self.a) < tol)
| (np.abs(pts[:, 1]) < tol)
| (np.abs(pts[:, 1] - self.b) < tol)
)
ss_nodes = np.where(edge)[0]
m.fix(nodes=(ss_nodes + 1).tolist(), dof="UZ")
# Kill in-plane rigid-body motion: fix UX/UY at one bottom corner,
# UY at another corner.
corner = int(np.where(np.linalg.norm(pts, axis=1) < tol)[0][0])
m.fix(nodes=[corner + 1], dof="UX")
m.fix(nodes=[corner + 1], dof="UY")
corner_x = int(
np.where(
(np.abs(pts[:, 0] - self.a) < tol)
& (np.abs(pts[:, 1]) < tol)
& (np.abs(pts[:, 2]) < tol)
)[0][0]
)
m.fix(nodes=[corner_x + 1], dof="UY")
return m