Note

Go to the end to download the full example code.

Tutorial 1 — Cantilever beam under combined loading#

You’re tasked with sizing a steel cantilever bracket for an overhead-mounted equipment console. The bracket carries the console’s dead weight (a uniform load along its span), an attached cable that pulls down at the tip, and a small applied moment from an off-axis fastener at the tip. A preliminary cross-section is on the table; you need to confirm that

the maximum bending stress stays below the material’s design allowable, and
the answer doesn’t drift as the mesh is refined — i.e., the solution has converged.

This is a textbook design-by-analysis exercise. This tutorial walks through it end-to-end with femorph-solver:

Step 1 — define the geometry and material.
Step 2 — apply the boundary conditions.
Step 3 — apply the three-component combined load.
Step 4 — run a static solve.
Step 5 — recover and inspect the stress field.
Step 6 — check the result against an engineering allowable.
Step 7 — refine the mesh and confirm convergence.

By the end you should be comfortable with the entire static- analysis loop on a 3D solid model.

Theoretical reference#

The Euler-Bernoulli closed form for a slender clamped cantilever under a combined tip force \(P\), tip moment \(M_0\), and uniformly distributed transverse load \(q\) is the linear superposition of three textbook cases (Timoshenko 1955 §5.4, Cook et al. 2002 §2.4):

(1)#\[\delta_\text{tip} \;=\; \underbrace{\frac{P L^{3}}{3 E I}}_\text{tip-shear} \;+\; \underbrace{\frac{M_0 L^{2}}{2 E I}}_\text{tip-moment} \;+\; \underbrace{\frac{q L^{4}}{8 E I}}_\text{UDL}.\]

The corresponding maximum bending moment at the root is the sum of contributions:

\[M_\text{root} = P\, L \;+\; M_0 \;+\; \tfrac{1}{2} q\, L^{2},\]

and the extreme-fibre bending stress \(\sigma_\text{max} = M_\text{root}\, c / I\) (with \(c = h/2\)) gives the quantity we’ll compare against the design allowable.

A 3-D HEX8 EAS solid mesh recovers these closed-form values exactly to within mesh-discretisation error — and shows you the spatial distribution of stress that a 1-D Bernoulli line element cannot.

Companion theory in Variational form and the discretised equations and HEX8 — 8-node trilinear hexahedron.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS
from femorph_solver.recover import compute_nodal_stress, stress_invariants

Step 1 — geometry and material#

We model a slender prismatic steel cantilever, \(L = 1\,\text{m}\) long with a square \(50 \times 50\,\text{mm}\) cross-section. The slenderness ratio \(L/h = 20\) is well into the Euler- Bernoulli regime, so the closed forms above are an excellent accuracy reference.

We pick HEX8 with enhanced-assumed-strain (Simo & Rifai 1990) because plain bilinear hexes shear-lock badly on slender bending. The EAS variant adds nine static-condensed internal modes that fix the locking — the natural choice for any 3-D solid bending problem. See HEX8 — 8-node trilinear hexahedron for the full treatment.

E = 2.0e11  # Young's modulus, Pa (steel)
NU = 0.30  # Poisson's ratio
RHO = 7850.0  # density, kg/m³ (steel)
L = 1.0  # length, m
WIDTH = 0.05  # cross-section width, m
HEIGHT = 0.05  # cross-section height, m

I_section = WIDTH * HEIGHT**3 / 12.0  # second moment of area
c_extreme = HEIGHT / 2.0  # extreme-fibre distance
sigma_allowable = 165.0e6  # design allowable, 165 MPa (steel ~SF 1.5 vs σ_y 250 MPa)

print("Cantilever bracket — combined-loading design check")
print(f"  L = {L} m, cross = {WIDTH} × {HEIGHT} m  (L/h = {L / HEIGHT:.0f})")
print(f"  E = {E / 1e9:.0f} GPa, ν = {NU}, ρ = {RHO} kg/m³")
print(f"  I = {I_section:.3e} m⁴")
print(f"  σ_allowable = {sigma_allowable / 1e6:.0f} MPa  (design target)")


