企业官网建设流程全解析

简介：面向结构工程与振动分析学习者，这份MATLAB脚本包围绕多自由度体系动力响应计算，重点演示振型叠加法在求解固有频率、振型及地震、风荷载响应中的应用。资源共5个文件，全部为.m源码脚本，压缩包仅3KB，涵盖Duhamel积分、矩阵迭代求解、特征多项式法以及多自由度主计算程序等模块，结构简洁清晰。已有216人学习下载。脚本代码量小但逻辑完整，适合对照理论公式逐行理解振型叠加法的实现流程，既可用于课堂教学演示与课程作业，也可作为小型科研项目验证工具。通过运行这些脚本，读者能够快速掌握从建立质量、刚度矩阵，到求解特征值问题，再到叠加各阶振型获得动力响应的完整计算链条。

1. 为什么说振型叠加是“多自由度动力响应”的必选项

接手的模型只要超过三个自由度，我第一反应不会直接写Newmark–β做全自由度时程积分。步长被最高固有频率卡死，模型稍大一点，几秒钟的瞬态就要等半天。而振型叠加法的核心是先把N维耦合方程投影到特征空间，得到N个独立的单自由度方程，再挑出若干参与度高的模态做积分，最后乘回振型矩阵。它不只是加速手段，更能告诉你哪些模态真正参与了响应。本文按“运动方程→特征值解耦→截断准则→响应重构→误差修正”这条路径，给出可直接复制的Python实现，重点讲清楚阻尼比怎么给、保留几阶模态、以及为什么加速度响应不能直接用模态位移法叠加。适合结构动力分析、机械振动仿真和想把算法落到脚本里的工程师。

2. 运动方程到模态解耦：振型叠加的前提与截断原则

2.1 多自由度体系的运动方程与矩阵组装

多自由度体系在时域内的运动方程通常写成：

M ü + C u̇ + K u = F(t)

M、C、K 分别是质量、阻尼和刚度矩阵，F(t) 是荷载向量。矩阵的组装方式直接决定后续特征值问题的规模。常见做法是先对所有自由度编号，再按单元装配。质量矩阵有两种典型选择：集中质量矩阵把质量放在对角线上，形式简单，低阶频率会略高；一致质量矩阵由形函数积分得到，更接近真实惯性分布，但计算代价更高。下表是两者在动力分析中的取舍。

一个简单的两自由度串联质量块模型，可以用下面代码组装出 M 和 K：

import numpy as np # 两自由度系统，k1连接固定端与m1，k2连接m1与m2 k1, k2 = 1000.0, 800.0 m1, m2 = 2.0, 1.0 K = np.array([[k1 + k2, -k2], [-k2, k2]]) M = np.diag([m1, m2]) print("K =\n", K) print("M =\n", M)

这里 K[0,0] 同时包含 k1 和 k2 贡献，K[0,1] 表示两个自由度之间的弹性耦合。集中质量下 M 是对角的；如果使用一致质量，M 的非对角元会被填充，后续特征值求解仍不变。

2.2 广义特征值问题与模态正交性

振型叠加的数学基础是广义特征值问题：

K φ = ω² M φ

其中 ω 是固有圆频率，φ 是振型向量。在M正定、K对称的前提下，特征值都是实数。实际计算中，我优先使用scipy.linalg.eigh而不是np.linalg.eig，因为它内部先做Cholesky分解，把广义问题转成标准特征值问题，数值稳定性更好，对接近奇异的M也能给出警告而不是直接报错。

from scipy.linalg import eigh # 求解广义特征值问题，特征值升序排列 omega2, phi = eigh(K, M) omega = np.sqrt(omega2) print("固有圆频率(rad/s):", omega) print("振型矩阵(按列):\n", phi)

eigh返回的振型矩阵已经相对质量矩阵归一化，即 φ.T @ M @ φ ≈ I。这个性质是模态坐标解耦的起点。验证方式是打印phi.T @ M @ phi，如果非对角元素接近 1e-12 量级，就可以放心往下做。若出现明显非对角项，通常是M或K组装错误，需要回头检查节点编号顺序。

当M包含零质量自由度时，M奇异，eigh无法直接求解。工程上的处理办法是先做静力凝聚，把无质量自由度消去，再求解凝聚后的模型。

2.3 模态截断：保留“被激励起来”的模态

特征值解出后，N阶模态都可用，但振型叠加的意义在于只保留少数模态。截断不是简单按频率排序，而是要结合激励的频率成分和空间方向。对于没有外荷载位置的局部响应，高频模态可能贡献显著。常用判据是累计有效质量参与系数。

在质量归一化振型下，第 i 阶模态参与系数为 Γ_i = φ_iᵀ M r，其中 r 是荷载空间分布向量。累计有效质量比的计算公式为：

ratio_j = Σ_{i≤j} Γ_i² / (rᵀ M r)

# r 表示荷载方向，例如所有自由度水平移动 r = np.array([1.0, 1.0]) gamma = phi.T @ M @ r total_mass = r @ M @ r cum_ratio = np.cumsum(gamma**2) / total_mass print("模态参与系数:", gamma) print("累计有效质量比:", cum_ratio)

