中心差分法_matlab力学_中心差分法_源码

function FiniteDifferenceMethod_Example4_3 %中心差分法求解动力响应 clc;clear all; %%%%%%%第1步：输入K，M，P，C；初始条件计算 M=[2, 0, 0; 0, 1.5, 0; 0, 0, 1.0];%3x3的质量矩阵 K=[3000, -1200, 0; -1200, 1800, -600; 0, -600, 600]; %3x3的刚度矩阵 a0=0.4417; a1=0.0007; C=a0*M+a1*K;%采用Rayleigh阻尼 dt=0.01;%时间步长 N=1000;%计算时刻数 t=((0:1:N-1)*dt)';%时间系列 %假设在第2个自由度作用半周荷载，进行荷载离散化 nDof=length(M(:,1));%自由度 P=zeros(nDof,N);%荷载矩阵 for i=1:1:N if (i-1)*dt < 2 P(2, i) = 5*sin(10*(i-1)*dt); %假设在第2个自由度作用半周荷载 end end [V D]=eigs(K,M); %模态分析函数 DD=diag(D);%取对角阵构成特征值向量 [sorted_D sorted_loc]=sort(DD,1,'ascend'); %特征值由小到大排序 D=diag(sorted_D); V=V; sorted_V=V(:, sorted_loc);%振型向量的顺序要与特征值相对应 Omiga=sqrt(sorted_D) ;%圆频率 scaled_V=sorted_V; for i=1:1:length(Omiga) %振型对质量矩阵归一化 scaled_V(:,i)=scaled_V(:,i)/sqrt( scaled_V(:,i)'*M*scaled_V(:,i) ); end %正交性验证 Mnn=zeros(length(Omiga), length(Omiga));%定义关于质量矩阵的相关性矩阵 Knn=zeros(length(Omiga), length(Omiga));%定义关于刚度矩阵的相关性矩阵 for i=1:1:length(Omiga) for j=1:1:length(Omiga) Mnn(i,j)=sorted_V(:,i)'*M*sorted_V(:,j); Knn(i,j)=sorted_V(:,i)'*K*sorted_V(:,j); end end u0=[0.5; 0; 1.0]; v0=[2; 0.25; 0]; %初始位移、速度 a0=inv(M)*(P(:,1)-C*v0-K*u0);%初始加速度 u_1=u0-dt*v0+dt^2*a0/2;%-1个时刻的位移 %%%%%%%第2步：等效刚度、中心差分计算公式 K_eq=M/dt/dt+C/2/dt; a=K-2*M/dt/dt; b=M/dt/dt-C/dt; %%%%%%%第3、4步：计算当前时刻的运动 u=zeros(nDof,N);%位移反应矩阵 u(:,1)=u0; for i=2:1:N if i==2 P_eq=P(:,i-1)-a*u0-b*u_1; u(:,i)=inv(K_eq)*P_eq; else P_eq=P(:,i-1)-a*u(:,i-1)-b*u(:,i-2); u(:,i)=inv(K_eq)*P_eq; end end for i=1:1:nDof plot(t,u(i,:)) hold on end return; save('Knn_Mnn','Mnn','Knn'); Mnn [diag(Knn) sorted_D]; %比较振型刚度和特征值的大小 u=[0.5,1,1.5; 0.1,1.2, 0.75; 2.0, 4.2, 1.5]';%u=[1;1;1]; u是动力响应时程 q=inv(scaled_V)*u %%%%%%%%%%%%由两阶模态的频率和阻尼比一起确定Rayleigh阻尼矩阵 zeta=zeros(length(Omiga),1);%由第1、3阶来求解 zeta(1)=0.02; zeta(3)=0.02; Coef=[Omiga(3), -Omiga(1); -1/Omiga(3), 1/Omiga(1)]; %系数矩阵 a01=2*Omiga(1)*Omiga(3)/(Omiga(3)^2-Omiga(1)^2)*Coef*[zeta(1); zeta(3)] % a01=2*zeta(1)/(Omiga(1)+Omiga(3))*[Omiga(1)*Omiga(3);1] C_rayleigh=a01(1)*M+a01(2)*K Cn_cal=0*zeta;%振型阻尼 for i=1:1:length(Omiga) Cn_cal(i)=scaled_V(:,i)'*C_rayleigh*scaled_V(:,i); end Cn_cal zeta_cal=0.5*Cn_cal./Omiga./diag(Mnn) % inv(0.5*[1/Omiga(1), Omiga(1); 1/Omiga(3), Omiga(3)]) % 2*Omiga(1)*Omiga(3)/(Omiga(3)^2-Omiga(1)^2)*Coef % a01(1)*Mnn+a01(2)*Knn