no6mainchimneymatriaxiteration_waterboxbuffmotivation.m

主要用于高层建筑的减振设计,根据线性振动理论,用matlab高级语言编写而成,由一个主程序和四个子程序组成

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%*******************************大作业六***********************************%
%*********MainChimneyMatriaxIteration_WaterBoxBuffMotivation.m*************%
%***在本程序中将调用FlexM，MassM，CM三个子程序求固有频率，主振型，阻尼阵***%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;
%输入梁的初值
n=20;           %将杆分为n段
L=60;           %L--为杆长
d1=3.9;         %钢管的内半径m 
d2=4.0;         %钢管的外半径m
 
Density=7800;   %钢管的密度
E=2.1*10^11;
c=2*10^3;         %烟囱阻尼系数
cwater=4*10^4;    %水箱阻尼系数
%调用子程序MassM与FlexM求质量阵与柔度阵
FlexM=SUBFunFlexibilityMatrix_WaterBoxCantilever(L,n,d1,d2,E);                %flexibility
MassM=SUBFunMassMatrix_Cantilever(L,n,d1,d2,Density);
CM=SUBFunCMatrix_Cantilever(n,c,E);                                          %阻尼阵

%质量阵的扩充
r1=d1/2;      %钢管的内半径
r2=d2/2;      %钢管的外半径
A=pi*(r2^2-r1^2); 
MM=Density*A*L;
Linem=zeros(1,21);Rowm=zeros(20,1);
Linem(1,21)=MM*0.1;
MassM=[MassM Rowm;Linem];

%阻尼阵的扩充
CM(20,20)=cwater+c;
Linec=zeros(1,21);Rowc=zeros(20,1);
Rowc(20,1)=-cwater;
Linec(1,20)=-cwater;Linec(1,21)=cwater;
CM=[CM Rowc;Linec];

alpha=0;
% D=KM+alpha*MassM;
% delta=inv(D);
% Dstart0=delta*MassM;
Dstart0=FlexM*MassM;
n=n+1;

A0=ones(n,1);
B1=Dstart0*A0; AA1=B1/B1(n,1); 
B2=Dstart0*AA1;AA2=B2/B2(n,1); 
while  max(abs(AA2-AA1))>1e-6
    AA1=AA2;    
    B=Dstart0*AA1;
    AA2=B/B(n,1);
end
PstartSqure=1/B(n,1);
P(1,1)=sqrt(PstartSqure-alpha);
A(1:n,1)=AA2;

AP=8;
Dstart(1:n,1:n,1)=Dstart0;
for i=1:AP-1
  Mp=A(1:n,i)'*MassM*A(1:n,i);
  Dstart(1:n,1:n,i+1)=Dstart(1:n,1:n,i)-A(1:n,i)*A(1:n,i)'*MassM/(P(i,1)^2*Mp);
  
  A0=ones(n,1);
  B1=Dstart(1:n,1:n,i+1)*A0; AA1=B1/B1(n,1);
  B2=Dstart(1:n,1:n,i+1)*AA1;AA2=B2/B2(n,1);
  false=1;
  while  max(abs(AA2-AA1))>1e-6
      AA1=AA2;    
      B=Dstart(1:n,1:n,i+1)*AA1;
      AA2=B/B(n,1);
      if false==1e6
          break;
      end
      false=false+1;
  end
  PstartSqure=1/B(n,1);
  P(i+1,1)=sqrt(PstartSqure-alpha);
  A(1:n,i+1)=AA2;  
end
%********************************************
%求响应
%*******************************************
%----------------------------------------------
%激励条件
%---------------------------------------------
f=2;
w=2*pi*f; 
T=1/f;
F0=zeros(n,1);
F0(20,1)=5000; %N
%----------------------------------------
%解耦
Mp=A'*MassM*A;
for j=1:8
AN(:,j)=A(:,j)/sqrt(Mp(j,j));
end
FN0=AN'*F0;          %正则激励
CMN=AN'*CM*AN;       %正则阻尼系数
for j=1:8
    radm(j,1)=CMN(j,j)/(2*P(j,1));%阻尼比ratio of damping
    Pd(j,1)=P(j,1)*sqrt(1-radm(j,1)^2);