gamma的正负号取决于振型方向，对有效质量计算没有影响。累计有效质量比超过 0.9 后，再增加模态对整体位移的改进就非常有限。但若是计算层间剪力或构件内力，还需要检查响应点的模态振型分量。

下表是三种常用截断准则的适用场景：

对剪切型结构，前三阶模态通常已经贡献 90% 以上的基底剪力，但扭转模态或局部弯曲模态需要单独检查。很多工程脚本只按频率截止，容易把参与质量很低但局部影响大的模态错误截掉。

3. 振型叠加实现路径：从模态坐标到物理响应

3.1 坐标变换与模态力计算

把物理位移写成模态坐标线性组合：

u = Φ q

代入运动方程并左乘 Φᵀ，利用正交性后，M_modal 和 K_modal 都变成对角阵。阻尼若采用 Rayleigh 阻尼 C = αM + βK，则 C_modal 也为对角阵，方程组解耦为 N 个单自由度方程。Rayleigh 阻尼的两个系数通过两个参考频率的期望阻尼比确定。

# 第1阶和第3阶目标阻尼比设为0.02 zeta1, zeta3 = 0.02, 0.02 A = np.array([[1/omega[0], omega[0]], [1/omega[2], omega[2]]]) alpha, beta = np.linalg.solve(A, [2*zeta1, 2*zeta3]) C = alpha * M + beta * K M_modal = phi.T @ M @ phi C_modal = phi.T @ C @ phi K_modal = phi.T @ K @ phi print("M_modal:\n", np.round(M_modal, 4)) print("C_modal:\n", np.round(C_modal, 4))

上面代码里，omega[2]是第三阶圆频率，选取第一阶和第三阶作为两个控制点，是因为通常一阶是结构主频，第三阶能覆盖目标激励频率范围。如果直接给模态阻尼比，就不需要组装 C，每阶单独指定 zeta_i 即可。下表对比了两种阻尼建模的特点：

非比例阻尼在建筑结构里较少见，如果要处理，就不能按单自由度逐模态积分，得转而求解完整状态空间方程，计算量会明显增加。

3.2 各模态单自由度响应的时域积分

解耦后的第 i 阶方程为：

q̈_i + 2ζ_i ω_i q̇_i + ω_i² q_i = F_i(t) / M_i

工程上最稳妥的求解方式是调用scipy.integrate.odeint，把二阶方程转化为状态空间形式：

from scipy.integrate import odeint def modal_sdof(omega_i, zeta_i, F_i, t, dt): """求解单阶模态的位移时程。F_i 是等时间步长的离散模态力。""" n_steps = len(t) def rhs(y, t_float): # 用最近邻时间点索引对应力值 idx = int(round(t_float / dt)) idx = np.clip(idx, 0, n_steps - 1) x, v = y a = F_i[idx] - 2 * zeta_i * omega_i * v - omega_i**2 * x return [v, a] y0 = [0.0, 0.0] sol = odeint(rhs, y0, t, hmax=dt) return sol[:, 0]

hmax=dt强制内部积分器不跳过用户时间点，避免模态力索引错位。这里用最近邻索引而不是插值，是因为在时程分析中 F_i 的时间间距已经足够细；如果荷载很光滑，插值也只会带来微小差别。对每一阶模态重复调用此函数，就能得到独立的 q_i(t)，天然适合并行计算。

3.3 物理响应重构与速度/加速度叠加

得到模态位移时程后，物理位移直接做矩阵乘法：

Phi_keep = phi[:, :n_modes_keep] q_keep = q_all[:n_modes_keep, :] u_physical = Phi_keep @ q_keep

其中n_modes_keep是预留模态数，u_physical的每一行对应一个自由度，每一列对应一个时间步。矩阵乘法的代价可忽略，真正的计算量集中在单自由度积分器上。

速度和加速度也可以用同样方式叠加，但截断误差随导数阶数递增。用模态位移法直接叠加加速度，误差会放大到一个不可接受的程度，因为被截断的高频模态在物理位移坐标中比重小，在加速度坐标中比重大。下表给出不同响应类型的误差随保留模态数增加的趋势：

所以工程项目里，位移、速度常用模态位移法，加速度或内力时最好切换到模态加速度法。

4. 三层剪切结构算例：模态截断对动力响应的影响

4.1 模型参数与模态分析

用一个三层剪切型框架验证完整流程。每层质量 m = 1000 kg，每层层间刚度 k = 1e6 N/m，底部固定。基底激励为 0.2g、频率 1.5 Hz 的正弦加速度，持时 2 秒。Rayleigh 阻尼用第1阶和第3阶目标阻尼比 0.02 反算系数。

m = 1000.0 k = 1.0e6 M = np.eye(3) * m K = np.array([[2*k, -k, 0.0], [-k, 2*k, -k], [0.0, -k, k]]) omega2, phi = eigh(K, M) omega = np.sqrt(omega2) freq = omega / (2*np.pi) print("频率 Hz:", freq)

