资源描述：
《结构的动力方程.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、结构的动力学方程clear;clc;n=4;II=sqrt(-1);%复数i%主结构质量、阻尼、刚度矩阵M=eye(n)*1.0e+4;%对角阵K=eye(n)*1.6*1.0e+7;%主结构刚度矩阵聚合zk=zeros(n);forj=1:(n-1)zk(j,j)=K(j,j)+K(j+1,j+1);zk(j,j+1)=-K(j+1,j+1);zk(j+1,j)=-K(j+1,j+1);endzk(n,n)=K(n,n);k=zk;m=M;%此处我的结构模型的质量、刚度矩阵应该从其它有限元软件导出来。%求解各阶模态频率
2、[tzxl,tzz]=eig(k,m);d=diag(sqrt(tzz));%振型规一化fori=1:n[tzz1(i),j]=min(d);tzxl1(:,i)=tzxl(:,j);d(j)=max(d)+1;end%振型归一化取第一层振型为1forj=1:ntzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);endw0=tzz1;w=tzz1/(2*pi);zhx=tzxl1;广义阻尼矩阵各阶模态阻尼比都取%阻尼比ks0=0.05;ks=ones(n,1)*ks0;第n阶广义质量：%求广义质量Mn=zh
3、x'*m*zhx;阻尼矩阵为：%求阻尼矩阵C=zeros(n);fori=1:nC(i,i)=2*ks(i)*w0(i)*Mn(i,i);endc=(zhx')C/zhx;参数eg即%过滤白噪声参数eg=0.6;wg=15.708;S0=0.001574;%按照书上的要求，取频率和时间的最大值和步长%频率间隔dw=0.3;%最大频率范围maxw=45;%最大时间值maxt=40;%时间间隔dt=0.2;%各层各时间点频率点的功率谱密度,循环变量：层数，时间点，频率点Pwt=zeros(n,maxt/dt,maxw/dw
4、);%频率点数循环变量wnwn=1;%对频率进行循环，求解各频率点的时间历程forw=0:dw:maxwx1=1+4*eg^2*(w/wg)^2;x2=(1-(w/wg)^2)^2+4*eg^2*(w/wg)^2;Sgw=x1*S0/x2;s=sqrt(Sgw);%采用精细积分法进行求解时间历程，得到位移和速度时程[disp,velp]=JINGXI67(M,zk,c,dt,maxt,w,s,n);Ywt=disp;forkkk=1:maxt/dt%求确定频率下各时间点的功率谱Yw=Ywt(:,kkk);每一时刻和频率点
5、的位移向量，对其进行求共轭和装置得到协方差矩阵，对角上的元素即是每一时刻的各层的功率谱y1=conj(Yw);%共轭y2=transpose(Yw);%转置矩阵%确定时间点确定频率下的功率谱Yw,取对角线元素Syyw=y1*y2;forkk=1:nPwt(kk,kkk,wn)=Syyw(kk,kk);endendwn=wn+1;end%求解完成后实际上wn为maxw/dw+2，实际频率点个数为maxw/dw+1%各层的时变方差，循环变量为：层数，时间点Fangcha=zeros(n,maxt/dt);fortn=1:ma
6、xt/dt%求解各层的时变方差forkk=1:nxx1=zeros(wn-1,1);%每一个时刻的方差对各频率点进行积分，频率点数取maxw/dw+1，即wn-1forwn0=1:wn-1xx1(wn0)=Pwt(kk,tn,wn0);end%采用复合梯形求积公式对功率谱进行积分得到方差Fangcha(kk,tn)=(xx1(1)+xx1(wn-1)+2*sum(xx1(2:wn-1-1)))*dw;endend%画图c1=(1:maxt/dt)*dt;d1=Fangcha(1,:)/S0;d2=Fangcha(2,:)
7、/S0;d3=Fangcha(3,:)/S0;d4=Fangcha(4,:)/S0;figure(3)plot(c1,d1,'k',c1,d2,'r',c1,d3,'m',c1,d4,'r-')精细积分的程序function[disp,velp]=JINGXI67(m,k,c,dt,maxt,w,s,n)%虚数单位II=sqrt(-1);%中的IIW=II*w;I=eye(n);Z=zeros(n);离散化n自由度结构受均匀调制演变随机激励时的运动微分方程可表示为：其中为平稳高斯白噪声随机过程向量，为调制函数。该方程可表
8、示为其中I为单位矩阵，记，，程序中的A即上面的矩阵A，AF即GA=[Z,I;-mk,-mc];AF=-1*[zeros(n,1);ones(n,1)];对线性体系来说，方程为线性向量常微分方程其解为对上式进行离散化为时间步长矩阵T可用数值方法精细算得采用函数EXPHDT计算，H即输入的矩阵A，dt即为时间间隔%采
