program1.m

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

    end
    if size(evectors,1) ~= Ndofs 
        errordlg('Variable ''evectors'' not properly sized (Ndofs x RNmodes)')
        error   ('Variable ''evectors'' not properly sized (Ndofs x RNmodes)')
    end
    Fai      = evectors; % Assign internal variable to 'evectors'
    faix     = Fai(1:(Ndofs/Fdofs),1:RNmodes);                               % Assign internal variable for d.o.f. 1
    if RNmodes > 1
    faiy     = Fai((Ndofs/Fdofs)+1:(Fdofs-1)*(Ndofs/Fdofs),1:RNmodes);       % Assign internal variable for d.o.f. 2
    end
    if RNmodes > 2
    fairz    = Fai((Fdofs-1)*(Ndofs/Fdofs)+1:Ndofs,1:RNmodes);               % Assign internal variable for d.o.f. 3
    end
end    % loop back on AnL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% WIND EFFECT CALCULATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate peak wind effect of interest for each wind direction (WD) and wind speed (WS)
w = waitbar(0,'Running Program1.m   Please wait...'); c=0;
a               = 0;                            % Counting variable for WD
for ddd = 1:length(WD)
    dir         = WD(ddd);                      % Internal variable for direction (WD)
    a           = a+1;
    b           = 0;                            % Counting variable for WS
    dir_list(a) = dir;                          % Variable to store wind directions (WD) used
    Cp_tr       = F_Cp_tr(dir,Ntaps,Npoints); % 'F_Cp_tr' loads the pressure coefficient time series 
                                                % and stores them in a matrix Cp_tr of size Ntaps x Npoints 
    for vvv = 1:length(WS)
        b           = b+1;
        V           = WS(vvv);                    % Reassign variable
        V_list(b)   = V;                          % Variable to store wind speeds (WS) used
        dt          = (1/freq)*(Vm/V)*(ms);       % Time increment such that (nD/V)model = (nD/V)prototype
        F_t         = F_Ft(Cp_tr,V);              % 'F_Ft' calculates the wind load matrix, F_t, of size
                                                  % Ntaps x Npoints, for the dynamic MDOF system
        WS_WD       = [V dir]
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Run a loop to integrate equations of motion if inertial loads are to be included
        if AnL_A == 1
            acc_y = zeros(Nfloors,Npoints);
            acc_rz = zeros(Nfloors,Npoints);
            % This function calculates the generalized displacements, velocities and acceleration
            [gen_d,gen_v,gen_a,Npointsout] = F_Load_Vector(M,Fai,F_t,dt,h,omega);
            disp_x      = faix * gen_d;               % Mass displacements, d.o.f. 1 (Nfloors x Npoints)
            vel_x       = faix * gen_v;               % Mass velocities, d.o.f. 1    (Nfloors x Npoints)
            acc_x       = faix * gen_a;               % Mass accelerations, d.o.f. 1 (Nfloors x Npoints)
            if RNmodes > 1
            disp_y      = faiy * gen_d;               % Mass displacements, d.o.f. 2 (Nfloors x Npoints)
            vel_y       = faiy * gen_v;               % Mass velocities, d.o.f. 2    (Nfloors x Npoints)
            acc_y       = faiy * gen_a;               % Mass accelerations, d.o.f. 2 (Nfloors x Npoints)
            end
            if RNmodes > 2
            disp_rz     = fairz * gen_d;              % Mass displacements, d.o.f. 3 (Nfloors x Npoints)
            vel_rz      = fairz * gen_v;              % Mass velocities, d.o.f. 3    (Nfloors x Npoints)
            acc_rz      = fairz * gen_a;              % Mass accelerations, d.o.f. 3 (Nfloors x Npoints)
            end
        else
            Npointsout  = Npoints;                    % Assign internal variable
        end    % loop back on AnL
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MEMBER EFFECTS
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for q=1:length(mem_list)
            mem         = mem_list(q);                 % Internal variable
            %%%%%%%%%% Dynamic influence coefficients : dif, Page 3, [ 3*(Ndofs + Ntaps) x 6 ]%%%%%%%%
            % NOTE: The dif's are dimensionless
            % Load flDif and check for errors
            clear fid; fid=fopen([flDif  '\dif_'  num2str(mem) '.mat']);
            if fid == -1
                errordlg(['Unable to open ''flDif  dif_' num2str(mem) ' - Page 3'])
                error   (['Unable to open ''flDif  dif_' num2str(mem) ' - Page 3'])
            end
            fclose(fid); clear fid;
            clear dif
            dif=[];                                    % Initialize 'dif'
            load ([flDif  '\dif_'  num2str(mem) ])
            if isempty(dif)==1
                errordlg(['Variable ''dif'' not defined for dif_' num2str(mem_list(q)) '. '...
                    'See flDif on Page 3'])
                error   (['Variable ''dif'' not defined for dif_' num2str(mem_list(q)) '. '...
                    'See flDif on Page 3'])
            end
            if size(dif,1) ~= 3*(Ndofs+Ntaps) | size(dif,2) ~= 6
                errordlg(['Variable ''dif'' not properly sized (3*[Ndofs+Ntaps] x 6)'])
                error   (['Variable ''dif'' not properly sized (3*[Ndofs+Ntaps] x 6)'])
            end
            P_col_no    = 1;                           % Assign local variable, axial load
            M2_col_no   = 5;                           % Assign local variable, moment M2
            M3_col_no   = 6;                           % Assign local variable, moment M3
            for axpos_index = 1:3  
                P_inertial  = zeros(1,Npointsout);     % Initialize variables
                P_applied   = zeros(1,Npoints);
                M2_inertial = zeros(1,Npointsout);
                M2_applied  = zeros(1,Npoints);
                M3_inertial = zeros(1,Npointsout);
                M3_applied  = zeros(1,Npoints);          
                P_applied   = dif([axpos_index*Ndofs+(axpos_index-1)*Ntaps+1:axpos_index*(Ndofs+Ntaps)], P_col_no)'  * F_t;
                M2_applied  = dif([axpos_index*Ndofs+(axpos_index-1)*Ntaps+1:axpos_index*(Ndofs+Ntaps)], M2_col_no)' * F_t;
                M3_applied  = dif([axpos_index*Ndofs+(axpos_index-1)*Ntaps+1:axpos_index*(Ndofs+Ntaps)], M3_col_no)'* F_t;

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                
                % Run a loop to calculate internal forces if inertial loads are to be included
                if AnL_A == 1
                    for i=1:Nfloors                     % Matrix multiply dif* mass*acceleration
                        P_inertial  = P_inertial+ ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+i,P_col_no)'*mx(i)*acc_x(i,:) + ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+Nfloors+i,P_col_no)'*my(i)*acc_y(i,:) + ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+2*Nfloors+i,P_col_no)'*mrz(i)*acc_rz(i,:);
                        M2_inertial = M2_inertial+ ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+i,M2_col_no)'*mx(i)*acc_x(i,:) + ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+Nfloors+i,M2_col_no)'*my(i)*acc_y(i,:) + ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+2*Nfloors+i,M2_col_no)'*mrz(i)*acc_rz(i,:);
                        M3_inertial = M3_inertial+ ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+i,M3_col_no)'*mx(i)*acc_x(i,:) + ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+Nfloors+i,M3_col_no)'*my(i)*acc_y (i,:)+ ...
                            dif((axpos_index-1)*(Ndofs+Ntaps)+2*Nfloors+i,M3_col_no)'*mrz(i)*acc_rz(i,:);
                    end
                end     % loop on AnL
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Otherwise, assign internal forces to applied loads if inertial loads are not included
                P_t         = P_applied (1:Npointsout) - P_inertial (1:Npointsout) ;
                M2_t        = M2_applied(1:Npointsout) - M2_inertial(1:Npointsout);
                M3_t        = M3_applied(1:Npointsout) - M3_inertial(1:Npointsout);
                % Sample the time history using tstep = (rf:rf:Npoints)
                tstep        = rf:rf:Npoints;
                P_tr         = P_t(tstep);
                M2_tr        = M2_t(tstep);
                M3_tr        = M3_t(tstep);
                % Factor the internal forces using (internal) load factors
                Pu           = DL_fact*(P_DL(q)  + P_SDL(q))  + LL_fact*(P_LL(q))  + WL_fact*P_tr;
                M2u          = DL_fact*(M2_DL(q) + M2_SDL(q)) + LL_fact*(M2_LL(q)) + WL_fact*M2_tr;
                M3u          = DL_fact*(M3_DL(q) + M3_SDL(q)) + LL_fact*(M3_LL(q)) + WL_fact*M3_tr;
                % Calculate wind effect of interest (LRFD Spec. Interaction Equations)
                aa1          = find( Pu > 0   &   Pu/(0.9*Pnt(q))< 0.2  );
                Bij(aa1)     = abs(Pu(aa1))/(2*0.9*Pnt(q)) +        abs(M2u(aa1))/(0.9*Mn2(q)) + ...
                    abs(M3u(aa1))/(0.9*Mn3(q));
                aa2          = find( Pu> 0   &   Pu/(0.9*Pnt(q))>=0.2  );
                Bij(aa2)     = abs(Pu(aa2))/(  0.9*Pnt(q)) + (8/9)*(abs(M2u(aa2))/(0.9*Mn2(q)) + ...
                    abs(M3u(aa2))/(0.9*Mn3(q)));
                bb1          = find( Pu<=0   &   abs(Pu/(0.9*Pnc(q)))<0.2  );
                Bij(bb1)     = abs(Pu(bb1))/(2*0.9*Pnc(q)) +        abs(M2u(bb1))/(0.9*Mn2(q)) + ...
                    abs(M3u(bb1))/(0.9*Mn3(q));
                bb2          = find( Pu<=0   &   abs(Pu/(0.9*Pnc(q)))>=0.2  );
                Bij(bb2)     = abs(Pu(bb2))/(  0.9*Pnc(q)) + (8/9)*(abs(M2u(bb2))/(0.9*Mn2(q)) + ...
                    abs(M3u(bb2))/(0.9*Mn3(q)));
                % Save the wind effect time series for this member => creates very large files
                K1 = findstr(flnSaveBij,'.mat');
                if ~isempty(K1)
                    flnSaveBij = flnSaveBij(1:K1-1);
                end
                if saveBij==1 
                    save ([flnSaveBij '_' num2str(mem) '_' num2str(dir) '_' num2str(V)],...
                        'Bij','mem','P_DL','P_SDL','P_LL','M2_DL','M2_SDL','M2_LL','M3_DL',...
                        'M3_SDL','M3_LL','DL_fact','LL_fact','WL_fact', 'disp_x', 'disp_y');
                end
                % Estimate peaks (if selected) using function 'maxminest'
                if Qs == 1
                    [max_est,min_est] = maxminest (Bij);
                    eval ([ 'Bij_obs_max_' num2str(mem) '(a,b,axpos_index) = max(Bij);']);
                    eval ([ 'Bij_obs_min_' num2str(mem) '(a,b,axpos_index) = min(Bij);']);
                    eval ([ 'Bij_est_max_' num2str(mem) '(a,b,axpos_index) = max_est;']);
                    eval ([ 'Bij_est_min_' num2str(mem) '(a,b,axpos_index) = min_est; ']);
                else
                    eval ([ 'Bij_obs_max_' num2str(mem) '(a,b,axpos_index) = max(Bij);']);
                    eval ([ 'Bij_obs_min_' num2str(mem) '(a,b,axpos_index) = min(Bij);']);
                    eval ([ 'Bij_est_max_' num2str(mem) '(a,b,axpos_index) = 0; ']);
                    eval ([ 'Bij_est_min_' num2str(mem) '(a,b,axpos_index) = 0; ']);
                end
            end  % loop on member axial position index (axpos_index)
        end   % loop on member (q)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        c           = c+1;
        counter     = c/(length(WD)*length(WS));
        waitbar(counter,w);
    end   % loop on WS
end   % loop on WD
close(w)

%%%%%%%%%%%%%%%%%%%%%%%%%%% SAVE PEAK WIND EFFECTS FOR EACH MEMBER %%%%%%%%%%%%%%%%%%%%%%%%%
% Save the four matrices using 'flnSavemaxBij' and the member number
K2 = findstr(flnSavemaxBij,'.mat');
if ~isempty(K2)
    flnSavemaxBij = flnSavemaxBij(1:K2-1);
end
for i=1:length(mem_list)
    mem     = mem_list(i);
    save ( [flnSavemaxBij '_' num2str(mem)],...
        ['Bij_obs_max_' num2str(mem)], ['Bij_obs_min_' num2str(mem)], ...
        ['Bij_est_max_' num2str(mem)], ['Bij_est_min_' num2str(mem)], ...
        'V_list','dir_list','P_DL','P_SDL','P_LL','M2_DL','M2_SDL','M2_LL',...
        'M3_DL','M3_SDL','M3_LL','DL_fact','LL_fact','WL_fact','intmeth','anal');
end