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