地震作用下房屋的碰撞编程求助

胡适之520

有没有会地震作用下房屋碰撞编程的大神，求助多谢。

胡适之520

求助大神，帮忙检查一下地震作用下碰撞部分的代码是否正确    [K,Kxo]=stiffnessTorsion(kx,ky,x,y,xm,ym,floor);  %结构的刚度矩阵       newmark法求解
    [C,T,z]=dampR(K,M,E,es,1);                        %结构的阻尼矩阵
    K1=K+(1/dt/dt/beta).*M+(1/2/dt/beta).*C;
    dF1=dF(:,i)+(M.*1/beta/dt+1/2/beta.*C)*vel(:,i)+(1/2/beta.*M-dt*(1-1/4/beta).*C)*acce(:,i);
    dx=inv(K1)*dF1;
    dv=(1/2/beta/dt).*dx-(1/2/beta).*vel(:,i)+((1-1/4/beta)*dt).*acce(:,i);
    dis(:,i+1)=dis(:,i)+dx;      
    vel(:,i+1)=vel(:,i)+dv;
    acce(:,i+1)=inv(M)*(F(:,i+1)-C*vel(:,i+1)-K*dis(:,i+1));                        
    %求解结束
   
    %第二栋newmark求解
    [K_,Kxo_]=stiffnessTorsion(kx_,ky_,x_,y_,xm_,ym_,floor_);  %结构的刚度矩阵
    [C_,T_,z_]=dampR(K_,M_,E_,es_,1);                          %结构的阻尼矩阵
    K2=K_+(1/dt_/dt_/beta_).*M_+(1/2/dt_/beta_).*C_;
    dF2=dF_(:,i)+(M_.*1/beta_/dt_+1/2/beta_.*C_)*vel_(:,i_)+(1/2/beta_.*M_-dt_*(1-1/4/beta_).*C_)*acce_(:,i);
    dx_=inv(K2)*dF2;
    dv_=(1/2/beta_/dt_).*dx_-(1/2/beta_).*vel_(:,i)+((1-1/4/beta_)*dt_).*acce_(:,i);
    dis_(:,i_+1)=dis_(:,i_)+dx_;        
    vel_(:,i_+1)=vel_(:,i_)+dv_;
    acce_(:,i_+1)=inv(M_)*(F_(:,i_+1)-C_*vel_(:,i_+1)-K_*dis_(:,i_+1));
    %第二栋参数结束
   
    %碰撞求解hertz-damp碰撞模型
    e1=zeros(cn1,1);
    e_=zeros(cn1,1);
    for z=1:cn1                                                  %同高度相碰撞
    e1(z)=-Kxo(z,z)/kx(z);                                             
    e_(z)=-Kxo_(z,z)/kx_(z);                                          
    end
    % 碰撞参数求解完毕
       e=0.65;                                                                   %hertz-damp碰撞模型
       kh=3e5;                                                                   %取值乱取
       gap=[0.001,0.001,0.001]';
%      jianxi=[dx(1:3,:)+10*abs(dx(7:9,:))-dx_(1:3,:)-10*abs(dx_(7:9,:))]-gap;   %dx第一栋 dx_第二栋
       jianxi=dx(1:3,:)-dx_(1:3,:)-gap;   %dx第一栋 dx_第二栋
       if max(jianxi(:))>0
       lc=find(jianxi>0);                                                     %找到间隙小于0.001,碰撞位置
       N=length(lc);                                                          %判断有几个发生碰撞
       jianxi(jianxi<0)=0;                                                    %未碰撞则为0
       Fc=kh*(jianxi.^(1.5)).*[1+3*kh*(1-e*e)/4./(vel(1:3,i+1)-vel_(1:3,i+1)).*(dis(1:3,i+1)-dis_(1:3,i+1))]; %全部X方向碰撞力求解
%      jianxi(jianxi<0)=0;
%      jj=find(jianxi<0);              %找到不碰撞的位置
%      N1=length(jj);
%      for ii=1:N1
%          Fc(jj(ii))=0;               %无碰撞为0，X方向碰撞力
%      end
      
       cell=mat2cell(K,seg,seg);         %矩阵划分小块
       Kxx=cell(1,1);                    %提取平移、扭转刚度
       Kyy=cell(2,2);
       Koo=cell(3,3);
       cell_=mat2cell(K_,seg_,seg_);
       Kxx_=cell_(1,1);                  
       Kyy_=cell_(2,2);
       Koo_=cell_(3,3);
      
      
       c=mat2cell(C,seg,seg);
       c_=mat2cell(C_,seg,seg);
       CC=blkdiag(cell2mat(c(1,1)),cell2mat(c(2,2)),cell2mat(c(3,3)),cell2mat(c_(1,1)),cell2mat(c_(2,2)),cell2mat(c_(3,3)));   
       KK=blkdiag(cell2mat(Kxx),cell2mat(Kyy),cell2mat(Koo),cell2mat(Kxx_),cell2mat(Kyy_),cell2mat(Koo_));







%      KK=blkdiag(K,K_);                                                 %刚度矩阵组合
%      CC=blkdiag(C,C_);                                                 %阻尼矩阵组合
       FFp=[Fc;zeros(cn1,1);Fc.*e1;-Fc;zeros(cn1_,1);-Fc.*e_];             %碰撞力
      
       KK_=KK+(1/dt/dt/beta).*MM+(1/2/dt/beta).*CC;
       dFF_=-dFp(:,i)+(MM.*1/beta/dt+1/2/beta.*CC)*velp(:,i)+(1/2/beta.*MM-dt*(1-1/4/beta).*CC)*accep(:,i)-FFp;   %根据求解公式，地震部分为负
       dxp=inv(KK_)*dFF_;
       dvp=(1/2/beta/dt).*dxp-(1/2/beta).*velp(:,i)+((1-1/4/beta)*dt).*accep(:,i);
       disp(:,i+1)=disp(:,i)+dxp;
       velp(:,i+1)=velp(:,i)+dvp;
       accep(:,i+1)=inv(MM)*(-FF(:,i+1)-CC*velp(:,i+1)-KK*disp(:,i+1)-FFp);         

      dis(:,i+1)=dis(:,i)+disp(1:3*cn1,i+1);                               %dxp位移                                 
      vel(:,i+1)=vel(:,i)+velp(1:3*cn1,i+1);                               %dvp速度
      acce(:,i+1)=acce(:,i)+accep(1:3*cn1,i+1);                            %accep加速度
      
      dis_(:,i_+1)=dis_(:,i_)+disp(3*cn1+1:3*(cn1+cn1_),i+1);        
      vel_(:,i_+1)=vel_(:,i_)+velp(3*cn1+1:3*(cn1+cn1_),i+1);
      acce_(:,i_+1)=acce_(:,i)+accep(3*cn1+1:3*(cn1+cn1_),i+1);
考虑X,Y方向以及转角的三方向。希望大神帮忙看看，在此感激不尽