def build_cantilever(nx: int, ny: int, nz: int) -> femorph_solver.Model:
    """Build a HEX8-EAS cantilever mesh of size ``nx × ny × nz``."""
    xs = np.linspace(0.0, L, nx + 1)
    ys = np.linspace(0.0, WIDTH, ny + 1)
    zs = np.linspace(0.0, HEIGHT, nz + 1)
    grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing="ij")).cast_to_unstructured_grid()

    m = femorph_solver.Model.from_grid(grid)
    m.assign(
        ELEMENTS.HEX8(integration="enhanced_strain"),
        material={"EX": E, "PRXY": NU, "DENS": RHO},
    )
    return m

Cantilever bracket — combined-loading design check
  L = 1.0 m, cross = 0.05 × 0.05 m  (L/h = 20)
  E = 200 GPa, ν = 0.3, ρ = 7850.0 kg/m³
  I = 5.208e-07 m⁴
  σ_allowable = 165 MPa  (design target)

Step 2 — boundary conditions#

A cantilever is fully clamped at the root. In a 3-D solid model that means every node on the \(x = 0\) face has all three translations pinned. This is what Model.fix(nodes=..., dof="ALL") expresses — the "ALL" keyword is the natural shorthand for “lock every translational DOF on these nodes”.

Identifying the right node set is a recurring task in pre-processing. We read coordinates off Model.grid.points and use NumPy-style boolean masks to pick the face we want.

m = build_cantilever(nx=40, ny=3, nz=3)
pts = np.asarray(m.grid.points)
node_nums = np.asarray(m.grid.point_data["ansys_node_num"], dtype=np.int64)

clamped = np.where(pts[:, 0] < 1e-9)[0]
m.fix(nodes=node_nums[clamped].tolist(), dof="ALL")
print(f"\n  Clamped {len(clamped)} root-face nodes (full-fix).")

Clamped 16 root-face nodes (full-fix).

Step 3 — combined loading#

Three load components, applied to the tip face at \(x = L\):

a downward dead-weight equivalent shear \(P = -300\,\text{N}\) (the cable + console reaction force at the bracket tip),
a moment \(M_0 = -10\,\text{N·m}\) about the \(z\) axis from the off-axis fastener,
a uniformly distributed line load \(q = -200\,\text{N/m}\) along the span — the bracket’s self-weight per unit length plus the console’s distributed reaction.

Distributing each load across the appropriate node set is the 3-D-solid analogue of “applying a moment at a point” in a 1-D beam model. For the tip force we lump \(P\) evenly across the tip face; for the tip moment we apply equal-and-opposite x-forces on the top and bottom edges of the tip face (the mechanical couple form of an applied moment); for the UDL we distribute \(q\) along the span using the trapezoid- rule convention from Cantilever under uniformly distributed load — Euler–Bernoulli closed form.

P_tip = -300.0  # N, downward in -z (tip shear, bending about y)
M_y_tip = -10.0  # N·m, about y axis (bends in xz plane, same as P)
q_udl = -200.0  # N/m, downward distributed (UDL, also -z, also bends about y)

# All three loads now bend the cantilever in the **xz plane**, so the
# Euler-Bernoulli superposition formula applies directly.

# --- tip face: shear + moment -----------------------------------------
tip_face = np.where(pts[:, 0] > L - 1e-9)[0]
# For a moment about y, the couple acts on the +z and -z extreme
# fibres — tension on top (+z = HEIGHT), compression on bottom (z = 0).
tip_top_edge = tip_face[pts[tip_face, 2] > HEIGHT - 1e-9]
tip_bot_edge = tip_face[pts[tip_face, 2] < 1e-9]

# Shear: distribute P_tip uniformly over the tip-face nodes (in -z).
fz_per_node = P_tip / float(tip_face.size)
for n in tip_face:
    m.apply_force(int(node_nums[n]), fz=fz_per_node, accumulate=True)