end
%---------------------------------------
FN00=FN0;
PP=P;
Pdd=Pd;
%------------------------------------------
% t<T/2时，烟囱顶点的响应
%------------------------------------------
for cc=1:8
    FN0=FN00(cc,1);
    rp=radm(cc,1)*PP(cc,1); %衰减系数
    Pd=Pdd(cc,1);    
    for j=1:25
        t=j/100;
        %符号积分，求杜哈美积分
        XliT(j,cc)=FN0*(exp(-rp*t)*w^3*sin(Pd*t)-sin(w*t)*w^2*Pd+(-2*rp*cos(w*t)*Pd+exp(-rp*t)*sin(Pd*t)*rp^2+2*exp(-rp*t)*rp*cos(Pd*t)*Pd-...
                     exp(-rp*t)*sin(Pd*t)*Pd^2)*w+sin(w*t)*Pd*rp^2+sin(w*t)*Pd^3)/Pd/(rp^2+w^2+2*w*Pd+Pd^2)/(rp^2+w^2-2*w*Pd+Pd^2);
     end
end

Txli=0;
for cc=1:8
    Txli=Txli + AN(20,cc)* XliT(:,cc) %0-0.25s（即0-T/2）之间烟囱顶点的响应
end 
%----------------------------------------------
% t>T/2时，烟囱顶点的响应
%----------------------------------------------
for cc=1:8
    FN0=FN00(cc,1);
    rp=radm(cc,1)*PP(cc,1); %衰减系数
    Pd=Pdd(cc,1); 
    jj=0;
    for j=25:10000
        t=j/100;
        t1=0.25;
        jj=jj+1;
        %符号积分，求杜哈美积分
        XlaTW(jj,cc)=-1/2*FN0*((exp(-rp*t+rp*t1)*sin(-t1*w-t1*Pd+Pd*t)-2*exp(-rp*t)*sin(Pd*t)+exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t))*...
                   w^3+(-exp(-rp*t+rp*t1)*rp*cos(-t1*w-t1*Pd+Pd*t)+exp(-rp*t+rp*t1)*rp*cos(t1*w-t1*Pd+Pd*t)-exp(-rp*t+rp*t1)*...
                   sin(-t1*w-t1*Pd+Pd*t)*Pd+exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t)*Pd)*w^2+(-exp(-rp*t+rp*t1)*sin(-t1*w-t1*Pd+Pd*t)*...
                   Pd^2+exp(-rp*t+rp*t1)*sin(-t1*w-t1*Pd+Pd*t)*rp^2+2*exp(-rp*t+rp*t1)*rp*cos(-t1*w-t1*Pd+Pd*t)*Pd+2*exp(-rp*t)*sin(Pd*t)*...
                   Pd^2+2*exp(-rp*t+rp*t1)*rp*cos(t1*w-t1*Pd+Pd*t)*Pd+exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t)*rp^2-4*exp(-rp*t)*rp*cos(Pd*t)*...
                   Pd-exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t)*Pd^2-2*exp(-rp*t)*sin(Pd*t)*rp^2)*w-exp(-rp*t+rp*t1)*rp^3*cos(-t1*w-t1*Pd+Pd*t)-...
                   exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t)*Pd*rp^2-exp(-rp*t+rp*t1)*sin(t1*w-t1*Pd+Pd*t)*Pd^3-exp(-rp*t+rp*t1)*rp*...
                   cos(-t1*w-t1*Pd+Pd*t)*Pd^2+exp(-rp*t+rp*t1)*sin(-t1*w-t1*Pd+Pd*t)*Pd*rp^2+exp(-rp*t+rp*t1)*sin(-t1*w-t1*Pd+Pd*t)*...
                   Pd^3+exp(-rp*t+rp*t1)*rp^3*cos(t1*w-t1*Pd+Pd*t)+exp(-rp*t+rp*t1)*rp*cos(t1*w-t1*Pd+Pd*t)*Pd^2)/Pd/(rp^2+w^2+2*w*Pd+Pd^2)...
                   /(rp^2+w^2-2*w*Pd+Pd^2); 
    end
end
TxlaW=0;
for cc=1:8
    TxlaW=TxlaW + AN(20,cc)* XlaTW(:,cc) %0-0.25s（即0-T/2）之间烟囱顶点的响应
end
%--------------------------------------------
%可视化
%-------------------------------------------
figure(1)
 plot(Txli(:,1));
 ylabel('塔高(60米)');
 title('0-0.25s顶点的响应图');
 set(gca,'ytick',[0]);grid on;
 
 figure(2) 
 plot(TxlaW(:,1))
 ylabel('塔高(60米)');
 title('0.25-200s顶点的响应图');
 set(gca,'ytick',[0]);grid on;