悬臂梁模态程序

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# ====================== 1. 悬臂梁模态基础参数 ======================
# 悬臂梁：x∈[0, L]，x=0固支，x=L自由端
L = 4.0  # 梁总长
num_node = 100  # 梁离散节点数
x_beam = np.linspace(0, L, num_node)

# 悬臂梁无量纲频率系数（前两阶）
lambda1 = 1.8751
lambda2 = 4.6941

# 模态振型函数 悬臂梁解析振型
def cantilever_mode_shape(x, lam):
    """
    悬臂梁归一化振型函数
    w(x) = cosh(lam x/L) - cos(lam x/L) - sigma*(sinh(lam x/L)-sin(lam x/L))
    sigma=(coshλ+cosλ)/(sinhλ+sinλ)
    """
    xi = x / L
    sigma = (np.cosh(lam) + np.cos(lam)) / (np.sinh(lam) + np.sin(lam))
    w = np.cosh(lam * xi) - np.cos(lam * xi) - sigma * (np.sinh(lam * xi) - np.sin(lam * xi))
    return w

# 两阶振型空间分布
phi1 = cantilever_mode_shape(x_beam, lambda1)
phi2 = cantilever_mode_shape(x_beam, lambda2)
# 归一化，自由端幅值统一为1
phi1 = phi1 / phi1[-1]
phi2 = phi2 / phi2[-1]

# 固有角频率
omega1 = lambda1**2
omega2 = lambda2**2

# 时间轴
tau_full = np.linspace(0, 15, 800)

# 各阶模态振动时程（自由端位移）
amp = 0.6
disp_mode1 = amp * np.cos(omega1 * tau_full)
disp_mode2 = amp * np.cos(omega2 * tau_full)
disp_total = disp_mode1 + disp_mode2

# 全梁任意时刻挠度：模态叠加 w(x,t)=phi1(x)*q1(t)+phi2(x)*q2(t)
def get_beam_deflection(t):
    q1 = amp * np.cos(omega1 * t)
    q2 = amp * np.cos(omega2 * t)
    w1 = phi1 * q1
    w2 = phi2 * q2
    w_tot = w1 + w2
    return w1, w2, w_tot

# ====================== 2. 画布布局 2行3列 ======================
fig = plt.figure(figsize=(18, 9))
gs = fig.add_gridspec(2, 3, height_ratios=[2, 1])

# 上排：自由端位移时程曲线
ax_time = fig.add_subplot(gs[0, :])  # 整行放时程
ax_mod1 = fig.add_subplot(gs[1, 0])
ax_mod2 = fig.add_subplot(gs[1, 1])
ax_total = fig.add_subplot(gs[1, 2])

# 时程曲线初始化
l_tot, = ax_time.plot([], [], "k-", lw=2, label="总振动（模态叠加）")
l_m1, = ax_time.plot([], [], "r-", lw=1.2, label="一阶弯曲模态")
l_m2, = ax_time.plot([], [], "b--", lw=1.2, label="二阶弯曲模态")

ax_time.set_xlim(0, 15)
ax_time.set_ylim(-amp*2.2, amp*2.2)
ax_time.set_xlabel(r"无量纲时间 $\tau$")
ax_time.set_ylabel("梁自由端挠度 w")
ax_time.grid(alpha=0.3)
ax_time.legend(loc="upper right")
ax_time.set_title("悬臂梁自由端振动时程曲线")

# 统一绘制悬臂梁画布函数
def setup_beam_axes(ax, title):
    ax.set_xlim(-0.5, L+0.5)
    ax.set_ylim(-1.0, 1.0)
    ax.set_aspect("equal")
    ax.set_title(title)
    ax.axis("off")
    # 左侧固支墙体
    wall = plt.Rectangle((-0.5, -1), 0.5, 2, color="dimgray", hatch="//")
    ax.add_patch(wall)
    # 梁主体线
    beam_line, = ax.plot([], [], color="#222222", lw=3)
    # 节点散点
    beam_scat = ax.scatter([], [], c="darkorange", s=12)
    return beam_line, beam_scat

# 三个视图初始化
line_m1, scat_m1 = setup_beam_axes(ax_mod1, "悬臂梁 第一阶模态")
line_m2, scat_m2 = setup_beam_axes(ax_mod2, "悬臂梁 第二阶模态")
line_tot, scat_tot = setup_beam_axes(ax_total, "两阶模态叠加总振动")

# ====================== 3. 动画更新逻辑 ======================
def update(frame):
    t_now = tau_full[frame]
    tau_data = tau_full[:frame]

    # 更新时程图
    l_tot.set_data(tau_data, disp_total[:frame])
    l_m1.set_data(tau_data, disp_mode1[:frame])
    l_m2.set_data(tau_data, disp_mode2[:frame])

    # 获取当前时刻全梁挠度
    w1, w2, wt = get_beam_deflection(t_now)

    # 一阶模态画面
    line_m1.set_data(x_beam, w1)
    scat_m1.set_offsets(np.column_stack([x_beam, w1]))

    # 二阶模态画面
    line_m2.set_data(x_beam, w2)
    scat_m2.set_offsets(np.column_stack([x_beam, w2]))

    # 叠加总振动画面
    line_tot.set_data(x_beam, wt)
    scat_tot.set_offsets(np.column_stack([x_beam, wt]))

    return l_tot, l_m1, l_m2, line_m1, scat_m1, line_m2, scat_m2, line_tot, scat_tot

# ========== 关键修改：interval从18改为9，帧率翻倍，速度快一倍 ==========
ani = FuncAnimation(
    fig, update, frames=len(tau_full),
    interval=9, blit=True, repeat=True
)

fig.suptitle("悬臂梁一阶、二阶弯曲模态与模态叠加同步动态仿真", fontsize=16)
plt.tight_layout()

# 导出GIF配套修改：fps改成60（原30），保证导出视频同样快一倍
# ani.save("cantilever_2modes_fast.gif", writer="pillow", fps=60)
plt.show()