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