可直接运行 基于Matlab的DQM方法代码，用于求解输液管道的失稳问题，在边界弹簧刚度变化的情况下，给出了临界失稳流速.rar

%微分求积法求解直管的流致振动 clc; clear; format long; a=0.002; tao=0; c=100; kt=10; u=4; N=10;%截点数目 beta=0.2; zx=ones(1,N); % 定义delat，截点坐标 delta = 0.0001; coord_x(1,1) = 0; coord_x(1,2) = delta; coord_x(1,N-1) = 1-delta; coord_x(1,N) = 1; %中间截点 for i=3:(N-2) coord_x(1,i)=(1/2)*(1-cos((i-2)/(N-3)*pi)); end %计算一阶权系数 for i=1:N for k=1:N if i~=k zx(1,i)=zx(1,i)*(coord_x(1,i)-coord_x(1,k)); end end end for i=1:N for j=1:N if i~=j wei_x_1(i,j)=zx(1,i)/zx(1,j)/(coord_x(1,i)-coord_x(1,j)); end end end for i=1:N for k=1:N if i~=k wei_x_1(i,i)=wei_x_1(i,i)-wei_x_1(i,k); end end end %计算二阶权系数 for i=1:N for j=1:N if i~=j wei_x_2(i,j)=2*(wei_x_1(i,i)*wei_x_1(i,j)-wei_x_1(i,j)/(coord_x(1,i)-coord_x(1,j))); end end end for i=1:N for k=1:N if i~=k wei_x_2(i,i)=wei_x_2(i,i)-wei_x_2(i,k); end end end %计算三阶权系数 for i=1:N for j=1:N if i~=j wei_x_3(i,j)=3*(wei_x_2(i,i)*wei_x_1(i,j)-wei_x_2(i,j)/(coord_x(1,i)-coord_x(1,j))); end end end for i=1:N for k=1:N if i~=k wei_x_3(i,i)=wei_x_3(i,i)-wei_x_3(i,k); end end end %计算四阶权系数 for i=1:N for j=1:N if i~=j wei_x_4(i,j)=4*(wei_x_3(i,i)*wei_x_1(i,j)-wei_x_3(i,j)/(coord_x(1,i)-coord_x(1,j))); end end end for i=1:N for k=1:N if i~=k wei_x_4(i,i)=wei_x_4(i,i)-wei_x_4(i,k); end end end %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ii=1:100 c=ii rm=-1 im=0 cc=0 for cc=1:100 u=0.1*cc K=zeros(N-4,N-4); G=zeros(N-4,N-4); M=zeros(N-4,N-4); %kbb kbb=[1,0,0,0;wei_x_1(2,1),wei_x_1(2,2),wei_x_1(2,N-1),wei_x_1(2,N);wei_x_2(N-1,1)+kt*wei_x_1(N-1,1),wei_x_2(N-1,2)+kt*wei_x_1(N-1,2),wei_x_2(N-1,N-1)+kt*wei_x_1(N-1,N-1),wei_x_2(N-1,N)+kt*wei_x_1(N-1,N);wei_x_3(N,1),wei_x_3(N,2),wei_x_3(N,N-1),(wei_x_3(N,N)-c)]; %kbb=[wei_x_3(1,1),wei_x_3(1,2),wei_x_3(1,N-1),(wei_x_3(1,N)-r);wei_x_2(2,1),wei_x_2(2,2),wei_x_2(2,N-1),wei_x_2(2,N);wei_x_2(N-1,1),wei_x_2(N-1,2),wei_x_2(N-1,N-1),wei_x_2(N-1,N);wei_x_3(N,1),wei_x_3(N,2),wei_x_3(N,N-1),(wei_x_3(N,N)-r)]; %kbd for j=1:N-4 kbd(1,j)=0; end for j=1:N-4 kbd(2,j)=wei_x_1(2,j+2); end for j=1:N-4 kbd(3,j)=wei_x_2(N-1,j+2)+kt*wei_x_1(N-1,j+2); end for j=1:N-4 kbd(4,j)=wei_x_3(N,j+2); end %kdb for i=1:N-4 kdb(i,1)=wei_x_4(i+2,1)+u^2*wei_x_2(i+2,1); end for i=1:N-4 kdb(i,2)=wei_x_4(i+2,2)+u^2*wei_x_2(i+2,2); end for i=1:N-4 kdb(i,3)=wei_x_4(i+2,N-1)+u^2*wei_x_2(i+2,N-1); end for i=1:N-4 kdb(i,4)=wei_x_4(i+2,N)+u^2*wei_x_2(i+2,N); end %kdd for i=1:N-4 for j=1:N-4 kdd(i,j)=wei_x_4(i+2,j+2)+u^2*wei_x_2(i+2,j+2); end end %Gdb for i=1:N-4 Gdb(i,1)=2*u*sqrt(beta)*wei_x_1(i+2,1)+a*wei_x_4(i+2,1); end for i=1:N-4 Gdb(i,2)=2*u*sqrt(beta)*wei_x_1(i+2,2)+a*wei_x_4(i+2,2); end for i=1:N-4 Gdb(i,3)=2*u*sqrt(beta)*wei_x_1(i+2,N-1)+a*wei_x_4(i+2,N-1); end for i=1:N-4 Gdb(i,4)=2*u*sqrt(beta)*wei_x_1(i+2,N)+a*wei_x_4(i+2,N); end %组装GDD矩阵 for i=1:N-4 for j=1:N-4 GDD(i,j)=2*u*sqrt(beta)*wei_x_1(i+2,j+2)+a*wei_x_4(i+2,j+2); end end %G G=-Gdb*inv(kbb)*kbd+GDD; %K K=-kdb*(inv(kbb))*kbd+kdd; %组装M矩阵 for i=1:N-4 M(i,i)=1; end MI=inv(M); MIK=MI*K; MIG=MI*G; R=zeros(2*(N-4),2*(N-4)); %组 装R for i=1:N-4 R(i+N-4,i)=1; end for i=1:N-4 for j=1:N-4 R(i,j)=(-1)*MIG(i,j); R(i,N-4+j)=(-1)*MIK(i,j); end end eigr=eig(R); [PIM,sor]=sort(real(eigr)); for j=1:2*(N-4) numb=sor(j); mmr(j,1)=eigr(numb,1); end mmr rm=real(mmr(2*(N-4),1)) im=imag(mmr(2*(N-4),1)) cc=cc+1 if rm>0&&im~=0 plot(c,u,'or') hold on elseif rm>0&&im==0 plot(c,u,'*b') hold on else plot(c,u,'+k') end end end