计算有阻尼的斜拉桥振动固有频率和固有模态的matlab代码

计算有阻尼的斜拉桥振动固有频率和固有模态需要用到特征值问题的求解方法，下面是一个MATLAB代码示例，仅供参考。

% 计算有阻尼的斜拉桥振动固有频率和固有模态
clear all; clc;

% 桥梁参数
L = 200; % 桥梁长度
H = 50; % 桥塔高度
h = 5; % 斜拉索高度
theta = pi/6; % 斜拉角度
E = 200e9; % 弹性模量
rho = 7800; % 密度
A = pi*(0.05)^2; % 斜拉索截面积
m = rho*A*L/2; % 桥梁质量
k = E*A/L; % 斜拉索刚度
zeta = 0.02; % 阻尼比

% 构造系统矩阵
K = [2*k*cos(theta), -k*cos(theta), 0, 0; ...
    -k*cos(theta), 2*k*cos(theta), -k*cos(theta), 0; ...
    0, -k*cos(theta), 2*k*cos(theta), -k*cos(theta); ...
    0, 0, -k*cos(theta), k*cos(theta)];
M = [m/2, 0, 0, 0; 0, m/2, 0, 0; 0, 0, m/2, 0; 0, 0, 0, m/2];
C = zeta*M + 0.5*K; % 阻尼矩阵

% 求解特征值问题
[phi, omega2] = eig(K, M);

% 将特征向量单位化
for i = 1:4
    phi(:, i) = phi(:, i)/sqrt(phi(:, i)'*M*phi(:, i));
end

% 计算有阻尼的固有频率
omega = sqrt(diag(omega2));
damping_ratio = diag(phi'*C*phi)./(2*omega); % 阻尼比

% 绘制固有模态图
figure;
x = [0, L/2*cos(theta), L/2, L/2*cos(theta), 0, -L/2*cos(theta), -L/2, -L/2*cos(theta), 0];
y = [0, -L/2*sin(theta), -H, L/2*sin(theta), 0, L/2*sin(theta), H, -L/2*sin(theta), 0];
for i = 1:4
    subplot(2, 2, i);
    plot(x, y, 'r-', 'LineWidth', 2);
    hold on;
    mode_shape = [phi(1:2, i)', 0, phi(3:4, i)', 0, -phi(1:2, i)', 0, -phi(3:4, i)', 0];
    fill(x, y, mode_shape);
    axis equal;
    xlabel('x (m)');
    ylabel('y (m)');
    title(['Mode ', num2str(i), ', \omega = ', num2str(omega(i)), ', \zeta = ', num2str(damping_ratio(i))]);
end

上述代码中，构造了系统矩阵K、质量矩阵M和阻尼矩阵C，并使用eig函数求解特征值问题。通过归一化特征向量，可以得到固有模态，并通过计算特征向量对应的阻尼比，得到有阻尼的固有频率。最后，通过绘制固有模态图，可以直观地了解斜拉桥的振动模态和固有频率。需要注意的是，在实际应用中，还需要考虑更多的因素，如桥梁的几何形状、材料参数等，需要更加复杂的模型和算法。