运行后的固有频率与参与系数如下：

第一阶累计有效质量比只有 0.785，不足以精确描述基底剪力；加上第二阶后达到 0.987，此时位移响应已经能用两阶模态覆盖。第三阶对位移贡献很小，但对层加速度可能仍有明显作用。

4.2 基底激励下的模态力与响应求解

基底加速度激励可以转化为等效荷载 F_eff = -M r a_g(t)，r 是单位向量。模态力为 F_i(t) = -Γ_i a_g(t)，负号来自运动方程的惯性项。

Amp = 0.2 * 9.8 dt = 0.005 t = np.arange(0, 2.0 + dt, dt) def ag(t_val): return Amp * np.sin(2 * np.pi * 1.5 * t_val) r = np.ones(3) gamma = phi.T @ M @ r # 模态力矩阵，shape (n_modes, n_steps) F_modal = -np.outer(gamma, ag(t))

这里的gamma与振型符号一致，因此模态力符号会自动匹配。接着计算每阶实际阻尼比，并逐阶积分。

# Rayleigh 系数 alpha, beta zeta1, zeta3 = 0.02, 0.02 A = np.array([[1/omega[0], omega[0]], [1/omega[2], omega[2]]]) alpha, beta = np.linalg.solve(A, [2*zeta1, 2*zeta3]) zeta_i = 0.5 * (alpha / omega + beta * omega) q_all = np.zeros((3, len(t))) for i in range(3): q_all[i, :] = modal_sdof(omega[i], zeta_i[i], F_modal[i, :], t, dt)

注意zeta_i中第2阶的实际阻尼比并非 0.02，这是 Rayleigh 阻尼的正常表现；只要两阶目标频率阻尼比正确，所有中间阶阻尼比都在合理范围内。

4.3 截断误差量化与直接积分对比

把两阶截断结果与全模态结果对比：

# 只保留前两阶 u3_two = phi[2, :2] @ q_all[:2, :] u3_three = phi[2, :] @ q_all error = np.max(np.abs(u3_two - u3_three)) / np.max(np.abs(u3_three)) * 100 print("顶层位移峰值误差: %.2f%%" % error)

不同保留模态数的顶层位移峰值误差如下：

加入第二阶后位移误差降到 0.4%，第三阶对位移几乎没有影响。但如果把输出换成层间加速度，第三阶贡献会升到 10% 以上。因此，截断阶数必须按输出量决定，而不是笼统取“前三阶”。

5. 模态加速度修正与数值验证的“最后一公里”

5.1 先验证，再谈优化

写完振型叠加脚本后，我习惯先用完整的 Newmark 直接积分结果做一次对拍。如果两条位移时程曲线在激励持续期间几乎重合，说明截断和积分误差都控制在合理范围；如果只有开始阶段重合，后面逐渐漂移，优先怀疑阻尼矩阵符号、模态力符号或Rayleigh系数算错。

在代码里加一个正交性检查非常简单：

print(phi.T @ M @ phi) print(phi.T @ K @ phi)

只要非对角元在 1e-10 量级以下，就可以继续排查其他环节。参与曲线的累计有效质量比是否达到 90% 只是第一道门槛，换输出量时必须重新评估。

5.2 模态加速度法：用静力修正补偿截断部分

需要计算加速度或内力时，直接用模态位移叠加会放大高频误差。常见做法是在位移叠加结果上增加静力修正项：

u_corr = K^{-1} (F_eff - Σ_{保留} φ_i Γ_i F_i(t))

修正项把截断模态的静力解补回来，忽略其动态惯性影响，对被截断的高频模态是合理近似。代码实现如下：

K_inv = np.linalg.inv(K) F_eff = -M @ r * ag(t) # shape (3, n_steps) F_modal_all = phi.T @ F_eff # 全部模态力 # 保留前两阶，修正因第3阶截断产生的静力误差 u_corr = K_inv @ (F_eff - phi[:, :2] @ F_modal_all[:2, :]) u_corrected = phi[:, :2] @ q_all[:2, :] + u_corr print("顶层加速度修正后峰值:", np.max(np.abs(u_corrected[2, :])))

u_corr在每个时间步都被重新计算，等效荷载较小的时候修正项也小，等效荷载达到峰值时修正项同步增大。它不会改变动力相位，但会把截断模态的准静态贡献补回来，是最划算的误差修正手段。

5.3 工程判据与最终建议

给保留模态数设定双阈值：累计有效质量比≥90%，且保留模态的最高频率不低于激励最高频率的2倍。两者都满足时，振型叠加结果才值得信任。时间步长建议取保留模态最高周期的1/50，比只按信号最高频率加密一倍，加速度峰值的稳定性会明显提升。把K逆的静力修正量折算到每个时步，你会发现截断模态的贡献并没有消失，而是变成了一张随荷载幅值变化的静力云图——这正好解释了为什么模态加速度法能显著改善加速度误差。

本文还有配套的精品资源，点击获取