虚拟激励法程序求助.docVIP

虚拟激励法程序求助 hustdd 发表于: 2008-12-14 11:34 来源: 振动资讯 学习林家浩教授和张亚辉教授的专著《随机振动的虚拟激励法》时，编制的程序和书上的结果对不上，请大家帮忙看看。clear;clc;n=4;II=sqrt(-1);%主结构质量、阻尼、刚度矩阵M=eye(n)*1.0e+4;K=eye(n)*1.6*1.0e+7;%主结构刚度矩阵聚合zk=zeros(n);for j=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;%求解各阶模态频率[tzxl,tzz]=eig(k,m);d=diag(sqrt(tzz));%振型规一化for i=1:n? ?[tzz1(i),j]=min(d);? ?tzxl1(:,i)=tzxl(:,j);? ?d(j)=max(d)+1;end%振型归一化取第一层振型为1for j=1:n? ?tzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);endw0=tzz1;w=tzz1/(2*pi);zhx=tzxl1;广义阻尼矩阵? ?各阶模态阻尼比都取 %阻尼比ks0=0.05;ks=ones(n,1)*ks0;第n阶广义质量：? ? %求广义质量Mn=zhx*m*zhx;阻尼矩阵为： %求阻尼矩阵C=zeros(n);for i=1:n? ?C(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);%频率点数循环变量wnwn=1;%对频率进行循环，求解各频率点的时间历程for w=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;for kkk=1:maxt/dt%求确定频率下各时间点的功率谱Yw=Ywt(:,kkk);??每一时刻和频率点的位移向量，对其进行求共轭和装置得到协方差矩阵，对角上的元素即是每一时刻的各层的功率谱y1=conj(Yw);y2=transpose(Yw);%确定时间点确定频率下的功率谱Yw,取对角线元素Syyw=y1*y2;? ?for kk=1:n? ?Pwt(kk,kkk,wn)=Syyw(kk,kk);? ?endendwn=wn+1;end%求解完成后实际上wn为maxw/dw+2，实际频率点个数为maxw/dw+1%各层的时变方差，循环变量为：层数，时间点Fangcha=zeros(n,maxt/dt);for tn=1:maxt/dt%求解各层的时变方差for kk=1:nxx1=zeros(wn-1,1);%每一个时刻的方差对各频率点进行积分，频率点数取maxw/dw+1，即wn-1? ?for wn0=1:wn-1? ?xx1(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,:)/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