考虑扭转耦联的振型分解反应谱法代码实现

%% 一层框架扭转惯量计算，单位制为(N, m, kg, s)clear;close all;clc;%% 1 参数设置Lx = 6.0;      % 框架X方向跨度  Ly = 4.0;      % 框架Y方向跨度  h = 3.0;       % 柱高  b_beam = 0.3;   % 梁宽度  h_beam = 0.5;   % 梁高度  b_column = 0.4; % 柱宽度  h_column = 0.4; % 柱高度  t_slab = 0.12;  % 楼板厚度  rho_concrete = 2500; % 混凝土密度 %% 2 计算各部分质量和转动惯量% 2.1 楼板质量和转动惯量slab_area = Lx * Ly;           % 楼板面积 slab_volume = slab_area * t_slab; % 楼板体积 slab_mass = slab_volume * rho_concrete; % 楼板质量 % 楼板绕质心的转动惯量 slab_Iz = slab_mass * (Lx^2 + Ly^2) / 12; % 绕z轴转动惯量 % 2.2 梁的质量和转动惯量% X方向梁x_beam_volume = b_beam * h_beam * Ly;x_beam_mass_per = x_beam_volume * rho_concrete; x_beam_count = 2;           % Y方向梁       y_beam_volume = b_beam * h_beam * Lx; y_beam_mass_per = y_beam_volume * rho_concrete; y_beam_count = 2;            % 梁的转动惯量计算x_beam_Iz_per = x_beam_mass_per * (Ly^2) / 12; y_beam_Iz_per = y_beam_mass_per * (Lx^2) / 12; % 2.3 柱的质量和转动惯量column_volume = b_column * h_column * h; column_mass_per = column_volume * rho_concrete; column_count = 4;% 柱的转动惯量计算column_positions = [0, 0; Lx, 0; 0, Ly; Lx, Ly]; % 四根柱的位置坐标center = [Lx/2, Ly/2]; % 质心坐标column_Iz_total = 0;for i = 1:column_count    % 计算第i根柱到质心的距离    dx = column_positions(i, 1) - center(1);    dy = column_positions(i, 2) - center(2);    r = sqrt(dx^2 + dy^2);       column_Iz_total = column_Iz_total + column_mass_per * r^2;end%% 3 总扭转惯量计算% 总转动惯量（绕质心）total_Iz = slab_Iz;% 添加梁的转动惯量for i = 1:x_beam_count    beam_x = (i-1) * Lx;    dx = beam_x - center(1);    total_Iz = total_Iz + x_beam_Iz_per + x_beam_mass_per * dx^2;endfor i = 1:y_beam_count    beam_y = (i-1) * Ly;    dy = beam_y - center(2);    total_Iz = total_Iz + y_beam_Iz_per + y_beam_mass_per * dy^2;end% 添加柱的转动惯量total_Iz = total_Iz + column_Iz_total;%% 4 结果输出fprintf('该层框架的扭转惯量为: %.2f kg·m²\n', total_Iz);

