SOLID185 — cantilever plate static analysis

Static analysis of a cantilever plate under a distributed tip load. This example walks through the full femorph_solver.Model.solve() → StaticResult → femorph_solver.io.static_result_to_grid() → pyvista rendering loop and checks static equilibrium via StaticResult.reaction.

Euler–Bernoulli beam theory is included as a back-of-envelope reference; SOLID185 is a first-order hex with full 2 × 2 × 2 Gauss integration, which exhibits well-known shear locking in thin-bending problems unless many elements are used through the thickness. The difference between the two is a feature of the element, not a bug — swap in SOLID186 (quadratic) to recover EB to a few percent with the same mesh.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Geometry + material

Steel cantilever, 1 m × 0.1 m × 0.05 m, meshed 40 × 4 × 4 hex (640 SOLID185 elements, 1 025 nodes).

E = 2.0e11  # Pa
NU = 0.30
RHO = 7850.0
LX, LY, LZ = 1.0, 0.1, 0.05
NX, NY, NZ = 40, 4, 4
F_TIP = -5.0e3  # N (downward)

xs = np.linspace(0.0, LX, NX + 1)
ys = np.linspace(0.0, LY, NY + 1)
zs = np.linspace(0.0, LZ, NZ + 1)
xx, yy, zz = np.meshgrid(xs, ys, zs, indexing="ij")
points = np.stack([xx.ravel(), yy.ravel(), zz.ravel()], axis=1)


def _node_id(i: int, j: int, k: int) -> int:
    return (i * (NY + 1) + j) * (NZ + 1) + k + 1


cells = []
for i in range(NX):
    for j in range(NY):
        for k in range(NZ):
            cells.append(
                [
                    _node_id(i, j, k),
                    _node_id(i + 1, j, k),
                    _node_id(i + 1, j + 1, k),
                    _node_id(i, j + 1, k),
                    _node_id(i, j, k + 1),
                    _node_id(i + 1, j, k + 1),
                    _node_id(i + 1, j + 1, k + 1),
                    _node_id(i, j + 1, k + 1),
                ]
            )

Build the model

m = femorph_solver.Model()
m.et(1, "SOLID185")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
m.mp("DENS", 1, RHO)
for nn, (x, y, z) in enumerate(points, start=1):
    m.n(nn, x, y, z)
for conn in cells:
    m.e(*conn)

# Clamp the ``x = 0`` face in all 3 DOFs.
x0_nodes = [nn for nn, p in enumerate(points, start=1) if p[0] < 1e-9]
for nn in x0_nodes:
    m.d(nn, "UX")
    m.d(nn, "UY")
    m.d(nn, "UZ")

# Distributed downward tip load.
tip_nodes = [nn for nn, p in enumerate(points, start=1) if p[0] > LX - 1e-9]
for nn in tip_nodes:
    m.f(nn, "FZ", F_TIP / len(tip_nodes))

/home/runner/_work/solver/solver/examples/analyses/static/example_cantilever_plate.py:73: 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, "SOLID185")
/home/runner/_work/solver/solver/examples/analyses/static/example_cantilever_plate.py:74: 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/analyses/static/example_cantilever_plate.py:75: 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/analyses/static/example_cantilever_plate.py:76: 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/analyses/static/example_cantilever_plate.py:78: 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(nn, x, y, z)
/home/runner/_work/solver/solver/examples/analyses/static/example_cantilever_plate.py:80: 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(*conn)
/home/runner/_work/solver/solver/examples/analyses/static/example_cantilever_plate.py:85: 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/analyses/static/example_cantilever_plate.py:86: 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/analyses/static/example_cantilever_plate.py:87: 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/analyses/static/example_cantilever_plate.py:92: 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, "FZ", F_TIP / len(tip_nodes))

Solve + reaction check

Model.solve() returns a StaticResult with displacement, reaction, and free_mask. Reactions are nonzero only at constrained DOFs; summing FZ at the clamp must equal -F_TIP to machine precision for a well-posed static solve.

res = m.solve()

dof = m.dof_map()
fz_clamp = 0.0
for nn in x0_nodes:
    rows = np.where((dof[:, 0] == nn) & (dof[:, 1] == 2))[0]
    for r in rows:
        fz_clamp += float(res.reaction[r])
print(f"Σ FZ reaction at clamp = {fz_clamp:.4e} N   (expected {-F_TIP:.4e})")

Σ FZ reaction at clamp = 5.0000e+03 N   (expected 5.0000e+03)

Tip deflection vs Euler–Bernoulli

\(\\delta_\\mathrm{EB} = F L^3 / (3 E I)\) with \(I = b h^3 / 12\) is the slender-beam estimate. With 4 elements through the thickness, SOLID185’s shear locking gives a few percent error — swap in SOLID186 (see SOLID186 — uniaxial tension on a 20-node hex) to remove it entirely.

I_y = LY * LZ**3 / 12.0
delta_eb = F_TIP * LX**3 / (3.0 * E * I_y)

grid = femorph_solver.io.static_result_to_grid(m, res)
tip_mask = grid.points[:, 0] > LX - 1e-9
w_tip_femorph_solver = grid.point_data["displacement"][tip_mask, 2].min()
print(f"Euler-Bernoulli tip deflection = {delta_eb:.3e} m")
print(f"femorph-solver tip deflection (min UZ)   = {w_tip_femorph_solver:.3e} m")
print(
    f"relative error                  = {abs(w_tip_femorph_solver - delta_eb) / abs(delta_eb):.2%}"
)

Euler-Bernoulli tip deflection = -8.000e-03 m
femorph-solver tip deflection (min UZ)   = -7.348e-03 m
relative error                  = 8.15%

Render the deformed plate, coloured by displacement magnitude

warped = grid.warp_by_vector("displacement", factor=20.0)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    m.grid,
    style="wireframe",
    color="gray",
    opacity=0.35,
    label="undeformed",
)
plotter.add_mesh(
    warped,
    scalars="displacement_magnitude",
    show_edges=True,
    cmap="viridis",
    scalar_bar_args={"title": "|u| [m]"},
    label="deformed ×20",
)
plotter.add_legend()
plotter.add_axes()
plotter.camera_position = [(2.4, -1.6, 1.0), (0.5, 0.05, 0.0), (0.0, 0.0, 1.0)]
plotter.show()