Plates (thin-plate bending)

Overview

Thin-plate (Kirchhoff–Love) bending: the deflection w of a plate of flexural rigidity D = E t³ / (12(1−ν²)) under a transverse load q solves the biharmonic equation D ∇⁴w = q. For the simply-supported (“Navier”) plate the problem decouples exactly into two Poisson problems —

\[ \Delta m = q/D \quad (m = 0 \text{ on } \partial\Omega), \qquad \Delta w = m \quad (w = 0 \text{ on } \partial\Omega) \]

— so PlateBending() reuses the scalar Laplacian already used by the heat solver. This is the FEA kernel behind plate/shell topology optimization.

using Rasmah
using LinearAlgebra

# a 16×16 triangular mesh of the unit square [0,1]² (z = 0)
function grid_tri(nx, ny)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    for j in 0:ny, i in 0:nx
        idx = i + j * (nx + 1) + 1
        verts[1, idx] = i / nx; verts[2, idx] = j / ny; verts[3, idx] = 0.0
    end
    faces = Matrix{Int}(undef, 3, 2nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

plate = grid_tri(16, 16)

# simply-supported: w = 0 on the boundary (box + filter)
edge = nodes_in_box(plate, [0.0, 0.0], [1.0, 1.0];
                    filter = p -> p[1] < 1e-9 || p[1] > 1 - 1e-9 ||
                                  p[2] < 1e-9 || p[2] > 1 - 1e-9)

w = solve(PlateBending(), plate, steel, 0.01; fixed = edge, pressure = 1e3)
println("max deflection = ", round(maximum(w); digits = 6), " m")

max deflection = 0.000221 m m

Navier-series verification

For a simply-supported square plate of side a under uniform pressure q the center deflection is w = 16 q a⁴ / (π⁶ D) · Σ_{m,n odd} 1/(mn(m²+n²)²):

E = youngs_modulus(steel); ν = poisson_ratio(steel); t = 0.01; q = 1e3
D = E * t^3 / (12 * (1 - ν^2))
s = sum(1 / (m * n * (m^2 + n^2)^2) for m in 1:2:199, n in 1:2:199)
w_navier = 16 * q / (π^6 * D) * s

ci = argmin([norm(plate.vertices[1:2, i] .- 0.5) for i in 1:size(plate.vertices, 2)])
println("FEM center   = ", round(w[ci]; digits = 6))
println("Navier center = ", round(w_navier; digits = 6))

FEM center   = 0.000221 m
Navier center = 0.000236

Notes

The field is :w (field_names(PlateBending()) == (:w,)).
pressure_load(plate, q) builds the consistent nodal source ∫ q N; a raw nodal loads vector is also accepted.
The two-Poisson splitting is exact only for the simply-supported (Navier) boundary conditions (w = 0, Δw = 0); clamped or free edges couple the two Poisson problems.
solve_plate(m, D, fixed, source) exposes the raw two-solve kernel directly.