勘探震源激发平板模态分析
import numpy as np
import scipy.linalg as la
from scipy.sparse import coo_matrix, csr_matrix, lil_matrix
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# ======================
# 有限元建模模块
# ======================
class PlateStructure:
def __init__(self, Lx=2.0, Ly=1.5, h=0.05, E=210e9, rho=7850):
"""震源激发平板结构参数
Args:
Lx: 平板长度(m)
Ly: 平板宽度(m)
h: 厚度(m)
E: 弹性模量(Pa)
rho: 密度(kg/m³)
"""
self.geometry = (Lx, Ly, h)
self.material = (E, rho)
# 网格参数
self.nx = 8 # X方向单元数
self.ny = 6 # Y方向单元数
self.nodes = None
self.elements = []
def generate_mesh(self):
"""生成结构化四边形网格"""
x = np.linspace(0, self.geometry[0], self.nx + 1)
y = np.linspace(0, self.geometry[1], self.ny + 1)
# 生成节点坐标 (Nx3数组: [x, y, z])
self.nodes = np.array([[xi, yi, 0] for yi in y for xi in x])
# 生成单元连接关系 (四节点单元)
for j in range(self.ny):
for i in range(self.nx):
n1 = i + j * (self.nx + 1)
n2 = n1 + 1
n3 = n2 + (self.nx + 1)
n4 = n3 - 1
self.elements.append([n1, n2, n3, n4])
def get_stiffness_matrix(self):
"""生成四节点板单元刚度矩阵(简化版本)"""
E, nu = self.material[0], 0.3
D = E * self.geometry[2] ** 3 / (12 * (1 - nu ** 2))
# 简化的12x12刚度矩阵（弯曲主导）
Ke = np.array([
[12, 3, 0, -6, 3, 0, -12, -3, 0, 6, -3, 0],
[3, 12, 0, -3, 6, 0, -3, -12, 0, 3, -6, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-6, -3, 0, 12, -3, 0, 6, 3, 0, -12, 3, 0],
[3, 6, 0, -3, 12, 0, -3, -6, 0, 3, -12, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-12, -3, 0, 6, -3, 0, 12, 3, 0, -6, 3, 0],
[-3, -12, 0, 3, -6, 0, 3, 12, 0, -3, 6, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[6, 3, 0, -12, 3, 0, -6, -3, 0, 12, -3, 0],
[-3, -6, 0, 3, -12, 0, 3, 6, 0, -3, 12, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
]) * D / (self.geometry[0] ** 3)
return csr_matrix(Ke)
def get_mass_matrix(self):
"""生成一致质量矩阵"""
rho = self.material[1]
h = self.geometry[2]
area = (self.geometry[0] / self.nx) * (self.geometry[1] / self.ny)
# 单元质量均分到各节点
Me = np.eye(12) * (rho * h * area) / 4
return csr_matrix(Me)
# ======================
# 模态分析模块
# ======================
class ModalAnalyzer:
def __init__(self, structure):
self.struct = structure
self.K = None # 全局刚度矩阵
self.M = None # 全局质量矩阵
self.freqs = None # 固有频率(Hz)
self.modes = None # 模态振型
def assemble_matrices(self):
"""组装全局刚度矩阵和质量矩阵"""
n_nodes = len(self.struct.nodes)
ndof = 3 * n_nodes # 总自由度
# 初始化稀疏矩阵
K = lil_matrix((ndof, ndof), dtype=np.float64)
M = lil_matrix((ndof, ndof), dtype=np.float64)
# 遍历所有单元进行组装
for elem in self.struct.elements:
Ke = self.struct.get_stiffness_matrix()
Me = self.struct.get_mass_matrix()
# 获取单元自由度索引
dof_indices = []
for node in elem:
dof_indices.extend([3 * node, 3 * node + 1, 3 * node + 2])
# 将单元矩阵添加到全局矩阵
for i in range(12):
for j in range(12):
K[dof_indices[i], dof_indices[j]] += Ke[i, j]
M[dof_indices[i], dof_indices[j]] += Me[i, j]
self.K = K.tocsr()
self.M = M.tocsr()
def apply_boundary_conditions(self):
"""应用固定边界条件（四角节点全约束）"""
fix_nodes = [0, self.struct.nx, -1, -self.struct.nx - 1]
fixed_dofs = []
for n in fix_nodes:
fixed_dofs.extend([3 * n, 3 * n + 1, 3 * n + 2]) # 约束x,y,z方向
# 创建自由度的掩码
free_dofs = np.setdiff1d(np.arange(self.K.shape[0]), fixed_dofs)
return free_dofs
def solve_modes(self, num_modes=6):
"""求解特征值问题"""
free_dofs = self.apply_boundary_conditions()
# 提取缩减矩阵
K_red = self.K[free_dofs, :][:, free_dofs]
M_red = self.M[free_dofs, :][:, free_dofs]
# 求解广义特征值问题
evals, evecs = la.eigh(K_red.toarray(), M_red.toarray())
# 重构全局模态振型
self.freqs = np.sqrt(np.abs(evals[:num_modes])) / (2 * np.pi)
self.modes = np.zeros((self.K.shape[0], num_modes))
self.modes[free_dofs, :] = evecs[:, :num_modes]
# ======================
# 可视化模块
# ======================
class Visualizer:
@staticmethod
def plot_mode_shape(struct, mode, freq, ax=None):
"""绘制三维模态振型
Args:
mode: 模态振型向量
freq: 对应频率
ax: 可选的绘图轴对象
"""
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
nodes = struct.nodes.copy()
scale = 0.2 * np.max(nodes[:, :2]) / np.max(np.abs(mode))
# 将振型位移叠加到坐标
nodes[:, 2] += mode[::3] * scale # 取Z方向位移
# 绘制变形后的网格
x = nodes[:, 0]
y = nodes[:, 1]
z = nodes[:, 2]
# 绘制单元面
for elem in struct.elements:
verts = nodes[elem, :]
ax.plot_trisurf(verts[:, 0], verts[:, 1], verts[:, 2], alpha=0.6)
ax.set_title(f'Mode Shape @ {freq:.1f}Hz')
ax.set_xlabel('X')
ax.set_ylabel('Y')
return ax
@staticmethod
def compare_frequencies(freq_orig, freq_opt):
"""绘制频率对比柱状图"""
plt.figure(figsize=(10, 6))
x = np.arange(len(freq_orig))
plt.bar(x - 0.2, freq_orig, 0.4, label='Original Design')
plt.bar(x + 0.2, freq_opt, 0.4, label='Optimized Design')
plt.xticks(x, [f'Mode {i + 1}' for i in x])
plt.xlabel('Mode Order')
plt.ylabel('Frequency (Hz)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# ======================
# 主程序
# ======================
if __name__ == "__main__":
# 原始设计建模
print("正在建立原始设计模型...")
orig = PlateStructure(Lx=2.0, Ly=1.5, h=0.05)
orig.generate_mesh()
# 模态分析
analyzer_orig = ModalAnalyzer(orig)
analyzer_orig.assemble_matrices()
analyzer_orig.solve_modes(num_modes=6)
print("原始设计前6阶频率:", analyzer_orig.freqs)
# 优化设计（增加厚度）
print("\n正在建立优化设计模型...")
optimized = PlateStructure(Lx=2.0, Ly=1.5, h=0.07)
optimized.generate_mesh()
analyzer_opt = ModalAnalyzer(optimized)
analyzer_opt.assemble_matrices()
analyzer_opt.solve_modes(num_modes=6)
print("优化设计前6阶频率:", analyzer_opt.freqs)
# 结果可视化
print("\n生成可视化结果...")
fig = plt.figure(figsize=(15, 10))
# 原始设计前三阶振型
for i in range(3):
ax = fig.add_subplot(2, 3, i + 1, projection='3d')
Visualizer.plot_mode_shape(orig, analyzer_orig.modes[:, i],
analyzer_orig.freqs[i], ax)
# 优化设计前三阶振型
for i in range(3):
ax = fig.add_subplot(2, 3, i + 4, projection='3d')
Visualizer.plot_mode_shape(optimized, analyzer_opt.modes[:, i],
analyzer_opt.freqs[i], ax)
plt.tight_layout()
# 频率对比
Visualizer.compare_frequencies(analyzer_orig.freqs, analyzer_opt.freqs)
plt.show()