简介：面向结构动力学与非线性地震响应分析，这是一份基于Matlab编写的非线性Newmark算法小型工具，适用于需要处理几何非线性、材料非线性或强震作用下动力时程分析的工程师与科研人员。压缩包共3个文件，包含线性和非线性二自由度体系两个计算脚本，以及一条地震加速度记录数据文件，包体仅11KB，结构紧凑、便于快速上手。目前已有601人学习/下载，说明其在相关领域具备一定参考价值。通过运行和修改源码，可以直观理解Newmark方法中β、γ参数对积分稳定性与数值耗散的影响，掌握逐时间步迭代求解非线性平衡方程、更新刚度与阻尼矩阵、以及提取位移/速度/加速度响应曲线的完整流程；同时程序的时间步选择与参数设置保留了清晰的输入接口，方便替换为其他地震波或自定义荷载，用于教学演示或初步验证。
1. 非线性动力响应计算为什么绕不开 Newmark-Beta 法
做结构地震响应或机械冲击分析的人，大概率遇到过这样的情况：从网上下载一个Newmark beta的压缩包，解压后是一个线性单自由度体系的 MATLAB 脚本，换一个非线性恢复力模型就不知道怎么改了。Newmark-Beta 法本身不复杂，复杂的是它默认的线性假设——刚度恒定、恢复力与位移成正比。实际结构进入弹塑性后，恢复力随加载历史变化，每一步的等效刚度都在变，这时如果不理解算法内部如何更新刚度、如何迭代消除不平衡力，写出来的非线性动力响应程序很容易发散或精度失控。
这篇把 Newmark-Beta 法从线性到非线性的 MATLAB 实现路径拆开讲：先解决加速度假设和稳定条件这两个理论问题，再给可复现的最小代码，最后落到滞回模型、迭代收敛和结果验证。适合土木、机械方向做动力分析的工程师和研究生阅读；已经写过线性求解器的读者，可以直接跳到第 4 章看非线性迭代的实现细节。
2. Newmark-Beta 法的数学骨架：加速度假设、稳定性边界与增量方程
2.1 加速度在时间步内的变化假设决定了 β 和 γ 的物理含义
单自由度体系在动力荷载下的运动方程可以写成：
m·ü(t) + c·ů(t) + f_s(u, ů) = p(t)其中 m 是质量，c 是阻尼系数，f_s 是恢复力。线性问题里 f_s = k·u，非线性问题里 f_s 是位移和速度的函数，甚至带有滞回记忆。Newmark-Beta 法的核心思想，是在一个时间步 [t_n, t_{n+1}] 内，对加速度的变化方式做某种假设，从而把 t_{n+1} 时刻的位移、速度用当前时刻的状态量表达出来。
具体假设写成数学形式是：
u_{n+1} = u_n + Δt·ů_n + Δt²·[(1/2 - β)·ü_n + β·ü_{n+1}] ů_{n+1} = ů_n + Δt·[(1 - γ)·ü_n + γ·ü_{n+1}]这里 β 决定位移更新式中终点加速度的权重，γ 决定速度更新式中终点加速度的权重。β = 0 时相当于用起点加速度外推，γ = 0 时速度更新完全不参考终点加速度。
实际工程中默认的组合是 β = 1/4、γ = 1/2，称为平均加速度法，它假设加速度在一个时间步内取起点和终点的平均值。另一种常见的 β = 1/6、γ = 1/2 是线性加速度法，假设加速度在步内线性变化。两者的差别在高频响应上：平均加速度法对高频分量有轻微的能量耗散，线性加速度法对周期偏短的振型更敏感，容易在非线性计算中激发出伪振荡。
2.2 无条件稳定边界与工程默认参数
Newmark-Beta 法稳定性的结论，是基于线性单自由度体系特征值分析得到的。当 γ ≥ 1/2 且 β ≥ (γ + 1/2)²/4 时，算法对任意时间步长都是无条件稳定的；不满足这个条件时，只有 Δt 小于体系自振周期的一定比例才能保稳定。
| 方法名称 | β | γ | 稳定性 | 精度 |
|---|---|---|---|---|
| 平均加速度法 | 1/4 | 1/2 | 无条件稳定 | 二阶 |
| 线性加速度法 | 1/6 | 1/2 | 条件稳定（Δt ≤ 0.551T） | 二阶 |
| Fox-Goodwin 法 | 1/12 | 1/2 | 条件稳定（Δt ≤ 0.389T） | 二阶 |
| 中心差分法（前身） | 0 | 1/2 | 条件稳定（Δt ≤ T/π） | 二阶 |
注意「无条件稳定」不等于「无条件准确」。无条件稳定只说明解不会随着时间步的推进指数式发散，但数值阻尼会衰减真实的高频响应。非线性问题中刚度不断变化，等效自振周期也在变，起始无条件稳定的参数组合在某个步内完全可能表现为条件稳定行为，所以时间步长的选取仍然必要。
2.3 非线性动力响应方程的增量改写
线性问题可以直接用全量位移求解，非线性问题最好写成增量方程。原因是恢复力 f_s(u) 的切线刚度 k_T = ∂f_s/∂u 在每个时间步内都要重新计算，而增量位移 Δu 是每次迭代的直接未知量。
把速度、加速度更新式写成增量形式：
Δü = (1/(β·Δt²))·Δu - (1/(β·Δt))·ů_n - (1/(2β))·ü_n Δů = (γ/(β·Δt))·Δu - (γ/β)·ů_n + Δt·(1 - γ/(2β))·ü_n代入运动方程后，得到等效刚度：
k̂ = k_T + (γ/(β·Δt))·c + (1/(β·Δt²))·m每一步要解的方程是：
k̂·Δu = Δp + (m/(β·Δt) + (γ/β)·c)·ů_n + (m/(2β) - Δt·(1 - γ/(2β))·c)·ü_n这个形式里 Δp = p_{n+1} - p_n，右端全部是已知量，只要 k_T 已知就能直接解出 Δu。线性问题 k_T 恒等于初始刚度 k，所以 k̂ 全时段不变；非线性问题每次迭代都要更新 k̂，这就是第 3、4 章代码实现的分水岭。
3. 用 MATLAB 实现线性 Newmark-Beta 动力响应计算
3.1 增量位移解法的 MATLAB 最小实现
线性单自由度版本的代码不复杂，重点是循环结构的清晰性。下面是一个可直接复用的函数，输入质量和刚度、阻尼系数、荷载时程，输出位移、速度、加速度时程。
function [u, v, a] = newmark_linear(m, c, k, p, dt, beta, gamma) % 输入: m 质量, c 阻尼, k 刚度, p 荷载时程, dt 时间步长, beta, gamma 为 Newmark 参数 % 输出: u 位移时程, v 速度时程, a 加速度时程 nt = length(p); u = zeros(1, nt); v = zeros(1, nt); a = zeros(1, nt); % 初始加速度由平衡方程求出，假设初始位移和速度为零 a(1) = (p(1) - c * v(1) - k * u(1)) / m; % 线性问题中等效刚度只计算一次 kHat = k + gamma / (beta * dt) * c + 1 / (beta * dt^2) * m; % 增量方程右端两个常系数矩阵 A = m / (beta * dt) + gamma / beta * c; B = m / (2 * beta) - dt * (1 - gamma / (2 * beta)) * c; for i = 1:nt-1 dp = p(i+1) - p(i); % 计算等效荷载增量 dP = dp + A * v(i) + B * a(i); % 解增量位移 du = dP / kHat; % 由增量位移反推速度和加速度增量 dv = gamma / (beta * dt) * du - gamma / beta * v(i) + dt * (1 - gamma / (2 * beta)) * a(i); da = 1 / (beta * dt^2) * du - 1 / (beta * dt) * v(i) - 1 / (2 * beta) * a(i); % 累加到当前状态 u(i+1) = u(i) + du; v(i+1) = v(i) + dv; a(i+1) = a(i) + da; end end这段代码里 kHat 是 2.3 节推导的等效刚度，A 和 B 对应增量方程右端两个状态系数项。调用时如果用平均加速度法，传入 beta = 0.25、gamma = 0.5。需要注意初始加速度不能设为 0，必须从平衡方程求，否则前几步会出现非物理的抖动。
3.2 β、γ 与 Δt 的取值参照表
参数取值的经验值，我整理成下方的表。Δt 的选取对结果的影响往往比 β、γ 更大，特别是在荷载时程是地震波或冲击脉冲时，必须保证一个荷载峰值周期内至少有 10 到 20 个计算步。
| 计算场景 | β | γ | Δt 建议 | 备注 |
|---|---|---|---|---|
| 地震响应（弹塑性） | 1/4 | 1/2 | 最短荷载周期的 1/20 以内 | 无条件稳定，推荐默认 |
| 冲击荷载 | 1/6 | 1/2 | 荷载脉宽的 1/20 以内 | 精度高于平均加速度，但需检查稳定性 |
| 正弦稳态激励 | 1/4 | 1/2 | 激励周期的 1/50 | 减小数值阻尼对幅值的影响 |
| 非线性滞回 | 1/4 | 1/2 | 最高有效频率周期的 1/20 | 建议再减半，因为切线刚度变化 |
提示：无条件稳定只保证不发散，不保证不衰减。平均加速度法在低频段幅值误差小，但在高频段会引入人为阻尼，如果关心高频响应细节，应把 γ 设为 0.5 并减小 Δt，而不是调 β。
3.3 用自由振动解析解验证 newmark 求解器
写完求解器后第一件事不是直接算地震波，而是用无阻尼自由振动检查代码。取 m = 1，k = (2π)²，初始位移 u(0) = 1，初始速度 0，理论解是 u(t) = cos(2πt)。用下面的脚本对比数值解与解析解：
m = 1; c = 0; k = (2*pi)^2; dt = 0.01; tEnd = 5; t = 0:dt:tEnd; p = zeros(size(t)); % 无外荷载 [u, v, a] = newmark_linear(m, c, k, p, dt, 0.25, 0.5); uExact = cos(2*pi*t); plot(t, u, 'b-', t, uExact, 'r--'); xlabel('时间 t'); ylabel('位移 u'); legend('Newmark 数值解', '解析解');判断标准有两条：一是 5 秒内峰值衰减不超过百分之几，二是数值解的周期不能和解析解有明显偏移。如果发现幅值持续增长，先检查 kHat 公式里的系数是不是写成 betadt^2 而不是 (betadt)^2，这是最容易写错的地方。确认线性版本无误后，再进入非线性版本。
4. 非线性 Newmark-Beta：滞回模型、Newton-Raphson 迭代与 MATLAB 代码
4.1 恢复力模型的差异决定非线性计算的成色
非线性动力响应里，恢复力 f_s(u, ů) 的形式决定了系统行为。常见的模型有三种：双线性滞回模型，屈服前刚度为 k0，屈服后刚度为 α·k0，适合钢材等弹塑性材料；Bouc-Wen 模型，用带记忆的微分方程描述光滑滞回曲线，适合混凝土和土体；还有刚度退化模型，用于循环荷载下刚度逐渐下降的结构。
双线性模型简单，但进入屈服段后切线刚度发生突变，给迭代收敛带来压力。Bouc-Wen 模型的优势是滞回曲线连续可导，切线刚度表达式是解析的，迭代稳定性更好。工程上做参数敏感性和多工况扫描时，我一般先用 Bouc-Wen 模型把算法调通，再替换成具体的材料本构。
4.2 在每一个时间步内做 Newton-Raphson 迭代
非线性问题的核心变化是：2.3 节的增量方程不再能一步解出精确的 Δu，因为 k_T 本身依赖于 u_{n+1} 的最终值。常见做法是在每个时间步内做 Newton-Raphson 迭代，把非线性平衡方程逐步线性化。
迭代思路是：先假设一个位移增量 Δu，用当前恢复力模型算出对应的恢复力 f_s，再检查 t_{n+1} 时刻的平衡残差：
R = p_{n+1} - (m·ü_{n+1} + c·ů_{n+1} + f_s(u_{n+1}))如果 R 不为零，就利用切线刚度修正 Δu，直到 R 足够小。每一步修正量 δ(Δu) = R / k̂，其中 k̂ 用当前状态的切线刚度计算，重新形成等效刚度矩阵。
下面是带 Bouc-Wen 滞回模型的单自由度非线性求解关键循环：
function [u, v, a, fs, z] = newmark_boucwen(m, c, k0, alpha, A, betaBw, gammaBw, nBw, p, dt, tEnd) nt = length(p); u = zeros(1, nt); v = zeros(1, nt); a = zeros(1, nt); fs = zeros(1, nt); z = zeros(1, nt); % z 为滞回位移 beta = 0.25; gamma = 0.5; % Newmark 参数，固定为平均加速度法 tol = 1e-6; maxIter = 200; for i = 1:nt-1 du = 0; % 时间步内累计位移增量，初始猜测为 0 zIter = z(i); fsIter = fs(i); for iter = 1:maxIter % 计算当前猜测位移下的速度、加速度 du = du; % 注意：第一轮迭代 du=0，后续在下方更新 dv = gamma/(beta*dt)*du - gamma/beta*v(i) + dt*(1-gamma/(2*beta))*a(i); da = 1/(beta*dt^2)*du - 1/(beta*dt)*v(i) - 1/(2*beta)*a(i); uGuess = u(i) + du; % Bouc-Wen 恢复力与当前切线刚度 signDu = sign(du + eps); dzdu = A - betaBw*signDu*abs(zIter)^(nBw-1)*zIter - gammaBw*abs(zIter)^nBw; fsGuess = alpha*k0*uGuess + (1-alpha)*k0*zIter; kT = alpha*k0 + (1-alpha)*k0*dzdu; % 平衡残差 R = p(i+1) - (m*da + c*dv + fsGuess); % 求解位移增量修正量 kHat = kT + gamma/(beta*dt)*c + 1/(beta*dt^2)*m; ddu = R / kHat; du = du + ddu; if abs(ddu) < tol * max(1, abs(du)) break; end end % 时间步末更新状态量 z(i+1) = zIter + dzdu * du; fs(i+1) = alpha*k0*u(i+1) + (1-alpha)*k0*z(i+1); u(i+1) = u(i) + du; v(i+1) = v(i) + dv; a(i+1) = a(i) + da; end end这段代码在每次迭代里同时做两件事：用当前 du 计算恢复力和切线刚度，再用平衡残差修正 du。其中 signDu 计算加载方向，dzdu 是 Bouc-Wen 模型对位移的导数，它出现在切线刚度中，保证 Newton-Raphson 迭代具有二阶收敛速度。收敛判据用的是位移增量修正量的相对值，而不是力的残差——在刚度很小（比如接近零切线刚度）时，力残差判据容易误判收敛。
4.3 Bouc-Wen 滞回模型的参数说明
Bouc-Wen 模型有两个 β、γ，和 Newmark 的 β、γ 同名但完全无关，这是代码移植时最容易出错的点。为避免混淆，Bouc-Wen 参数在代码里命名为 betaBw 和 gammaBw。四个核心参数对滞回形状的影响如下：
| 参数 | 典型取值 | 对滞回曲线的作用 |
|---|---|---|
| A | 1.0 | 控制滞回环整体幅值，对应屈服位移的倒数 |
| betaBw | 0.5 | 控制滞回环的胖瘦，增大则耗能增加 |
| gammaBw | 0.5 | 控制软化或硬化趋势，与 betaBw 共同决定滞回形状 |
| nBw | 1 ~ 2 | 控制屈服过渡的平滑程度，越大越接近双线性 |
当 nBw = 1 且 A = 1、betaBw = gammaBw = 0.5 时，Bouc-Wen 模型近似等于 Masing 规则下的滞回行为，适合做初始验证。注意 z 的初值必须是 0，否则滞回曲线起点会偏移，后续所有圈都跟着偏。
4.4 不收敛与发散的三种典型场景
非线性迭代最常见的失败有三种。第一种是时间步跨越滞回拐点。双线性模型屈服时刻切线刚度突变，Newton-Raphson 迭代在拐点附近震荡。解决办法是减小 Δt，或者在检测到拐点后对该步做二分加密。
第二种是负切线刚度段的发散。软化结构在峰值荷载后切线刚度为负，k̂ 可能接近零甚至为负，直接除会出现巨大增量。常见做法是改用 Modified Newton-Raphson，迭代过程中保持起始刚度不变，虽然收敛变慢但稳定。也可以对 kHat 设置下限，比如防止它小于初始刚度的 10%。
第三种是收敛判据选错。力残差判据在切线刚度极小时会把很小的位移误差误判为收敛，导致滞回曲线出现平台锯齿。我一般用位移修正量的相对值abs(ddu) < tol * max(1,abs(du))，同时保留一个最大迭代次数上限作为兜底。
5. Newmark-Beta 计算结果的验证技巧：能量平衡与高频振荡排查
5.1 能量平衡检查：判定非线性 newmark 有没有算错
非线性时程序写完后，比对着滞回曲线更重要的一件事是能量平衡。对单自由度体系，输入能应当等于动能、弹性变形能、滞回耗能和阻尼耗能之和。检查能量平衡能同时暴露两类错误：迭代不收敛导致的额外能量注入，以及状态更新时的丢步误差。
Wkin = 0.5 * m * v.^2; Wel = 0.5 * alpha * k0 * u.^2; % 弹性部分变形能 Whys = zeros(1, nt); % 滞回耗能累积 for i = 1:nt-1 du = u(i+1) - u(i); Whys(i+1) = Whys(i) + (fs(i) + fs(i+1)) / 2 * du - ... 0.5 * alpha * k0 * (u(i+1)^2 - u(i)^2); end Wdis = cumtrapz(t, c * v.^2); % 阻尼耗能 Wtotal = Wkin + Wel + Whys + Wdis; max(abs(Wtotal - cumtrapz(t, p .* v))) / max(abs(Wtotal))正常结果中最后一行给出的相对误差应小于 1%。如果误差明显偏大，优先检查位移时程末尾是否还有残余振荡，并回看滞回曲线是否出现不闭合的圈。能量检查比单纯看位移时程更严格，因为位移曲线看起来合理时能量守恒可能已经在慢慢失守。
5.2 高频振型的伪振荡与排查顺序
多自由度系统扩展之前，先确认单自由度版本没有隐藏的数值振荡。判断方法很简单：把加速度时程画出来，如果加载段之外出现周期近似等于几个 Δt 的高频抖动，就说明算法在往模型里注入虚假能量。排查顺序是先确认阻尼设置，过小的 Rayleigh 阻尼系数会让高阶振型的能量无法耗散；再检查 Δt 是不是大于最高有效频率周期的 1/10；最后看滞回模型的状态变量 z 是否在每个时间步末尾都被正确刷新。
如果高频抖动仍然存在，把 γ 从 0.5 略微上调到 0.55 会引入少量数值阻尼，能有效抑制振荡，但代价是精度降为一阶。我通常会先用两组 Δt（比如 dt 和 dt/2）各算一遍，结果差异小于 5% 才认为高频分量处理得当。在非线性滞回模型中，还要同时输出 u 和 f_s 的滞回图，若滞回圈出现锯齿或负斜率异常尖峰，优先怀疑迭代收敛容差放太宽或状态变量更新顺序写反，而不是急着调模型参数。
本文还有配套的精品资源，点击获取
