MATLAB 地震动非线性时程分析程序

MATLAB 地震动非线性时程分析程序，面向单自由度（SDOF）弹塑性结构，采用Newmark‑β 隐式积分 + Bouc–Wen 光滑滞回模型，这是结构抗震分析中最常用、数值稳定性最好的组合之一。

---

一、理论模型

1、运动方程（非线性）

\[m\ddot{u}+c\dot{u}+f_s(u,\dot{u})= -m\ddot{u}_g(t) \]

2、Bouc–Wen 滞回恢复力

\[f_s = \alpha k u + (1-\alpha)k z \]

\[\dot{z}=A\dot{u}-\beta|\dot{u}||z|^{n-1}z-\gamma\dot{u}|z|^n \]

• (\(z\))：滞回位移
• (\(\alpha\))：屈服后刚度比
• (\(A,\beta,\gamma,n\))：控制滞回形状

---

二、主程序：非线性时程分析（SDOF）

%% Earthquake Nonlinear Time History Analysis (SDOF)
% Bouc-Wen hysteretic model + Newmark-beta
clear; clc; close all;

%% -------------------- 结构参数 --------------------
m  = 1000;          % 质量 (kg)
k  = 2e5;           % 初始刚度 (N/m)
zeta = 0.05;        % 阻尼比
c  = 2*zeta*sqrt(m*k);  % 阻尼 (Ns/m)

%% -------------------- Bouc-Wen 参数 ----------------
alpha = 0.1;        % 屈服后刚度比
A = 1.0;
beta = 0.5;
gamma = 0.5;
n = 2;              % 控制平滑程度

%% -------------------- 地震动 -----------------------
% 读取地震波（加速度，单位：m/s^2）
load('elcentro.mat');   % 变量名应为 acc, dt
% 若没有 elcentro.mat，可用内置示例
if ~exist('acc','var')
    dt = 0.02;
    t = 0:dt:30;
    acc = 0.3*sin(2*pi*1.5*t).*exp(-0.1*t);
else
    t = (0:length(acc)-1)*dt;
end

N = length(acc);

%% -------------------- Newmark-beta 参数 -----------
beta_n = 0.25;      % Newmark-β
gamma_n = 0.5;      % 无条件稳定

%% -------------------- 初始化 -----------------------
u   = zeros(N,1);   % 位移
ud  = zeros(N,1);   % 速度
udd = zeros(N,1);   % 加速度
z   = zeros(N,1);   % Bouc-Wen 滞回变量

% 初始加速度
udd(1) = -acc(1);

%% -------------------- 时程积分 ---------------------
for i = 1:N-1
    % ---- 当前恢复力 ----
    fs = alpha*k*u(i) + (1-alpha)*k*z(i);

    % ---- Newmark 预测 ----
    u_pred = u(i) + dt*ud(i) + ...
             dt^2*(0.5-beta_n)*udd(i);
    ud_pred = ud(i) + dt*(1-gamma_n)*udd(i);

    % ---- 等效刚度 ----
    dfdz = (1-alpha)*k;
    K_eq = k + dfdz * ...
           (A - beta*abs(ud_pred)*abs(z(i))^(n-1)*z(i) ...
              - gamma*ud_pred*abs(z(i))^n);

    % ---- 有效刚度 ----
    K_eff = m/(beta_n*dt^2) + gamma_n*c/(beta_n*dt) + K_eq;

    % ---- 有效荷载 ----
    P_eff = -m*acc(i+1) ...
            - c*(ud_pred + gamma_n*dt*udd(i)/(beta_n)) ...
            - m*(u_pred/dt^2 + ud_pred/dt + (0.5-beta_n)*udd(i)) ...
            + fs;

    % ---- 位移修正 ----
    du = P_eff / K_eff;
    u(i+1) = u_pred + du;

    % ---- 速度、加速度 ----
    ud(i+1) = ud_pred + gamma_n*du/(beta_n*dt);
    udd(i+1) = (du/(beta_n*dt^2)) - ud(i)/dt ...
               - (0.5-beta_n)*udd(i)/beta_n;

    % ---- Bouc-Wen 滞回变量 ----
    dz = A*ud(i+1) ...
         - beta*abs(ud(i+1))*abs(z(i))^(n-1)*z(i) ...
         - gamma*ud(i+1)*abs(z(i))^n;
    z(i+1) = z(i) + dz*dt;
end

%% -------------------- 结果计算 ---------------------
fs_all = alpha*k*u + (1-alpha)*k*z;
base_shear = fs_all;
drift_ratio = u / (k/m * 9.81);   % 近似层间位移角

%% -------------------- 绘图 ------------------------
figure('Color','white','Position',[100 100 1000 700])

subplot(3,2,1)
plot(t, acc,'k')
xlabel('Time (s)'); ylabel('Acc (m/s^2)')
title('Ground Motion')
grid on

subplot(3,2,2)
plot(t, u*100,'b')
xlabel('Time (s)'); ylabel('Disp (cm)')
title('Displacement')
grid on

subplot(3,2,3)
plot(t, ud,'r')
xlabel('Time (s)'); ylabel('Vel (m/s)')
title('Velocity')
grid on

subplot(3,2,4)
plot(t, base_shear,'m')
xlabel('Time (s)'); ylabel('Base Shear (N)')
title('Restoring Force')
grid on

subplot(3,2,5)
plot(u, fs_all,'k','LineWidth',1.5)
xlabel('Displacement (m)')
ylabel('Restoring Force (N)')
title('Hysteresis Loop')
grid on
axis equal

subplot(3,2,6)
plot(t, drift_ratio,'g')
xlabel('Time (s)')
ylabel('Drift Ratio')
title('Interstory Drift Ratio')
grid on

sgtitle('Nonlinear Seismic Time History Analysis (Bouc-Wen Model)')

---

三、双线性滞回模型

如果你只想用理想弹塑性，替换 Bouc–Wen 为：

fy = 5000;          % 屈服力
k = 2e5;            % 弹性刚度
ke = 0.1*k;         % 屈服后刚度

for i = 1:N-1
    fs = min(fy, max(-fy, k*u(i))) + ke*(u(i)-sign(u(i))*fy/k);
    ...
end

---

四、多自由度（MDOF）扩展思路

运动方程：

\[\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{F}_s=-\mathbf{M}\mathbf{r}\ddot{u}_g \]

• 集中质量矩阵
• Rayleigh 阻尼
• 层间恢复力：每层一个 Bouc–Wen / 双线性模型
• 积分：Newmark‑β 对全系统联立求解

参考代码 地震动非线性时程分析程序 www.youwenfan.com/contentcnw/82474.html

五、工程实用建议

地震波

• PEER NGA-West2
• 中国规范《建筑抗震设计规范》人工波
• 统一单位：m / s / s²

数值稳定性

• beta = 0.25, gamma = 0.5（平均加速度法）
• 时间步 dt ≤ 0.02 s

验证

• 弹性阶段 → 与反应谱对比
• 能量平衡检查：

\[\int f_s du = E_{damp}+E_{kin} \]