%% 考虑扭转耦联的振型分解反应谱法，单位制为(N,m,kg,s)clear;clc;close all;%% 1输入%% 1.1 输入结构信息% 每层质量（kg）m = 1e7 * [1, 1.1, 1.2, 1.1, 1.3, 0.9]; % 每层转动惯量（kg·m²）J = 1e6 * [1, 1.1, 1.2, 1.1, 1.3, 0.9];% 每层X方向刚度（N/m）k_xx = 1e9 * [1.1, 1.2, 1.3, 1.1, 1.2, 0.85]; % 每层Y方向刚度（N/m）k_yy = 1e9 * [1.0, 1.1, 1.2, 1.0, 1.1, 0.8]; % 每层扭转刚度（N·m/rad）k_xy = 1e9 * [0.5, 0.6, 0.7, 0.6, 0.7, 0.4]; % 每层高度（m）h = [3.3, 3.0, 3.0, 3.0, 3.0, 2.8];% 每层X方向偏心距（m）e_x = [0.1, 0.15, 0.15, 0.15, 0.15, 0.1];% 每层Y方向偏心距（m）e_y = [0.1, 0.15, 0.15, 0.15, 0.15, 0.1];% 层数n = length(m);%% 1.2 输入地震相关信息% 特征周期(s)Tg = 0.40;% 地震影响系数最大值a_max = 0.08;% 阻尼比znb = 0.05;%% 2计算%% 2.1 求解质量矩阵和刚度矩阵% 质量矩阵M = zeros(3*n, 3*n);for i = 1:n    M(3*i-2, 3*i-2) = m(i);        M(3*i-1, 3*i-1) = m(i);       M(3*i, 3*i) = J(i);         end% 刚度矩阵K = zeros(3*n, 3*n);for i = 1:n    if i == 1        K(1,1) = k_xx(1) + k_xx(2);        K(1,4) = -k_xx(2);    elseif i == n        K(3*n-2, 3*n-5) = -k_xx(n);        K(3*n-2, 3*n-2) = k_xx(n);    else        K(3*i-2, 3*i-5) = -k_xx(i);        K(3*i-2, 3*i-2) = k_xx(i) + k_xx(i+1);        K(3*i-2, 3*i+1) = -k_xx(i+1);    endendfor i = 1:n    if i == 1        K(2,2) = k_yy(1) + k_yy(2);        K(2,5) = -k_yy(2);    elseif i == n        K(3*n-1, 3*n-4) = -k_yy(n);        K(3*n-1, 3*n-1) = k_yy(n);    else        K(3*i-1, 3*i-4) = -k_yy(i);        K(3*i-1, 3*i-1) = k_yy(i) + k_yy(i+1);        K(3*i-1, 3*i+2) = -k_yy(i+1);    endendfor i = 1:n    k_xx_i = k_xx(i);    k_yy_i = k_yy(i);    k_xy_i = k_xy(i);    e_x_i = e_x(i);    e_y_i = e_y(i);    if i == 1        K(3,3) = k_xy(1) + k_xy(2) + k_xx_i*e_y_i^2 + k_yy_i*e_x_i^2;        K(3,6) = -k_xy(2);    elseif i == n        K(3*n, 3*n-3) = -k_xy(n);        K(3*n, 3*n) = k_xy(n) + k_xx_i*e_y_i^2 + k_yy_i*e_x_i^2;    else        K(3*i, 3*i-3) = -k_xy(i);        K(3*i, 3*i) = k_xy(i) + k_xy(i+1) + k_xx_i*e_y_i^2 + k_yy_i*e_x_i^2;        K(3*i, 3*i+3) = -k_xy(i+1);    end    K(3*i-2, 3*i) = -k_xx_i*e_y_i;    K(3*i, 3*i-2) = -k_xx_i*e_y_i;    K(3*i-1, 3*i) = k_yy_i*e_x_i;    K(3*i, 3*i-1) = k_yy_i*e_x_i;end%% 2.2 求解特征值和特征向量[V, D] = eig(K, M); % 得到圆频率pinlv = sqrt(abs(diag(D)));  % 得到周期T = 2*pi./pinlv;% 按周期从大到小排序[T, index] = sort(T, 'descend'); % 根据周期，排序频率和振型pinlv = pinlv(index);V = V(:, index);% 取前6个振型num_modes = min(6, length(T));T = T(1:num_modes);pinlv = pinlv(1:num_modes);V = V(:, 1:num_modes); %% 2.3计算各振型地震影响系数% aa：地震影响系数列矩阵,《高规》4.3.8aa = zeros(num_modes, 1);for i = 1:num_modes    if T(i) < 0.1        aa(i) = (0.45 + 5.5 * T(i)) * a_max;    elseif T(i) <= Tg        aa(i) = a_max;    elseif T(i) <= 5*Tg        aa(i) = (Tg / T(i))^0.9 * a_max;    else        aa(i) = (0.2^0.9 - 0.02*(T(i) - 5*Tg)) * a_max;    endend%% 2.4 计算振型参与系数% GM:振型参与系数GM = zeros(3, num_modes);for i = 1:num_modes    V_x = V(1:3:end, i); % x方向振型    V_y = V(2:3:end, i); % y方向振型      V_theta = V(3:3:end, i); % 扭转振型    % 计算各方向振型参与系数    GM(1, i) = (V_x' * m(:)) / (V(:, i)' * M * V(:, i));    GM(2, i) = (V_y' * m(:)) / (V(:, i)' * M * V(:, i));    GM(3, i) = (V_theta' * J(:)) / (V(:, i)' * M * V(:, i));end%% 2.5 计算CQC耦联系数function rho = calculate_cqc_coefficients(T, znb)    n_modes = length(T);    rho = zeros(n_modes, n_modes);        for i = 1:n_modes        for j = 1:n_modes            if T(i) == 0 || T(j) == 0                beta_ij = 0;            else                beta_ij = T(j)/T(i);            end                        rho(i,j) = (8*znb^2*(1+beta_ij)*beta_ij^(1.5)) / ...                ((1-beta_ij^2)^2 + 4*znb^2*beta_ij*(1+beta_ij)^2);        end    endend%% 2.6各振型地震作用标准值计算% F：各振型各层地震作用标准值矩阵F = zeros(3*n, num_modes);          for i = 1:num_modes    % 计算各方向地震作用    F_x = aa(i) * GM(1, i) * V(1:3:end, i) .* m(:);    F_y = aa(i) * GM(2, i) * V(2:3:end, i) .* m(:);    F_theta = aa(i) * GM(3, i) * V(3:3:end, i) .* J(:);    % 组合地震作用矩阵    F(1:3:end, i) = F_x;    F(2:3:end, i) = F_y;    F(3:3:end, i) = F_theta;end%% 2.7 CQC组合计算function response = cqc_combination(F, rho)    n_modes = size(F, 2);    n_points = size(F, 1);    response = zeros(n_points, 1);    for k = 1:n_points        temp = 0;        for i = 1:n_modes            for j = 1:n_modes                temp = temp + F(k,i) * rho(i,j) * F(k,j);            end        end        response(k) = sqrt(temp);    endend% 计算CQC耦联系数rho = calculate_cqc_coefficients(T, znb);% 采用CQC组合来计算每层地震作用标准值V_CQC_x = cqc_combination(F(1:3:end, :), rho);V_CQC_y = cqc_combination(F(2:3:end, :), rho);V_CQC_theta = cqc_combination(F(3:3:end, :), rho);%% 3输出%% 3.1周期与频率fprintf('各振型周期与频率\n');for i = 1:num_modes    fprintf('第%d阶: 周期=%.2fs, 频率=%.2fHz\n', i, T(i), pinlv(i)/(2*pi));end%% 3.2地震作用标准值fprintf('\n各层地震作用标准值（CQC组合）\n');fprintf('楼层\tX方向(kN)\tY方向(kN)\t扭转(kN·m)\n');for i = 1:n    fprintf('%d\t%.2f\t\t%.2f\t\t%.2f\n', i, V_CQC_x(i)/1000, ...        V_CQC_y(i)/1000, V_CQC_theta(i)/1000);end