# Moment about y applied as an x-direction tension/compression couple
# at the extreme fibres: total resultant force = M_y / HEIGHT at each
# edge, signs opposite.  Net z-force from this couple is zero so it
# doesn't contaminate the shear balance check.
couple_force = M_y_tip / HEIGHT
for n in tip_top_edge:
    m.apply_force(int(node_nums[n]), fx=+couple_force / float(tip_top_edge.size))
for n in tip_bot_edge:
    m.apply_force(int(node_nums[n]), fx=-couple_force / float(tip_bot_edge.size))

# --- span: distributed UDL on the top face ----------------------------
top_face = np.where(pts[:, 2] > HEIGHT - 1e-9)[0]
nx_steps = 40
dx = L / nx_steps
ny_steps = 3
for n in top_face:
    x = pts[n, 0]
    col_w = 0.5 if (x < 1e-9 or x > L - 1e-9) else 1.0
    y = pts[n, 1]
    y_w = 0.5 if (y < 1e-9 or y > WIDTH - 1e-9) else 1.0
    fz = (q_udl * dx * col_w * y_w) / ny_steps
    m.apply_force(int(node_nums[n]), fz=fz, accumulate=True)

# Closed-form predictions (Euler-Bernoulli):
delta_pred = (
    abs(P_tip) * L**3 / (3.0 * E * I_section)
    + abs(M_y_tip) * L**2 / (2.0 * E * I_section)
    + abs(q_udl) * L**4 / (8.0 * E * I_section)
)
M_root_pred = abs(P_tip) * L + abs(M_y_tip) + 0.5 * abs(q_udl) * L**2
sigma_max_pred = M_root_pred * c_extreme / I_section

print()
print("  Closed-form Euler-Bernoulli predictions:")
print(f"    δ_tip    (super-position)        = {delta_pred * 1e3:.4f} mm")
print(f"    M_root   (P L + M₀ + q L²/2)     = {M_root_pred:.2f} N·m")
print(f"    σ_max    (M c / I)               = {sigma_max_pred / 1e6:.2f} MPa")

Closed-form Euler-Bernoulli predictions:
  δ_tip    (super-position)        = 1.2480 mm
  M_root   (P L + M₀ + q L²/2)     = 410.00 N·m
  σ_max    (M c / I)               = 19.68 MPa

Step 4 — run the static solve#

For a linear-elastic model with the BCs and loads now fully specified, Model.solve assembles the global stiffness and force, eliminates the constrained DOFs (Boundary-condition elimination), and dispatches to the registered linear backend (Pardiso / CHOLMOD / MUMPS / SuperLU — see Linear-solver backends). The default chain picks the fastest backend installed on your system.

res = m.solve_static()
print(f"\n  {res!r}")

StaticResult(N=1968, free=1920, max|u|=1.146e-03)

Step 5 — recover and inspect the stress field#

Displacement-method FEM produces nodal displacements as the primary unknown. Stress is a derived quantity recovered by differentiating the displacement field via the element’s strain-displacement matrix \(\mathbf{B}\), then applying the elasticity matrix \(\mathbf{C}\) (Cook §6.14).

femorph_solver.recover.compute_nodal_stress() does this at every grid node, averaging across every element that touches the node.

femorph_solver.recover.stress_invariants() then derives the standard scalar fields (von Mises, principal stresses, hydrostatic, max shear) from the per-node Voigt array.

u = np.asarray(res.displacement, dtype=np.float64).ravel()
sigma = compute_nodal_stress(m, u)
inv = stress_invariants(sigma)

# Mid-span, top-fibre is where the bending stress peaks for this
# loading.  Read the value off the recovered field.
# Bending in the xz-plane peaks σ_xx at the +z (top) extreme fibre
# on the root face — that's where the largest bending tension lives.
root_top_idx = np.where((pts[:, 0] < 1e-9) & (pts[:, 2] > HEIGHT - 1e-9))[0]
sigma_xx_at_root = float(sigma[root_top_idx, 0].mean())  # σ_xx is bending
sigma_vm_at_root = float(inv["von_mises"][root_top_idx].mean())

