Note
Go to the end to download the full example code.
Elas6#
A bi-fixed beam undergoing bending deformation.
Beam: section uses linear TRI3 elements — _Get_shear_correction_factor converges at O(h²). Use TRI6 / QUAD8 (or finer mesh) for accurate k.
err beam model : 0.00 %
err 3d model : 0.42 %
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Models, Mesher, ElemType, Simulations, PyVista
17 from EasyFEA.Geoms import Point, Points, Line
18
19 if __name__ == "__main__":
20 Display.Clear()
21
22 # ----------------------------------------------
23 # Configuration
24 # ----------------------------------------------
25
26 # geom
27 L = 500 # mm
28 hx = 13
29 hy = 20
30 e = 2
31
32 # model
33 E = 210000 # MPa (Young's modulus)
34 v = 0.3 # Poisson's ratio
35 coef = 1
36
37 # load
38 load = 800 # N
39
40 # ----------------------------------------------
41 # Section
42 # ----------------------------------------------
43
44 elemType = ElemType.TETRA4
45
46 def DoSym(p: Point, n: np.ndarray) -> Point:
47 pc = p.copy()
48 pc.Symmetry(n=n)
49 return pc
50
51 p1 = Point(-hx / 2, -hy / 2)
52 p2 = Point(hx / 2, -hy / 2)
53 p3 = Point(hx / 2, -hy / 2 + e)
54 p4 = Point(e / 2, -hy / 2 + e, r=e)
55 p5 = DoSym(p4, (0, 1))
56 p6 = DoSym(p3, (0, 1))
57 p7 = DoSym(p2, (0, 1))
58 p8 = DoSym(p1, (0, 1))
59 p9 = DoSym(p6, (1, 0))
60 p10 = DoSym(p5, (1, 0))
61 p11 = DoSym(p4, (1, 0))
62 p12 = DoSym(p3, (1, 0))
63 section = Points([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], e / 1)
64 meshSection = section.Mesh_2D()
65
66 section.Get_Contour().Plot()
67
68 section.Rotate(-90, direction=(0, 1))
69
70 # ----------------------------------------------
71 # Mesh
72 # ----------------------------------------------
73 mesh1 = section.Mesh_Extrude([], [L / 2, 0, 0], [50], elemType)
74 mesh2 = mesh1.copy()
75 mesh2.Translate(-L / 2)
76 mesh = mesh1.Merge([mesh1, mesh2])
77
78 nodes_fixed = mesh.Nodes_Conditions(lambda x, y, z: (x == -L / 2) | (x == L / 2))
79 nodes_load = mesh.Nodes_Line(Line(p7, p8))
80
81 # ----------------------------------------------
82 # Beam simulation
83 # ----------------------------------------------
84
85 beam = Models.Beam.Isotropic(
86 2, Line(Point(-L / 2), Point(L / 2), L / 10), meshSection, E, v
87 )
88
89 mesh_beam = Mesher().Mesh_Beams([beam], ElemType.SEG3)
90
91 beams = Models.Beam.BeamStructure([beam])
92 simu_beam = Simulations.Beam(mesh_beam, beams)
93
94 simu_beam.add_dirichlet(
95 mesh_beam.Nodes_Conditions(lambda x, y, z: (x == -L / 2) | (x == L / 2)),
96 [0] * simu_beam.Get_dof_n(),
97 simu_beam.Get_unknowns(),
98 )
99 simu_beam.add_neumann(
100 mesh_beam.Nodes_Conditions(lambda x, y, z: (x == 0)), [-load], ["y"]
101 )
102 simu_beam.Solve()
103
104 Display.Plot_Result(simu_beam, "uy")
105
106 u_an = load * L**3 / (192 * E * beam.Iz)
107
108 uy_1d = np.abs(simu_beam.Result("uy").min())
109
110 print(f"err beam model : {np.abs(u_an - uy_1d) / u_an * 100:.2f} %")
111
112 # ----------------------------------------------
113 # Elastic simulation
114 # ----------------------------------------------
115
116 material = Models.Elastic.Isotropic(3, E, v)
117 simu = Simulations.Elastic(mesh, material)
118
119 simu.add_dirichlet(nodes_fixed, [0] * 3, simu.Get_unknowns())
120 simu.add_lineLoad(nodes_load, [-load / hx], ["y"])
121 sol = simu.Solve()
122
123 uy_3d = np.abs(simu.Result("uy").min())
124
125 print(f"err 3d model : {np.abs(u_an - uy_3d) / u_an * 100:.2f} %")
126
127 # ----------------------------------------------
128 # Results
129 # ----------------------------------------------
130
131 plotter = PyVista._Plotter(shape=(2, 1))
132 PyVista.Plot(
133 simu,
134 "ux",
135 coef=1 / coef,
136 nColors=20,
137 plotMesh=True,
138 plotter=plotter,
139 verticalColobar=False,
140 )
141 PyVista.Plot_BoundaryConditions(simu, plotter=plotter)
142 plotter.subplot(1, 0)
143 PyVista.Plot(
144 simu, "uy", coef=1 / coef, nColors=20, plotter=plotter, verticalColobar=False
145 )
146 plotter.show()
147
148 plt.show()
Total running time of the script: (0 minutes 2.596 seconds)