大型土木工程结构风致位移响应与湍流估计（振动理论 + MATLAB）

原创于 2026-06-22 11:08:34 发布 · 26 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#matlab #开发语言

核心逻辑：把风看成随机激励，用模态叠加法（振型分解）把大型结构降维，再基于随机振动理论计算位移 RMS，并反过来用响应估计湍流强度。

一、物理模型与假设

结构：大型土木工程（高层/桥梁/塔架）
→ 多自由度线性系统，质量矩阵 M，刚度矩阵 K，阻尼矩阵 C（Rayleigh 阻尼）。
风荷载：
- 平均风：静力荷载 $Fmean=12ρCDAUˉ2F_{\text{mean}} = \frac{1}{2}\rho C_D A \bar{U}^2$
- 脉动风（湍流）：零均值平稳随机过程，功率谱用 Davenport 谱。
响应：
位移响应 $y (x, t)$ 是随机过程 → 计算 RMS 位移 $σy\sigma_y$ （工程设计最关心）。
湍流估计：
已知响应 RMS → 反推湍流强度 $I_u$ 。

二、结构建模（多自由度剪切层模型）

2.1 质量、刚度矩阵（高层建筑示例）

%% 结构参数
H = 300;          % 建筑高度 (m)
N = 10;           % 层数（模态截断）
m = 5e6;          % 每层质量 (kg)
k = 2e9;          % 层间刚度 (N/m)

% 质量矩阵 M
M = m * eye(N);

% 刚度矩阵 K（剪切层）
K = zeros(N);
for i = 1:N
    if i == 1
        K(i,i) = k;
    else
        K(i,i) = 2*k;
    end
    if i < N
        K(i,i+1) = -k;
        K(i+1,i) = -k;
    end
end

% Rayleigh 阻尼 C = αM + βK
xi = 0.02;        % 阻尼比 2%
[V, D] = eig(K, M);
fn = sqrt(diag(D))/(2*pi);   % 频率 (Hz)
alpha = zeros(N,1); beta = zeros(N,1);
% 取前两阶模态计算 α,β
fn1 = fn(1); fn2 = fn(2);
alpha = (2*xi*(fn1*fn2)/(fn1+fn2));
beta  = (2*xi)/(fn1+fn2);
C = alpha*M + beta*K;

三、风荷载建模（Davenport 谱 + 相干性）

3.1 Davenport 脉动风速谱

$S_u(f) = \frac{4 K \bar{U}^2 x_1^2}{(1+x_1^2)^{4/3}}, \quad x_1 = \frac{1200 f}{\bar{U}}$

function Su = davenport_spectrum(f, U_bar, K)
% f: 频率 (Hz)
% U_bar: 平均风速 (m/s)
% K: 地面粗糙度系数 (0.003~0.03)
x1 = 1200*f/U_bar;
Su = 4*K*U_bar^2 * x1.^2 ./ (1+x1.^2).^(4/3);
end

3.2 空间相干性（考虑不同楼层风的相关性）

function coh = spatial_coherence(f, dz, U_bar, Lc)
% dz: 两点竖向距离
% Lc: 积分尺度 (~50~100m)
coh = exp(-0.5 * f * dz / U_bar);
end

四、模态叠加法

4.1 模态坐标下的风荷载

%% 模态分析
[Phi, Lambda] = eig(K, M);
fn = sqrt(diag(Lambda))/(2*pi);   % 频率 (Hz)
Phi = Phi ./ vecnorm(Phi);        % 振型归一化

