经过几天的程序处理程序的编写,昨天完成了平均风压系数和等效静力风荷载程序的编写。
  与师兄之前的项目相比,由一千多行的程序缩减到两百行;
  从数据统计来说,避免了很多重复性工作,取而代之的是用了一个interp1的函数完成了线性插值的过程;
  从程序效率上来说,在单次循环过程中,内部循环次数减少为一次,为不可避免的循环过程。
  程式如下:

%% 程序信息
% 程序内容:测压试验数据处理程序
% 参考资料:《建筑结构荷载规范》
% 编译单位:湖南大学建筑安全与节能教育部重点实验室
% 创建时间:2014年12月15日 
%% 初始化
clear,clc,close all
%% 计算总程式
for i=35:-1:0;
    %% 计算平均风压系数
    tic,disp(['当前计算风向角为',num2str(10*i),'°']);
    PP=load ([pwd,'\风压\',num2str(10*i),'.txt']);    %加载测点风压时程
    P=[PP(:,498:513),PP(:,194:497),PP(:,2:49)];       %提取主体测点层数据
    PX=[PP(:,50:65),PP(:,130:135)]-PP(:,136:157);
    PY=PP(:,158:169)-PP(:,170:181);
    % 测压坏点修正
    P(:,353)=0.95*P(:,354);              %W9
    P(:,392)=0.5*(P(:,391)+P(:,393));    %Y2
    % 风压系数计算
    [fid,message]=fopen([pwd,'\风速\',num2str(10*i),' (Ve).asA']);
    fseek(fid,523,'bof');
    v=fscanf(fid,'%f',1);               %实验室参考高度风速,300m;
    fclose(fid);
    h0=281.7;                           %参考高度
    v=v*(h0/300)^0.22;                  %建筑物顶部采样风速
    disp(['当前计算风速为',num2str(v),'m/s']);
    pr=0.5*1.25*v^2;                    %参考风压
    P0=8;                               %参考静压
    Pmr=[P-P0,PX,PY];                   %测点净风压
    
    Cpall=Pmr./pr;                      %测点的风压系数时程
    Cpmean=mean(Cpall);                 %测点的平均风压系数
    Cpstd=std(Cpall);                   %测点的脉动风压系数
    % 汇总各风向角下风压系数
    CpmeanALL(:,i+1)=Cpmean';
    CpstdALL(:,i+1)=Cpstd';
    %% 计算等效静力风荷载
    %采样环境信息
    w0=400;                             %建筑地区50年基本风压,1.1倍50年风压
    wr=0.544*w0*(h0/10)^0.44;           %参考高度地区风压
    vp=(2*wr/1.25)^0.5;                 %参考高度地区风速
    fn1=312.5;                          %扫描阀采样频率
    fn2=fn1*vp/(v*300);                 %实际结构采样频率,300为缩尺比
    %建筑结构信息(建筑信息由上往下统计!)
    nx=0;                               %X向结构频率
    ny=0;                               %Y向结构频率
    nt=0;                               %扭转结构频率
    damptio=0.035;                      %阻尼比,分析时取值0.035,0.04,0.05
    zx=load('zx,txt');                  %X向结构振型
    zy=load('zy,txt');                  %Y向结构振型
    zt=load('zt,txt');                  %扭转结构振型
    mass=load('mass.txt');              %结构质量,kg
    zdgl=load('zdgl.txt');              %转动惯量,kg
    dH=load('dH.txt');                  %各层层高
    H=load('H.txt');                    %各层标高
    %测点信息(测点信息由下往上统计!)
    lx=load('lx.txt');                  %测点积分长度
    ly=load('ly.txt');                  %测点积分长度
    dx=load('dx.txt');                  %测点力臂
    dy=load('dy.txt');                  %测点力臂
    h=load('h.txt');                    %测点标高
    %计算测点层风荷载
    fxa=wr.*Cpall.*repmat(lx',size(Cpall,1),1);
    fya=wr.*Cpall.*repmat(ly',size(Cpall,1),1);
    fta=wr.*Cpall.*repmat(lx'.*dx'+ly'.*dy',size(Cpall,1),1);
    fxall=[
        sum(fxa(:,1:20),2),...          %C层(01)
        sum(fxa(:,21:40),2),...         %D层(02)
        sum(fxa(:,41:60),2),...         %E层(03)
        sum(fxa(:,61:80),2),...         %F层(04)
        sum(fxa(:,81:100),2),...        %G层(05)
        sum(fxa(:,101:120),2),...       %H层(06)
        sum(fxa(:,121:140),2),...       %J层(07)
        sum(fxa(:,141:160),2),...       %K层(08)
        sum(fxa(:,161:180),2),...         %L层(09)
        sum(fxa(:,181:200),2),...       %M层(10)
        sum(fxa(:,201:220),2),...        %N层(11)
        sum(fxa(:,221:240),2),...       %P层(12)
        sum(fxa(:,241:260),2),...       %Q层(13)
        sum(fxa(:,261:280),2),...       %R层(14)
        sum(fxa(:,281:300),2),...       %S层(15)
        sum(fxa(:,301:320),2),...       %T层(16)
        sum(fxa(:,321:344),2),...        %U层(17)
        sum(fxa(:,345:368),2),...       %W层(18)
        sum(fxa(:,369:390),2),...         %X层(19)
        sum(fxa(:,391:402),2)...        %Y层(20)
        ];
    fyall=[
        sum(fya(:,1:20),2),...          %C层
        sum(fya(:,21:40),2),...         %D层
        sum(fya(:,41:60),2),...         %E层
        sum(fya(:,61:80),2),...         %F层
        sum(fya(:,81:100),2),...        %G层
        sum(fya(:,101:120),2),...       %H层
        sum(fya(:,121:140),2),...       %J层
        sum(fya(:,141:160),2),...       %K层
        sum(fya(:,161:180),2),...         %L层
        sum(fya(:,181:200),2),...       %M层
        sum(fya(:,201:220),2),...        %N层
        sum(fya(:,221:240),2),...       %P层
        sum(fya(:,241:260),2),...       %Q层
        sum(fya(:,261:280),2),...       %R层
        sum(fya(:,281:300),2),...       %S层
        sum(fya(:,301:320),2),...       %T层
        sum(fya(:,321:344),2),...        %U层
        sum(fya(:,345:368),2),...       %W层
        sum(fya(:,369:390),2),...         %X层
        sum(fya(:,391:402),2)...        %Y层
        ];
    ftall=[        
        sum(fta(:,1:20),2),...          %C层
        sum(fta(:,21:40),2),...         %D层
        sum(fta(:,41:60),2),...         %E层
        sum(fta(:,61:80),2),...         %F层
        sum(fta(:,81:100),2),...        %G层
        sum(fta(:,101:120),2),...       %H层
        sum(fta(:,121:140),2),...       %J层
        sum(fta(:,141:160),2),...       %K层
        sum(fta(:,161:180),2),...         %L层
        sum(fta(:,181:200),2),...       %M层
        sum(fta(:,201:220),2),...        %N层
        sum(fta(:,221:240),2),...       %P层
        sum(fta(:,241:260),2),...       %Q层
        sum(fta(:,261:280),2),...       %R层
        sum(fta(:,281:300),2),...       %S层
        sum(fta(:,301:320),2),...       %T层
        sum(fta(:,321:344),2),...        %U层
        sum(fta(:,345:368),2),...       %W层
        sum(fta(:,369:390),2),...         %X层
        sum(fta(:,391:402),2)...        %Y层
        ];
    fpxall=mean(fxall);                          %测点层平均压力
    fpyall=mean(fyall);                          %测点层平均压力
    fptall=mean(ftall);                          %测点层平均压力
    fmxall=fxall-repmat(fpxall,size(fpxall,1),1);%测点层脉动时程
    fmyall=fyall-repmat(fpyall,size(fpyall,1),1);%测点层脉动时程
    fmtall=ftall-repmat(fptall,size(fptall,1),1);%测点层脉动时程
    %计算各结构楼层风荷载
    Fpx=interp1([0,h'],[fpxall(1),fpxall],H','linear')'.*dH;
    Fpy=interp1([0,h'],[fpyall(1),fpyall],H','linear')'.*dH;
    Fpt=interp1([0,h'],[fptall(1),fptall],H','linear')'.*dH;
    for j=size(Cpall):-1:1
        Fmx(:,j)=interp1([0,h'],[fmxall(j,1),fmxall(j,:)],H','linear')'.*dH;
        Fmy(:,j)=interp1([0,h'],[fmyall(j,1),fmyall(j,:)],H','linear')'.*dH;
        Fmt(:,j)=interp1([0,h'],[fmtall(j,1),fmtall(j,:)],H','linear')'.*dH;
    end
    %施加非结构层荷载(最顶层楼板上)
    Fpx(1)=Fpx(1)+fpxall(end-1)*8.4+fpxall(end)*8.23;
    Fpy(1)=Fpy(1)+fpyall(end-1)*8.4+fpyall(end)*8.23;
    Fpt(1)=Fpt(1)+fptall(end-1)*8.4+fptall(end)*8.23;
    Fmx(1,:)=Fmx(1,:)+fmx(:,end-1)'*8.4+fmx(:,end)'*8.23;
    Fmy(1,:)=Fmy(1,:)+fmy(:,end-1)'*8.4+fmy(:,end)'*8.23;
    Fmt(1,:)=Fmt(1,:)+fmt(:,end-1)'*8.4+fmt(:,end)'*8.23;
    %广义荷载
    pfx=zx'*Fmx;
    pfy=zy'*Fmy;
    pft=zt'*Fmx;
    %广义质量
    mx=mass'*zx.^2;
    my=mass'*zy.^2;
    mt=zdgl'*zt.^2;
    %广义荷载谱以及频响函数的平方
    nft=1024*8;sf=fn2;window=boxcar(nft/2);noverlap=0;
    [Spx,fx]=pmtm(pfx,4,nft,sf); 
    [Spy,fy]=pmtm(pfy,4,nft,sf);
    [Spt,ft]=pmtm(pft,4,nft,sf);
    sHx=(1/(mx^2*(2*pi*nx)^4))*(1./((1-fx.^2/nx^2).^2+(2*damptio*fx/nx).^2));%X向频响
    sHy=(1/(my^2*(2*pi*ny)^4))*(1./((1-fy.^2/ny^2).^2+(2*damptio*fy/ny).^2));%Y向频响
    sHt=(1/(mt^2*(2*pi*nt)^4))*(1./((1-ft.^2/nt^2).^2+(2*damptio*ft/nt).^2));%扭转频响
    nnn1=length(fx);nnn2=length(fy);nnn3=length(ft);
    sHSpx=Spx*sHx'.*fx(nnn1)/nnn1;
    sHSpy=Spy*sHy'.*fy(nnn2)/nnn2;
    sHSpt=Spt*sHt'.*ft(nnn3)/nnn3;
    sigmax=abs(zx*sHSpx^0.5);   %X向位移响应均方根
    sigmay=abs(zy*sHSpy^0.5);   %Y向位移响应均方根
    sigmat=abs(zt*sHSpt^0.5);   %扭转位移响应均方根
    %保证因子
    %grx=(2*log(nx*600))^0.5+0.5772/(2*log(nx*600))^0.5;
    %gry=(2*log(ny*600))^0.5+0.5772/(2*log(ny*600))^0.5;
    %grt=(2*log(nt*600))^0.5+0.5772/(2*log(nt*600))^0.5;
    grx=3.5;
    gry=3.5;
    grt=3.5;
    %脉动部分等效静力风荷载
    Fefx=grx*(2*pi*nx)^2*mass.*sigmax;
    Fefy=gry*(2*pi*ny)^2*mass.*sigmay;
    Feft=grt*(2*pi*nt)^2*zdgl.*sigmat;
    %总的等效静力风荷载
    Feswlx1=Fpx+Fefx.*sign(sum(Fpx));
    Feswly1=Fpy+Fefy.*sign(sum(Fpy));
    Feswlt1=Fpt+Feft.*sign(sum(Fpt));
    %汇总各风向等效静力风荷载
    AAFeswl1(:,1+3*i)=Feswlx1;
    AAFeswl1(:,2+3*i)=Feswly1;
    AAFeswl1(:,3+3*i)=Feswlt1;
    %% 基础等效静力弯矩
    jPX1=sum(Feswlx1);
    jPY1=sum(Feswly1);
    Meswlx1=Feswly1'.*H;
    Meswly1=Feswlx1'.*H;
    Meswlt1=sum(Feswlt1);
    AAPeqwl1(i,:)=[jPX1,jPY1,Meswlx1,Meswly1,Meswlt1]; 
    %% 计算顶部加速度
    aHSpx=sum(Spx.*(2.*pi.*fx).^4.*sHx.*fx(nnn1)./nnn1);
    aHSpy=sum(Spy.*(2.*pi.*fy).^4.*sHy.*fy(nnn2)./nnn2);
    aHSpt=sum(Spt.*(2.*pi.*ft).^4.*sHt.*ft(nnn3)./nnn3);
    accx=grx*zx*aHSpx^0.5;
    accy=gry*zy*aHSpy^0.5;
    acct=gry*zt*aHSpt^0.5;
    %汇总各风向顶部加速度
    AACC(i,:)=[accx(1),accy(1),acct(1)];
    disp(['风向角为',num2str(10*i),'°用时:',num2str(round(toc)),'秒'])
end

标签: 山东临沂

评论已关闭