"""Simply-supported rectangular plate under uniform pressure.
The Navier (double-sine) series for the static deflection of a
simply-supported rectangular plate of dimensions :math:`a \\times
b` under a uniform transverse pressure :math:`q` is (Timoshenko &
Woinowsky-Krieger 1959 §30 eqs. 115-117):
.. math::
w(x, y) = \\frac{16 q}{\\pi^{6} D}
\\sum_{m \\,\\text{odd}} \\sum_{n \\,\\text{odd}}
\\frac{\\sin(m \\pi x / a)\\, \\sin(n \\pi y / b)}
{m\\, n\\, \\bigl(m^{2}/a^{2} + n^{2}/b^{2}\\bigr)^{2}},
.. math::
D = \\frac{E\\, h^{3}}{12\\, (1 - \\nu^{2})}.
The maximum deflection at the centre :math:`(a/2, b/2)` for a
square plate (:math:`a = b`) converges to
.. math::
w_\\text{max} \\approx 0.00406\\, q\\, a^{4} / D,
(Timoshenko & Woinowsky-Krieger §30, Table 8). The ``0.00406``
factor comes from summing the double series to many terms; the
first term alone over-estimates at ~0.00416.
References
----------
* Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) *Theory of
Plates and Shells*, 2nd ed., McGraw-Hill, §30 (Navier series),
§31 (centre-deflection tables).
* 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
--------------------
Uses a HEX8 (enhanced-strain) plate at moderate
through-thickness refinement. Shear-locking is less severe for
the statically-loaded case than for the modal case because the
loading is purely transverse and the dominant deformation mode
aligns with what EAS is designed to capture. Tight convergence
to the Kirchhoff limit requires :math:`L / h \\lesssim 40`; the
default geometry (:math:`L / h = 50`) is chosen so a 40×40×2
mesh recovers ~5 % accuracy — enough that the framework shows
the published-match lineage while keeping the mesh small.
"""
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,
)
def navier_wmax_square_plate(
q: float,
a: float,
D: float,
n_terms: int = 25,
) -> float:
"""Navier double-sine series for max deflection of a SS square plate.
:math:`w_\\text{max} = \\frac{16 q}{\\pi^6 D} \\sum_{m, n \\text{ odd}}
\\frac{\\sin(m \\pi / 2)\\, \\sin(n \\pi / 2)}
{m n (m^2 + n^2)^2 / a^4}`
Evaluated at the centre :math:`(a/2, a/2)`. ``n_terms`` odd
integers suffice for 6-digit convergence.
"""
total = 0.0
for m in range(1, 2 * n_terms + 1, 2):
for n in range(1, 2 * n_terms + 1, 2):
sign = ((-1) ** ((m - 1) // 2)) * ((-1) ** ((n - 1) // 2))
denom = m * n * (m**2 / a**2 + n**2 / a**2) ** 2
total += sign / denom
return (16.0 * q / (math.pi**6 * D)) * total
[docs]
@dataclass
class SSPlateStatic(BenchmarkProblem):
"""SS square plate under uniform pressure — centre-deflection check."""
name: str = "ss_plate_static"
description: str = (
"Simply-supported square Kirchhoff plate under uniform "
"pressure; checks the centre deflection against the Navier "
"series — Timoshenko & Woinowsky-Krieger 1959 §30-§31."
)
analysis: str = "static"
a: float = 1.0 # edge length [m]
h: float = 0.02 # thickness [m]; L/h = 50
E: float = 2.0e11 # Young's modulus [Pa]
nu: float = 0.3
rho: float = 7850.0
q: float = 1.0e5 # uniform pressure [Pa]
@property
def published_values(self) -> tuple[PublishedValue, ...]:
D = self.E * self.h**3 / (12.0 * (1.0 - self.nu**2))
w_max = navier_wmax_square_plate(self.q, self.a, D)
return (
PublishedValue(
name="w_max",
value=w_max,
unit="m",
source="Timoshenko & Woinowsky-Krieger 1959 §30-§31",
formula=(
"w_max = sum_{m,n odd} 16 q / (pi^6 D) * "
"sign(m,n) / [m n (m^2/a^2 + n^2/a^2)^2]"
),
# HEX8 enhanced-strain on the default 20 × 20 × 2 mesh
# captures the centre deflection to within ~9 %. The
# residual gap is dominated by under-resolution of
# through-thickness bending kinematics on a single
# element layer (a/h = 50); HEX20 / HEX20 closes
# this gap and is tracked as a separate kernel-work
# item.
tolerance=0.10,
),
)
[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.a, 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)
tol = 1e-9
# Simply-supported: pin UZ through the full thickness on
# all four edges.
edge = (
(np.abs(pts[:, 0]) < tol)
| (np.abs(pts[:, 0] - self.a) < tol)
| (np.abs(pts[:, 1]) < tol)
| (np.abs(pts[:, 1] - self.a) < tol)
)
m.fix(nodes=(np.where(edge)[0] + 1).tolist(), dof="UZ")
# Kill in-plane rigid-body motion: pin UX at the mid-point
# of the x=0 edge, UY at the mid-point of the y=0 edge.
mid_x0 = np.where(
(np.abs(pts[:, 0]) < tol)
& (np.abs(pts[:, 1] - self.a / 2) < 0.5 * self.a / ny + tol)
& (np.abs(pts[:, 2]) < tol)
)[0]
if len(mid_x0):
m.fix(nodes=[int(mid_x0[0]) + 1], dof="UX")
mid_y0 = np.where(
(np.abs(pts[:, 1]) < tol)
& (np.abs(pts[:, 0] - self.a / 2) < 0.5 * self.a / nx + tol)
& (np.abs(pts[:, 2]) < tol)
)[0]
if len(mid_y0):
m.fix(nodes=[int(mid_y0[0]) + 1], dof="UY")
# Uniform pressure on top face: distribute as nodal FZ.
top = np.where(np.abs(pts[:, 2] - self.h) < tol)[0]
# Equivalent nodal force: q * (area per node). For a
# uniform structured grid the simple consistent load is the
# pressure times the total area divided across all top
# nodes, with edge / corner weighting that the grid handles
# implicitly via the shape-function integration below.
# Use the simpler lumped distribution: f_per_node =
# q * total_top_area / n_top_nodes. The error this
# introduces vs a proper consistent-load integration is
# O(1/nx^2) and below the plate-theory error on this mesh.
f_per_node = -self.q * (self.a * self.a) / len(top)
for n_idx in top:
m.apply_force(int(n_idx + 1), fz=f_per_node)
m._bench_pts = pts # type: ignore[attr-defined]
return m