本程序模拟人在板上行走时，板产生的动力响应。采用移动荷载模型 + 模态叠加法，这是建筑结构中人行振动分析的经典方法。

一、理论模型

1.1 板的振动方程

$$m\ddot{w}+c\dot{w}+kw=p(x,t)$$

其中：

• $w$ = 板挠度
• $m$ = 模态质量
• $c$ = 模态阻尼
• $k$ = 模态刚度
• $p(x,t)$ = 人行走产生的移动荷载

1.2 人行荷载模型

$$p(t)=G\left[1+\sum _{i=1}^3{\alpha }_i\sin (2\pi i{f}_{s}t+{\varphi }_i)\right]$$

• $G$ = 人的重量（≈700 N）
• $f_s$ = 步行频率（≈2 Hz）
• ${\alpha }_i$ = 谐波系数

---

二、完整 MATLAB 程序

2.1 主程序 (human_walking_plate_vibration.m)

%% 人行走作用下板的振动响应仿真
% 作者：结构动力学工程师
% 日期：2024年

clear; clc; close all;

fprintf('=== 人行走作用下板的振动仿真 ===\n');

%% 1. 板参数（混凝土楼板）
plate = struct();
plate.E = 3.0e10;          % 弹性模量 (Pa)
plate.nu = 0.2;            % 泊松比
plate.rho = 2500;          % 密度 (kg/m³)
plate.thickness = 0.15;     % 板厚 (m)
plate.length = 6.0;        % 板长 (m)
plate.width = 4.0;         % 板宽 (m)

% 计算弯曲刚度
plate.D = plate.E * plate.thickness^3 / (12*(1-plate.nu^2));

fprintf('板参数:\n');
fprintf('  尺寸: %.1f m × %.1f m × %.3f m\n', ...
        plate.length, plate.width, plate.thickness);
fprintf('  弯曲刚度: %.2e N·m\n', plate.D);

%% 2. 人行走参数
human = struct();
human.weight = 700;         % 体重 (N) ≈ 70 kg
human.step_freq = 2.0;      % 步行频率 (Hz)
human.step_length = 0.75;   % 步长 (m)
human.walking_speed = human.step_freq * human.step_length; % 行走速度 (m/s)

% 动力荷载系数（ISO 10137标准）
human.alpha1 = 0.4;        % 第一谐波系数
human.alpha2 = 0.1;        % 第二谐波系数
human.alpha3 = 0.05;       % 第三谐波系数

fprintf('\n人行走参数:\n');
fprintf('  体重: %.1f N (%.1f kg)\n', human.weight, human.weight/9.81);
fprintf('  步行频率: %.1f Hz\n', human.step_freq);
fprintf('  行走速度: %.2f m/s\n', human.walking_speed);

%% 3. 模态分析（四边简支板）
fprintf('\n计算板的前几阶模态...\n');
modes = plate_modal_analysis(plate);

%% 4. 时间参数
T_total = plate.length / human.walking_speed + 2;  % 总时长
dt = 0.001;                                        % 时间步长
t = 0:dt:T_total;
N_steps = length(t);

fprintf('仿真时长: %.1f s, 时间步长: %.3f s\n', T_total, dt);

%% 5. 计算移动荷载时程
fprintf('计算人行荷载时程...\n');
[P_total, P_positions] = walking_load_time_history(human, plate, t);

%% 6. 模态叠加法计算响应
fprintf('模态叠加法计算板振动响应...\n');
displacement = modal_superposition(plate, modes, human, P_total, P_positions, t);

%% 7. 结果可视化
fprintf('绘制结果...\n');
plot_results(plate, human, modes, displacement, t, P_total);

%% 8. 舒适度评估
comfort_assessment(displacement, human.step_freq, plate);

2.2 板模态分析 (plate_modal_analysis.m)

function modes = plate_modal_analysis(plate)
% 四边简支矩形板的模态分析
% 模态频率: ωmn = π²√(D/ρh) * √((m/a)² + (n/b)²)

a = plate.length;       % 板长
b = plate.width;        % 板宽
h = plate.thickness;    % 板厚
rho = plate.rho;        % 密度
D = plate.D;            % 弯曲刚度

% 计算前9阶模态 (m,n = 1,2,3)
m_vals = [1 1 1 2 2 2 3 3 3];
n_vals = [1 2 3 1 2 3 1 2 3];

frequencies = zeros(9,1);
modal_mass = zeros(9,1);
modal_stiffness = zeros(9,1);
modal_damping = zeros(9,1);

