zhongxinchafen.rar_中心差分法_中心差分法matlab_振动 微分_振动特性_结构振动

K=[4 -2 0 0 0 0; -2 4 -2 0 0 0; 0 -2 4 -2 0 0; 0 0 -2 4 -2 0; 0 0 0 -2 4 -2; 0 0 0 0 -2 2]*1e+8; %刚度矩阵 M=4e5*eye(6); %质量矩阵 [fi,o2]=eig(M^-1*K) %自振频率 w=[sqrt(o2(1,1));sqrt(o2(2,2));sqrt(o2(3,3));sqrt(o2(4,4));sqrt(o2(5,5));sqrt(o2(6,6))] A01=2*0.05/(sqrt(o2(1,1))+sqrt(o2(2,2)))*[sqrt(o2(1,1)*o2(2,2));1]; %Rayleigh阻尼的常数a0,a1 C=A01(1,1)*M+A01(2,1)*K; %Rayleigh阻尼矩阵 t=0.02; %步长 kk=M/(t^2)+C/(2*t); a=K-2*M/(t^2); %中心差分法计算 b=M/(t^2)-C/(2*t); p=zeros(6,2503); pp=zeros(6,2503); u=zeros(6,2504); %位移计算 for j=1:2503 p(:,j)=-4e5*e(2,j); if j==1 pp(:,j)=p(:,j); else pp(:,j)=p(:,j)-a*u(:,j)-b*u(:,j-1); end u(:,j+1)=kk^-1*pp(:,j); end umax=max(abs(u'))' %各层最大位移输出 %绘制位移图 for i=1:6 subplot(3,2,i) plot(e(1,1:2502),u(i,1:2502),'m') xlabel('时间 (s)'); ylabel('位移 (m)'); title(['第' num2str(i) '层时间-位移图']); grid on end % 速度计算 sd=zeros(6,2503); for j=1:2503 if j==1 sd(:,j)=(u(:,j+1))/2/t; else sd(:,j)=(u(:,j+1)-u(:,j-1))/2/t; end end sdmax=max(abs(sd'))' %各层最大速度输出 %绘制速度图 for i=1:6 subplot(3,2,i) plot(e(1,1:2502),sd(i,1:2502),'m') xlabel('时间 (s)'); ylabel('速度 (m/s)'); title(['第' num2str(i) '层时间-速度图']); grid on end %加速度计算 jsd=zeros(6,2503); for j=1:2502 if j==1 jsd(:,j)=(u(:,j+1))/t/t; else jsd(:,j)=(u(:,j+1)-2*u(:,j)+u(:,j-1))/t/t; end end jsdmax=max(abs(jsd'))' %各层最大加速度输出 %绘制加速度 for i=1:6 subplot(3,2,i) plot(e(1,1:2502),jsd(i,1:2502),'m') xlabel('时间 (s)'); ylabel('加速度 (m/s2)'); title(['第' num2str(i) '层时间-加速度图']); grid on end