Note
Go to the end to download the full example code.
Modal2#
Modal analysis of a structure.
12 import matplotlib.pyplot as plt
13 import numpy as np
14 from scipy.sparse import linalg, eye
15
16 from EasyFEA import Display, Models, Mesher, ElemType, Mesh, Simulations, PyVista
17
18
19 def Construct_struct(
20 L: float, e: float, t: float, meshSize: float = 0.0, openGmsh=False, verbosity=False
21 ) -> Mesh:
22 mesher = Mesher()
23
24 h = L - e - t
25
26 factory = mesher._factory
27
28 # create the pilars
29 pilar1 = [(3, factory.addBox(0, 0, 0, e, e, h))]
30 pilar2 = factory.copy(pilar1)
31 factory.translate(pilar2, L - e, 0, 0)
32 pilar3 = factory.copy(pilar1)
33 factory.translate(pilar3, L - e, L - e, 0)
34 pilar4 = factory.copy(pilar1)
35 factory.translate(pilar4, 0, L - e, 0)
36 pilars = factory.getEntities(3)
37
38 # creates the plate
39 plate = [(3, factory.addBox(0, 0, h, L, L, t))]
40
41 # creates the table (pilars + plate)
42 table, __ = factory.fragment(plate, pilars)
43
44 # creates the cuve (Empty Box)
45 box = [(3, factory.addBox(0, 0, L - e, L, L, L))]
46 inc = [(3, factory.addBox(e, e, L, L - 2 * e, L - 2 * e, L - 2 * e))]
47 cuve, __ = factory.cut(box, inc)
48
49 # creates the structure (table + cuve)
50 struct, __ = factory.fragment(table, cuve)
51
52 if meshSize > 0:
53 mesher.Set_meshSize(meshSize)
54
55 mesher._Set_PhysicalGroups()
56
57 mesher._Mesh_Generate(3, ElemType.TETRA10)
58
59 mesh = mesher._Mesh_Get_Mesh()
60
61 return mesh
62
63
64 if __name__ == "__main__":
65 Display.Clear()
66
67 # ----------------------------------------------
68 # Configuration
69 # ----------------------------------------------
70
71 # geom
72 L = 21 # m
73 e = 1
74 t = 0.5
75
76 # outputs
77 isFixed = True
78 Nmode = 3
79
80 # ----------------------------------------------
81 # Mesh
82 # ----------------------------------------------
83
84 mesh = Construct_struct(L, e, t, 0, False, False)
85
86 nodes_pilars = mesh.Nodes_Tags(["V0", "V1", "V2", "V3"])
87 elems_pilars = mesh.Elements_Tags(["V0", "V1", "V2", "V3"])
88
89 nodes_plate = mesh.Nodes_Tags(["V4"])
90 nodes_cuve = mesh.Nodes_Tags(["V5"])
91 nodesZ0 = mesh.Nodes_Conditions(lambda x, y, z: z == 0)
92 nodesSupZ0 = mesh.Nodes_Conditions(lambda x, y, z: z > 0)
93
94 plotter = PyVista.Plot_Elements(
95 mesh, nodes_pilars, color="red", alpha=0.5, label="Pilars"
96 )
97 PyVista.Plot_Elements(
98 mesh, nodes_plate, color="blue", alpha=0.5, plotter=plotter, label="Plate"
99 )
100 PyVista.Plot_Elements(
101 mesh, nodes_cuve, color="green", alpha=0.5, plotter=plotter, label="Cuve"
102 )
103 PyVista._setCameraPosition(plotter, 3, "yz", roll=-90)
104 plotter.zoom_camera(0.8)
105 plotter.add_legend()
106 plotter.show()
107
108 plotter = PyVista.Plot_Mesh(mesh, alpha=0.5)
109 PyVista._setCameraPosition(plotter, 3, "yz", roll=-90)
110 plotter.zoom_camera(0.8)
111 plotter.show()
112
113 # ----------------------------------------------
114 # Material
115 # ----------------------------------------------
116
117 E_pilars = 2000 * 1e9 # GPa
118 E_cuve = 20 * 1e9
119 E_plate = E_cuve
120
121 E = np.ones(mesh.Ne) * E_cuve
122 E[elems_pilars] = E_pilars
123
124 material = Models.Elastic.Isotropic(3, E, 0.3)
125
126 # ----------------------------------------------
127 # Simulation
128 # ----------------------------------------------
129 simu = Simulations.Elastic(mesh, material)
130 simu.rho = 7860 # kg/m3
131
132 simu.add_dirichlet(nodesZ0, [0] * 3, simu.Get_unknowns())
133 known, unknown = simu.Bc_dofs_known_unknown(simu.problemType)
134
135 plotter = PyVista.Plot_BoundaryConditions(simu)
136 PyVista._setCameraPosition(plotter, 3, "yz", roll=-90)
137 plotter.zoom_camera(0.8)
138 plotter.show()
139
140 K, C, M, F = simu.Get_K_C_M_F()
141
142 if isFixed:
143 K_t = K[unknown, :].tocsc()[:, unknown].tocsr()
144 M_t = M[unknown, :].tocsc()[:, unknown].tocsr()
145 else:
146 K_t = K + K.min() * eye(K.shape[0]) * 1e-12
147 M_t = M
148
149 eigenValues, eigenVectors = linalg.eigs(K_t, Nmode, M_t, sigma=0, which="LR")
150
151 eigenValues = eigenValues.real
152 eigenVectors = eigenVectors.real
153 freq_t = np.sqrt(eigenValues) / 2 / np.pi
154
155 # ----------------------------------------------
156 # Plot modes
157 # ----------------------------------------------
158 for n, eigenValue in enumerate(eigenValues):
159 if isFixed:
160 mode = np.zeros((mesh.Nn, 3))
161 mode[nodesSupZ0, :] = np.reshape(eigenVectors[:, n], (-1, 3))
162 else:
163 mode = np.reshape(eigenVectors[:, n], (-1, 3))
164
165 simu._Set_solutions(simu.problemType, mode.ravel())
166 simu.Save_Iter()
167
168 sol = np.linalg.norm(mode, axis=1)
169 deformFactor = L / 5 / np.abs(sol).max()
170
171 plotter = PyVista.Plot(simu, alpha=0.5)
172 PyVista.Plot(simu, None, deformFactor, alpha=0.8, color="r", plotter=plotter)
173 plotter.add_title(f"mode {n + 1}")
174 PyVista._setCameraPosition(plotter, 3, "yz", roll=-90)
175 plotter.zoom_camera(0.8)
176 plotter.show()
177
178 axModes = Display.Init_Axes()
179 axModes.plot(np.arange(eigenValues.size), freq_t, ls="", marker=".")
180 axModes.set_xticks(np.arange(eigenValues.size))
181 axModes.set_xlabel("modes")
182 axModes.set_ylabel("freq [Hz]")
183 axModes.grid()
184
185 plt.show()
Total running time of the script: (0 minutes 4.952 seconds)