pem.rar_振型分解_振型分解法_虚拟激励_虚拟激励法

%-------------------------stationary poisson process----------- %-------10story frame structure unisolation--------% %****************************nonisolation*************************% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc close all; clear warning off damp=0.05; n=10;% m=2e4;%kg k=4e7;%n/m %地震动参数 w=0.0:0.1:30;%compute frequency dw=0.1; % s0=141.0438*1e-4;%8 度0.3g s0=1.5365*1e-4;%8 度0.035g duoyu wg=12.57; cg=0.72; S=(1+4*(cg*w/wg).^2)./((1-w.^2/wg^2).^2+4*(cg*w/wg).^2)*s0;%guo lv bai zhao sheng %%%-----------------非隔震--虚拟激励法---------------------%% M0=M0_made(m,n); K0=K0_made(k,n); [v0,d]=eig(K0,M0); w0=sort(diag(sqrt(d))); disp('非隔震结构周期') T=2*pi/w0(1) a0=2*damp*w0(1)*w0(2)/(w0(1)+w0(2));%瑞雷阻尼 b0=2*damp/(w0(1)+w0(2)); C0=a0*M0+b0*K0;%主体结构阻尼 A0=[zeros(n) eye(n);-inv(M0)*K0 -inv(M0)*C0]; B0=[zeros(n,1);-ones(n,1)]; CC0=[eye(n) zeros(n);zeros(n) eye(n);-inv(M0)*K0 -inv(M0)*C0]; D0=[zeros(n,1);zeros(n,1);-ones(n,1)]; t=0:0.02:10; gt=1; for j=1:length(w) ug1=sqrt(S(j))*gt.*cos(w(j).*t); ug2=sqrt(S(j))*gt.*sin(w(j).*t); %%%%%%%%抗震结构%%%%%%%%%%%5 [Y01]=lsim(A0,B0,CC0,D0,ug1,t); [Y02]=lsim(A0,B0,CC0,D0,ug2,t); dis0=(Y01(:,30)+ug1')+i*(Y02(:,30)+ug2');%顶层绝对加速度 Syy0=conj(dis0).*dis0; %形成Syy(w,t)矩阵 Sywt0(j,:)=Syy0; end %形成时变方差 disp('非隔震顶层绝对加速度方差-虚拟激励法') cigam0=2*sum(Sywt0(:,end))*dw ;%对列进行相加