附结构动力学计算程序,结构力学结构动力计算

%********************************************************************%子空间迭代%********************************************************************FP1=fopen('input.txt','rt');%打开数据输入文本niter=fscanf(FP1,'%f',1);%读入最多迭代次数nitern=fscanf(FP1,'%f',1);%读入自由度数目ncolumn=fscanf(FP1,'%f',1);%读入假设振型阶数columnerr=fscanf(FP1,'%f',1);%读入误差限errA0=xlsread('假设振型');%从excle读入假设振型M=xlsread('质量矩阵');%从excle读入质量矩阵K=xlsread('刚度矩阵');%从excle读入刚度矩阵forj=1:columnA0(:,j)=A0(:,j)/max(abs(A0(:,j)));%将假设振型A0基准化end%----开始迭代foriter=1:niter%----第iter次迭代F0=A0;F1=K\M*F0;forj=1:columnF1(:,j)=F1(:,j)/max(abs(F1(:,j)));endM1=F1'*M*F1;K1=F1'*K*F1;[V1,W1]=eig(M1\K1);forj=1:columna1(:,j)=V1(:,j)/V1(j,j);endA1=F1*a1;forj=1:columnA1(:,j)=A1(:,j)/max(abs(A1(:,j)));end%----精度判断fori=1:nforj=1:columnifabs(A1(i,j)-A0(i,j))>=errA0=A1;elseifabs(A1(i,j)-A0(i,j))<errbreakendendendendendW=sqrt(diag(W1));A=A1;%----输出自振特性xlswrite('自振特性.xlsx',A,1,'A2');%将数据写入xlsx文件中，数据的开始位置为sheet1，A2xlswrite('自振特性.xlsx',W,2,'A2');%将数据写入xlsx文件中，数据的开始位置为sheet2，A2%**************************************************************************%振型叠加法%**************************************************************************FP1=fopen('input.txt','rt');%打开数据输入文本dt=fscanf(FP1,'%f',1);%读入时间步长dtnt=fscanf(FP1,'%f',1);%读入计算次数ntn=fscanf(FP1,'%f',1);%读入自由度数目nnmode=fscanf(FP1,'%f',1);%读入选取振型数目nmodep=fscanf(FP1,'%f',1);%读入外力荷载参数wbar=fscanf(FP1,'%f',1);%读入外力荷载频率rdamp=xlsread('阻尼比');%从excle读入阻尼比M=xlsread('质量矩阵');%从excle读入质量矩阵K=xlsread(刚度矩阵');%从excle读入刚度矩阵A=xlsread('振型');%从excle读入前nmode阶振型W=xlsread('频率');%从excle读入前nmode阶频率%----确定与正则坐标对应的广义荷载向量fori=1:nt;forj=1:n;ifi>=0&i<=nt;Q(j,i)=p*sin(i*dt*wbar);end;end;end;%----将振型进行正则化处理Mgnrl=diag(A'*M*A);fori=1:nmodeforj=1:nAbar(j,i)=A(j,i)/sqrt(Mgnrl(i));endend%----本程序以例4-5为实例，基于其正则坐标响应的解析式，写出其稳态响应的程序，对于无解析式的荷载工况，此处需改为杜哈美数值积分程序fori=1:nmodeD(i)=1/sqrt((1-wbar^2/W(i)^2)^2+(2*rdamp(i)*wbar/W(i))^2);cta(i)=atan((2*rdamp(i)*wbar/W(i))/(1-wbar^2/W(i)^2));endfori=1:nmodeforj=1:ntTbar(i,j)=p*Abar(:,i)'*ones(n,1)*D(i)*sin(j*wbar*dt-cta(i))/W(i)^2;endend%----体系的稳态响应fori=1:nforj=1:ntd(i,j)=Abar(i,:)*Tbar(:,j);endend%----输出稳态响应i=1:nt;xlswrite('稳态响应.xlsx',d(:,i),1,'A1');%将数据写入xlsx文件中，数据的开始位置为sheet1，A1%**************************************************************************%newmark-β法%**************************************************************************FP1=fopen('input.txt','rt');%打开数据输入文件dt=fscanf(FP1,'%f',1);%读入时间步长dtnt=fscanf(FP1,'%f',1);%读入计算次数ntn=fscanf(FP1,'%f',1);%读入自由度数目nalfa=fscanf(FP1,'%f',1);%读入积分控制参数alfadta=fscanf(FP1,'%f',1);%读入积分控制参数dtaRayleigha0=fscanf(FP1,'%f',1);%读入瑞利阻尼待定常数Rayleigha0Rayleigha1=fscanf(FP1,'%f',1);%读入瑞利阻尼待定常数Rayleigha1p=fscanf(FP1,'%f',1);%读入外力荷载参数wbar=fscanf(FP1,'%f',1);%读入外力荷载频率M=xlsread('质量矩阵');%从excle上读入质量矩阵K=xlsread('刚度矩阵');%从excle上读入刚度矩阵C=Rayleigha0*M+Rayleigha1*K;%定义阻尼矩阵C%----外力荷载fori=1:ntforj=1:nifi>0&&i<=ntQ(j,i)=p*sin(i*dt*wbar);endendend%----计算动力响应b0=1/(alfa*dt^2);b1=dta/(alfa*dt);b2=1/(alfa*dt);b3=1/(2*alfa)-1;b4=dta/alfa-1;b5=dt/2*(dta/alfa-2);d(:,1)=zeros(n,1);v(:,1)=zeros(n,1);a(:,1)=M\(Q(:,1)-C*v(:,1)+K*d(:,1));Keff=b0*M+b1*C+K;fori=2:ntQeff=Q(:,i)+M*(b0*d(:,i-1)+b2*v(:,i-1)+b3*a(:,i-1))+C*(b1*d(:,i-1)+b4*v(:,i-1)+b5*a(:,i-1));d(:,i)=Keff\Qeff;a(:,i)=b0*(d(:,i)-d(:,i-1))-b2*v(:,i-1)-b3*a(:,i-1);v(:,i)=b1*(d(:,i)-d(:,i-1))-b4*v(:,i-1)-b5*a(:,i-1);end%----输出动力响应i=1:nt;xlswrite('动力响应.xlsx',d(:,i),1,'A1');%将数据写入xlsx文件中，数据的开始位置为sheet1，A1xlswrite('动力响应.xlsx',v(:,i),2,'A1');%将数据写入xlsx文件中，数据的开始位置为sheet2，A1xlswrite('动力响应.xlsx',a(:,i),3,'A1');%将数据写入xlsx文件中，数据的开始位置为sheet3，A1%**************************************************************************%wilson-θ法%**************************************************************************FP1=fopen('input.txt','rt');%打开数据输入文件dt=fscanf(FP1,'%f',1);%读入时间步长dtnt=fscanf(FP1,'%f',1);%读入计算次数ntn=fscanf(FP1,'%f',1);%读入自由度数目ncta=fscanf(FP1,'%f',1);%读入参数ctaRayleigha0=fscanf(FP1,'%f',1);%读入瑞利阻尼待定常数Rayleigha0Rayleigha1=fscanf(FP1,'%f',1);%读入瑞利阻尼待定常数Rayleigha1p=fscanf(FP1,'%f',1);%读入外力荷载参数wbar=fscanf(FP1,'%f',1);%读入外力荷载频率M=xlsread('质量矩阵');%从excle上读入质量矩阵K=xlsread('刚度矩阵');%从excle上读入刚度矩阵C=Rayleigha0*M+Rayleigha1*K;%定义阻尼矩阵C%----外力荷载fori=1:ntforj=1:nifi>0&&i<=ntQ(j,i)=p*sin(i*dt*wbar);endendend%----计算动力响应b0=6/(cta^2*dt^2);b1=3/(cta*dt);b2=6/(cta*dt);b3=cta*dt/2;b4=6/(cta^3*dt^2);b5=-6/(cta^2*dt);b6=1-3/cta;b7=dt/2;b8=dt^2/6;d(:,1)=zeros(n,1);v(:,