華科流固耦聯(lián)懸臂輸液管大作業(yè)微分求積法詳細(xì)推導(dǎo)與源代碼

1、華中科技大學(xué)流固耦聯(lián) 懸臂輸液管動態(tài)特性 微分求積法 附源代碼研究一端固支一端懸臂輸液管道的動態(tài)特性。管道的流固質(zhì)量比為，無量綱流速為，最大值為14，懸臂輸液管道的無量綱控制方程為為管道的橫向位移，為無量綱的軸向位置，范圍為01，為時(shí)間項(xiàng) 圖1 一端固支一端懸臂輸液管模型示意圖1、問題分析與方法選擇對于一邊固支，一邊懸臂形式的輸液管道，邊界條件為當(dāng)時(shí)撓度（橫向位移對軸向位移零階導(dǎo)數(shù)）與轉(zhuǎn)角（橫向位移對軸向位移一階導(dǎo)數(shù)）均為0，即 當(dāng)時(shí)，曲率（橫向位移對軸向位移二階導(dǎo)數(shù)）為0，而且曲率的空間變化率（橫向位移對軸向位移三階導(dǎo)數(shù)）為0，即最近 30 年，微分求積法（DQM, Differential

2、 Quadrature Method）逐漸成為一種求解邊值和初值問題的有效方法，逐漸被人們認(rèn)可，廣泛地應(yīng)用于工程結(jié)構(gòu)的振動分析之中。該方法的本質(zhì)是把函數(shù)在給定網(wǎng)點(diǎn)處的導(dǎo)數(shù)值近似地用域上全部網(wǎng)點(diǎn)處函數(shù)值的加權(quán)和表示對于輸液管道模型,傳統(tǒng)的作法是采用Galekrin法將控制方程進(jìn)行離散化,僅提取有限的低階振型函數(shù)作近似處理；DQM則克服了Galekrin法的這種局限性,本質(zhì)上考慮了所有振型的綜合貢獻(xiàn)。而且DQM的實(shí)施過程避免了Galekrin的大量積分運(yùn)算,其計(jì)算精度和效率較好。所以本文最終選取DQM方法。2、結(jié)構(gòu)離散、計(jì)算權(quán)系數(shù)矩陣微分求積法的核心是將給定網(wǎng)格節(jié)點(diǎn)的導(dǎo)數(shù)值由由網(wǎng)格域上所有節(jié)點(diǎn)的函

3、數(shù)值的加權(quán)和近似地表示表示。為此，首先需要對結(jié)構(gòu)進(jìn)行網(wǎng)格離散，對于一維的輸液管結(jié)構(gòu)，可以劃分成軸向分布的一系列網(wǎng)點(diǎn)。假定待求函數(shù)為，自變量的范圍為，網(wǎng)點(diǎn)的劃分總個(gè)數(shù)為N，如下圖所示，各網(wǎng)點(diǎn)的為待求函數(shù)。需要指出的是，在實(shí)用傳統(tǒng)的DQM求解超過二階的微分方程時(shí)，一般在邊界點(diǎn)需引入微小量（），以便于添加系統(tǒng)的邊界條件。邊界點(diǎn)之外的網(wǎng)點(diǎn)稱為內(nèi)部網(wǎng)點(diǎn)，可以采取不同的劃分方式：均勻劃分與非均勻劃分。，（邊界網(wǎng)點(diǎn)預(yù)設(shè)微小量）， （內(nèi)部網(wǎng)點(diǎn)均勻劃分）， （內(nèi)部網(wǎng)點(diǎn)非均勻劃分）圖2 輸液管網(wǎng)點(diǎn)劃分示意圖加權(quán)系數(shù)的計(jì)算方法一般由拉格朗日插值多項(xiàng)式推導(dǎo)得到。在區(qū)間內(nèi)，假設(shè)函數(shù)的n個(gè)節(jié)點(diǎn)的函數(shù)值利用Lagrange

4、插值近似為：其中為拉格朗日插值多項(xiàng)式則一階導(dǎo)數(shù)為將多項(xiàng)式形式轉(zhuǎn)化為權(quán)系數(shù)矩陣形式，權(quán)系數(shù)矩陣有如下規(guī)則：通過矩陣連乘得到更高階的權(quán)系數(shù)矩陣3、將DQM格式應(yīng)用于輸液管的運(yùn)動方程和邊界條件，并組裝矩陣，計(jì)算求解 at at 可得將上述關(guān)系合并，寫成矩陣形式為了方便求解方程，一般將邊界網(wǎng)點(diǎn)（包括網(wǎng)點(diǎn)1,2,N-1,N，用表示）放在一起，剩余為內(nèi)部場點(diǎn)（用表示），并相應(yīng)地對矩陣內(nèi)部元素的位置進(jìn)行調(diào)整，所以其中為4X4矩陣：其中為4X(N-4)矩陣：由第一行的關(guān)系，可以得到等式由此可以得到邊界點(diǎn)位移到內(nèi)部點(diǎn)位移的轉(zhuǎn)移矩陣T將矩陣T帶入到針對內(nèi)部網(wǎng)點(diǎn)的式子，消除邊界網(wǎng)點(diǎn)即有對于，從而得到假設(shè)解的形式，代

5、入上面的動力學(xué)方程組可得： 有非零解的條件是： 對于這種廣義特征值問題，可以令，從而得到利用MATLAB得到矩陣特征值。由此得到復(fù)數(shù)形式的運(yùn)動頻率。再據(jù)此分析懸臂輸液管的運(yùn)動特性3、分析所得結(jié)果，繪制運(yùn)動響應(yīng)圖形（略）分析前三階結(jié)果，繪制分叉圖與Argand圖圖3 懸臂輸液管的無量綱頻率的實(shí)部、虛部隨無量綱流速的變化曲線圖4 懸臂輸液管的Argand圖時(shí)，一階模態(tài)的實(shí)頻率變?yōu)?，反映的是失穩(wěn)。（First Mode Divergency）時(shí)，輸液管的運(yùn)動不穩(wěn)定，可能同時(shí)存在第二階或第三階運(yùn)動時(shí)，輸液管的運(yùn)動不穩(wěn)定，可能同時(shí)存在第一階或第二階運(yùn)動時(shí)，一階模態(tài)的實(shí)頻率變?yōu)?，反映的是失穩(wěn)。時(shí)，一階

6、模態(tài)的實(shí)頻率變?yōu)?，可能同時(shí)存在第二階或第三階運(yùn)動附：程序源代碼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