Clamped square plate under uniform pressure (NAFEMS LE5)

Classical fully-clamped thin-plate verification: a square plate of side \(a\) and thickness \(h\) is clamped on all four edges (every DOF pinned, full thickness) and loaded by a uniform transverse pressure \(q\). The maximum centre deflection has no elementary closed form but is published from the converged double-Fourier solution (NAFEMS R0015 §2.5 LE5; Timoshenko & Woinowsky-Krieger 1959 §32 Table 12) as

\[w_\mathrm{max} \;\approx\; 0.00126 \, \frac{q\, a^{4}}{D}, \qquad D = \frac{E\, h^{3}}{12 (1 - \nu^{2})}.\]

The 0.00126 coefficient (the NAFEMS LE5 reference value) is a factor of ~3.2 stiffer than the simply-supported counterpart (0.00406 in Simply-supported plate under uniform pressure (Navier series)) — clamping roughly triples the bending stiffness of the slab.

Implementation

Drives the existing ClampedPlateStatic problem class on a default 20×20×2 HEX8 enhanced-strain mesh. The validation suite tolerates ~10–15 % deviation on this mesh (NAFEMS LE5 specifies a 5 % engineering tolerance for the benchmark itself, and a solid-element analogue at \(L/h = 50\) lands within ~10 % of the Kirchhoff limit).

References

• Timoshenko, S. P. and Woinowsky-Krieger, S. (1959) Theory of Plates and Shells, 2nd ed., McGraw-Hill, §32 Table 12 (clamped-plate deflection coefficients).

• NAFEMS, Standard Benchmark Tests for Linear Elastic Analysis, R0015 (1990), §2.5 LE5 “Clamped square plate”.

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

• (Industry-neutral reference — NAFEMS R0015 §2.5 LE5); plate-bending companion sweeps in Abaqus AVM 1.4.1 clamped-plate-pressure family (problem-id citations only).

from __future__ import annotations

import numpy as np
import pyvista as pv

from femorph_solver.validation.problems import ClampedPlateStatic

Build the model from the validation problem class

problem = ClampedPlateStatic()
m = problem.build_model()
print(
    f"clamped plate mesh: {m.grid.n_points} nodes, {m.grid.n_cells} HEX8 cells; "
    f"a = {problem.a} m, h = {problem.h} m, L/h = {problem.a / problem.h:.0f}"
)
print(f"E = {problem.E / 1e9:.0f} GPa, ν = {problem.nu}, q = {problem.q / 1e3:.1f} kPa")

clamped plate mesh: 11163 nodes, 7200 HEX8 cells; a = 1.0 m, h = 0.02 m, L/h = 50
E = 200 GPa, ν = 0.3, q = 100.0 kPa

Closed-form reference + flexural rigidity

D = problem.E * problem.h**3 / (12.0 * (1.0 - problem.nu**2))
w_published = 0.00126 * problem.q * problem.a**4 / D
print(f"D            = {D:.3e} N m  (flexural rigidity)")
print(f"w_max  (NAFEMS LE5 / Timoshenko Table 12 ≈ 0.00126 q a^4/D) = {w_published * 1e6:.3f} µm")

D            = 1.465e+05 N m  (flexural rigidity)
w_max  (NAFEMS LE5 / Timoshenko Table 12 ≈ 0.00126 q a^4/D) = 859.950 µm

Static solve + centre-deflection extraction

res = m.solve()
w_computed = problem.extract(m, res, "w_max")
err_pct = (w_computed - w_published) / w_published * 100.0
print()
print(f"w_max computed (HEX8 EAS, default mesh) = {w_computed * 1e6:+.3f} µm")
print(f"w_max published (NAFEMS LE5)            = {w_published * 1e6:+.3f} µm")
print(f"relative error                          = {err_pct:+.2f} %")

# 5 % tolerance — the 60 × 60 × 2 HEX8-EAS default mesh in
# ``ClampedPlateStatic.build_model`` lands within ~3 % of the
# Kirchhoff closed form.  Tighter convergence (down to ~1 %) needs a
# dedicated SHELL kernel (tracked separately).
assert abs(err_pct) < 5.0, f"w_max deviation {err_pct:.2f}% exceeds 5 % tolerance"

w_max computed (HEX8 EAS, default mesh) = +836.751 µm
w_max published (NAFEMS LE5)            = +859.950 µm
relative error                          = -2.70 %

Render the deflected plate

grid = m.grid.copy()
u = np.asarray(res.displacement, dtype=np.float64).reshape(-1, 3)
grid.point_data["displacement"] = u
grid.point_data["UZ"] = u[:, 2]

warp = grid.warp_by_vector("displacement", factor=5.0e2)

plotter = pv.Plotter(off_screen=True, window_size=(720, 480))
plotter.add_mesh(
    grid,
    style="wireframe",
    color="grey",
    opacity=0.35,
    line_width=1,
    label="undeformed",
)
plotter.add_mesh(
    warp,
    scalars="UZ",
    cmap="viridis",
    show_edges=False,
    scalar_bar_args={"title": "UZ [m]  (deformed ×500)"},
)
plotter.view_isometric()
plotter.camera.zoom(1.05)
plotter.show()