Note
Go to the end to download the full example code.
Beam2#
A cantilever beam undergoing bending deformation.
err uy: 1.21e-10 %err rz: 1.12e-10 %
==================== Mesh ====================
Element type: SEG2
Ne = 10, Nn = 11
==================== Model ====================
<EasyFEA.Models.Beam._beam.BeamStructure object at 0x700eaf5ce410>
solver : scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
Ux max = 0.00e+00
Ux min = 0.00e+00
Uy max = 0.00e+00
Uy min = -9.22e-01
=================== TicTac ===================
Mesh : 15.353 ms
Boundary Conditions : 21.935 µs
Matrix : 8.173 ms
Solver : 1.539 ms
Display : 51.848 ms
12 import matplotlib.pyplot as plt
13 import numpy as np
14
15 from EasyFEA import Display, Models, Mesher, ElemType, Simulations
16 from EasyFEA.Geoms import Domain, Line
17
18 if __name__ == "__main__":
19 Display.Clear()
20
21 # ----------------------------------------------
22 # Configuration
23 # ----------------------------------------------
24
25 # geom
26 L = 120
27 nL = 10
28 h = 13
29 b = 13
30
31 # model
32 E = 210000
33 v = 0.3
34
35 # load
36 load = 800
37
38 # ----------------------------------------------
39 # Mesh
40 # ----------------------------------------------
41
42 elemType = ElemType.SEG2
43 beamDim = 2 # must be >= 2
44
45 # Create a section object for the beam mesh
46 mesher = Mesher()
47 section = mesher.Mesh_2D(Domain((0, 0), (b, h)))
48
49 p1 = (0, 0)
50 p2 = (L, 0)
51 line = Line(p1, p2, L / nL)
52 beam = Models.Beam.Isotropic(beamDim, line, section, E, v)
53
54 mesh = mesher.Mesh_Beams([beam], elemType=elemType)
55
56 # ----------------------------------------------
57 # Simulation
58 # ----------------------------------------------
59
60 Iy = beam.Iy
61 Iz = beam.Iz
62
63 # Initialize the beam structure with the defined beam segments
64 beamStructure = Models.Beam.BeamStructure([beam])
65
66 # Create the beam simulation
67 simu = Simulations.Beam(mesh, beamStructure)
68 dof_n = simu.Get_dof_n()
69
70 # Apply boundary conditions
71 simu.add_dirichlet(mesh.Nodes_Point(p1), [0] * dof_n, simu.Get_unknowns())
72 simu.add_neumann(mesh.Nodes_Point(p2), [-load], ["y"])
73
74 # Solve the beam problem and get displacement results
75 sol = simu.Solve()
76 simu.Save_Iter()
77
78 # ----------------------------------------------
79 # Results
80 # ----------------------------------------------
81
82 Display.Plot_Mesh(simu, deformFactor=-L / 10 / sol.min())
83 Display.Plot_BoundaryConditions(simu)
84 Display.Plot_Result(simu, "uy")
85
86 rz = simu.Result("rz")
87 v = simu.Result("uy")
88
89 x = np.linspace(0, L, 100)
90 uy_x = load / (E * Iz) * (x**3 / 6 - (L * x**2) / 2)
91
92 flecheanalytique = load * L**3 / (3 * E * Iz)
93 err_uy = np.abs(flecheanalytique + v.min()) / flecheanalytique
94 Display.MyPrint(f"err uy: {err_uy * 100:.2e} %")
95
96 # Plot the analytical and finite element solutions for vertical displacement (v)
97 axUy = Display.Init_Axes()
98 axUy.plot(x, uy_x, label="Analytical", c="blue")
99 axUy.scatter(mesh.coord[:, 0], v, label="FE", c="red", marker="x", zorder=2)
100 axUy.set_title("$u_y(x)$")
101 axUy.legend()
102
103 rz_x = load / E / Iz * (x**2 / 2 - L * x)
104 rotalytique = load * L**2 / (2 * E * Iz)
105 err_rz = np.abs(rotalytique + rz.min()) / rotalytique
106 Display.MyPrint(f"err rz: {err_rz * 100:.2e} %")
107
108 # Plot the analytical and finite element solutions for rotation (rz)
109 axRz = Display.Init_Axes()
110 axRz.plot(x, rz_x, label="Analytical", c="blue")
111 axRz.scatter(mesh.coord[:, 0], rz, label="FE", c="red", marker="x", zorder=2)
112 axRz.set_title("$r_z(x)$")
113 axRz.legend()
114
115 print(simu)
116
117 plt.show()
Total running time of the script: (0 minutes 0.432 seconds)