PLANE182 — plane-stress uniaxial tension

A square piece of plate is fixed on its left edge and pulled on its right edge. We check the tip displacement against Hooke’s law u_x = σ L / E and the Poisson contraction u_y = −ν u_x in plane stress, then render the deformed mesh.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem set-up

Unit square of steel, 1 mm thick, pulled with 1 MN total along +x.

E = 2.0e11  # Pa
NU = 0.3
RHO = 7850.0  # kg/m³
THK = 1.0e-3  # m
F_TOTAL = 1.0e6  # N

m = femorph_solver.Model()
m.et(1, "PLANE182")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
# Plane-stress KEYOPT(3)=0 is the default; annotate it explicitly:
m.materials[1]["_PLANE_MODE"] = "stress"
m.r(1, THK)

for i, (x, y) in enumerate([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)], start=1):
    m.n(i, x, y, 0.0)
m.e(1, 2, 3, 4)

# Clamp the x=0 edge in UX, and kill rigid-body UY without locking Poisson.
for nn in (1, 4):
    m.d(nn, "UX")
for nn in (1, 2):
    m.d(nn, "UY")

F_each = F_TOTAL / 2.0
for nn in (2, 3):
    m.f(nn, "FX", F_each)

/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:32: DeprecationWarning: Model.et(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).et(et_id, name)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material)` for new code.
  m.et(1, "PLANE182")
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:33: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("EX", 1, E)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:34: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("PRXY", 1, NU)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:35: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("DENS", 1, RHO)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:38: DeprecationWarning: Model.r(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).r(real_id, *values)` for line-by-line APDL deck porting, or the native `Model.assign(element, material, real=[...])` for new code.
  m.r(1, THK)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:41: DeprecationWarning: Model.n(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).n(num, x, y, z)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.n(i, x, y, 0.0)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:42: DeprecationWarning: Model.e(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).e(*node_nums)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.e(1, 2, 3, 4)
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:46: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(nn, "UX")
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:48: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(nn, "UY")
/home/runner/_work/solver/solver/examples/elements/plane182/example_plane182.py:52: DeprecationWarning: Model.f(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).f(node, label, value)` for line-by-line APDL deck porting, or the native `Model.apply_force(node, fx=..., fy=..., fz=...)` for new code.
  m.f(nn, "FX", F_each)

Static solve + analytical comparison

Out-of-plane stress is zero (plane stress), so the effective axial stress on the x=1 edge is σ = F_TOTAL / (A · thickness) with edge length 1 m, giving strain σ / E and tip displacement σ · L / E.

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

sigma = F_TOTAL / (1.0 * THK)
ux_expected = sigma / E  # L = 1 m
uy_expected = -NU * ux_expected

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

print(f"Expected u_y (+y edge) = {uy_expected:.6e} m")
for nn in (3, 4):
    row = np.where((dof[:, 0] == nn) & (dof[:, 1] == 1))[0][0]
    print(f"  node {nn} UY         = {res.displacement[row]:.6e}")
    assert np.isclose(res.displacement[row], uy_expected, rtol=1e-8)

Expected u_x (+x edge) = 5.000000e-03 m
  node 2 UX         = 5.000000e-03
  node 3 UX         = 5.000000e-03
Expected u_y (+y edge) = -1.500000e-03 m
  node 3 UY         = -1.500000e-03
  node 4 UY         = -1.500000e-03

Plot the deformed quad, coloured by displacement magnitude

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:
        disp[i, int(dof[r, 1])] = res.displacement[r]
grid.point_data["displacement"] = disp

plotter = pv.Plotter(off_screen=True)
warped = grid.warp_by_vector("displacement", factor=50.0)
plotter.add_mesh(
    warped,
    scalars=np.linalg.norm(disp, axis=1),
    scalar_bar_args={"title": "disp [m]"},
    show_edges=True,
    cmap="viridis",
)
plotter.add_mesh(grid, style="wireframe", color="gray", opacity=0.4)
plotter.view_xy()
plotter.show()