华科流固耦联悬臂输液管大作业微分求积法详细推导与源代码

1、华中科技大学流固耦联 悬臂输液管动态特性 微分求积法 附源代码研究一端固支一端悬臂输液管道的动态特性。管道的流固质量比为，无量纲流速为，最大值为14，悬臂输液管道的无量纲控制方程为为管道的横向位移，为无量纲的轴向位置，范围为01，为时间项 图1 一端固支一端悬臂输液管模型示意图1、问题分析与方法选择对于一边固支，一边悬臂形式的输液管道，边界条件为当时挠度（横向位移对轴向位移零阶导数）与转角（横向位移对轴向位移一阶导数）均为0，即 当时，曲率（横向位移对轴向位移二阶导数）为0，而且曲率的空间变化率（横向位移对轴向位移三阶导数）为0，即最近 30 年，微分求积法（DQM, Differential

2、 Quadrature Method）逐渐成为一种求解边值和初值问题的有效方法，逐渐被人们认可，广泛地应用于工程结构的振动分析之中。该方法的本质是把函数在给定网点处的导数值近似地用域上全部网点处函数值的加权和表示对于输液管道模型,传统的作法是采用Galekrin法将控制方程进行离散化,仅提取有限的低阶振型函数作近似处理；DQM则克服了Galekrin法的这种局限性,本质上考虑了所有振型的综合贡献。而且DQM的实施过程避免了Galekrin的大量积分运算,其计算精度和效率较好。所以本文最终选取DQM方法。2、结构离散、计算权系数矩阵微分求积法的核心是将给定网格节点的导数值由由网格域上所有节点的函

3、数值的加权和近似地表示表示。为此，首先需要对结构进行网格离散，对于一维的输液管结构，可以划分成轴向分布的一系列网点。假定待求函数为，自变量的范围为，网点的划分总个数为N，如下图所示，各网点的为待求函数。需要指出的是，在实用传统的DQM求解超过二阶的微分方程时，一般在边界点需引入微小量（），以便于添加系统的边界条件。边界点之外的网点称为内部网点，可以采取不同的划分方式：均匀划分与非均匀划分。，（边界网点预设微小量）， （内部网点均匀划分）， （内部网点非均匀划分）图2 输液管网点划分示意图加权系数的计算方法一般由拉格朗日插值多项式推导得到。在区间内，假设函数的n个节点的函数值利用Lagrange

4、插值近似为：其中为拉格朗日插值多项式则一阶导数为将多项式形式转化为权系数矩阵形式，权系数矩阵有如下规则：通过矩阵连乘得到更高阶的权系数矩阵3、将DQM格式应用于输液管的运动方程和边界条件，并组装矩阵，计算求解 at at 可得将上述关系合并，写成矩阵形式为了方便求解方程，一般将边界网点（包括网点1,2,N-1,N，用表示）放在一起，剩余为内部场点（用表示），并相应地对矩阵内部元素的位置进行调整，所以其中为4X4矩阵：其中为4X(N-4)矩阵：由第一行的关系，可以得到等式由此可以得到边界点位移到内部点位移的转移矩阵T将矩阵T带入到针对内部网点的式子，消除边界网点即有对于，从而得到假设解的形式，代

5、入上面的动力学方程组可得： 有非零解的条件是： 对于这种广义特征值问题，可以令，从而得到利用MATLAB得到矩阵特征值。由此得到复数形式的运动频率。再据此分析悬臂输液管的运动特性3、分析所得结果，绘制运动响应图形（略）分析前三阶结果，绘制分叉图与Argand图图3 悬臂输液管的无量纲频率的实部、虚部随无量纲流速的变化曲线图4 悬臂输液管的Argand图时，一阶模态的实频率变为0，反映的是失稳。（First Mode Divergency）时，输液管的运动不稳定，可能同时存在第二阶或第三阶运动时，输液管的运动不稳定，可能同时存在第一阶或第二阶运动时，一阶模态的实频率变为0，反映的是失稳。时，一阶

6、模态的实频率变为0，可能同时存在第二阶或第三阶运动附：程序源代码clc;clear;%»®·ÖÍø¸ñ %¹ÜµÀ²ÎÁ¿beta=0.2;%Á÷¹ÌÖÊÁ¿±ÈΪ0.2umax=14;uu=0:0.1:umax;%Á÷ËÙ·¶Î§ll=

7、length(uu);omega1=1:ll;omega2=1:ll;omega3=1:ll; for mm=1:ll u=uu(mm); %u=2;%Á÷ËÙ %µÈ¼ä¸ô»®·ÖΪ30·ÝN=35;delta=1e-4;%¾ùÔÈ»®·Ö% xx=1:N;% for ii=3:N-2% xx(ii)=(ii-2)/(N-3);%

8、end% xx(1)=0;% xx(2)=delta;% xx(N-1)=1-delta;% xx(N)=1; %·Ç¾ùÔÈ»®·Öxx=1:N;for ii=3:N-2 xx(ii)=0.5*(1-cos(ii-2)/(N-3)*pi);endxx(1)=0;xx(2)=delta;xx(N-1)=1-delta;xx(N)=1; %Çóµ¼¼ÓȨA¾ØÕóµ&

9、#196;±àдAA1=zeros(N,N);AA2=zeros(N,N);AA3=zeros(N,N);AA4=zeros(N,N);KK=zeros(N,N);CC=zeros(N,N);MM=zeros(N,N);%A¾ØÕóµÄ±àд %»»Ò»ÖÖÕýÈ·µÄA¾ØÕóÌ&#

