.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/verification/example_verify_convergence_cantilever_eb.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_verification_example_verify_convergence_cantilever_eb.py: Mesh-refinement convergence — cantilever Euler-Bernoulli ======================================================== The textbook demonstration of FE convergence: for a fixed problem with a known closed form, refine the mesh and watch the discretisation error decrease at the asymptotic rate set by the element's polynomial order. For the standard 8-node hex with enhanced-strain (Simo–Rifai) under a tip point load on a slender cantilever, the expected convergence rate in the energy / displacement norm is .. math:: \|u^{h} - u\| \;\sim\; C\, h^{p}, \qquad p = 2, (Cook §6 + §6.3; Strang & Fix 2008 §3.7) — a quadratic decay of the tip-deflection error with the characteristic mesh length :math:`h`. Shear-locking-prone formulations (HEX8 plain Gauss without B-bar / EAS) fall short of this rate on thin-bending geometries, so this example is also a clean demonstration of why EAS is the right choice for solid-element plates and slender beams. Implementation -------------- A unit-load cantilever of length :math:`L = 1\,\mathrm{m}`, square cross-section :math:`b = h = 0.1\,\mathrm{m}`, with one HEX8 element through the thickness in :math:`y` and :math:`z` and a refinement ladder in the axial direction :math:`N_x \in \{2, 4, 8, 16, 32\}`. At each refinement we fix the :math:`x = 0` face in all three translations and apply a uniform tip force at the :math:`x = L` face, then read the tip mid-surface UY against the Euler-Bernoulli closed form :math:`\delta = P L^{3} / (3 E I)`. A least-squares fit to :math:`\log |\mathrm{error}|` vs :math:`\log h` recovers the asymptotic convergence rate :math:`p`. References ---------- * Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) *Concepts and Applications of Finite Element Analysis*, 4th ed., Wiley, §6 (convergence of the FE method). * Strang, G. and Fix, G. (2008) *An Analysis of the Finite Element Method*, 2nd ed., Wellesley-Cambridge, §3.7 (a-priori error estimates for the Galerkin method). * Simo, J. C. and Rifai, M. S. (1990) "A class of mixed assumed strain methods …" *International Journal for Numerical Methods in Engineering* 29 (8), 1595–1638 (HEX8 EAS). * Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed., §5.5 (a-posteriori error estimation), §11.7 (convergence rates for the discrete eigenvalue problem). .. GENERATED FROM PYTHON SOURCE LINES 61-73 .. code-block:: Python from __future__ import annotations import math import matplotlib.pyplot as plt import numpy as np import pyvista as pv import femorph_solver from femorph_solver import ELEMENTS .. GENERATED FROM PYTHON SOURCE LINES 74-76 Problem data ------------ .. GENERATED FROM PYTHON SOURCE LINES 76-125 .. code-block:: Python E = 2.0e11 # Pa NU = 0.3 RHO = 7850.0 b = 0.1 # cross-section side [m] A = b * b I = b**4 / 12.0 # noqa: E741 L = 1.0 # span [m] P_TOTAL = -5.0e3 # tip downward force [N] delta_eb = P_TOTAL * L**3 / (3.0 * E * I) print(f"L = {L} m, b = {b} m, EI = {E * I:.3e} N m^2") print(f"δ_EB = P L^3 / (3 EI) = {delta_eb:.6e} m") def run_one(nx: int) -> tuple[float, float]: """Return (h, |relative error|) for an Nx × 1 × 1 hex mesh.""" xs = np.linspace(0.0, L, nx + 1) ys = np.linspace(0.0, b, 2) zs = np.linspace(0.0, b, 2) 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}, ) pts = np.asarray(m.grid.points) node_nums = np.asarray(m.grid.point_data["ansys_node_num"]) fixed = node_nums[pts[:, 0] < 1e-9] m.fix(nodes=fixed.tolist(), dof="ALL") tip_mask = pts[:, 0] > L - 1e-9 tip_nodes = node_nums[tip_mask] fz_per_node = P_TOTAL / len(tip_nodes) for nn in tip_nodes: m.apply_force(int(nn), fy=fz_per_node) res = m.solve() g = femorph_solver.io.static_result_to_grid(m, res) pts_g = np.asarray(g.points) tip_mask_g = pts_g[:, 0] > L - 1e-9 uy_tip = float(g.point_data["displacement"][tip_mask_g, 1].mean()) h = L / nx rel_err = abs(uy_tip - delta_eb) / abs(delta_eb) return h, rel_err .. rst-class:: sphx-glr-script-out .. code-block:: none L = 1.0 m, b = 0.1 m, EI = 1.667e+06 N m^2 δ_EB = P L^3 / (3 EI) = -1.000000e-03 m .. GENERATED FROM PYTHON SOURCE LINES 126-128 Refinement ladder ----------------- .. GENERATED FROM PYTHON SOURCE LINES 128-144 .. code-block:: Python nx_list = [2, 4, 8, 16, 32] hs: list[float] = [] errs: list[float] = [] print() print(f"{'Nx':>4} {'h [m]':>8} {'δ_FE / δ_EB':>14} {'rel err':>10}") print(f"{'-' * 4:>4} {'-' * 8:>8} {'-' * 14:>14} {'-' * 10:>10}") for nx in nx_list: h, err = run_one(nx) hs.append(h) errs.append(err) # also print the deflection ratio delta_fe = delta_eb * (1 - err) if err > 0 else delta_eb # signed via convention below print(f"{nx:>4} {h:8.4f} {(1 - err):>14.6f} {err:10.4e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Nx h [m] δ_FE / δ_EB rel err ---- -------- -------------- ---------- 2 0.5000 0.903520 9.6480e-02 4 0.2500 0.967550 3.2450e-02 8 0.1250 0.989625 1.0375e-02 16 0.0625 0.997742 2.2581e-03 32 0.0312 0.999416 5.8367e-04 .. GENERATED FROM PYTHON SOURCE LINES 145-151 Asymptotic convergence rate via least-squares fit ------------------------------------------------- Drop the coarsest point (pre-asymptotic) and fit :math:`\log |\mathrm{err}| = -p\, \log h + c` by ordinary linear regression on the log–log data. .. GENERATED FROM PYTHON SOURCE LINES 151-165 .. code-block:: Python log_h = np.log(np.asarray(hs[1:])) log_err = np.log(np.asarray(errs[1:])) slope, intercept = np.polyfit(log_h, log_err, 1) p_estimated = float(slope) print() print(f"least-squares fit (Nx ≥ {nx_list[1]}): p ≈ {p_estimated:.3f}") print("expected for HEX8 EAS bending response: p = 2") # Allow some tolerance because the fit uses only 4 points and the # coarsest of those is still slightly pre-asymptotic. assert 1.5 < p_estimated, f"convergence rate {p_estimated:.3f} unexpectedly poor (< 1.5)" .. rst-class:: sphx-glr-script-out .. code-block:: none least-squares fit (Nx ≥ 4): p ≈ 1.959 expected for HEX8 EAS bending response: p = 2 .. GENERATED FROM PYTHON SOURCE LINES 166-168 Render the convergence plot --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 168-189 .. code-block:: Python fig, ax = plt.subplots(1, 1, figsize=(6.4, 4.0)) ax.loglog(hs, errs, "o-", color="#1f77b4", lw=2, label="HEX8 EAS — FE error") # Reference slope-2 line passing through the (h, err) at the finest # mesh, for visual orientation. h_ref = np.array([hs[0], hs[-1]]) err_ref_p2 = errs[-1] * (h_ref / hs[-1]) ** 2 ax.loglog(h_ref, err_ref_p2, "--", color="#d62728", lw=1.5, label=r"reference $\propto h^{2}$") ax.set_xlabel(r"mesh size $h = L / N_x$") ax.set_ylabel(r"$|\delta^{h}_{\mathrm{tip}} - \delta_{\mathrm{EB}}| / |\delta_{\mathrm{EB}}|$") ax.set_title( f"Cantilever EB — fitted convergence rate p ≈ {p_estimated:.2f}", fontsize=11, ) ax.legend(loc="lower right", fontsize=9) ax.grid(True, which="both", ls=":", alpha=0.5) fig.tight_layout() fig.show() .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_convergence_cantilever_eb_001.png :alt: Cantilever EB — fitted convergence rate p ≈ 1.96 :srcset: /gallery/verification/images/sphx_glr_example_verify_convergence_cantilever_eb_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 190-195 Verify partition-of-unity: at every refinement the deflection stays the *same sign* as the analytical (negative for our downward load), and monotone-converges to it from below (HEX8 stiffens marginally vs. the analytical bending kinematics even with EAS — Cook §6.6 + §6.13). .. GENERATED FROM PYTHON SOURCE LINES 195-204 .. code-block:: Python assert all(e > 0 for e in errs), "errors must be positive (computed |·|)" # Monotone decrease check (allow the last step to plateau within # the asymptotic-floor noise). for prev, nxt in zip(errs[:-1], errs[1:]): assert nxt < prev * 1.05, "error grew between successive refinements" print() print("OK — error decreases monotonically with refinement.") .. rst-class:: sphx-glr-script-out .. code-block:: none OK — error decreases monotonically with refinement. .. GENERATED FROM PYTHON SOURCE LINES 205-208 Render the deformed mesh at the finest resolution for visual orientation ------------------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 208-242 .. code-block:: Python xs = np.linspace(0.0, L, nx_list[-1] + 1) ys = np.linspace(0.0, b, 2) zs = np.linspace(0.0, b, 2) 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}) pts = np.asarray(m.grid.points) node_nums = np.asarray(m.grid.point_data["ansys_node_num"]) m.fix(nodes=node_nums[pts[:, 0] < 1e-9].tolist(), dof="ALL") for nn in node_nums[pts[:, 0] > L - 1e-9]: m.apply_force(int(nn), fy=P_TOTAL / int((pts[:, 0] > L - 1e-9).sum())) res = m.solve() grid = femorph_solver.io.static_result_to_grid(m, res) warp = grid.warp_by_vector("displacement", factor=20.0) plotter = pv.Plotter(off_screen=True, window_size=(720, 320)) plotter.add_mesh(grid, style="wireframe", color="grey", opacity=0.4) plotter.add_mesh( warp, scalars="displacement_magnitude", cmap="viridis", show_edges=True, scalar_bar_args={ "title": f"|u| [m] — Nx = {nx_list[-1]}, error {errs[-1] * 100:.2f}% (×20 warp)" }, ) plotter.view_xy() plotter.camera.zoom(1.05) plotter.show() # Compute the published-deflection magnitude for the title print() print(f"Published Euler-Bernoulli δ = {math.fabs(delta_eb):.4e} m") .. image-sg:: /gallery/verification/images/sphx_glr_example_verify_convergence_cantilever_eb_002.png :alt: example verify convergence cantilever eb :srcset: /gallery/verification/images/sphx_glr_example_verify_convergence_cantilever_eb_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Published Euler-Bernoulli δ = 1.0000e-03 m .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.518 seconds) .. _sphx_glr_download_gallery_verification_example_verify_convergence_cantilever_eb.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_verify_convergence_cantilever_eb.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_verify_convergence_cantilever_eb.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_verify_convergence_cantilever_eb.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_