.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/Beam/Beam3.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_examples_Beam_Beam3.py: Beam3 ===== A bi-fixed beam undergoing bending deformation. .. GENERATED FROM PYTHON SOURCE LINES 12-181 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_001.png :alt: SEG3: Ne = 10, Nn = 21 :srcset: /examples/Beam/images/sphx_glr_Beam3_001.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_002.png :alt: Section :srcset: /examples/Beam/images/sphx_glr_Beam3_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_003.png :alt: Boundary conditions :srcset: /examples/Beam/images/sphx_glr_Beam3_003.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_004.png :alt: $uy^{e}$ :srcset: /examples/Beam/images/sphx_glr_Beam3_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_005.png :alt: $u_y(x)$ :srcset: /examples/Beam/images/sphx_glr_Beam3_005.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_006.png :alt: $M_z(x)$ :srcset: /examples/Beam/images/sphx_glr_Beam3_006.png :class: sphx-glr-multi-img * .. image-sg:: /examples/Beam/images/sphx_glr_Beam3_007.png :alt: $T_y(x)$ :srcset: /examples/Beam/images/sphx_glr_Beam3_007.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none err uy: 1.99e-11% err Mz : 1.17e-11 % err Ty : 1.43e-11 % ==================== Mesh ==================== Element type: SEG3 Ne = 10, Nn = 21 ==================== Model ==================== solver:scipy ============= Boundary Conditions ============= Unspecified. =================== Results =================== Ux max = -0.00e+00 Ux min = -0.00e+00 Uy max = 0.00e+00 Uy min = -3.61e-03 =================== TicTac =================== Mesh: 44.024 ms Boundary Conditions: 26.226 µs Matrix: 24.618 ms Solver: 3.196 ms Display: 50.361 ms | .. code-block:: Python :lineno-start: 13 import matplotlib.pyplot as plt import numpy as np from EasyFEA import Display, Models, Mesher, ElemType, Simulations from EasyFEA.Geoms import Line, Point, Points if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Configuration # ---------------------------------------------- # geom L = 120 h = 20 b = 13 e = 2 # model E = 210000 v = 0.3 beamDim = 2 # must be >= 2 # load F = -800 # ---------------------------------------------- # Section # ---------------------------------------------- def DoSym(p: Point, n: np.ndarray) -> Point: pc = p.copy() pc.Symmetry(n=n) return pc p1 = Point(-b / 2, -h / 2) p2 = Point(b / 2, -h / 2) p3 = Point(b / 2, -h / 2 + e) p4 = Point(e / 2, -h / 2 + e, r=e) p5 = DoSym(p4, (0, 1)) p6 = DoSym(p3, (0, 1)) p7 = DoSym(p2, (0, 1)) p8 = DoSym(p1, (0, 1)) p9 = DoSym(p6, (1, 0)) p10 = DoSym(p5, (1, 0)) p11 = DoSym(p4, (1, 0)) p12 = DoSym(p3, (1, 0)) contour = Points([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], e / 3) section = Mesher().Mesh_2D(contour, elemType=ElemType.TRI3) # ---------------------------------------------- # Mesh # ---------------------------------------------- p1 = Point() pL = Point(x=L * 0.75) p2 = Point(x=L) line = Line(p1, p2, L / 9) beam = Models.Beam.Isotropic(beamDim, line, section, E, v) mesh = Mesher().Mesh_Beams([beam], additionalPoints=[pL]) # ---------------------------------------------- # Simulation # ---------------------------------------------- # Initialize the beam structure with the defined beam segments beamStructure = Models.Beam.BeamStructure([beam]) # Create the beam simulation simu = Simulations.Beam(mesh, beamStructure) dof_n = simu.Get_dof_n() # Apply boundary conditions simu.add_dirichlet(mesh.Nodes_Point(p1), [0] * dof_n, simu.Get_unknowns()) simu.add_dirichlet(mesh.Nodes_Point(p2), [0] * dof_n, simu.Get_unknowns()) simu.add_neumann(mesh.Nodes_Point(pL), [F], ["y"]) # Solve the beam problem and get displacement results sol = simu.Solve() simu.Save_Iter() # ---------------------------------------------- # Results # ---------------------------------------------- Display.Plot_Mesh(simu, L / 20 / sol.min()) ax = Display.Plot_Mesh(section) ax.set_title("Section") Display.Plot_BoundaryConditions(simu) Display.Plot_Result(simu, "uy", L / 20 / sol.min()) # ------------------------ # uy # ------------------------ # beam properties Iz = beam.Iz G = beam.mu A = section.area # general reactions and fixed-end moment (valid for arbitrary pL.x) a = pL.x b = L - a Ra = -F * b**2 * (3 * a + b) / L**3 # upward reaction at x=0 Rb = -F * a**2 * (a + 3 * b) / L**3 # upward reaction at x=L Ma = F * a * b**2 / L**2 # fixed-end moment at x=0 x = np.linspace(0, L, 100) x_n = mesh.coord[:, 0] x_e = x_n[mesh.connect].mean(1) # element centroid x-coords def uy_x(x): x = np.asarray(x) macaulay = np.where(x > a, (x - a) ** 3 / 6, 0.0) return (Ma / 2 * x**2 + Ra / 6 * x**3 + F * macaulay) / (E * Iz) uy = simu.Result("uy") err_uy = np.abs(uy_x(x_n) - uy).max() / np.abs(uy_x(a)) Display.MyPrint(f"\nerr uy: {err_uy * 100:.2e}%") ax_uy = Display.Init_Axes() ax_uy.plot(x, uy_x(x), label="Analytical", c="blue") ax_uy.scatter(x_n, uy, label="FE", c="red", marker="x", zorder=2) ax_uy.set_title("$u_y(x)$") ax_uy.legend() # ------------------------ # Mz # ------------------------ def Mz_x(x): x = np.asarray(x) Ma_plus_Ra_x = Ma + Ra * x return np.where(x <= a, Ma_plus_Ra_x, Ma_plus_Ra_x + F * (x - a)) Mz = simu.Result("Mz", nodeValues=False) err_Mz = np.abs(Mz_x(x_e) - Mz).max() / np.abs(Mz_x(x_e)).max() Display.MyPrint(f"\nerr Mz : {err_Mz * 100:.2e} %") axMz = Display.Init_Axes() axMz.plot(x, Mz_x(x), label="Analytical", c="blue") axMz.scatter(x_e, Mz, label="FE", c="red", marker="x", zorder=2) axMz.set_title("$M_z(x)$") axMz.legend() # ------------------------ # Ty # ------------------------ def Ty_x(x): x = np.asarray(x) return np.where(x < a, -Ra, Rb) Ty = simu.Result("Ty", nodeValues=False) err_Ty = np.abs(Ty_x(x_e) - Ty).max() / max(Ra, Rb) Display.MyPrint(f"\nerr Ty : {err_Ty * 100:.2e} %") ax_Ty = Display.Init_Axes() ax_Ty.step(x, Ty_x(x), label="Analytical", c="blue", where="mid") ax_Ty.scatter(x_e, Ty, label="FE", c="red", marker="x", zorder=2) ax_Ty.set_title("$T_y(x)$") ax_Ty.legend() print(simu) plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.532 seconds) .. _sphx_glr_download_examples_Beam_Beam3.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Beam3.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Beam3.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Beam3.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_