fprintf('模态序号   m   n   频率(Hz)  阻尼比\n');
fprintf('-----------------------------------\n');

for i = 1:9
    m = m_vals(i);
    n = n_vals(i);
    
    % 固有频率 (rad/s)
    omega = (pi^2) * sqrt(D/(rho*h)) * sqrt((m/a)^2 + (n/b)^2);
    frequencies(i) = omega / (2*pi);  % Hz
    
    % 模态质量（单位模态质量归一化）
    modal_mass(i) = rho * h * a * b / 4;
    
    % 模态刚度
    modal_stiffness(i) = omega^2 * modal_mass(i);
    
    % 模态阻尼（结构阻尼比2%）
    modal_damping(i) = 0.02;
    
    fprintf('  %d        %d   %d    %6.2f    %5.1f%%\n', ...
            i, m, n, frequencies(i), modal_damping(i)*100);
end

% 返回模态参数
modes.frequencies = frequencies;
modes.modal_mass = modal_mass;
modes.modal_stiffness = modal_stiffness;
modes.modal_damping = modal_damping;
modes.m_vals = m_vals;
modes.n_vals = n_vals;
modes.a = a;
modes.b = b;
end

2.3 人行荷载时程 (walking_load_time_history.m)

function [P_total, P_positions] = walking_load_time_history(human, plate, t)
% 计算人行走产生的移动荷载时程

N_steps = length(t);
P_total = zeros(N_steps,1);
P_positions = zeros(N_steps,2);  % [x, y] 位置

% 起始位置（板的一端）
x_start = 0.1;  % 距边缘0.1m
y_center = plate.width/2;

for i = 1:N_steps
    current_time = t(i);
    
    % 当前位置
    x_pos = x_start + human.walking_speed * current_time;
    y_pos = y_center;
    
    % 如果走到板末端，停止
    if x_pos > plate.length - 0.1
        x_pos = plate.length - 0.1;
    end
    
    P_positions(i,:) = [x_pos, y_pos];
    
    % 计算当前时刻的荷载值
    fs = human.step_freq;
    
    % 傅里叶级数形式的步行荷载
    P = human.weight * (1 + ...
        human.alpha1 * sin(2*pi*fs*current_time) + ...
        human.alpha2 * sin(2*pi*2*fs*current_time) + ...
        human.alpha3 * sin(2*pi*3*fs*current_time));
    
    P_total(i) = P;
end
end

2.4 模态叠加法 (modal_superposition.m)

function displacement = modal_superposition(plate, modes, human, P_total, P_positions, t)
% 使用模态叠加法计算板振动响应

N_steps = length(t);
dt = t(2) - t(1);
N_modes = length(modes.frequencies);

% 初始化模态坐标响应
q_response = zeros(N_modes, N_steps);

% 模态形状函数（简支板）
m_vals = modes.m_vals;
n_vals = modes.n_vals;
a = plate.length;
b = plate.width;

fprintf('模态叠加法计算进度:\n');
progress_step = ceil(N_modes/10);

for mode = 1:N_modes
    if mod(mode, progress_step) == 0
        fprintf('  计算第 %d/%d 阶模态...\n', mode, N_modes);
    end
    
    m = m_vals(mode);
    n = n_vals(mode);
    
    % 模态阻尼比
    zeta = modes.modal_damping(mode);
    
    % 模态圆频率
    omega_n = 2*pi*modes.frequencies(mode);
    
    % 模态刚度
    k_m = modes.modal_stiffness(mode);
    
    % 模态质量
    m_m = modes.modal_mass(mode);
    
    % 模态阻尼
    c_m = 2*zeta*sqrt(k_m*m_m);
    
    % 计算该模态的广义荷载
    Q_generalized = zeros(N_steps,1);
    
    for i = 1:N_steps
        % 荷载位置
        x_load = P_positions(i,1);
        y_load = P_positions(i,2);
        
        % 模态形状在该点的值
        phi = sin(m*pi*x_load/a) * sin(n*pi*y_load/b);
        
        % 广义荷载
        Q_generalized(i) = P_total(i) * phi;
    end
    
    % 单自由度系统响应（Duhamel积分）
    q_mode = sdof_response(Q_generalized, m_m, c_m, k_m, dt);
    
    q_response(mode,:) = q_mode;
end

% 合成总响应
fprintf('合成各阶模态响应...\n');
displacement = zeros(N_steps,1);