% 模态质量、刚度、阻尼
M_n = diag(Phi'*M*Phi);
K_n = diag(Phi'*K*Phi);
C_n = diag(Phi'*C*Phi);

%% 风荷载（模态投影）
Nf = 512; df = 0.01;
f = (0:Nf-1)*df;
U_bar = 30;        % 10m 高度平均风速 (m/s)
K_rough = 0.01;    % 粗糙度
A = 50;            % 迎风面积 (m^2)
Cd = 1.2;          % 阻力系数
rho = 1.25;        % 空气密度

% 脉动风谱
Su = davenport_spectrum(f, U_bar, K_rough);

% 模态风荷载谱矩阵 S_Fn
Nm = 5;            % 取前5阶模态
S_Fn = zeros(Nm, Nm, Nf);
for i = 1:Nm
    for j = 1:Nm
        % 振型投影
        phi_i = Phi(:,i);
        phi_j = Phi(:,j);
        % 简化：假设各层风荷载完全相关
        F_mag = 0.5*rho*Cd*A*U_bar^2;
        Fi = F_mag * phi_i(1);   % 基底模态力
        Fj = F_mag * phi_j(1);
        for ff = 1:Nf
            S_Fn(i,j,ff) = Fi * Fj * Su(ff);
        end
    end
end

五、位移响应 RMS 计算

5.1 模态响应谱

$S_{q_i}(f) = \frac{S_{F_i}(f)}{| - (2\pi f)^2 M_i + K_i + j 2\pi f C_i |^2}$

S_qn = zeros(Nm, Nf);
for i = 1:Nm
    for ff = 1:Nf
        w = 2*pi*f(ff);
        H = -w^2*M_n(i) + K_n(i) + 1i*w*C_n(i);
        S_qn(i,ff) = S_Fn(i,i,ff) / abs(H)^2;
    end
end

% 模态位移 RMS
sigma_q = sqrt(sum(S_qn, 2) * df);

5.2 物理坐标位移 RMS

% 顶层位移 RMS
sigma_y_top = sqrt(sum((Phi(end,:).*sigma_q').^2));
fprintf('顶层位移 RMS = %.4f m\n', sigma_y_top);

六、湍流强度估计

已知响应 RMS $σy\sigma_y$ ，反推湍流强度 $I_u$ ：

$\sigma_y^2 = \sum_{i=1}^{N_m} \phi_i^2 \int_{-\infty}^{\infty} |H_i(j\omega)|^2 S_{F_i}(\omega) d\omega$

其中 $SFi∝Iu2S_{F_i} \propto I_u^2$ ，因此：

$I_u = \sqrt{\frac{\sigma_y^2}{\text{常数}}}$

%% 湍流强度反演（简化）
% 假设：响应主要由第一阶模态控制
phi1 = Phi(end,1);
H1 = 1 / (K_n(1) - (2*pi*fn(1))^2 + 1i*2*pi*fn(1)*C_n(1));
% 谱积分近似
integral_Su = trapz(f, davenport_spectrum(f, U_bar, K_rough));
% 风荷载方差
sigma_F1_sq = (0.5*rho*Cd*A*U_bar^2)^2 * integral_Su * Iu^2;
% 响应方差
sigma_y_sq = phi1^2 * abs(H1)^2 * sigma_F1_sq * df;
% 反推 Iu
Iu_est = sqrt(sigma_y_top^2 / (phi1^2 * abs(H1)^2 * sigma_F1_sq * df));
fprintf('估计湍流强度 Iu = %.4f\n', Iu_est);

七、完整仿真主程序

%% main_wind_response.m
clear; clc; close all;

% 1. 结构建模
H = 300; N = 10; m = 5e6; k = 2e9;
M = m*eye(N);
K = zeros(N);
for i = 1:N
    K(i,i) = 2*k*(i<N) + k*(i==1);
    if i < N
        K(i,i+1) = -k; K(i+1,i) = -k;
    end
end
xi = 0.02; [V,D]=eig(K,M); fn=sqrt(diag(D))/(2*pi);
alpha = 2*xi*fn(1)*fn(2)/(fn(1)+fn(2));
beta = 2*xi/(fn(1)+fn(2));
C = alpha*M + beta*K;

% 2. 模态分析
[Phi, Lambda] = eig(K, M);
fn = sqrt(diag(Lambda))/(2*pi);
[fn, idx] = sort(fn);
Phi = Phi(:,idx);
M_n = diag(Phi'*M*Phi);
K_n = diag(Phi'*K*Phi);
C_n = diag(Phi'*C*Phi);

% 3. 风荷载
U_bar = 30; rho = 1.25; Cd = 1.2; A = 50;
K_rough = 0.01;
f = (0:511)*0.01;
Su = davenport_spectrum(f, U_bar, K_rough);

% 4. 模态响应
Nm = 5; S_qn = zeros(Nm, length(f));
for i = 1:Nm
    for ff = 1:length(f)
        w = 2*pi*f(ff);
        H = -w^2*M_n(i) + K_n(i) + 1i*w*C_n(i);
        F_mag = 0.5*rho*Cd*A*U_bar^2 * Phi(1,i);
        S_F = F_mag^2 * Su(ff);
        S_qn(i,ff) = S_F / abs(H)^2;
    end
end

% 5. 位移 RMS
sigma_q = sqrt(sum(S_qn,2)*0.01);
sigma_y_top = sqrt(sum((Phi(end,1:Nm).*sigma_q').^2));
fprintf('顶层位移 RMS = %.4f m\n', sigma_y_top);

% 6. 湍流估计
Iu = sqrt(sigma_y_top^2 / (Phi(end,1)^2 * (1/K_n(1)^2) * trapz(f,Su)*0.01));
fprintf('估计湍流强度 Iu = %.4f\n', Iu);

% 7. 绘图
figure('Color','w');
plot(fn(1:Nm), sigma_q, 'o-', 'LineWidth', 1.5);
xlabel('模态频率 (Hz)'); ylabel('模态位移 RMS (m)');
title('各阶模态位移响应'); grid on;