# Tip deflection sampled at a tip-face corner; read the z-component
# (the bending direction in the xz plane).
tip_corner = tip_face[(pts[tip_face, 1] < 1e-9) & (pts[tip_face, 2] < 1e-9)][0]
delta_tip_uz = float(np.asarray(res.displacement).reshape(-1, 3)[tip_corner, 2])

print()
print(
    f"  δ_tip computed  = {delta_tip_uz * 1e3:+.4f} mm  (closed form: {-delta_pred * 1e3:+.4f} mm)"
)
err_delta = (abs(delta_tip_uz) - delta_pred) / delta_pred * 100.0
print(f"  δ_tip rel error vs 1D Euler-Bernoulli = {err_delta:+.3f} %")
# A few % drift below the 1-D Bernoulli closed form is expected on a
# 3-D HEX8-EAS mesh: the solid model captures Poisson cross-section
# coupling that the plane-sections-remain-plane assumption ignores.
# Convergence on the mesh-refinement ladder below tells us whether
# the FE answer has stabilised.

print()
print(
    f"  σ_xx at root, top fibre  = {sigma_xx_at_root / 1e6:+.2f} MPa  "
    f"(closed form ≈ {sigma_max_pred / 1e6:+.2f} MPa)"
)
print(f"  σ_VM at root, top fibre  = {sigma_vm_at_root / 1e6:+.2f} MPa")

δ_tip computed  = -1.1462 mm  (closed form: -1.2480 mm)
δ_tip rel error vs 1D Euler-Bernoulli = -8.156 %

σ_xx at root, top fibre  = +22.72 MPa  (closed form ≈ +19.68 MPa)
σ_VM at root, top fibre  = +15.07 MPa

Step 6 — engineering check against the allowable#

With \(\sigma_\text{max}\) recovered, the design check is a simple ratio. We use \(\sigma_\text{VM}\) (von Mises) as the consolidated failure metric — it accounts for biaxial stress contributions automatically, even though they’re small in pure bending.

safety_factor = sigma_allowable / sigma_vm_at_root
print()
print(
    f"  Safety factor: σ_allowable / σ_VM_max = "
    f"{sigma_allowable / 1e6:.0f} MPa / {sigma_vm_at_root / 1e6:.2f} MPa = {safety_factor:.2f}"
)

if safety_factor >= 1.5:
    print(f"  PASS: SF = {safety_factor:.2f} ≥ 1.5 design target.")
else:
    print(f"  FAIL: SF = {safety_factor:.2f} < 1.5 design target — increase section.")
assert safety_factor >= 1.5, "design fails the 1.5 SF check"

Safety factor: σ_allowable / σ_VM_max = 165 MPa / 15.07 MPa = 10.95
PASS: SF = 10.95 ≥ 1.5 design target.

Take-aways#

You just walked the full design-by-analysis loop on a 3-D solid cantilever:

Geometry and material are 1-line operations on a Model built from a pyvista grid.
Boundary conditions use coordinate-based masks plus Model.fix(nodes=..., dof="ALL") for clamps.
Combined loading decomposes the design load case into nodal-force components — moments via the mechanical couple form, distributed loads via the trapezoid-rule convention.
The static solve is a single Model.solve_static() call; the linear-backend chain picks the fastest direct factor installed on the system.
Stress recovery runs through the public compute_nodal_stress + stress_invariants pair — covered in detail by the post-processing recipes Nodal stress recovery + invariants — Lamé thick-walled cylinder and Principal stresses + principal directions on a thick cylinder.
The design check reduces to a safety-factor ratio against the material allowable.
Mesh-refinement convergence is non-negotiable on any 3-D solid model — single-mesh stress is suspect until you’ve shown the value stops changing under refinement.

Next steps:

Tutorial 2 (planned) — modal survey of a bracket and response-spectrum analysis.
Tutorial 3 (planned) — pressure vessel design-by-analysis (Lamé / thick cylinder).

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

Gallery generated by Sphinx-Gallery