.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tutorials/tutorial_01_cantilever_combined.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_tutorials_tutorial_01_cantilever_combined.py: 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 1. the maximum bending stress stays below the material's design allowable, and 2. 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 :math:`P`, tip moment :math:`M_0`, and uniformly distributed transverse load :math:`q` is the linear superposition of three textbook cases (Timoshenko 1955 §5.4, Cook et al. 2002 §2.4): .. math:: :label: combined-tip-deflection \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: .. math:: M_\text{root} = P\, L \;+\; M_0 \;+\; \tfrac{1}{2} q\, L^{2}, and the extreme-fibre bending stress :math:`\sigma_\text{max} = M_\text{root}\, c / I` (with :math:`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 :doc:`/reference/theory/variational_form` and :doc:`/reference/elements/hex8`. .. GENERATED FROM PYTHON SOURCE LINES 71-81 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 82-97 Step 1 — geometry and material ------------------------------ We model a slender prismatic steel cantilever, :math:`L = 1\,\text{m}` long with a square :math:`50 \times 50\,\text{mm}` cross-section. The slenderness ratio :math:`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 :doc:`/reference/elements/hex8` for the full treatment. .. GENERATED FROM PYTHON SOURCE LINES 97-131 .. code-block:: Python 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 .. rst-class:: sphx-glr-script-out .. code-block:: none 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) .. GENERATED FROM PYTHON SOURCE LINES 132-145 Step 2 — boundary conditions ---------------------------- A cantilever is *fully clamped* at the root. In a 3-D solid model that means every node on the :math:`x = 0` face has all three translations pinned. This is what :meth:`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. .. GENERATED FROM PYTHON SOURCE LINES 145-154 .. code-block:: Python 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).") .. rst-class:: sphx-glr-script-out .. code-block:: none Clamped 16 root-face nodes (full-fix). .. GENERATED FROM PYTHON SOURCE LINES 155-177 Step 3 — combined loading ------------------------- Three load components, applied to the tip face at :math:`x = L`: * a downward dead-weight equivalent shear :math:`P = -300\,\text{N}` (the cable + console reaction force at the bracket tip), * a moment :math:`M_0 = -10\,\text{N·m}` about the :math:`z` axis from the off-axis fastener, * a uniformly distributed line load :math:`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 :math:`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 :math:`q` along the span using the trapezoid- rule convention from :ref:`sphx_glr_gallery_verification_example_verify_cantilever_udl.py`. .. GENERATED FROM PYTHON SOURCE LINES 177-235 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 236-247 Step 4 — run the static solve ----------------------------- For a linear-elastic model with the BCs and loads now fully specified, :meth:`Model.solve ` assembles the global stiffness and force, eliminates the constrained DOFs (:doc:`/reference/theory/bc_elimination`), and dispatches to the registered linear backend (Pardiso / CHOLMOD / MUMPS / SuperLU — see :doc:`/reference/solvers/linear_backends`). The default chain picks the fastest backend installed on your system. .. GENERATED FROM PYTHON SOURCE LINES 247-251 .. code-block:: Python res = m.solve_static() print(f"\n {res!r}") .. rst-class:: sphx-glr-script-out .. code-block:: none StaticResult(N=1968, free=1920, max|u|=1.146e-03) .. GENERATED FROM PYTHON SOURCE LINES 252-268 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 :math:`\mathbf{B}`, then applying the elasticity matrix :math:`\mathbf{C}` (Cook §6.14). :func:`femorph_solver.recover.compute_nodal_stress` does this at every grid node, averaging across every element that touches the node. :func:`femorph_solver.recover.stress_invariants` then derives the standard scalar fields (von Mises, principal stresses, hydrostatic, max shear) from the per-node Voigt array. .. GENERATED FROM PYTHON SOURCE LINES 268-305 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none δ_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 .. GENERATED FROM PYTHON SOURCE LINES 306-314 Step 6 — engineering check against the allowable ------------------------------------------------ With :math:`\sigma_\text{max}` recovered, the design check is a simple ratio. We use :math:`\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. .. GENERATED FROM PYTHON SOURCE LINES 314-328 .. code-block:: Python 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" .. rst-class:: sphx-glr-script-out .. code-block:: none Safety factor: σ_allowable / σ_VM_max = 165 MPa / 15.07 MPa = 10.95 PASS: SF = 10.95 ≥ 1.5 design target. .. GENERATED FROM PYTHON SOURCE LINES 329-341 Step 7 — mesh-refinement convergence study ------------------------------------------ A single-mesh result is never enough. Even on a problem with a clean closed form, the recovered stress depends on the mesh — the 3-D solid captures Poisson contributions that the 1-D closed form ignores, and the stress recovery itself smooths across element boundaries. A converged answer is one whose value stops changing as we refine. We sweep three meshes (coarse / medium / fine) and tabulate the recovered :math:`\sigma_\text{xx}` at the root-top fibre. .. GENERATED FROM PYTHON SOURCE LINES 341-406 .. code-block:: Python print() print(f" {'mesh':>16} {'σ_xx at root [MPa]':>20} {'δ tip [mm]':>11} {'safety factor':>14}") print(" " + "-" * 70) def recover_at_root(model: femorph_solver.Model, result): pts_local = np.asarray(model.grid.points) u_local = np.asarray(result.displacement, dtype=np.float64).ravel() sigma_local = compute_nodal_stress(model, u_local) inv_local = stress_invariants(sigma_local) root_top_local = np.where((pts_local[:, 0] < 1e-9) & (pts_local[:, 2] > HEIGHT - 1e-9))[0] sxx = float(sigma_local[root_top_local, 0].mean()) vm = float(inv_local["von_mises"][root_top_local].mean()) tip_face_local = np.where(pts_local[:, 0] > L - 1e-9)[0] tip_corner_local = tip_face_local[ (pts_local[tip_face_local, 1] < 1e-9) & (pts_local[tip_face_local, 2] < 1e-9) ][0] delta_local = float(np.asarray(result.displacement).reshape(-1, 3)[tip_corner_local, 2]) return sxx, vm, delta_local for tag, (nx_r, ny_r, nz_r) in ( ("coarse 20×3×3", (20, 3, 3)), ("medium 40×3×3", (40, 3, 3)), ("fine 80×3×3", (80, 3, 3)), ): m_r = build_cantilever(nx=nx_r, ny=ny_r, nz=nz_r) pts_r = np.asarray(m_r.grid.points) nn_r = np.asarray(m_r.grid.point_data["ansys_node_num"], dtype=np.int64) clamped_r = np.where(pts_r[:, 0] < 1e-9)[0] m_r.fix(nodes=nn_r[clamped_r].tolist(), dof="ALL") tip_face_r = np.where(pts_r[:, 0] > L - 1e-9)[0] tip_top_r = tip_face_r[pts_r[tip_face_r, 2] > HEIGHT - 1e-9] tip_bot_r = tip_face_r[pts_r[tip_face_r, 2] < 1e-9] fz_r = P_tip / float(tip_face_r.size) for n in tip_face_r: m_r.apply_force(int(nn_r[n]), fz=fz_r, accumulate=True) couple_r = M_y_tip / HEIGHT for n in tip_top_r: m_r.apply_force(int(nn_r[n]), fx=+couple_r / float(tip_top_r.size)) for n in tip_bot_r: m_r.apply_force(int(nn_r[n]), fx=-couple_r / float(tip_bot_r.size)) top_r = np.where(pts_r[:, 2] > HEIGHT - 1e-9)[0] dx_r = L / nx_r for n in top_r: x = pts_r[n, 0] col_w = 0.5 if (x < 1e-9 or x > L - 1e-9) else 1.0 y = pts_r[n, 1] y_w = 0.5 if (y < 1e-9 or y > WIDTH - 1e-9) else 1.0 fz = (q_udl * dx_r * col_w * y_w) / ny_r m_r.apply_force(int(nn_r[n]), fz=fz, accumulate=True) res_r = m_r.solve_static() sxx_r, vm_r, delta_r = recover_at_root(m_r, res_r) sf_r = sigma_allowable / vm_r print(f" {tag:>16} {sxx_r / 1e6:>+18.3f} {delta_r * 1e3:>+9.3f} {sf_r:>13.2f}") print() print( f" Closed form (EB): σ_xx ≈ {-sigma_max_pred / 1e6:+.2f} MPa, " f"δ ≈ {-delta_pred * 1e3:+.3f} mm" ) .. rst-class:: sphx-glr-script-out .. code-block:: none mesh σ_xx at root [MPa] δ tip [mm] safety factor ---------------------------------------------------------------------- coarse 20×3×3 +22.271 -1.140 9.44 medium 40×3×3 +22.716 -1.146 10.95 fine 80×3×3 +22.903 -1.149 11.08 Closed form (EB): σ_xx ≈ -19.68 MPa, δ ≈ -1.248 mm .. GENERATED FROM PYTHON SOURCE LINES 407-441 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 :ref:`sphx_glr_gallery_post-processing_example_nodal_stress_recovery.py` and :ref:`sphx_glr_gallery_post-processing_example_principal_stress.py`. - **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). .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.063 seconds) .. _sphx_glr_download_gallery_tutorials_tutorial_01_cantilever_combined.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tutorial_01_cantilever_combined.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tutorial_01_cantilever_combined.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tutorial_01_cantilever_combined.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_