流_固耦合問(wèn)題邊界元_有限元耦合方法分析

1、第13卷第2期石家莊鐵道學(xué)院學(xué)報(bào)V o l.13N o.2 2000年6月JOU RNAL O F SH I J I A ZHUAN G RA I L WA Y I N ST ITU TE Jun.2000流-固耦合問(wèn)題邊界元-有限元耦合方法分析張效松1葉天麒2(石家莊鐵道學(xué)院基礎(chǔ)部1石家莊050043西北工業(yè)大學(xué)飛機(jī)系2西安710072【摘要】利用邊界元-有限元耦合方法對(duì)流-固耦合振動(dòng)問(wèn)題進(jìn)行了分析,假設(shè)流體控制方程為L(zhǎng) ap lace方程,利用非連續(xù)邊界元對(duì)流體域進(jìn)行離散,從而有效地解決了邊界元分析中的“角點(diǎn)效應(yīng)”問(wèn)題。固體以平面梁為模型,采用有限元進(jìn)行離散,對(duì)非連續(xù)邊界元和有限元的耦合問(wèn)

2、題進(jìn)行了分析,通過(guò)對(duì)懸臂梁在一側(cè)受液體作用時(shí)的瞬態(tài)響應(yīng)分析的數(shù)值解同解析解的比較,表明了本文所給方法的有效性,同時(shí)為利用邊界元-有限元耦合方法對(duì)流-固耦合分析奠定了基礎(chǔ)?！娟P(guān)鍵詞】邊界元-有限元耦合流-固耦合非連續(xù)元角點(diǎn)效應(yīng)【分類(lèi)號(hào)】TB115用B E M2FE M耦合方法分析流-固耦合問(wèn)題引起了一些研究人員的關(guān)視。E sto rff和A n tes1將FE M2B E M耦合方法用于流-固耦合問(wèn)題瞬態(tài)響應(yīng)分析,得到了同邊界元分析一致的計(jì)算結(jié)果,但在耦合方法分析中,他們沒(méi)有考慮對(duì)邊界元分析精度有較大影響的“角點(diǎn)效應(yīng)”問(wèn)題。T sai等2應(yīng)用FE M-B E M耦合方法對(duì)流-固耦合問(wèn)題的固有振動(dòng)

3、進(jìn)行了分析,他們的研究也沒(méi)有對(duì)邊界元離散的“角點(diǎn)效應(yīng)”給予足夠的重視?；谌缟戏治?本文利用FE M2B E M耦合方法對(duì)流-固耦合問(wèn)題進(jìn)行了分析,為了有效地解決邊界元分析中的“角點(diǎn)效應(yīng)”問(wèn)題,采用非連續(xù)元離散邊界積分方程,推導(dǎo)了有限元同非連續(xù)邊界元耦合的公式,并對(duì)梁在一側(cè)有流體作用的流-固耦合問(wèn)題進(jìn)行了數(shù)值實(shí)施,同時(shí)為了驗(yàn)證數(shù)值算法的精度和有效性,對(duì)梁在一側(cè)有流體作用時(shí)的解析解進(jìn)行了推導(dǎo),通過(guò)數(shù)值結(jié)果同解析解的比較,表明了本文所給方法的有效性。1流一固耦合問(wèn)題基本方程及解析解分析如圖1所示,懸臂梁在一側(cè)受有無(wú)粘、無(wú)旋、不可壓縮的理想液體作用,則流體的控制方程為52P 5x2+52P5z2=0

4、(x,z85P5x x=0=-52w5t2;P z=H=05P5z z=0=0,5Px x=L=0(1式中,P為流體液動(dòng)壓力,w為梁的撓度,為流體的密度。利用分離變量技術(shù)可得流體液動(dòng)壓力的表達(dá)式為收稿日期1999206229張效松男1966年10月出生副教授圖1流-固耦合系統(tǒng)P (x ,z ,t =6n =0-2n H (E n -1(E n e -n x+e-n xh52w5t 2co s n z d z co s n z (2式中,n =(2n +12H;E n =e -2n L.計(jì)及液體的液動(dòng)壓力后,梁振動(dòng)的控制方程為m 52w 5t 2+E I 54w5z4=q (z ,t -P (x

5、 ,z x =0(3(3式為一微分積分方程,其求解比較困難,下面利用文獻(xiàn)3提出的方法進(jìn)行求解,即利用無(wú)水梁的振型展開(kāi)有水梁的撓度設(shè)w (z ,t =6n =1A i (t Y i (z (4其中,Y i (z 為無(wú)水懸臂梁的第i 階固有振型,其具體表達(dá)式為Y i (z =1(co s k i z -ch k i z +(sin k i z -sh k i z (5式中,i =-sin k i H +sh k i Hco s k i H +ch k i H。將(4式代入(3式,經(jīng)過(guò)一系列推導(dǎo)可得如下控制方程m 20s D s A s (t +m D s A s (t =Q s (t -6i =1

6、6n =0-2(B n +1n H (B n -1D i B in D s B sn A i (t (6式中,20s為無(wú)水梁的第s 階固有頻率,D s 和B sn 的具體表達(dá)式可參見(jiàn)文獻(xiàn)3,Q s (t =Hq (z ,t Y s(z d z .2非連續(xù)邊界元一有限元耦合分析對(duì)結(jié)構(gòu)域采用有限元方法進(jìn)行離散,相應(yīng)的有限元方程為M +K =F P -F W (7式中,F P 為等效外載荷向量;F W 為等效液動(dòng)壓力向量。對(duì)流體控制微分方程采用非連續(xù)邊界元進(jìn)行離散,具體離散方式可參見(jiàn)文獻(xiàn)4,將離散后的線(xiàn)性代數(shù)方程組根據(jù)是否同結(jié)構(gòu)相接觸對(duì)自由度進(jìn)行重新排序后可得H N N H N C HN CHCCP

7、 N P C=G N N G N C G CNG CC5P N 5n 5P C 5n(8式中,含下標(biāo)C 的量為同固體接觸自由度,含下標(biāo)N 的量為不同結(jié)構(gòu)接觸自由度。對(duì)(8進(jìn)行邊界條件處理,并將不同結(jié)構(gòu)相接觸的自由度進(jìn)行縮聚后可得V 1P C =V 25P C 5n +F C (9式中,V 1,V 2和F C 的具體表達(dá)式和計(jì)算方法可參閱文獻(xiàn)3。將(9左乘V 1-1后,可得P C =V125P C 5n +F 12(10式中,V12=V 1-1V 2;F 12=V 1-1F C 。將液動(dòng)壓力分布列陣P C 轉(zhuǎn)換為等效結(jié)點(diǎn)力向量F i ,假定轉(zhuǎn)換矩陣為T(mén) ,則有F W =T P C =T V125

8、P C 5n +T F 12(11因?yàn)椴捎梅沁B續(xù)邊界元離散流體域控制方程,所以配位點(diǎn)不是元素的端點(diǎn),為了實(shí)現(xiàn)B E M 和FE M 耦合,需進(jìn)行相應(yīng)的變換,變換的具體形式如下5P C 5n =R 5P E 5n (12式中,5P C 5n 和5P E5n分別為配位點(diǎn)和端點(diǎn)處液動(dòng)壓力的法向?qū)?shù),R 為轉(zhuǎn)換矩陣,其具體形式為7第2期張效松等:流-固耦合問(wèn)題邊界元-有限元耦合方法分析R =R 1R 20+00R 1R 2+0+000+R 2(13對(duì)于平面線(xiàn)性非連續(xù)邊界元R 1=12(1+12(1-T ;R 2=12(1-12(1+T 。式中,為非連續(xù)元的配位因子。考慮流-固交界面處的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)條

9、件后可得流-固耦合系統(tǒng)的振動(dòng)方程為(M a +M +K =F p -F F (14式中,M a =T V 12R ;F F =T F 12。由(14式可以看出,當(dāng)考慮流體作用后,梁有限元質(zhì)量矩陣增加了一項(xiàng)附加質(zhì)量矩陣M a 。借助(14式可對(duì)流-固耦合系統(tǒng)進(jìn)行固有特性和瞬態(tài)響應(yīng)分析。在實(shí)現(xiàn)B E M -FE M 耦合過(guò)程中,一個(gè)關(guān)鍵問(wèn)題是完成元素分布力向等效結(jié)點(diǎn)力轉(zhuǎn)換的轉(zhuǎn)換矩陣,以梁元素為例,對(duì)流-固交界面上的元素,利用虛功原理,可得元素分布矩陣的具體形式為T(mén) e =L 2a 1a 3a 2-a 4a 2a 4a 1-a 3(15式中,a 1=12+15,a 2=12-15,a 3=L 4(1

10、3-116,a 4=L 4(13+116;L 為梁元素的長(zhǎng)度。對(duì)于總體分布矩陣,需考慮邊界元液動(dòng)壓力的累積效應(yīng),疊加后可得總體分布矩陣的表達(dá)式為T(mén) =6N T E Fi =1T e i =A10+0A2A 1+00A2+00+A 2 (16式中,A 1=L2a 1a 3a 2a 4,A 2=L2a 3-a 4a 4-a 3,N T E F 為流-固耦合交界面單元數(shù)目.3數(shù)值算例及結(jié)果分析對(duì)流-固耦合系統(tǒng)瞬態(tài)響應(yīng)編制了相應(yīng)的Fo rtran 程序,下面給出具體算例和結(jié)果。如圖1所示,懸臂梁在一側(cè)有液體作用,梁的抗彎剛度E I =110×1011Pa m 2,質(zhì)量密度s =2450kg圖

11、2載荷(1情況下梁中跨撓度響應(yīng)曲線(xiàn)圖3載荷(2情況下梁中跨撓度響應(yīng)曲線(xiàn)8石家莊鐵道學(xué)院學(xué)報(bào)第13卷m 3,橫截面積A =110m 2,梁的高度H =L =40m ,流體密度f(wàn) =1000kg m 3,考慮兩種外載荷情況下耦合系統(tǒng)的瞬態(tài)響應(yīng):(1在梁端部受突加外載P =20N 作用;(2在梁端部受外載P =20sin (10t N 作用。對(duì)耦合系統(tǒng)的控制方程采用N ewm ark 2直接積分法進(jìn)行計(jì)算,在計(jì)算中取=015,=0125,圖2和圖3分別給出了兩種載荷作用情況下梁中點(diǎn)的撓度響應(yīng)曲線(xiàn),從圖中所示結(jié)果可以看出本文結(jié)果同解析解吻合良好。4結(jié)束語(yǔ)利用解析方法和B E M 2FE M 耦合方法,

12、以梁在一側(cè)受有液體作用為分析對(duì)象,對(duì)流-固耦合問(wèn)題進(jìn)行了分析,計(jì)算過(guò)程中采用非連續(xù)元離散流體域邊界積分方程,有效地解決了邊界元分析中由于存在“角點(diǎn)效應(yīng)”而產(chǎn)生的誤差,通過(guò)具體的數(shù)值算例和解析結(jié)果的比較可以看出本文所給方法的有效性。參考文獻(xiàn)1E sto rff O .V on ,A ntes H .O n FE M 2BE M coup ing fo r fluid 2structrue interacti on analysis in ti m e dom ain .Int .J .N um .M eth .Engng .1991,32:115111672T sai C .S .,L ee G

13、 .C .,Ketter R .L .So luti on of dam 2reservo ir interacti on p roblem using com bined FE M ,BE M w ithparticular integrals ,modal analysis and substructuring ,Engineering A nalysis w ith Boundary E lem ents .1992,9:219232Fluid -structure I n teraction Analysis by BE M -FE MCoupli ng ProcedureZhang

14、X i aosong 1Y e T i anq i2(D epartm ent of Basic Courses ,Sh ijiazhuang R ail w ay Institute 1,Sh ijiazhuang 050043D epartm ent of A ircraft Engineering ,N o rthw estern Po lytechnical U niversity 2,X i an ,710072【Abstract 】In th is p ap er ,B E M 2FE M coup ling p rocedu re is develop ed fo r flu i

15、d 2structu re in teracti onanalysis .T he flu id is suppo sed to be incom p ressib le ideal flow w ith the govern ing equati on as L ap lace equati on w h ich is discretized by B E M .Fo r si m p licity ,a beam as a structu re is show n in F ig .1.To circum ven t the freedom con strain t in FE M 2B

16、E M coup ling p rocedu re ,discon tinuou s boundary elem en ts are em p loyed to discretize boundary in tegral equati on fo r flu id dom ain .A p rob lem of a can tilever beam w ith the flu id on one side is analyzed .It is dem on strated ,by com paring w ith the analytical so lu ti on p resen ted in th is paper ,that the B