Note
Go to the end to download the full example code.
Beam1#
Beam subjected to pure tensile loading.
err ux: 1.07e-13%
err N: 1.32e-13%
==================== Mesh ====================
Element type: SEG2
Ne = 10, Nn = 11
==================== Model ====================
<EasyFEA.Models.Beam._beam.BeamStructure object at 0x742244715290>
solver:scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
Ux max = 4.45e-05
Ux min = 0.00e+00
=================== TicTac ===================
Mesh: 12.346 ms
Boundary Conditions: 27.657 µs
Matrix: 11.383 ms
Solver: 1.063 ms
Display: 38.842 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 ro = 7800
35 v = 0.3
36
37 # load
38 g = 10
39 q = ro * g * (h * b)
40 F = 5000
41
42 # ----------------------------------------------
43 # Mesh
44 # ----------------------------------------------
45
46 elemType = ElemType.SEG2
47 beamDim = 1
48
49 # Create a section for the beam
50 mesher = Mesher()
51 section = mesher.Mesh_2D(Domain((0, 0), (b, h)))
52
53 p1 = (0, 0)
54 p2 = (L, 0)
55 line = Line(p1, p2, L / nL)
56 beam = Models.Beam.Isotropic(beamDim, line, section, E, v)
57
58 mesh = mesher.Mesh_Beams([beam], elemType=elemType)
59
60 # ----------------------------------------------
61 # Simulation
62 # ----------------------------------------------
63
64 # Initialize the beam structure with the defined beam segments
65 beamStructure = Models.Beam.BeamStructure([beam])
66
67 # Create the beam simulation
68 simu = Simulations.Beam(mesh, beamStructure)
69 dof_n = simu.Get_dof_n()
70
71 # Apply boundary conditions
72 simu.add_dirichlet(mesh.Nodes_Point(p1), [0] * dof_n, simu.Get_unknowns())
73 simu.add_lineLoad(mesh.nodes, [q], ["x"])
74 simu.add_neumann(mesh.Nodes_Point(p2), [F], ["x"])
75
76 # Solve the beam problem and get displacement results
77 sol = simu.Solve()
78 simu.Save_Iter()
79
80 # ----------------------------------------------
81 # Results
82 # ----------------------------------------------
83
84 Display.Plot_Mesh(simu, deformFactor=L / 10 / sol.max())
85 Display.Plot_BoundaryConditions(simu)
86 Display.Plot_Result(simu, "ux")
87
88 # ------------------------
89 # ux
90 # ------------------------
91
92 A = section.area
93 x = np.linspace(0, L, 100)
94 ux_x = lambda x: (F * x / (E * A)) + (ro * g * x / 2 / E * (2 * L - x))
95
96 ux = simu.Result("ux")
97 x_n = mesh.coord[:, 0]
98 err_ux = np.abs(ux_x(x_n) - ux).max() / np.abs(ux_x(L))
99 Display.MyPrint(f"\nerr ux: {err_ux * 100:.2e}%")
100
101 axUx = Display.Init_Axes()
102 axUx.plot(x, ux_x(x), label="Analytical", c="blue")
103 axUx.scatter(x_n, ux, label="FE", c="red", marker="x", zorder=2)
104 axUx.set_title("$u_x(x)$")
105 axUx.legend()
106
107 # ------------------------
108 # N
109 # ------------------------
110
111 N_x = lambda x: F + q * (L - x)
112
113 x_e = x_n[mesh.connect].mean(1) # element centroid x-coords
114 N = simu.Result("N", nodeValues=False)
115 err_N = np.abs(N_x(x_e) - N).max() / np.abs(N_x(x_e)).max()
116 Display.MyPrint(f"\nerr N: {err_N * 100:.2e}%")
117
118 axN = Display.Init_Axes()
119 axN.plot(x, N_x(x), label="Analytical", c="blue")
120 axN.scatter(x_e, N, label="FE", c="red", marker="x", zorder=2)
121 axN.set_title("$N(x)$")
122 axN.legend()
123
124 print(simu)
125
126 plt.show()
Total running time of the script: (0 minutes 0.339 seconds)