计算有阻尼的斜拉桥在激励力作用下振动固有频率和固有模态的matlab代码

以下是一个简单的 MATLAB 代码，用于计算有阻尼的斜拉桥的固有频率和固有模态。该代码中使用了数值计算方法，即有限元分析方法。请注意，该代码仅供参考，具体实现可能会因实际情况而异。

% 定义斜拉桥的基本参数
L = 100;        % 桥梁长度
h = 10;         % 桥梁高度
d = 2;          % 桥塔距离
E = 2e11;       % 弹性模量
rho = 7850;     % 材料密度
A = 1;          % 截面面积
I = 1;          % 截面惯性矩
L1 = L/2;       % 桥梁长度的一半
L2 = L1;        % 桥梁长度的一半

% 定义有限元分析的参数
n = 10;         % 有限元数量
m = 3;          % 模态数量
dt = 0.01;      % 时间步长
t = 0:dt:10;    % 时间序列

% 定义初始条件
u0 = zeros(2*n, 1);     % 位移向量
v0 = zeros(2*n, 1);     % 速度向量
u0(2:2:end) = -h;       % 斜拉索的初始长度

% 构建刚度矩阵和质量矩阵
K = zeros(2*n, 2*n);
M = zeros(2*n, 2*n);

for i = 1:n
    % 计算每个有限元的长度和角度
    x1 = (i-1) * L1/n;
    y1 = h - (i-1) * h/n;
    x2 = i * L1/n;
    y2 = h - i * h/n;
    L_i = sqrt((x2-x1)^2 + (y2-y1)^2);
    theta_i = atan((y2-y1)/(x2-x1));
    
    % 计算每个有限元的刚度矩阵和质量矩阵
    k_i = E*I/L_i^3 * [12*sin(theta_i)^2, 6*L_i*sin(theta_i), -12*sin(theta_i)^2, 6*L_i*sin(theta_i);...
                       6*L_i*sin(theta_i), 4*L_i^2, -6*L_i*sin(theta_i), 2*L_i^2;...
                       -12*sin(theta_i)^2, -6*L_i*sin(theta_i), 12*sin(theta_i)^2, -6*L_i*sin(theta_i);...
                       6*L_i*sin(theta_i), 2*L_i^2, -6*L_i*sin(theta_i), 4*L_i^2];
    m_i = rho*A*L_i/6 * [2, 0, 1, 0;...
                         0, 2, 0, 1;...
                         1, 0, 2, 0;...
                         0, 1, 0, 2];
                     
    % 将每个有限元的刚度矩阵和质量矩阵加入总的刚度矩阵和质量矩阵中
    K(2*i-1:2*i+2, 2*i-1:2*i+2) = K(2*i-1:2*i+2, 2*i-1:2*i+2) + k_i;
    M(2*i-1:2*i+2, 2*i-1:2*i+2) = M(2*i-1:2*i+2, 2*i-1:2*i+2) + m_i;
end

% 添加约束条件
K_reduced = K(3:end-2, 3:end-2);
M_reduced = M(3:end-2, 3:end-2);

% 计算固有频率和固有模态
[V, D] = eigs(K_reduced, M_reduced, m, 'smallestabs');
f = sqrt(diag(D))/(2*pi);
u = zeros(2*n, m);
for i = 1:m
    u(:, i) = [0; V(:, i); 0];
end

% 添加阻尼
zeta = 0.05;    % 阻尼比
C = 2*zeta*sqrt(E*I/rho/A/L^3)*M;   % 阻尼矩阵
F = zeros(2*n, length(t));  % 激励力序列
F(2*n-1, :) = sin(2*pi*10*t);  % 在第一个斜拉索上施加振动力

% 使用 Newmark 方法进行时间历程分析
u(:, 1) = u0;
v = v0;
for i = 2:length(t)
    a = inv(M_reduced + (1/2)*dt*C + (1/4)*dt^2*K_reduced) * (F(3:end-2, i-1) - C*v - K_reduced*u(3:end-2, i-1));
    u(3:end-2, i) = u(3:end-2, i-1) + dt*v(3:end-2) + (1/2)*dt^2*a;
    v(3:end-2) = v(3:end-2) + (1/2)*dt*(a + inv(M_reduced + (1/2)*dt*C + (1/4)*dt^2*K_reduced)*(F(3:end-2, i) - C*v - K_reduced*u(3:end-2, i)));
end

% 绘制振动响应
figure
plot(t, u(2*n-1, :))
xlabel('Time (s)')
ylabel('Displacement (m)')
title('Response of the first cable to a sinusoidal force')

在上述代码中，斜拉桥的刚度矩阵和质量矩阵是通过有限元分析方法计算得到的，而固有频率和固有模态是通过对刚度矩阵和质量矩阵进行特征值分解得到的。在添加阻尼后，使用 Newmark 方法进行时间历程分析，从而得到斜拉桥在激励力作用下的振动响应。