Note
Go to the end to download the full example code.
Elas6#
A bi-fixed beam undergoing bending deformation.
err beam model : 0.00 %
err 3d model : 0.36 %
13 import matplotlib.pyplot as plt
14 import numpy as np
15
16 from EasyFEA import Display, Models, Mesher, ElemType, Simulations, PyVista, MeshIO
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 = MeshIO.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 1.303 seconds)