Simply-supported plate — first transverse bending mode

Classical thin-plate eigenvalue problem. A simply-supported square plate of side \(a\) and thickness \(h\) has Kirchhoff natural frequencies (Timoshenko & Woinowsky-Krieger 1959 §63 eq. 151)

\[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})},\]

so the (1, 1) fundamental of a square plate \(a = b\) is

\[f_{11} = \pi \frac{1}{a^{2}} \sqrt{\frac{D}{\rho\, h}}.\]

Implementation

A 30 × 30 × 2 HEX8 enhanced-strain (Simo–Rifai 1990) mesh on the unit square plate at \(L/h = 100\). Simply-supported boundary conditions pin \(u_z = 0\) along all four edges through the full thickness; in-plane rigid-body translation is removed at one corner. The mode-shape filter selects the first mode whose UZ root-mean-square dominates the in-plane components (i.e. the (1, 1) bending mode), in case any spurious in-plane modes precede it on a coarse mesh.

Modal analysis runs through Model.modal_solve() (ARPACK shift-invert via femorph_solver.solvers.modal).

References

• Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) Theory of Plates and Shells, 2nd ed., McGraw-Hill, §63 eq. 151 (rectangular plate eigenfrequencies).

• 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.

• Simo, J. C. and Rifai, M. S. (1990) “A class of mixed assumed-strain methods …” IJNME 29 (8), 1595–1638.

from __future__ import annotations

import math

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

Problem data

E = 2.0e11  # Pa (steel)
NU = 0.3
RHO = 7850.0  # kg/m^3
a = 1.0  # plate edge length [m]
h = 0.01  # thickness [m]; L/h = 100

NX = NY = 30
NZ = 2

D = E * h**3 / (12.0 * (1.0 - NU**2))
f11_published = (math.pi / 2.0) * (1.0 / a**2 + 1.0 / a**2) * math.sqrt(D / (RHO * h))
print(f"a = {a} m, h = {h} m, L/h = {a / h:.0f}")
print(f"E = {E / 1e9:.0f} GPa, ν = {NU}, ρ = {RHO} kg/m^3, D = {D:.3e} N m")
print(f"f_11 published (Timoshenko §63)   = {f11_published:.4f} Hz")

a = 1.0 m, h = 0.01 m, L/h = 100
E = 200 GPa, ν = 0.3, ρ = 7850.0 kg/m^3, D = 1.832e+04 N m
f_11 published (Timoshenko §63)   = 47.9865 Hz

Build a 30 × 30 × 2 HEX8 mesh

xs = np.linspace(0.0, a, NX + 1)
ys = np.linspace(0.0, a, NY + 1)
zs = np.linspace(0.0, h, NZ + 1)
grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing="ij")).cast_to_unstructured_grid()
print(f"mesh: {grid.n_points} nodes, {grid.n_cells} HEX8 cells")

m = femorph_solver.Model.from_grid(grid)
m.assign(
    ELEMENTS.HEX8(integration="enhanced_strain"),
    material={"EX": E, "PRXY": NU, "DENS": RHO},
)

mesh: 2883 nodes, 1800 HEX8 cells

SS BCs + rigid-body suppression

pts = np.asarray(m.grid.points)
tol = 1e-9
edge = (
    (np.abs(pts[:, 0]) < tol)
    | (np.abs(pts[:, 0] - a) < tol)
    | (np.abs(pts[:, 1]) < tol)
    | (np.abs(pts[:, 1] - a) < tol)
)
m.fix(nodes=(np.where(edge)[0] + 1).tolist(), dof="UZ")
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] - a) < tol) & (np.abs(pts[:, 1]) < tol) & (np.abs(pts[:, 2]) < tol))[
        0
    ][0]
)
m.fix(nodes=[corner_x + 1], dof="UY")

Modal solve + filter for (1, 1) bending mode

n_modes = 10
res = m.modal_solve(n_modes=n_modes)
freqs = np.asarray(res.frequency, dtype=np.float64)
shapes = np.asarray(res.mode_shapes).reshape(-1, 3, n_modes)

f11_idx = None
for k in range(n_modes):
    ux_rms = float(np.sqrt((shapes[:, 0, k] ** 2).mean()))
    uy_rms = float(np.sqrt((shapes[:, 1, k] ** 2).mean()))
    uz_rms = float(np.sqrt((shapes[:, 2, k] ** 2).mean()))
    if uz_rms > 3.0 * max(ux_rms, uy_rms):
        f11_idx = k
        break

assert f11_idx is not None, f"no UZ-dominant mode found in first {n_modes}"

f11_computed = float(freqs[f11_idx])
err_pct = (f11_computed - f11_published) / f11_published * 100.0
print()
print(f"first 5 frequencies (Hz)        = {[round(float(x), 3) for x in freqs[:5]]}")
print(f"selected (1,1) bending mode idx = {f11_idx}")
print(f"f_11 computed                   = {f11_computed:.4f} Hz")
print(f"f_11 published                  = {f11_published:.4f} Hz")
print(f"relative error                  = {err_pct:+.2f} %")

# 0.5 % tolerance — HEX8 enhanced-strain on the default 30 × 30 × 2
# mesh recovers the Kirchhoff fundamental to within ~0.06 % at the
# default thin geometry (L/h = 100).  Tighter convergence (sub-0.1 %)
# would need a dedicated SHELL kernel (tracked separately).
assert abs(err_pct) < 0.5, f"f_11 deviation {err_pct:.2f}% exceeds 0.5%"

first 5 frequencies (Hz)        = [48.016, 120.495, 120.495, 193.089, 242.576]
selected (1,1) bending mode idx = 0
f_11 computed                   = 48.0161 Hz
f_11 published                  = 47.9865 Hz
relative error                  = +0.06 %

Render the (1, 1) bending mode shape

phi = shapes[:, :, f11_idx]
warped = m.grid.copy()
warped.points = m.grid.points + (h / np.max(np.abs(phi)) * 0.4) * phi
warped["UZ"] = phi[:, 2]

plotter = pv.Plotter(off_screen=True, window_size=(720, 480))
plotter.add_mesh(
    m.grid,
    style="wireframe",
    color="grey",
    opacity=0.35,
    label="undeformed",
)
plotter.add_mesh(
    warped,
    scalars="UZ",
    cmap="coolwarm",
    show_edges=False,
    scalar_bar_args={"title": f"mode (1,1) — f = {f11_computed:.1f} Hz"},
)
plotter.view_isometric()
plotter.camera.zoom(1.05)
plotter.show()