QUAD4_SHELL — uniaxial stretch of a flat square plate

Single QUAD4_SHELL element (a 1 m × 1 m flat plate of thickness t) in pure in-plane tension. Out-of-plane DOFs are suppressed so the problem reduces to plane-stress membrane behaviour: u_x = σ L / E and the Poisson contraction u_y = −ν u_x.

from __future__ import annotations

import numpy as np
import pyvista as pv
from vtkmodules.util.vtkConstants import VTK_QUAD

import femorph_solver
from femorph_solver import ELEMENTS

Problem data

Steel plate, 10 mm thick, pulled with a total 100 kN along +x.

E = 2.1e11
NU = 0.3
RHO = 7850.0
THK = 0.010  # 10 mm
L = 1.0
F_TOTAL = 1.0e5  # N spread over the +x edge (two corner nodes)

Build the model

Corner ordering is VTK_QUAD: node 1 = (0,0), 2 = (1,0), 3 = (1,1), 4 = (0,1). We clamp UX on the -x edge, anchor UY on a single corner to kill the in-plane translation mode, and suppress every out-of-plane DOF (UZ, ROTX, ROTY, ROTZ) so the element behaves like a pure membrane.

coords = np.array(
    [
        [0.0, 0.0, 0.0],
        [L, 0.0, 0.0],
        [L, L, 0.0],
        [0.0, L, 0.0],
    ],
    dtype=np.float64,
)

cells = np.array([4, 0, 1, 2, 3], dtype=np.int64)
cell_types = np.array([VTK_QUAD], dtype=np.uint8)
grid = pv.UnstructuredGrid(cells, cell_types, coords)

m = femorph_solver.Model.from_grid(grid)
m.assign(
    ELEMENTS.QUAD4_SHELL,
    material={"EX": E, "PRXY": NU, "DENS": RHO},
    real=(THK,),
)

# -x edge clamped in UX
m.fix(nodes=[1], dof="UX")
m.fix(nodes=[4], dof="UX")
# Anchor UY on one node to stop rigid-body in-plane translation
m.fix(nodes=[1], dof="UY")
# Pure membrane: zero all out-of-plane DOFs everywhere
for nn in (1, 2, 3, 4):
    m.fix(nodes=[nn], dof="UZ")
    m.fix(nodes=[nn], dof="ROTX")
    m.fix(nodes=[nn], dof="ROTY")
    m.fix(nodes=[nn], dof="ROTZ")

F_each = F_TOTAL / 2.0
m.apply_force(2, fx=F_each)
m.apply_force(3, fx=F_each)

Static solve + analytical comparison

Stress on the +x face is σ = F / (width · t). On a 1 m width with 10 mm thickness: σ = 10 MPa. Expected u_x = σ L / E and u_y = −ν u_x on the free corner nodes.

res = m.solve()
dof = m.dof_map()

sigma = F_TOTAL / (L * THK)
ux_expected = sigma * L / E
uy_expected = -NU * ux_expected

print(f"In-plane stress σ    = {sigma:.3e} Pa")
print(f"Expected u_x         = {ux_expected:.6e} m")
for nn in (2, 3):
    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 0))[0][0]
    val = res.displacement[row]
    print(f"  node {nn} UX        = {val:.6e}")
    assert np.isclose(val, ux_expected, rtol=1e-8)

# Node 3 is free in UY (node 1 is anchored, so check there)
row_uy3 = np.where((dof[:, 0] == 3) & (dof[:, 1] == 1))[0][0]
uy3 = res.displacement[row_uy3]
print(f"Expected u_y (node 3) = {uy_expected:.6e}")
print(f"Computed u_y (node 3) = {uy3:.6e}")
assert np.isclose(uy3, uy_expected, rtol=1e-6)

In-plane stress σ    = 1.000e+07 Pa
Expected u_x         = 4.761905e-05 m
  node 2 UX        = 4.761905e-05
  node 3 UX        = 4.761905e-05
Expected u_y (node 3) = -1.428571e-05
Computed u_y (node 3) = -1.428571e-05

Plot the deformed plate

Scatter displacements onto point_data and warp the mesh. A modest magnification (~1000×) makes the Poisson contraction visible.

grid = m.grid.copy()
disp = np.zeros((grid.n_points, 3), dtype=np.float64)
for i, nn in enumerate(grid.point_data["ansys_node_num"]):
    rows = np.where(dof[:, 0] == int(nn))[0]
    for r in rows:
        d_idx = int(dof[r, 1])
        if d_idx < 3:
            disp[i, d_idx] = res.displacement[r]
grid.point_data["displacement"] = disp

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=2)
plotter.add_mesh(
    grid.warp_by_vector("displacement", factor=2.0e3),
    scalars=np.linalg.norm(disp, axis=1),
    show_edges=True,
    scalar_bar_args={"title": "|u| [m]"},
)
plotter.add_axes()
plotter.show()