for i = 1:N_steps
    total_disp = 0;
    for mode = 1:N_modes
        m = m_vals(mode);
        n = n_vals(mode);
        
        % 中心点模态位移
        x_center = a/2;
        y_center = b/2;
        phi_center = sin(m*pi*x_center/a) * sin(n*pi*y_center/b);
        
        total_disp = total_disp + q_response(mode,i) * phi_center;
    end
    displacement(i) = total_disp;
end

fprintf('模态叠加法计算完成\n');
end

function response = sdof_response(F, m, c, k, dt)
% 单自由度系统Duhamel积分
N = length(F);
response = zeros(N,1);

% 系统参数
omega_n = sqrt(k/m);
zeta = c/(2*sqrt(k*m));
omega_d = omega_n * sqrt(1-zeta^2);

% 初始条件
x0 = 0; v0 = 0;

% 递推计算
response(1) = x0;

for i = 2:N
    % 时间步
    dt_small = dt;
    
    % 状态转移矩阵
    A = [1, dt_small; 
         -omega_n^2*dt_small, 1-2*zeta*omega_n*dt_small];
    
    % 激励项
    B = [0; dt_small/m];
    
    % 状态向量
    state = [response(i-1); (response(i-1)-response(max(1,i-2)))/dt_small];
    
    % 下一时刻状态
    state_new = A*state + B*F(i);
    response(i) = state_new(1);
end
end

2.5 结果可视化 (plot_results.m)

function plot_results(plate, human, modes, displacement, t, P_total)
% 绘制分析结果

figure('Position', [100, 100, 1400, 900]);

