Note
Go to the end to download the full example code.
Modal1#
Modal analysis of a wall 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, ElemType, Simulations, PyVista
18 from EasyFEA.Geoms import Domain
19
20
21 if __name__ == "__main__":
22 Display.Clear()
23
24 # ----------------------------------------------
25 # Configuration
26 # ----------------------------------------------
27
28 dim = 3
29
30 # outputs
31 isFixed = True
32 Nmode = 3
33
34 # ----------------------------------------------
35 # Mesh
36 # ----------------------------------------------
37 contour = Domain((0, 0), (1, 1), 1 / 10)
38 thickness = 1 / 10
39
40 if dim == 2:
41 mesh = contour.Mesh_2D([], ElemType.QUAD9, isOrganised=True)
42 else:
43 mesh = contour.Mesh_Extrude(
44 [], [0, 0, -thickness], [2], ElemType.HEXA27, isOrganised=True
45 )
46 nodesY0 = mesh.Nodes_Conditions(lambda x, y, z: y == 0)
47 nodesSupY0 = mesh.Nodes_Conditions(lambda x, y, z: y > 0)
48
49 PyVista.Plot_Mesh(mesh).show()
50
51 # ----------------------------------------------
52 # Simulation
53 # ----------------------------------------------
54 material = Models.Elastic.Isotropic(dim, planeStress=True, thickness=thickness)
55
56 simu = Simulations.Elastic(mesh, material)
57
58 simu.Solver_Set_Hyperbolic_Algorithm(0.1)
59
60 K, C, M, F = simu.Get_K_C_M_F()
61
62 if isFixed:
63 simu.add_dirichlet(nodesY0, [0] * dim, simu.Get_unknowns())
64 known, unknown = simu.Bc_dofs_known_unknown(simu.problemType)
65 K_t = K[unknown, :].tocsc()[:, unknown].tocsr()
66 M_t = M[unknown, :].tocsc()[:, unknown].tocsr()
67
68 else:
69 K_t = K + K.min() * eye(K.shape[0]) * 1e-12
70 M_t = M
71
72 # eigenValues, eigenVectors = linalg.eigs(K_t, Nmode, M_t, which="SR")
73 eigenValues, eigenVectors = linalg.eigs(K_t, Nmode, M_t, sigma=0, which="LR")
74
75 eigenValues = eigenValues.real
76 eigenVectors = eigenVectors.real
77 freq_t = np.sqrt(eigenValues) / 2 / np.pi
78
79 # ----------------------------------------------
80 # Plot modes
81 # ----------------------------------------------
82 for n in range(eigenValues.size):
83 if isFixed:
84 mode = np.zeros((mesh.Nn, dim))
85 mode[nodesSupY0, :] = np.reshape(eigenVectors[:, n], (-1, dim))
86 else:
87 mode = np.reshape(eigenVectors[:, n], (-1, dim))
88
89 simu._Set_solutions(simu.problemType, mode.ravel())
90 simu.Save_Iter()
91
92 sol = np.linalg.norm(mode, axis=1)
93 deformFactor = 1 / 5 / np.abs(sol).max()
94
95 plotter = PyVista.Plot(simu, alpha=0.5)
96 PyVista.Plot(simu, None, deformFactor, alpha=0.8, color="r", plotter=plotter)
97 plotter.add_title(f"mode {n + 1}")
98 plotter.show()
99
100 axModes = Display.Init_Axes()
101 axModes.plot(np.arange(eigenValues.size), freq_t, ls="", marker=".")
102 axModes.set_xlabel("modes")
103 axModes.set_xticks(np.arange(eigenValues.size))
104 axModes.set_ylabel("freq [Hz]")
105 axModes.grid()
106
107 plt.show()
Total running time of the script: (0 minutes 1.390 seconds)