Note
Go to the end to download the full example code.
Beam1#
Beam subjected to pure tensile loading.
err ux: 1.52e-14%
err N: 4.40e-14%
==================== Mesh ====================
Element type: SEG2
Ne = 10, Nn = 11
==================== Model ====================
<EasyFEA.Models.Beam._beam.BeamStructure object at 0x7e08ddefee90>
solver:scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
Ux max = 4.45e-05
Ux min = 0.00e+00
=================== TicTac ===================
Mesh: 12.488 ms
Boundary Conditions: 38.862 µs
Matrix: 19.333 ms
Solver: 2.752 ms
Display: 38.957 ms
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Models, Mesher, ElemType, Simulations
17 from EasyFEA.Geoms import Domain, Line
18
19 if __name__ == "__main__":
20 Display.Clear()
21
22 # ----------------------------------------------
23 # Configuration
24 # ----------------------------------------------
25
26 # geom
27 L = 10
28 nL = 10
29 h = 0.1
30 b = 0.1
31
32 # model
33 E = 200000e6
34 v = 0.3
35 rho = 7800
36 beamDim = 1
37
38 # load
39 g = 10
40 q = rho * g * (h * b)
41 F = 5000
42
43 # ----------------------------------------------
44 # Mesh
45 # ----------------------------------------------
46
47 # Create a section for the beam
48 mesher = Mesher()
49 section = mesher.Mesh_2D(Domain((0, 0), (b, h)))
50
51 p1 = (0, 0)
52 p2 = (L, 0)
53 line = Line(p1, p2, L / nL)
54 beam = Models.Beam.Isotropic(beamDim, line, section, E, v)
55
56 mesh = mesher.Mesh_Beams([beam], elemType=ElemType.SEG2)
57
58 # ----------------------------------------------
59 # Simulation
60 # ----------------------------------------------
61
62 # Initialize the beam structure with the defined beam segments
63 beamStructure = Models.Beam.BeamStructure([beam])
64
65 # Create the beam simulation
66 simu = Simulations.Beam(mesh, beamStructure)
67 dof_n = simu.Get_dof_n()
68
69 # Apply boundary conditions
70 simu.add_dirichlet(mesh.Nodes_Point(p1), [0] * dof_n, simu.Get_unknowns())
71 simu.add_lineLoad(mesh.nodes, [q], ["x"])
72 simu.add_neumann(mesh.Nodes_Point(p2), [F], ["x"])
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.max())
83 Display.Plot_BoundaryConditions(simu)
84 Display.Plot_Result(simu, "ux")
85
86 # ------------------------
87 # ux
88 # ------------------------
89
90 A = section.area
91 x = np.linspace(0, L, 100)
92 ux_x = lambda x: (F * x / (E * A)) + (rho * g * x / 2 / E * (2 * L - x))
93
94 ux = simu.Result("ux")
95 x_n = mesh.coord[:, 0]
96 err_ux = np.abs(ux_x(x_n) - ux).max() / np.abs(ux_x(L))
97 Display.MyPrint(f"\nerr ux: {err_ux * 100:.2e}%")
98
99 axUx = Display.Init_Axes()
100 axUx.plot(x, ux_x(x), label="Analytical", c="blue")
101 axUx.scatter(x_n, ux, label="FE", c="red", marker="x", zorder=2)
102 axUx.set_title("$u_x(x)$")
103 axUx.legend()
104
105 # ------------------------
106 # N
107 # ------------------------
108
109 N_x = lambda x: F + q * (L - x)
110
111 x_e = x_n[mesh.connect].mean(1) # element centroid x-coords
112 N = simu.Result("N", nodeValues=False)
113 err_N = np.abs(N_x(x_e) - N).max() / np.abs(N_x(x_e)).max()
114 Display.MyPrint(f"\nerr N: {err_N * 100:.2e}%")
115
116 axN = Display.Init_Axes()
117 axN.plot(x, N_x(x), label="Analytical", c="blue")
118 axN.scatter(x_e, N, label="FE", c="red", marker="x", zorder=2)
119 axN.set_title("$N(x)$")
120 axN.legend()
121
122 print(simu)
123
124 plt.show()
Total running time of the script: (0 minutes 0.350 seconds)