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

联动态特性法码研究端固支端悬臂液管道动态特。管道流固质比为，量纲流速为，最大值为悬臂输液管的无量纲控方程为

的，的围，时项YηXζ

ζ=1βu=14m/smax图一端固支臂输液示意图、问题分析法选择对于边支，边悬臂形式输液道，边条件为挠（向位对轴向位零导数与转（横向位对向位移一阶数）为，即

时，曲（横向移对向位移阶导）为，而曲率的空变化率（横位移对轴向位移三阶导数）为，即最近30年，微分求积法（DQM,Method）渐成为一种求解值和初值问题的有效方法逐渐被人们可广泛应用于工程结构的振动分之中方法的质是把函在给定网点处的导数值近似地用域上全部网点处数值的加权和表示NnNn对于输液管道模型,传统的作法是采用法将控方程进离散化,仅提取有限低阶振函数作近处理；DQM则服了法这局性本上虑所振的合献而且的施程免了的大量积运,其计精度效率好。所以文最选取方。、、权矩微求积的核心是将给定网格节点导数值由由网域上所有节点的函数值的加权和近似地表示表示此首需要对结构进网格离散一维的输管结构，划分成轴向分布一系列点。假为f

，自变量的范围0

，点划总为N，如所示的函需要指的，实的i求解过二的微方程，一在边点需入微量

以便于添系统的边界条件。点之外的网称为内部网点，可以采取不同的划分方式均匀划分与非均划分。

，

（边界点预设微小量）i

iN

，

i2,3,...,

（内部点均匀分）

i

2N

，

i2,3,...,

（内部点非均划分）图2输管网划示图加系的算法一由格日值多式导到。在

内假函数节点的数值利用插值近似为：f

ljj

j其lj

kkj

kjk

为拉朗日插多项式则一导数为NrrNrAAANrrNrAAAf

j

ljj将多项式形式转化为权系数矩阵形式ddx

f

xi

k

Af

，1,2,...,：阵将DQM格式应运，并组装矩计

1/2

u

at

at可得将上述关系合并，写成阵形式为了便求解程，一将边网点（括网点，b在，剩内矩阵内部元素位d调，所以KGGηI0其

为矩：111122223333111111122223333111222ANAN......N333AN,34......AN,

1

AAAAAAA

ANA

N

其中矩阵：0......

2,4AN由一行的系，可得到等式bbbd由此可界点位点位移阵T

bb

将阵带入到对内网点式子消除界网即有dbdbddddb

dbbddddddd

对于d

，从而得到ddd

假设解的形式，代入上面的动力学方程组得：有非零解的条是：Det对这广特问题以令yxx从得到0I][]Min][]利用得矩特征值由得到复数式运动频率据分析悬臂输管运动特性、所结绘运动响应图形（略）前绘分叉与图图臂液管的量纲频率的实部、虚部随量纲流速的变曲线

图悬液的图时模态的实频率变为0映的ModeDivergencyuu

时，输管的运动不稳，可能同时存第二阶或第三运动时，输管的运动不稳，可能同时存第一阶或第二运动时一模的频率变为0，反映的是失。时一模的频率变为0，可能同时存第二阶或三运动附：序源代clc;clear;%»®·ÖÍø¸ñ%¹ÜµÀ²ÎÁ¿beta=0.2;%Á÷¹ÌÖÊÁ¿±ÈΪ0.2umax=14;uu=0:0.1:umax;%Á÷ËÙ·¶Î§ll=length(uu);omega1=1:ll;omega2=1:ll;omega3=1:ll;formm=1:ll%u=2;%Á÷ËÙ%µÈ¼ä¸ô»®·ÖΪ30·ÝN=35;delta=1e-4;%¾ùÔÈ»®·Ö%xx=1:N;%forii=3:N-2%%end%xx(1)=0;%xx(2)=delta;%xx(N-1)=1-delta;%xx(N)=1;%·Ç¾ùÔÈ»®·Öxx=1:N;forii=3:N-2endxx(1)=0;xx(2)=delta;xx(N-1)=1-delta;xx(N)=1;%Çóµ¼¼ÓȨA¾ØÕóµÄ±àд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¾ØÕóÌî³ä¼ÆËã·½·¨%forii=1:N%%%%%%%%%%%%%%%%%%%%%%%%%%%%endforii=1:1:Nforj=1:1:Nifii==jfork=1:1:Nifk~=jendendelsefork=1:Nifk~=ii&&k~=jendendendendend%¸ß´ÎÇóµ¼¼ÓȨ¾ØÕóµÄÇó·¨AA2=AA1*AA1;AA3=AA2*AA1;AA4=AA3*AA1;%Óɱ߽çÌõ¼þÇó½âT¾ØÕ󣬴ÓÄÚ²¿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),AA3(N,N)];forii=1:N-4endT=B^(-1)*H;%Çó½âÕûÌåÐÔµÄK¡¢C¡¢M¾ØÕó%1¡¢K¾ØÕóÏÂϽËĸöС·Ö¿é¾ØÕóKbb=B;Kbd=-H;Kdb=zeros(N-4,4);forii=1:N-4endKdd=zeros(N-4,N-4);forii=1:N-4forjj=1:N-4endend%2¡¢C¾ØÕóÏÂϽËĸö·Ö¿é¾ØÕóCbb=zeros(4,4);Cbd=zeros(4,N-4);Cdb=zeros(N-4,4);Cdd=zeros(N-4,N-4);forii=1:N-4endforii=1:N-4forjj=1:N-4endend%2¡¢M¾ØÕóΪһ¸öµ¥Î»¾ØÕóM=eye(N-4,N-4);%ת»»K¡¢C¡¢M¾ØÕ󣨽«±ß½çµãµÄλÒÆÈ«²¿´ú»»ÎªÄÚ²¿µãµÄλÒÆ£©£¬¼´³ËÉÏ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-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);endfiguresubplot(1,2,1),plot(uu,omega1r,':',uu,omega2r,'-.',uu,omega3r,);h=legend('1stMode','2ndMode','3rdMode')set(h,'interpreter','none')title('Re(f)~flowvelocity')('flowvelocity/(m/s)')('frequency/Hz')axis([0umax080])gridonsubplot(1,2,2),plot(uu,omega1i,':',uu,omega2i,'-.',uu,omega3i,);h=legend('1stMode','2ndMode','3rdMode')set(h,'interpreter','none')title('Im(f)~flowvelocity')('flowvelocity/(m/s)')('frequency/Hz')axis([0umax-6060])gridonfigurellx=length(omega1r);holdonplot(omega1r,omega1i,'--',omega2r,omega2i,'-.',omega3r,omega3i,);h1=legend('1stMode','2ndMode','3rdMode')set(h,'interpreter','none')title('Im(f)~Re(f)velocity')('Re(f)/Hz')('Im(f)/Hz')%axis([040-6060])gridonforii=1:llx/10str1=[num2str(uu(ii*10))];%str2=['('num2str(omega2r(ii))','num2str(omega2i(ii))')'];%str3=['('num2