10、238;³ä¼ÆËã·½·¨% for ii=1:N% for jj=1:N% if ii=jj % for kk=1:N% m1=0;% if kk=ii;% h=0;% else h=1/(xx(ii)-xx(kk);% end% m1=m1+h;% end% AA1(ii,jj)=m1;% else% for kk=1:N% m2=1;% if kk=ii% h=1;% else if kk=jj% h=1;% else% h=(xx(ii)-xx(kk)/(xx(jj)-xx(kk);%

11、 end % end% m2=m2*h;% end% AA1(ii,jj)=m2;% end% end% endfor ii=1:1:N for j=1:1:N if ii=j a=0; for k=1:1:N if k=j a=a+1/(xx(j)-xx(k); end end AA1(ii,j)=a; else b=1; for k=1:N if k=ii&&k=j b=b*(xx(ii)-xx(k)/(xx(j)-xx(k); end end AA1(ii,j)=b/(xx(j)-xx(ii); end endend %¸ß´Î&#

12、199;óµ¼¼ÓȨ¾ØÕóµÄÇó·¨AA2=AA1*AA1;AA3=AA2*AA1;AA4=AA3*AA1; %Óɱ߽çÌõ¼þÇó½âT¾ØÕ󣬴ÓÄÚ²¿

13、;N-4¸öµãµ½Íⲿ4¸öµãµÄ¾ØÕóת»¯B=zeros(4,4);H=zeros(4,N-4); B=1,0,0,0;AA1(2,1),AA1(2,2),AA1(2,N-1),AA1(2,N);AA2(N-1,1),AA2(N-1,2),AA2(N-1,N-1),AA2(N-1,N);AA3(N,1),AA3(N,2),AA3(N,N-1)

14、,AA3(N,N);for ii=1:N-4 H(1,ii)=0; H(2,ii)=-AA1(2,ii+2); H(3,ii)=-AA2(N-1,ii+2); H(4,ii)=-AA3(N,ii+2);endT=B(-1)*H; %Çó½âÕûÌåÐÔµÄK¡¢C¡¢M¾ØÕó%1¡¢K¾ØÕóÏÂÏ

15、9;ËĸöС·Ö¿é¾ØÕóKbb=B;Kbd=-H;Kdb=zeros(N-4,4);for ii=1:N-4 Kdb(ii,1)=AA4(ii+2,1)+u2*AA2(ii+2,1); Kdb(ii,2)=AA4(ii+2,2)+u2*AA2(ii+2,2); Kdb(ii,3)=AA4(ii+2,N-1)+u2*AA2(ii+2,N-1); Kdb(ii,4)=AA4(ii+2,N)+u2*AA2(ii+2,N);endKdd=zeros(N-4

16、,N-4);for ii=1:N-4 for jj=1:N-4 Kdd(ii,jj)=AA4(ii+2,jj+2)+u2*AA2(ii+2,jj+2); endend %2¡¢C¾ØÕóÏÂϽËĸö·Ö¿é¾ØÕóCbb=zeros(4,4);Cbd=zeros(4,N-4);Cdb=zeros(N-4,4);Cdd=zeros(N-4,N-4); for ii=1:N-4

17、 Cdb(ii,1)=2*beta0.5*u*AA1(ii+2,1); Cdb(ii,2)=2*beta0.5*u*AA1(ii+2,2); Cdb(ii,3)=2*beta0.5*u*AA1(ii+2,N-1); Cdb(ii,4)=2*beta0.5*u*AA1(ii+2,N);end for ii=1:N-4 for jj=1:N-4 Cdd(ii,jj)=2*beta0.5*u*AA1(ii+2,jj+2); endend %2¡¢M¾ØÕóΪһ¸öµ&

18、#165;λ¾ØÕóM=eye(N-4,N-4); %ת»»K¡¢C¡¢M¾ØÕ󣨽«±ß½çµãµÄλÒÆÈ«²¿´ú»»ÎªÄÚ

19、;²¿µãµÄλÒÆ£©£¬¼´³ËÉÏT¾ØÕóKT=Kdb*T+Kdd;CT=Cdb*T+Cdd;MT=M; det1=zeros(N-4,N-4);det2=eye(N-4);det3=-(MT-1)*KT;det4=-(MT-1)*CT; DET=det1,det2;det3,det4;V,D=eig(DET); omega1(mm)=D(2*N-

20、8,2*N-8)/1i;omega1r=abs(real(omega1);omega1i=imag(omega1); omega2(mm)=D(2*N-10,2*N-10)/1i;omega2r=abs(real(omega2);omega2i=imag(omega2); omega3(mm)=D(2*N-12,2*N-12)/1i;omega3r=abs(real(omega3);omega3i=imag(omega3);end figuresubplot(1,2,1),plot(uu,omega1r,':',uu,omega2r,'-.',uu,omega3

21、r,'-');h=legend('1st Mode','2nd Mode','3rd Mode')set(h,'interpreter','none')title('Re(f)flow velocity')xlabel('flow velocity/(m/s)')ylabel('frequency/Hz')axis(0 umax 0 80)grid on subplot(1,2,2),plot(uu,omega1i,':',uu,ome

22、ga2i,'-.',uu,omega3i,'-');h=legend('1st Mode','2nd Mode','3rd Mode')set(h,'interpreter','none')title('Im(f)flow velocity')xlabel('flow velocity/(m/s)')ylabel('frequency/Hz')axis(0 umax -60 60)grid on figurellx=length(omega1r);hold onplot(o