modal_eom_solve.m

关于高层建筑强风条件下风压系数计算

function [UUU,VVV,V_dot] = modal_eom_solve(h,omega,Modal_LoadM,time)

% The n x n matrix 'h' defines the modal damping ratios:
%   h = diag( [zeta_1 zeta_2 ... zeta_n] ) 
% The n x n matrix 'omega' defines the natural frequencies
%   omega = diag( [omega_1 omega_2 ... omega_n] )
% The n x s matrix 'Modal_LoadM' defines the time series of generalized modal forces
%   Modal_LoadM  = [ g_1(t_1) g_1(t_2) ... g_1(t_s); ...
%                    g_2(t_1) g_2(t_2) ... g_2(t_s); ...
%                                          ...                   ...
%                    g_n(t_1) g_n(t_2) ... g_n(t_s)]

n           = length(omega);            % number of modes
s           = length(Modal_LoadM);      % number of samples
T           = max(time)-min(time);      % record duration
dt          = T/(s-1);                  % time step -- MUST BE CONSTANT

% Initialize modal displacement, velocity, and acceleration matrices:
UUU         = zeros(n,s);
VVV         = zeros(n,s);
V_dot       = zeros(n,s);

for iii = 1:n
    wn      = omega(iii,iii);           % natural frequency
    zeta    = h(iii,iii);               % damping ratio
    wd      = wn*sqrt(1-zeta^2);        % damped frequency
    E       = exp(-zeta.*wn.*dt);
    K       = dt.*wd;
    C       = E.*cos(K);
    S       = E.*sin(K);
    Sp      = S./K;

    % Calculate filter coefficients
    B0      = (2*zeta*(C-1)+(2*zeta^2-1)*S/sqrt(1-zeta^2)+dt*wn)/(dt*wn);
    B1      = (-2*C*dt*wn+2*zeta*(1-E^2)-2*(2*zeta^2-1)*S/sqrt(1-zeta^2))/(dt*wn);
    B2      = (E^2*(dt*wn+2*zeta)-2*zeta*C+(2*zeta^2-1)*S/sqrt(1-zeta^2))/(dt*wn);
    A0      = 1;                        % A0 = constant = 1, see help filter
    A1      = -2.*C;
    A2      = E.^2;

    UUU(iii,:)   = filter([B0 B1 B2],[A0 A1 A2],Modal_LoadM(iii,:))/wn^2; % modal displacements
    
    % Calculate filter coefficients for modal velocities:
    B0      = (-dt*wn*Sp+2*zeta*(1-C)+(1-2*zeta^2)*S/sqrt(1-zeta^2))/(2*dt*zeta*wn);
    B1      = (2*dt*wn*Sp-2*zeta*(1-E^2)-2*(1-2*zeta^2)*S/sqrt(1-zeta^2))/(2*dt*zeta*wn);
    B2      = (-dt*wn*Sp-2*zeta*E^2+2*zeta*C+(1-2*zeta^2)*S/sqrt(1-zeta^2))/(2*dt*zeta*wn);
    
    % Compute modal velocities
    VVV(iii,:)   = filter([B0 B1 B2],[A0 A1 A2], Modal_LoadM(iii,:))/wn; 
    % Obtain modal accelerations from equation of motion
    V_dot(iii,:) = Modal_LoadM(iii,:) - 2*zeta*wn*VVV(iii,:) - wn^2*UUU(iii,:);
end