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