{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Mesh-refinement convergence \u2014 cantilever Euler-Bernoulli\n\nThe textbook demonstration of FE convergence: for a fixed\nproblem with a known closed form, refine the mesh and watch\nthe discretisation error decrease at the asymptotic rate set\nby the element's polynomial order. For the standard\n8-node hex with enhanced-strain (Simo\u2013Rifai) under a tip\npoint load on a slender cantilever, the expected convergence\nrate in the energy / displacement norm is\n\n\\begin{align}\\|u^{h} - u\\|\n \\;\\sim\\; C\\, h^{p},\n \\qquad p = 2,\\end{align}\n\n(Cook \u00a76 + \u00a76.3; Strang & Fix 2008 \u00a73.7) \u2014 a quadratic decay\nof the tip-deflection error with the characteristic mesh\nlength $h$. Shear-locking-prone formulations (HEX8\nplain Gauss without B-bar / EAS) fall short of this rate on\nthin-bending geometries, so this example is also a clean\ndemonstration of why EAS is the right choice for solid-element\nplates and slender beams.\n\n## Implementation\n\nA unit-load cantilever of length $L = 1\\,\\mathrm{m}$,\nsquare cross-section $b = h = 0.1\\,\\mathrm{m}$, with one\nHEX8 element through the thickness in $y$ and $z$\nand a refinement ladder in the axial direction\n$N_x \\in \\{2, 4, 8, 16, 32\\}$. At each refinement we\nfix the $x = 0$ face in all three translations and apply\na uniform tip force at the $x = L$ face, then read the\ntip mid-surface UY against the Euler-Bernoulli closed form\n$\\delta = P L^{3} / (3 E I)$.\n\nA least-squares fit to $\\log |\\mathrm{error}|$ vs\n$\\log h$ recovers the asymptotic convergence rate\n$p$.\n\n## References\n\n* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.\n (2002) *Concepts and Applications of Finite Element\n Analysis*, 4th ed., Wiley, \u00a76 (convergence of the FE method).\n* Strang, G. and Fix, G. (2008) *An Analysis of the Finite\n Element Method*, 2nd ed., Wellesley-Cambridge, \u00a73.7 (a-priori\n error estimates for the Galerkin method).\n* Simo, J. C. and Rifai, M. S. (1990) \"A class of mixed\n assumed strain methods \u2026\" *International Journal for\n Numerical Methods in Engineering* 29 (8), 1595\u20131638\n (HEX8 EAS).\n* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,\n \u00a75.5 (a-posteriori error estimation), \u00a711.7 (convergence\n rates for the discrete eigenvalue problem).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import annotations\n\nimport math\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pyvista as pv\n\nimport femorph_solver\nfrom femorph_solver import ELEMENTS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E = 2.0e11 # Pa\nNU = 0.3\nRHO = 7850.0\nb = 0.1 # cross-section side [m]\nA = b * b\nI = b**4 / 12.0 # noqa: E741\nL = 1.0 # span [m]\nP_TOTAL = -5.0e3 # tip downward force [N]\n\ndelta_eb = P_TOTAL * L**3 / (3.0 * E * I)\nprint(f\"L = {L} m, b = {b} m, EI = {E * I:.3e} N m^2\")\nprint(f\"\u03b4_EB = P L^3 / (3 EI) = {delta_eb:.6e} m\")\n\n\ndef run_one(nx: int) -> tuple[float, float]:\n \"\"\"Return (h, |relative error|) for an Nx \u00d7 1 \u00d7 1 hex mesh.\"\"\"\n xs = np.linspace(0.0, L, nx + 1)\n ys = np.linspace(0.0, b, 2)\n zs = np.linspace(0.0, b, 2)\n grid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\n\n m = femorph_solver.Model.from_grid(grid)\n m.assign(\n ELEMENTS.HEX8(integration=\"enhanced_strain\"),\n material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO},\n )\n\n pts = np.asarray(m.grid.points)\n node_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"])\n fixed = node_nums[pts[:, 0] < 1e-9]\n m.fix(nodes=fixed.tolist(), dof=\"ALL\")\n\n tip_mask = pts[:, 0] > L - 1e-9\n tip_nodes = node_nums[tip_mask]\n fz_per_node = P_TOTAL / len(tip_nodes)\n for nn in tip_nodes:\n m.apply_force(int(nn), fy=fz_per_node)\n\n res = m.solve_static()\n g = femorph_solver.io.static_result_to_grid(m, res)\n pts_g = np.asarray(g.points)\n tip_mask_g = pts_g[:, 0] > L - 1e-9\n uy_tip = float(g.point_data[\"displacement\"][tip_mask_g, 1].mean())\n h = L / nx\n rel_err = abs(uy_tip - delta_eb) / abs(delta_eb)\n return h, rel_err" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Refinement ladder\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nx_list = [2, 4, 8, 16, 32]\nhs: list[float] = []\nerrs: list[float] = []\n\nprint()\nprint(f\"{'Nx':>4} {'h [m]':>8} {'\u03b4_FE / \u03b4_EB':>14} {'rel err':>10}\")\nprint(f\"{'-' * 4:>4} {'-' * 8:>8} {'-' * 14:>14} {'-' * 10:>10}\")\nfor nx in nx_list:\n h, err = run_one(nx)\n hs.append(h)\n errs.append(err)\n # also print the deflection ratio\n delta_fe = delta_eb * (1 - err) if err > 0 else delta_eb # signed via convention below\n print(f\"{nx:>4} {h:8.4f} {(1 - err):>14.6f} {err:10.4e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Asymptotic convergence rate via least-squares fit\n\nDrop the coarsest point (pre-asymptotic) and fit\n$\\log |\\mathrm{err}| = -p\\, \\log h + c$ by ordinary\nlinear regression on the log\u2013log data.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "log_h = np.log(np.asarray(hs[1:]))\nlog_err = np.log(np.asarray(errs[1:]))\nslope, intercept = np.polyfit(log_h, log_err, 1)\np_estimated = float(slope)\n\nprint()\nprint(f\"least-squares fit (Nx \u2265 {nx_list[1]}): p \u2248 {p_estimated:.3f}\")\nprint(\"expected for HEX8 EAS bending response: p = 2\")\n\n# Allow some tolerance because the fit uses only 4 points and the\n# coarsest of those is still slightly pre-asymptotic.\nassert 1.5 < p_estimated, f\"convergence rate {p_estimated:.3f} unexpectedly poor (< 1.5)\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Render the convergence plot\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(6.4, 4.0))\nax.loglog(hs, errs, \"o-\", color=\"#1f77b4\", lw=2, label=\"HEX8 EAS \u2014 FE error\")\n\n# Reference slope-2 line passing through the (h, err) at the finest\n# mesh, for visual orientation.\nh_ref = np.array([hs[0], hs[-1]])\nerr_ref_p2 = errs[-1] * (h_ref / hs[-1]) ** 2\nax.loglog(h_ref, err_ref_p2, \"--\", color=\"#d62728\", lw=1.5, label=r\"reference $\\propto h^{2}$\")\n\nax.set_xlabel(r\"mesh size $h = L / N_x$\")\nax.set_ylabel(r\"$|\\delta^{h}_{\\mathrm{tip}} - \\delta_{\\mathrm{EB}}| / |\\delta_{\\mathrm{EB}}|$\")\nax.set_title(\n f\"Cantilever EB \u2014 fitted convergence rate p \u2248 {p_estimated:.2f}\",\n fontsize=11,\n)\nax.legend(loc=\"lower right\", fontsize=9)\nax.grid(True, which=\"both\", ls=\":\", alpha=0.5)\nfig.tight_layout()\nfig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify partition-of-unity: at every refinement the deflection\nstays the *same sign* as the analytical (negative for our\ndownward load), and monotone-converges to it from below\n(HEX8 stiffens marginally vs. the analytical bending kinematics\neven with EAS \u2014 Cook \u00a76.6 + \u00a76.13).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "assert all(e > 0 for e in errs), \"errors must be positive (computed |\u00b7|)\"\n# Monotone decrease check (allow the last step to plateau within\n# the asymptotic-floor noise).\nfor prev, nxt in zip(errs[:-1], errs[1:]):\n assert nxt < prev * 1.05, \"error grew between successive refinements\"\nprint()\nprint(\"OK \u2014 error decreases monotonically with refinement.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Render the deformed mesh at the finest resolution for visual\norientation\n------------------------------------------------------------\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xs = np.linspace(0.0, L, nx_list[-1] + 1)\nys = np.linspace(0.0, b, 2)\nzs = np.linspace(0.0, b, 2)\ngrid = pv.StructuredGrid(*np.meshgrid(xs, ys, zs, indexing=\"ij\")).cast_to_unstructured_grid()\nm = femorph_solver.Model.from_grid(grid)\nm.assign(ELEMENTS.HEX8(integration=\"enhanced_strain\"), material={\"EX\": E, \"PRXY\": NU, \"DENS\": RHO})\npts = np.asarray(m.grid.points)\nnode_nums = np.asarray(m.grid.point_data[\"ansys_node_num\"])\nm.fix(nodes=node_nums[pts[:, 0] < 1e-9].tolist(), dof=\"ALL\")\nfor nn in node_nums[pts[:, 0] > L - 1e-9]:\n m.apply_force(int(nn), fy=P_TOTAL / int((pts[:, 0] > L - 1e-9).sum()))\nres = m.solve_static()\ngrid = femorph_solver.io.static_result_to_grid(m, res)\n\nwarp = grid.warp_by_vector(\"displacement\", factor=20.0)\nplotter = pv.Plotter(off_screen=True, window_size=(720, 320))\nplotter.add_mesh(grid, style=\"wireframe\", color=\"grey\", opacity=0.4)\nplotter.add_mesh(\n warp,\n scalars=\"displacement_magnitude\",\n cmap=\"viridis\",\n show_edges=True,\n scalar_bar_args={\n \"title\": f\"|u| [m] \u2014 Nx = {nx_list[-1]}, error {errs[-1] * 100:.2f}% (\u00d720 warp)\"\n },\n)\nplotter.view_xy()\nplotter.camera.zoom(1.05)\nplotter.show()\n\n# Compute the published-deflection magnitude for the title\nprint()\nprint(f\"Published Euler-Bernoulli \u03b4 = {math.fabs(delta_eb):.4e} m\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 0 }