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