简介面向结构动力学与非线性地震响应分析这是一份基于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_{n1}] 内对加速度的变化方式做某种假设从而把 t_{n1} 时刻的位移、速度用当前时刻的状态量表达出来。具体假设写成数学形式是u_{n1} u_n Δt·ů_n Δt²·[(1/2 - β)·ü_n β·ü_{n1}] ů_{n1} ů_n Δt·[(1 - γ)·ü_n γ·ü_{n1}]这里 β 决定位移更新式中终点加速度的权重γ 决定速度更新式中终点加速度的权重。β 0 时相当于用起点加速度外推γ 0 时速度更新完全不参考终点加速度。实际工程中默认的组合是 β 1/4、γ 1/2称为平均加速度法它假设加速度在一个时间步内取起点和终点的平均值。另一种常见的 β 1/6、γ 1/2 是线性加速度法假设加速度在步内线性变化。两者的差别在高频响应上平均加速度法对高频分量有轻微的能量耗散线性加速度法对周期偏短的振型更敏感容易在非线性计算中激发出伪振荡。2.2 无条件稳定边界与工程默认参数Newmark-Beta 法稳定性的结论是基于线性单自由度体系特征值分析得到的。当 γ ≥ 1/2 且 β ≥ (γ 1/2)²/4 时算法对任意时间步长都是无条件稳定的不满足这个条件时只有 Δt 小于体系自振周期的一定比例才能保稳定。方法名称βγ稳定性精度平均加速度法1/41/2无条件稳定二阶线性加速度法1/61/2条件稳定Δt ≤ 0.551T二阶Fox-Goodwin 法1/121/2条件稳定Δt ≤ 0.389T二阶中心差分法前身01/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_{n1} - 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(i1) - 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(i1) u(i) du; v(i1) v(i) dv; a(i1) 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/41/2最短荷载周期的 1/20 以内无条件稳定推荐默认冲击荷载1/61/2荷载脉宽的 1/20 以内精度高于平均加速度但需检查稳定性正弦稳态激励1/41/2激励周期的 1/50减小数值阻尼对幅值的影响非线性滞回1/41/2最高有效频率周期的 1/20建议再减半因为切线刚度变化提示无条件稳定只保证不发散不保证不衰减。平均加速度法在低频段幅值误差小但在高频段会引入人为阻尼如果关心高频响应细节应把 γ 设为 0.5 并减小 Δt而不是调 β。3.3 用自由振动解析解验证 newmark 求解器写完求解器后第一件事不是直接算地震波而是用无阻尼自由振动检查代码。取 m 1k (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_{n1} 的最终值。常见做法是在每个时间步内做 Newton-Raphson 迭代把非线性平衡方程逐步线性化。迭代思路是先假设一个位移增量 Δu用当前恢复力模型算出对应的恢复力 f_s再检查 t_{n1} 时刻的平衡残差R p_{n1} - (m·ü_{n1} c·ů_{n1} f_s(u_{n1}))如果 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; % 注意第一轮迭代 du0后续在下方更新 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(i1) - (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(i1) zIter dzdu * du; fs(i1) alpha*k0*u(i1) (1-alpha)*k0*z(i1); u(i1) u(i) du; v(i1) v(i) dv; a(i1) a(i) da; end end这段代码在每次迭代里同时做两件事用当前 du 计算恢复力和切线刚度再用平衡残差修正 du。其中 signDu 计算加载方向dzdu 是 Bouc-Wen 模型对位移的导数它出现在切线刚度中保证 Newton-Raphson 迭代具有二阶收敛速度。收敛判据用的是位移增量修正量的相对值而不是力的残差——在刚度很小比如接近零切线刚度时力残差判据容易误判收敛。4.3 Bouc-Wen 滞回模型的参数说明Bouc-Wen 模型有两个 β、γ和 Newmark 的 β、γ 同名但完全无关这是代码移植时最容易出错的点。为避免混淆Bouc-Wen 参数在代码里命名为 betaBw 和 gammaBw。四个核心参数对滞回形状的影响如下参数典型取值对滞回曲线的作用A1.0控制滞回环整体幅值对应屈服位移的倒数betaBw0.5控制滞回环的胖瘦增大则耗能增加gammaBw0.5控制软化或硬化趋势与 betaBw 共同决定滞回形状nBw1 ~ 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(i1) - u(i); Whys(i1) Whys(i) (fs(i) fs(i1)) / 2 * du - ... 0.5 * alpha * k0 * (u(i1)^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 的滞回图若滞回圈出现锯齿或负斜率异常尖峰优先怀疑迭代收敛容差放太宽或状态变量更新顺序写反而不是急着调模型参数。本文还有配套的精品资源点击获取
非线性动力响应计算：Newmark-Beta法MATLAB实现与收敛策略
简介面向结构动力学与非线性地震响应分析这是一份基于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_{n1}] 内对加速度的变化方式做某种假设从而把 t_{n1} 时刻的位移、速度用当前时刻的状态量表达出来。具体假设写成数学形式是u_{n1} u_n Δt·ů_n Δt²·[(1/2 - β)·ü_n β·ü_{n1}] ů_{n1} ů_n Δt·[(1 - γ)·ü_n γ·ü_{n1}]这里 β 决定位移更新式中终点加速度的权重γ 决定速度更新式中终点加速度的权重。β 0 时相当于用起点加速度外推γ 0 时速度更新完全不参考终点加速度。实际工程中默认的组合是 β 1/4、γ 1/2称为平均加速度法它假设加速度在一个时间步内取起点和终点的平均值。另一种常见的 β 1/6、γ 1/2 是线性加速度法假设加速度在步内线性变化。两者的差别在高频响应上平均加速度法对高频分量有轻微的能量耗散线性加速度法对周期偏短的振型更敏感容易在非线性计算中激发出伪振荡。2.2 无条件稳定边界与工程默认参数Newmark-Beta 法稳定性的结论是基于线性单自由度体系特征值分析得到的。当 γ ≥ 1/2 且 β ≥ (γ 1/2)²/4 时算法对任意时间步长都是无条件稳定的不满足这个条件时只有 Δt 小于体系自振周期的一定比例才能保稳定。方法名称βγ稳定性精度平均加速度法1/41/2无条件稳定二阶线性加速度法1/61/2条件稳定Δt ≤ 0.551T二阶Fox-Goodwin 法1/121/2条件稳定Δt ≤ 0.389T二阶中心差分法前身01/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_{n1} - 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(i1) - 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(i1) u(i) du; v(i1) v(i) dv; a(i1) 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/41/2最短荷载周期的 1/20 以内无条件稳定推荐默认冲击荷载1/61/2荷载脉宽的 1/20 以内精度高于平均加速度但需检查稳定性正弦稳态激励1/41/2激励周期的 1/50减小数值阻尼对幅值的影响非线性滞回1/41/2最高有效频率周期的 1/20建议再减半因为切线刚度变化提示无条件稳定只保证不发散不保证不衰减。平均加速度法在低频段幅值误差小但在高频段会引入人为阻尼如果关心高频响应细节应把 γ 设为 0.5 并减小 Δt而不是调 β。3.3 用自由振动解析解验证 newmark 求解器写完求解器后第一件事不是直接算地震波而是用无阻尼自由振动检查代码。取 m 1k (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_{n1} 的最终值。常见做法是在每个时间步内做 Newton-Raphson 迭代把非线性平衡方程逐步线性化。迭代思路是先假设一个位移增量 Δu用当前恢复力模型算出对应的恢复力 f_s再检查 t_{n1} 时刻的平衡残差R p_{n1} - (m·ü_{n1} c·ů_{n1} f_s(u_{n1}))如果 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; % 注意第一轮迭代 du0后续在下方更新 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(i1) - (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(i1) zIter dzdu * du; fs(i1) alpha*k0*u(i1) (1-alpha)*k0*z(i1); u(i1) u(i) du; v(i1) v(i) dv; a(i1) a(i) da; end end这段代码在每次迭代里同时做两件事用当前 du 计算恢复力和切线刚度再用平衡残差修正 du。其中 signDu 计算加载方向dzdu 是 Bouc-Wen 模型对位移的导数它出现在切线刚度中保证 Newton-Raphson 迭代具有二阶收敛速度。收敛判据用的是位移增量修正量的相对值而不是力的残差——在刚度很小比如接近零切线刚度时力残差判据容易误判收敛。4.3 Bouc-Wen 滞回模型的参数说明Bouc-Wen 模型有两个 β、γ和 Newmark 的 β、γ 同名但完全无关这是代码移植时最容易出错的点。为避免混淆Bouc-Wen 参数在代码里命名为 betaBw 和 gammaBw。四个核心参数对滞回形状的影响如下参数典型取值对滞回曲线的作用A1.0控制滞回环整体幅值对应屈服位移的倒数betaBw0.5控制滞回环的胖瘦增大则耗能增加gammaBw0.5控制软化或硬化趋势与 betaBw 共同决定滞回形状nBw1 ~ 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(i1) - u(i); Whys(i1) Whys(i) (fs(i) fs(i1)) / 2 * du - ... 0.5 * alpha * k0 * (u(i1)^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 的滞回图若滞回圈出现锯齿或负斜率异常尖峰优先怀疑迭代收敛容差放太宽或状态变量更新顺序写反而不是急着调模型参数。本文还有配套的精品资源点击获取