% 1. 板中心位移时程
subplot(3,3,1);
plot(t, displacement*1000, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('位移 (mm)');
title('板中心振动位移时程');
grid on; box on;

% 标注峰值
[max_disp, max_idx] = max(abs(displacement));
hold on;
plot(t(max_idx), max_disp*1000, 'ro', 'MarkerSize', 8);
text(t(max_idx), max_disp*1000, sprintf('峰值: %.3f mm', max_disp*1000), ...
     'VerticalAlignment', 'bottom');

% 2. 人行荷载时程
subplot(3,3,2);
plot(t, P_total, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('荷载 (N)');
title('人行荷载时程');
grid on; box on;

% 3. 模态频率分布
subplot(3,3,3);
bar(modes.frequencies(1:6));
xlabel('模态序号');
ylabel('频率 (Hz)');
title('前6阶模态频率');
grid on; box on;

% 4. 位移频谱
subplot(3,3,4);
N = length(displacement);
fft_result = fft(displacement);
freq = (0:N-1)*(1/(t(2)-t(1)))/N;
magnitude = abs(fft_result)/N;

plot(freq(1:round(N/2)), magnitude(1:round(N/2))*1000, 'g-', 'LineWidth', 1.5);
xlabel('频率 (Hz)');
ylabel('位移幅值 (mm)');
title('位移频谱');
grid on; box on;

% 标注步行频率
hold on;
plot([human.step_freq, human.step_freq], [0, max(magnitude)*1000], 'r--', 'LineWidth', 1);
text(human.step_freq, max(magnitude)*1000*0.8, sprintf('步行频率\n%.1f Hz', human.step_freq));

% 5. 板振动动画（简化）
subplot(3,3,5);
[x_grid, y_grid] = meshgrid(linspace(0,plate.length,50), linspace(0,plate.width,30));
z_grid = zeros(size(x_grid));

% 显示当前时刻的变形（取中间时刻）
mid_idx = round(length(t)/2);
z_grid = displacement(mid_idx) * sin(pi*x_grid/plate.length) .* sin(pi*y_grid/plate.width);

surf(x_grid, y_grid, z_grid*1000);
xlabel('长度 (m)');
ylabel('宽度 (m)');
zlabel('位移 (mm)');
title(sprintf('板变形 (t=%.1f s)', t(mid_idx)));
view(45, 45);
colorbar;

% 6. 舒适度指标
subplot(3,3,6);
% 计算加速度
dt = t(2)-t(1);
velocity = gradient(displacement, dt);
acceleration = gradient(velocity, dt);

plot(t, acceleration, 'm-', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('加速度 (m/s²)');
title('板中心加速度时程');
grid on; box on;

% 7. RMS响应
subplot(3,3,7);
rms_disp = sqrt(mean(displacement.^2));
rms_acc = sqrt(mean(acceleration.^2));

bar([rms_disp*1000, rms_acc]);
set(gca, 'XTickLabel', {'RMS位移(mm)', 'RMS加速度(m/s²)'});
ylabel('数值');
title('RMS响应指标');
grid on; box on;

% 8. 峰值因子
subplot(3,3,8);
peak_factor_disp = max(abs(displacement)) / rms_disp;
peak_factor_acc = max(abs(acceleration)) / rms_acc;

bar([peak_factor_disp, peak_factor_acc]);
set(gca, 'XTickLabel', {'位移峰值因子', '加速度峰值因子'});
ylabel('峰值因子');
title('峰值因子分析');
grid on; box on;

% 9. 参数汇总
subplot(3,3,9);
axis off;
text(0.1, 0.8, sprintf('板尺寸: %.1f×%.1f×%.3f m', ...
    plate.length, plate.width, plate.thickness), 'FontSize', 10);
text(0.1, 0.7, sprintf('板频率: %.2f Hz', modes.frequencies(1)), 'FontSize', 10);
text(0.1, 0.6, sprintf('步行频率: %.1f Hz', human.step_freq), 'FontSize', 10);
text(0.1, 0.5, sprintf('峰值位移: %.3f mm', max(abs(displacement))*1000), 'FontSize', 10);
text(0.1, 0.4, sprintf('RMS加速度: %.3f m/s²', rms_acc), 'FontSize', 10);

% 舒适度评价
if rms_acc < 0.005
    comfort = '优秀';
elseif rms_acc < 0.015
    comfort = '良好';
elseif rms_acc < 0.05
    comfort = '可接受';
else
    comfort = '不舒适';
end
text(0.1, 0.3, sprintf('舒适度评价: %s', comfort), 'FontSize', 12, 'FontWeight', 'bold');

title('分析参数汇总');

sgtitle('人行走作用下板振动响应分析', 'FontSize', 16, 'FontWeight', 'bold');
end

2.6 舒适度评估 (comfort_assessment.m)

function comfort_assessment(displacement, step_freq, plate)
% 舒适度评估（基于国际标准）

fprintf('\n=== 舒适度评估 ===\n');

% 计算RMS加速度
dt = 0.001; % 假设时间步长
velocity = gradient(displacement, dt);
acceleration = gradient(velocity, dt);
rms_acc = sqrt(mean(acceleration.^2));

% ISO 2631-2 标准
fprintf('ISO 2631-2 标准评估:\n');
fprintf('  RMS加速度: %.4f m/s²\n', rms_acc);

if rms_acc < 0.005
    fprintf('  评价: 无感知 (优秀)\n');
elseif rms_acc < 0.015
    fprintf('  评价: 轻微感知 (良好)\n');
elseif rms_acc < 0.05
    fprintf('  评价: 明显感知 (可接受)\n');
else
    fprintf('  评价: 强烈感知 (需改进)\n');
end

% AISC 钢框架建筑标准
fprintf('\nAISC 标准评估:\n');
peak_disp = max(abs(displacement))*1000; % mm
if peak_disp < 0.5
    fprintf('  峰值位移: %.3f mm < 0.5 mm ✓ 满足\n', peak_disp);
else
    fprintf('  峰值位移: %.3f mm > 0.5 mm ✗ 不满足\n', peak_disp);
end

% 共振检查
fprintf('\n共振风险检查:\n');
fundamental_freq = 1/(2*pi)*sqrt(plate.E*plate.thickness^3/(12*plate.rho*plate.thickness));
fprintf('  板基频: %.2f Hz\n', fundamental_freq);
fprintf('  步行频率: %.2f Hz\n', step_freq);

freq_ratio = step_freq / fundamental_freq;
if abs(freq_ratio - 1) < 0.1
    fprintf('  ⚠ 警告: 步行频率接近板基频，存在共振风险！\n');
else
    fprintf('  ✓ 无共振风险\n');
end
end

参考代码 用MATLAB程序表示人行走下，板的振动情况 www.youwenfan.com/contentcsv/80860.html

三、使用说明

3.1 运行程序

% 直接运行主程序
human_walking_plate_vibration

3.2 参数调整

% 修改板参数
plate.thickness = 0.20;     % 增加板厚
plate.length = 8.0;        % 增加板长

% 修改人行走参数
human.weight = 800;        % 增加人体重
human.step_freq = 2.5;     % 增加步行频率

% 重新运行分析

3.3 结果解读

指标	标准	说明
RMS加速度	<0.005 m/s²	无感知
峰值位移	<0.5 mm	结构安全要求
共振比	≠1±0.1	避免共振

四、工程应用建议

1. 减振措施：当RMS加速度>0.015 m/s²时，考虑增加板厚或设置TMD
2. 频率调整：避免板基频与步行频率重合（2±0.2 Hz）
3. 人群效应：多人行走时，荷载应乘以1.85的人群折减系数