基于Newmark-β法的单自由度体系地震响应MATLAB实现

基于Newmark-β法的单自由度体系地震响应MATLAB实现

---

一、算法原理与公式

单自由度体系运动方程：

Newmark-β法递推公式：

1. 预测步：

   

2. 校正步：

   

---

二、MATLAB核心代码实现

function [time, disp, vel, acc] = Newmark_SDOF(m, c, k, ag, dt, beta, gamma, u0, v0)% 参数说明：% m: 质量 (kg)% c: 阻尼系数 (N·s/m)% k: 刚度 (N/m)% ag: 地震加速度时程 (m/s²)% dt: 时间步长 (s)% beta, gamma: Newmark参数% u0, v0: 初始位移和速度n = length(ag);time = 0:dt:(n-1)*dt;% 初始条件u = zeros(n,1); v = zeros(n,1); a = zeros(n,1);u(1) = u0; v(1) = v0;a(1) = (-m*ag(1) - c*v0 - k*u0)/m;% 迭代计算for i = 1:n-1% 预测步a_pred = a(i);v_pred = v(i) + dt*(1-gamma)*a_pred;u_pred = u(i) + dt*v_pred + 0.5*dt^2*(1-2*beta)*a_pred;% 校正步a(i+1) = (-m*ag(i+1) - c*v_pred - k*u_pred)/m;v(i+1) = v_pred + dt*( (1-gamma)*a_pred + gamma*a(i+1) );u(i+1) = u_pred + dt*v_pred + 0.5*dt^2*( (1-2*beta)*a_pred + 2*beta*a(i+1) );end% 输出结果disp = u;vel = v;acc = a;
end

---

三、参数设置与优化

1. 典型参数组合

参数类型	推荐值范围	物理意义
β	0.20-0.30	位移滞后因子
γ	0.50-0.60	速度滞后因子
时间步长	T/10 ≤ dt ≤ T/5	T为结构自振周期 (s)

2. 阻尼矩阵计算（瑞利阻尼）

function c = rayleigh_damping(m, k, xi, omega1, omega2)% 输入参数：% xi: 阻尼比% omega1, omega2: 特征频率 (rad/s)c = xi*m*2*omega1*omega2/(omega1+omega2);
end

---

四、测试案例（El Centro地震波）

% 参数设置
m = 1000;       % 质量 (kg)
k = 25000;      % 刚度 (N/m) → 自振周期 T=0.63s
xi = 0.05;      % 阻尼比
dt = 0.01;      % 时间步长
beta = 0.25;    % Newmark-β参数
gamma = 0.5;    % Newmark-γ参数% 加载地震波（示例：El Centro NS分量）
ag = load('elcentro.txt');  % 假设数据格式为时间加速度
ag = ag(:,2);               % 提取加速度列% 运行算法
[time, disp, vel, acc] = Newmark_SDOF(m, rayleigh_damping(m,k,xi,0.1,10), k, ag, dt, beta, gamma, 0, 0);% 可视化
figure;
subplot(3,1,1);
plot(time, ag*9.81, 'r', time, disp*1000, 'b--');
title('加速度时程对比 (m/s² vs mm)');
xlabel('时间 (s)'); ylabel('幅值');
legend('地震输入', '结构响应');subplot(3,1,2);
plot(time, vel*100, 'g', time, acc, 'm:');
title('速度与加速度响应');
xlabel('时间 (s)'); ylabel('幅值');subplot(3,1,3);
plot(disp, acc*1000, 'ko-');
title('位移-加速度滞回曲线');
xlabel('位移 (mm)'); ylabel('加速度 (m/s²)');

参考代码 采用Newmark法计算单自由度体系在地震作用下的响应 www.youwenfan.com/contentcnl/79193.html

五、结果分析要点

1. 频谱特性验证：

   [Pxx,f] = pwelch(disp,1024,512,1024,1/dt);
   figure;
   plot(f,10*log10(Pxx));
   title('位移响应功率谱密度');
   xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB/Hz)');
   

2. 关键指标计算： 最大位移：max(abs(disp)) 等效阻尼比：通过能量等效法计算 周期偏移：对比初始周期与响应周期