Note

Go to the end to download the full example code.

SHELL181 — uniaxial stretch of a flat square plate#

Single SHELL181 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

import femorph_solver

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,
)

m = femorph_solver.Model()
m.et(1, "SHELL181")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
m.r(1, THK)
for i, (x, y, z) in enumerate(coords, start=1):
    m.n(i, x, y, z)
m.e(1, 2, 3, 4)

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

F_each = F_TOTAL / 2.0
m.f(2, "FX", F_each)
m.f(3, "FX", F_each)

/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:50: 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, "SHELL181")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:51: 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/shell181/example_shell181.py:52: 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/shell181/example_shell181.py:53: 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/shell181/example_shell181.py:54: 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/shell181/example_shell181.py:56: 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, z)
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:57: 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/shell181/example_shell181.py:60: 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(1, "UX")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:61: 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(4, "UX")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:63: 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(1, "UY")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:66: 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, "UZ")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:67: 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, "ROTX")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:68: 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, "ROTY")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:69: 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, "ROTZ")
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:72: 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(2, "FX", F_each)
/home/runner/_work/solver/solver/examples/elements/shell181/example_shell181.py:73: 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(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()
example shell181

Total running time of the script: (0 minutes 0.271 seconds)

Gallery generated by Sphinx-Gallery