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