f_load_vector.m

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

function [gen_d,gen_v,gen_a,Npointsout] = F_Load_Vector(M,Fai,F_t,dt,h,omega)
global  flnTaps  Npoints  Ndofs  Ntaps  Nfloors  H_bldg  Fdofs
% This function calculates the generalized displacements, velocities and accelerations

time = linspace(0,dt*Npoints,Npoints);  % Generate a "time" vector, [Npoints x 1]

%%%%%%%%%%%%%%%%% Pressure tap information : flnTaps, Page 2, [Ntaps x 5] %%%%%%%%%%%%%%%%%%%%
% Load flnTaps and check for errors
clear fid; fid=fopen(flnTaps);
if fid == -1
    errordlg('Unable to open ''flnCp fnCp'' - Page 2')
    error   ('Unable to open ''flnCp fnCp'' - Page 2')
end
fclose(fid); clear fid;
taps=[]; load (flnTaps) 
if isempty(taps)==1
    errordlg('Variable ''taps'' not defined. See flnTaps on Page Two')
    error   ('Variable ''taps'' not defined. See flnTaps on Page Two')
end
if size(taps,1) ~= Ntaps | size(taps,2) ~= 5
    errordlg('Variable ''taps'' not properly sized (Ntaps x 5)')
    error   ('Variable ''taps'' not properly sized (Ntaps x 5)')
end
if min(taps(:,4))<1 | max(taps(:,4))>2
    errordlg('Column 4 of variable ''taps'' should be 1 or 2 only')
    error   ('Column 4 of variable ''taps'' should be 1 or 2 only')
end
X_loc           = taps(:,1);              % X coordinates
Y_loc           = taps(:,2);              % Y coordinates
Z_loc           = taps(:,3);              % Z coordinates
tap_dir         = taps(:,4);              % Tap directions (1-x direction, 2-y direction)
x_max           = max(X_loc);             % Internal variable
y_max           = max(Y_loc);             % Internal variable

%%%%%%%%%%%%%%% Change signs of pressure time series to match building axes %%%%%%%%%%%%%%
% F_t is the original matrix of forces at each tap calculated from the pressure
% coefficient files provided by the wind tunnel. The pressures are assumed
% positive when they act in the direction opposite to the outward normal to a 
% given face. It follows that on the face having outward normal coinciding with
% the positive x-axis, negative pressures act in the direction of the positive x-axis.
% The negative pressures and the negative forces calculated from them are multiplied by -1 to reflect 
% the fact that they acting in the direction of the positive x-axis. The same argument 
% applies to the face with outward normal coinciding the positive y-axis.  
k = find (X_loc==x_max | Y_loc==y_max);
F_t(k,:)        = F_t(k,:)*(-1);

%%%%%%%%%%%%%%%%%%%%%% Compile matrix: Load_Vector, [Ndofs x Npoints] %%%%%%%%%%%%%%%%%%%%%
% THE COMMENTED CODE BELOW WORKS ONLY WHEN TAP ELEVATIONS COINCIDE EXACTLY WITH FLOOR 
% ELEVATIONS
% NOTE: FLOORS ARE ASSUMED TO BE EQUALLY SPACED
% Load_Vector     = zeros(Ndofs,Npoints);   % Initial variable
% if Ndofs>2
%     for j=1:Ntaps
%         if tap_dir(j) == 1
%             Load_Vector(Z_loc(j)-2,:) = Load_Vector(Z_loc(j)-2,:) + F_t(j,:);                    
%             Load_Vector(Z_loc(j),:) = Load_Vector(Z_loc(j),:) - F_t(j,:)*Y_loc(j);        
%         else
%             Load_Vector(Z_loc(j)-1,:) = Load_Vector(Z_loc(j)-1,:) + F_t(j,:);     
%             Load_Vector(Z_loc(j),:) = Load_Vector(Z_loc(j),:) + F_t(j,:)*X_loc(j);         
%         end
%     end
% x               = zeros(Ndofs/3,Npoints); % Initialize variables
% y               = zeros(Ndofs/3,Npoints);
% z               = zeros(Ndofs/3,Npoints);
% for i=1:Ndofs/3
%     x(i,:)      = Load_Vector(3*i-2,:);
%     y(i,:)      = Load_Vector(3*i-1,:);
%     z(i,:)      = Load_Vector(3*i,:);
% end
% Load_Vector     = cat(1,x,y,z);           % [Ndofs x time]
% else
%     Load_Vector = F_t;                    % [Ndofs x time]
% end

%%%%%%%%%%%%%%%%%%%%%% Compile matrix: Load_Vector [Ndofs x Npoints] %%%%%%%%%%%%%%%%%%%%%% 
% THE CODE BELOW ASSUMES THAT THE TAP ELEVATIONS DO NOT NECESSARILLY
% COINCIDE EXACTLY WITH FLOOR ELEVATIONS; THE TAP IS "REASSIGNED" TO THE
% CLOSEST FLOOR
% NOTE: FLOORS ARE ASSUMED TO BE EQUALLY SPACED
FloorSeparation = H_bldg/Nfloors;
FloorHeightVec = FloorSeparation:FloorSeparation:H_bldg;
Load_Vector = zeros(Ndofs,Npoints);
for j = 1:Ntaps
    [Dist2NearestFloor,Index_NearestFloor] = min(abs(Z_loc(j)-FloorHeightVec));
    if tap_dir(j) == 1
        Load_Vector(Fdofs*Index_NearestFloor-2,:) = Load_Vector(Fdofs*Index_NearestFloor-2,:) + F_t(j,:);
        Load_Vector(Fdofs*Index_NearestFloor,:) = Load_Vector(Fdofs*Index_NearestFloor,:) - F_t(j,:)*Y_loc(j);
    else
        Load_Vector(Fdofs*Index_NearestFloor-1,:) = Load_Vector(Fdofs*Index_NearestFloor-1,:) + F_t(j,:);
        Load_Vector(Fdofs*Index_NearestFloor,:) = Load_Vector(Fdofs*Index_NearestFloor,:) + F_t(j,:)*X_loc(j);
    end
end
x               = zeros(Ndofs/Fdofs,Npoints);  % Initialize variables
y               = zeros(Ndofs/Fdofs,Npoints);
z               = zeros(Ndofs/Fdofs,Npoints);
for i = 1:Nfloors
    x(i,:) = Load_Vector(Fdofs*i-2,:);
    y(i,:) = Load_Vector(Fdofs*i-1,:);
    z(i,:) = Load_Vector(Fdofs*i,:);
end
Load_Vector = cat(1,x,y,z);  % [Ndofs x Npoints]

%%%%%%%%%%%%%%%%%%%%%%%% Integrate the Equation of Motion (EOM) %%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Call the function 'modal_eom_solve' to perform integrations
ModM            = Fai' * M * Fai;         % modal mass matrix     [ RNmodes x RNmodes ]
Modal_Load      = Fai' * Load_Vector;     % modal load vector     [ RNmodes x Npoints ]
Modal_LoadM     = inv(ModM)*Modal_Load;   % the right side of EOM [ RNmodes x Npoints ]

% Calculate generalized displacements (gen_d), velocities (gen_v), accelerations (gen_a)
% Each is [RNmodes x Npoints]
[gen_d,gen_v,gen_a] = modal_eom_solve(h,omega,Modal_LoadM,time);  
Npointsout      = length(gen_d);