SOLID186 — uniaxial tension on a 20-node hex

Single SOLID186 brick (a unit cube, 20 nodes including the mid-edge nodes of the serendipity family) loaded in uniaxial tension. The consistent-load vector for a serendipity 8-node face has corner weight \(-1/12\) and mid-edge weight \(+1/3\); applying these produces a uniform σxx field, so εxx = σ/E and εyy = εzz = −ν·εxx exactly at every node.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Reference 20-node unit cube

Corners 1-8 in MAPDL/VTK hex order, then mid-edge nodes 9-20 on the bottom, top, and vertical edges.

E = 2.1e11  # Pa
NU = 0.30
F_TOTAL = 4.2e4  # N

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

Build the model

m = femorph_solver.Model()
m.et(1, "SOLID186")
m.mp("EX", 1, E)
m.mp("PRXY", 1, NU)
for i, (x, y, z) in enumerate(coords, start=1):
    m.n(i, x, y, z)
m.e(*range(1, 21))

# Symmetry BC: x=0 face UX, y=0 face UY, z=0 face UZ.
x0 = [1, 4, 5, 8, 12, 16, 17, 20]
y0 = [1, 2, 5, 6, 9, 13, 17, 18]
z0 = [1, 2, 3, 4, 9, 10, 11, 12]
for n in x0:
    m.d(n, "UX", 0.0)
for n in y0:
    m.d(n, "UY", 0.0)
for n in z0:
    m.d(n, "UZ", 0.0)

# Consistent-load on x=1 face: 4 corners × −F/12 + 4 mid-edges × +F/3
# (integrates to F_TOTAL exactly for a uniform traction).
for n in (2, 3, 6, 7):
    m.f(n, "FX", -F_TOTAL / 12.0)
for n in (10, 14, 18, 19):
    m.f(n, "FX", F_TOTAL / 3.0)

/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:61: 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, "SOLID186")
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:62: 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/solid186/example_solid186.py:63: 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/solid186/example_solid186.py:65: 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/solid186/example_solid186.py:66: 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(*range(1, 21))
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:73: 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(n, "UX", 0.0)
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:75: 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(n, "UY", 0.0)
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:77: 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(n, "UZ", 0.0)
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:82: 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(n, "FX", -F_TOTAL / 12.0)
/home/runner/_work/solver/solver/examples/elements/solid186/example_solid186.py:84: 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(n, "FX", F_TOTAL / 3.0)

Static solve and strain recovery

res = m.solve()
eps = m.eel(res.displacement)

eps_xx_expected = F_TOTAL / E
eps_yy_expected = -NU * eps_xx_expected
print(f"εxx expected = {eps_xx_expected:.3e}  |  recovered (mean) = {eps[:, 0].mean():.3e}")
print(f"εyy expected = {eps_yy_expected:.3e}  |  recovered (mean) = {eps[:, 1].mean():.3e}")

εxx expected = 2.000e-07  |  recovered (mean) = 2.000e-07
εyy expected = -6.000e-08  |  recovered (mean) = -6.000e-08

Plot the deformed cube, coloured by εxx

grid = femorph_solver.io.static_result_to_grid(m, res)
node_nums = m.node_numbers
node_to_idx = {int(nn): i for i, nn in enumerate(node_nums)}
grid.point_data["eps_xx"] = np.array(
    [eps[node_to_idx[int(nn)], 0] for nn in grid.point_data["ansys_node_num"]]
)

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=1, opacity=0.4)
plotter.add_mesh(
    grid.warp_by_vector("displacement", factor=2.0e5),
    scalars="eps_xx",
    show_edges=True,
    cmap="viridis",
    scalar_bar_args={"title": "εxx"},
)
plotter.add_axes()
plotter.show()