function [SY2]=PEM(K,M,C,R,w,Sa);
%使用虚拟激励法求响应的自谱密度；
%输入变量为：
%K,M,C分别为结构的刚度、质量和阻尼矩阵；
%R为激励的位置和幅值向量；
%Sa为激励的自谱密度；
%输出变量为SY2表示响应的谱密度矩阵。
cn=length(K);
[mod,w2]=eig(K,M); %求解结构的振型
Kst=mod'*K*mod; %确定结构的广义质量、刚度阻尼矩阵
Mst=mod'*M*mod;
Cst=mod'*C*mod;
gama=mod'*R./diag(Mst); %求振型参与系数
damp=Cst./sqrt(Kst.*Mst)/2; %求各振型阻尼比
for j=1:cn %计算结构频率响应函数
r=w./sqrt(w2(j,j));
Hw(j,:)=1./w2(j,j)./((1-r.^2)+sqrt(-1)*2*damp(j)*r);
end
y=zeros(cn,length(w)); %求虚拟响应
for jj=1:cn
for kk=1:cn
y(jj,:)=y(jj,:)+gama(kk).*Hw(kk,:).*sqrt(Sa).*mod(jj,kk);
end
end
for jj=1:cn
for kk=1:cn
Sy(jj,kk,:)=conj(y(kk,:)).*y(jj,:);
end
end
for ii=1:cn
for jj=1:cn
SY2(ii,jj)=trapz(w,Sy(ii,jj,:));
end
end
end
评论0