流固耦合界面信息傳遞方法研究

流固耦合界面信息傳遞方法研究

隨著計(jì)算力學(xué)的發(fā)展和計(jì)算機(jī)計(jì)算計(jì)算的提高，多物理場的集成分析開始廣泛應(yīng)用于不同的領(lǐng)域，如航空、醫(yī)學(xué)、磁學(xué)等。在結(jié)構(gòu)工程領(lǐng)域，結(jié)構(gòu)風(fēng)梯度的破壞日益嚴(yán)重。為了揭示風(fēng)和結(jié)構(gòu)的耦合效應(yīng)，有必要根據(jù)流固耦合理論進(jìn)行分析，并從物理本質(zhì)上揭示風(fēng)和結(jié)構(gòu)的耦合效應(yīng)規(guī)律。此外，該項(xiàng)目結(jié)構(gòu)也日益與其他學(xué)科相結(jié)合。例如，關(guān)于太陽能熱態(tài)屋頂結(jié)構(gòu)的研究包括分析壓力場、位移場、磁體場和電爐等物理和化學(xué)性質(zhì)的耦合效應(yīng)。為能對復(fù)雜的耦合場進(jìn)行分析,需要在不同領(lǐng)域的計(jì)算程序間搭建信息傳遞平臺,因此,合理有效的界面信息傳遞方法是實(shí)現(xiàn)耦合分析的關(guān)鍵環(huán)節(jié)之一.本文針對流固耦合分析中的CFD、CSD間的界面信息傳遞進(jìn)行系統(tǒng)的闡述和研究,旨在為實(shí)現(xiàn)流固耦合分析提供參考,文中多數(shù)方法可拓展至其他領(lǐng)域.1cfd和csd流的基本原理流固耦合的分析過程一般可以描述為:在流場(如風(fēng)、水流等)作用下,結(jié)構(gòu)特別是柔性結(jié)構(gòu)(如大跨、高聳、膜結(jié)構(gòu)等)會產(chǎn)生較大的變形和振動,因此會對周圍流場會產(chǎn)生較大影響,而流場的改變會進(jìn)一步改變作用在結(jié)構(gòu)表面上的壓力大小,從而形成流體與結(jié)構(gòu)的相互耦合作用,圖1.要合理有效地在CFD和CSD流之間進(jìn)行數(shù)據(jù)交換,必須遵循一定的原理和法則.結(jié)合文獻(xiàn),概括為以下4個基本原理.1.1流體耦合邊界位移分量的uf、vf、wf運(yùn)動學(xué)連續(xù)性條件即在流固界面上流體、固體對應(yīng)點(diǎn)位移、速度的一致性.位移連續(xù)條件可表示為:df=ds(1)寫成分量的形式:式中,df為流體在耦合邊界的位移,dfx、dfy、dfz分別為df在x、y、z方向的位移分量;ds為結(jié)構(gòu)在耦合邊界的位移,dsx、dsy、dsz分別為ds在x、y、z方向的位移分量.速度連續(xù)條件可表示為:uf=us(3)寫成分量的形式:式中,uf為流體在耦合邊界的速度,uf、vf、wf為uf在x、y、z方向的速度分量;us為結(jié)構(gòu)在耦合邊界的速度,us、vs、wf為us在x、y、z方向的速度分量.1.2耦合邊界上的應(yīng)力動力連續(xù)性條件即接觸界面要滿足力的守恒.根據(jù)界面上任一點(diǎn)的力的平衡:σs·ns=σf·nf(5)式中,σs、σf分別為結(jié)構(gòu)、流體在耦合邊界上的柯西應(yīng)力張量;ns、nf分別為結(jié)構(gòu)、流體在耦合邊界上的外法線方向矢量.對固體域而言,耦合邊界作為結(jié)構(gòu)的紐曼邊界條件,式(5)可寫成分量的形式:式中,fx、fy、fz為固體在耦合邊界上任一點(diǎn)沿x、y、z方向的應(yīng)力分量;Px、Py、Pz為流體在耦合邊界上任一點(diǎn)沿x、y、z方向的壓強(qiáng)分量,其值由CFD計(jì)算給出.1.3耦合界面上固體、流體的變化耦合界面的能量守恒原理是指在耦合作用過程中,耦合界面上流體荷載(外力)、固體力(內(nèi)力)在界面位移上所做的虛功相等:δW=δuΤs·fs=δuΤf·ff(7)式中,δus、δuf分別為耦合界面上固體、流體虛位移;fs、ff分別為耦合界面上固體、流體表面力.耦合界面上流體、固體的虛位移之間關(guān)系可以用下式表示:δuf=Hδus(8)H為由不同耦合方法得到的傳遞矩陣.將式(8)代入式(7)可得:fs=HTff(9)根據(jù)式(7)～式(9),可知傳遞矩陣H的求解致關(guān)重要.文獻(xiàn)證明當(dāng)轉(zhuǎn)換矩陣H中行元素之和等于1時,可從能量守恒出發(fā),導(dǎo)出界面力的守恒.1.4能量守恒原理效率的高低可以根據(jù)求解精度和求解時間的比值進(jìn)行衡量.一般情況下,同時滿足上述四個條件的理想的數(shù)據(jù)傳輸比較困難,而能量守恒原理作為自然界的一般法則,近些年來得到廣泛重視和應(yīng)用.2網(wǎng)格密度函數(shù)問題如圖2所示,流體域和固體域在離散后,根據(jù)網(wǎng)格節(jié)點(diǎn)間的對應(yīng)關(guān)系可分為一致網(wǎng)格和非匹配網(wǎng)格.通常情況下CFD計(jì)算要求的網(wǎng)格密度比CSD要密得多,由此產(chǎn)生耦合界面上兩套非匹配網(wǎng)格之間的數(shù)據(jù)傳遞問題,在數(shù)學(xué)上這是一個雙向插值問題,對該問題的研究也成為流固耦合的熱點(diǎn)和難點(diǎn)之一.3局部插值法耦合界面信息傳遞方法從總體上可以分為局部插值法和整體插值法.局部插值法指在進(jìn)行信息傳遞時,流體或固體界面網(wǎng)格點(diǎn)的未知物理量可根據(jù)固體或流體界面上的部分單元或網(wǎng)格點(diǎn)的已知物理量求得,而整體插值法需要根據(jù)所有單元或網(wǎng)格點(diǎn)上的已知物理量求得.以下介紹局部插值法中的不同方法.3.1固體主單元節(jié)點(diǎn)數(shù)計(jì)算映射點(diǎn)插值法是最直接、最簡單的傳遞方法.其思想是:流固耦合邊界上流體節(jié)點(diǎn)的位移和其對應(yīng)的固體主單元中的映射點(diǎn)位移相等,而映射點(diǎn)位移由固體主單元中根據(jù)形函數(shù)插值得到,對力的傳遞則是一個相反的過程.關(guān)于主單元及映射點(diǎn)的概念和算法在3.4節(jié)給出.如圖2b所示,流體邊界單元中的節(jié)點(diǎn)A,其在結(jié)構(gòu)中的主單元為abc,映射點(diǎn)為A′,根據(jù)映射點(diǎn)插值法:dfA=dfA′=ns∑i=1Νsidsi=Νsa(xA′)dsa+Νsb(xA′)dsa+Νsc(xA′)dsc(10)式中,ns為固體主單元節(jié)點(diǎn)數(shù);dfA、dfA′分別為流體節(jié)點(diǎn)和映射點(diǎn)的位移矢量.寫成分量的形式:dfAx=dfA′x=ns∑i=1Νsixdsix=Νsa(xA′)dsax+Νsb(xA′)dsax+Νsc(xA′)dscx(11a)dfAy=dsA′y=ns∑i=1Νsiydsiy=Νsa(xA′)dsay+Νsb(xA′)dsay+Νsc(xA′)dscy(11b)dfAz=dsA′z=ns∑i=1Νsizdsiz=Νsa(xA′)dsaz+Νsb(xA′)dsaz+Νsc(xA′)dscz(11c)采用上述方法,流體界面所有節(jié)點(diǎn)位移可由固體界面節(jié)點(diǎn)位移表示出來,寫成矩陣的形式:式中,Hnf×ns為位移轉(zhuǎn)換矩陣,根據(jù)映射點(diǎn)插值法組裝得到,ns、nf分別為固體、流體耦合邊界上的節(jié)點(diǎn)數(shù).映射點(diǎn)插值法容易理解和應(yīng)用,但缺點(diǎn)是沒有考慮界面整體位移協(xié)調(diào)和力的平衡,使其應(yīng)用范圍受到很大的局限.3.2界面位移nisStein等根據(jù)任一點(diǎn)位移相等的概念,采用加權(quán)余量法進(jìn)行信息傳遞.界面上位移守恒,其連續(xù)形式:ds(x)=df(x)(13)ds、df可近似表示為:ds=ns∑i=1Νis(x)dsi,df=nf∑j=1Νif(x)dfi(14)式中,ds、df分別為界面上結(jié)構(gòu)、流體網(wǎng)格點(diǎn)上的位移矢量,Nis、Nif為結(jié)構(gòu)、流體的基函數(shù).根據(jù)迦遼金方法,對式(13)同乘以權(quán)函數(shù)?(x),得到:∫Γ?(x)ns∑i=1Νis(x)usidx=∫Γ?(x)nf∑j=1Νif(x)dfidx(15)對于?(x),可有兩種選擇,即可以為流體或結(jié)構(gòu)的基函數(shù).ns∑i=1[∫ΓΝkfΝisdx]dsi=nf∑j=1[∫ΓΝkfΝifdx]dfifork=1,...nf(16)可寫成矩陣的形式:Afsus=Affdf(17)式中,Afs為nf×ns的矩陣,Aff為nf×nf的矩陣.由式(17),流體界面位移可由界面結(jié)構(gòu)上已知位移求得:df=A-1ffAfsds(18)因此位移轉(zhuǎn)換矩陣寫為:H=A-1ffAfs(19)采用加權(quán)余量法時,會碰到形如∫ΓsΝfΝsdS的積分問題,即積分式里有關(guān)于流體和固體單元插值函數(shù),但積分域?yàn)楣腆w界面單元或流體界面單元.至今,這一混合積分未得到有效處理.3.3變密度插值Goura等人提出了一種常體積轉(zhuǎn)換方法(CVT),該方法是一種與結(jié)構(gòu)模態(tài)無關(guān)的局部插值方法.基本思想為:氣動點(diǎn)的正交投影在結(jié)構(gòu)三角形內(nèi)保持不變,與固體表面單元的線性彈性特性一致.通過確保由結(jié)構(gòu)點(diǎn)和氣動點(diǎn)組成的四面體的體積守恒,氣動點(diǎn)超出平板的分量就可以計(jì)算出來.徐敏、陳士櫓在CVT基礎(chǔ)上提出了一種改進(jìn)的常體積轉(zhuǎn)換法,即引入面積限制值來保證網(wǎng)格插值的質(zhì)量.改進(jìn)的CVT插值方法能避免原CVT方法可能出現(xiàn)的異常情況,且大大提高了插值精度.CVT法也需計(jì)算節(jié)點(diǎn)間的映射關(guān)系,且關(guān)于常體積的假定增加了數(shù)據(jù)傳遞的復(fù)雜性.3.4尋找主單元和映射點(diǎn)的算法應(yīng)用局部插值法不可避免地涉及到映射點(diǎn)的搜尋問題,這需要單獨(dú)的算法和計(jì)算程序.3.4.1單元內(nèi)部映射流體或固體界面上的某一節(jié)點(diǎn)向固體或流體界面上單元所在平面作正投影,當(dāng)映射點(diǎn)在單元內(nèi)部時,稱此單元為節(jié)點(diǎn)的主單元,其垂足為映射點(diǎn),圖3.3.4.2映射點(diǎn)位置的確定假設(shè)單元ijk為節(jié)點(diǎn)p(xp)的一個備選主單元,見圖4:節(jié)點(diǎn)p(xp)在單元ijk中映射點(diǎn)p′的坐標(biāo)xp′:xp′=2∑i=1[(xp-x0)?gi]gi(20)垂足:xn=2∑i=1[(xp-x0)?g3]g3(21)3.4.3單元特征特性考慮映射點(diǎn)p′(xp′),確定p′點(diǎn)是否在給定單元中的直接方法是根據(jù)p′點(diǎn)坐標(biāo)值計(jì)算的其形函數(shù)值.根據(jù)有限元概念:xp′=∑iΝixi(22)形函數(shù)之和:∑iΝi=1(23)由式(22)、式(23),對于三維四面體單元:式(24)可簡寫為:xp=xN(25)因此,可以求得N:N=x-1xp(26)對于三維四面單元:由此,判斷點(diǎn)xp′是否在單元中的條件為:min(Ni,1-Ni)≥0?i(28)在某些情況下,按式(28)得到的主單元可能不止一個,在此情況下,需要附加的判定標(biāo)準(zhǔn)為:dn=|xn|≤δ(29)式中,dn為法向距離;xn為依式(21)計(jì)算所得的垂足;δ為指定容差.3.4.4求解算法2:完全搜尋法根據(jù)上述原理,可采用完全搜尋法(bruteforcesearch)進(jìn)行主單元及映射點(diǎn)的搜尋.完全搜尋法首先假定映射點(diǎn)p′及其坐標(biāo)xp′,然后對邊界上所有的單元作一次循環(huán)計(jì)算,計(jì)算與xp′相關(guān)的形函數(shù).其計(jì)算步驟可描述為:(1)根據(jù)式(20)、(21)計(jì)算映射點(diǎn)和垂足的坐標(biāo).(2)根據(jù)式(27)計(jì)算Ni,并按式(28)進(jìn)行判斷,如符合條件的主單元不止一個,進(jìn)行步驟(3).(3)根據(jù)式(29)作進(jìn)一步判定.回到步驟1,直到完成對所有邊界單元節(jié)點(diǎn)掃描和計(jì)算.完全搜尋法算法,需要對所有單元進(jìn)行循環(huán),一般情況下要一半單元需要參與循環(huán).為提高計(jì)算效率,文獻(xiàn)還給出一些改進(jìn)方法:Octreesearch,NeighborttoNeighborsearch,Advancingsearch方法等.這些方法都是基于不同的假定,采用不同數(shù)據(jù)結(jié)構(gòu)類型,因此計(jì)算量比完全搜尋法算法有不同程度降低.如Octreesearch方法對計(jì)算域就有不同的數(shù)據(jù)結(jié)構(gòu)類型,圖5.需要指出的是Octreesearch方法通過分解域搜尋法可以提高計(jì)算效率,但由于單元組織相對復(fù)雜,可能在計(jì)算過程中失效.其他改進(jìn)的方法亦可能出現(xiàn)類似的問題,此時需重新啟用完全搜尋法.4一般插值方法4.1有限表面插值全域插值法最早來源于航空領(lǐng)域中對平面機(jī)翼變形的計(jì)算.Harder(1971)提出了無限平板樣條法(IPS),該方法基于一無限板微分平衡方程的疊加求解.Yu對IPS作了進(jìn)一步研究,并編制了相應(yīng)的計(jì)算程序.IPS是程序中非常普遍采用的一種插值方法,如MSC/NASTRAN.Appa(1989)認(rèn)為采用IPS法從結(jié)構(gòu)內(nèi)部網(wǎng)格點(diǎn)外插得到平板邊網(wǎng)格點(diǎn)的值不可靠,并提出有限表面插值(FPS).FPS引入了一系列的約束條件,使得變形后的板經(jīng)過給定數(shù)據(jù)點(diǎn).由此,建立了結(jié)構(gòu)和氣動網(wǎng)格點(diǎn)位移映射,得到傳遞矩陣.根據(jù)傳遞矩陣可以對流固界面間變形、壓強(qiáng)、溫度等插值計(jì)算.此外,還有InverseIsoparametricMapping(IIM),方法以有限元中的等參元為基礎(chǔ).但此方法只適用二維問題,且無法進(jìn)行外推.徐敏認(rèn)為上述樣條插值僅適合于薄板處于最小彎曲能(平衡位置)所確定的位置并且應(yīng)在滿足流體表面和結(jié)構(gòu)表面一致的條件下才能得到理想的結(jié)果.在IPS、FPS之外,出現(xiàn)了旨對不規(guī)則曲面的插值計(jì)算的樣條函數(shù):薄板樣條(TPS),Multiquadric-Biharmonic(MQ)等.Smith對上述方法作了技術(shù)上的全面闡述,并指出了每種方法適用范圍.4.2距離平均加權(quán)和距離平均加權(quán)Shepard方法首先由氣象學(xué)家與地質(zhì)工作者提出,后來由于Shepard的工作稱為Shepard方法,其基本思想是將插值函數(shù)定義為各數(shù)據(jù)點(diǎn)函數(shù)值的距離倒數(shù)加權(quán)和:s(x)=n∑j=0fjlj(x)=n∑j=0fj∏k≠jrk(x)∑nm=0∏k≠mrk(x)(30)式中,fj為x點(diǎn)處的函數(shù)值;lj(x)為權(quán)函數(shù);rk(x)=|x-xj|22.Shepard插值雖然是光滑的,但缺點(diǎn)是在數(shù)據(jù)點(diǎn)附近表現(xiàn)為平點(diǎn)(其各階導(dǎo)數(shù)為零),這樣限制了Shepard插值的精度.如果已經(jīng)知道了函數(shù)在數(shù)據(jù)點(diǎn)的高階導(dǎo)數(shù),可以通過構(gòu)造改進(jìn)的Shepard插值或帶導(dǎo)數(shù)條件的Shepard插值來解決此問題.4.3rbf插值函數(shù)近些年來,數(shù)學(xué)界對徑向基函數(shù)(radialbasisfunction,RBF)為基礎(chǔ)的近似方法做了大量的研究工作.其中,德國Justus-Liebig大學(xué)的Martin.D.Buhmann教授和國內(nèi)復(fù)旦大學(xué)數(shù)學(xué)系的吳宗敏教授對RBF函數(shù)的基礎(chǔ)理論和應(yīng)用研究起了重要推動作用.徑向基函數(shù)插值方法可用于大型散亂數(shù)據(jù)的處理而被廣泛應(yīng)用于地形學(xué)、醫(yī)學(xué)、計(jì)算機(jī)圖形學(xué)等領(lǐng)域.Frank曾做了大量的各種散亂數(shù)據(jù)插值方法的實(shí)例比較,得到的結(jié)論是徑向基函數(shù)插值的結(jié)果最能使人滿意.在d維歐幾里德空間里給定一組位置不同的點(diǎn)X={x1[WW(]...[WW)]xΝ}?Rd,稱之為“中心點(diǎn)”.進(jìn)一步,如果知道在這些中心點(diǎn)上的標(biāo)量值g1,g2,...,gN,可以由此確定一個連續(xù)函數(shù),使函數(shù)通過這些中心點(diǎn).當(dāng)連續(xù)函數(shù)采用式(1)表示的形式時,稱之為RBF插值函數(shù):s(x)=i=Ν∑i=1αi?(|x-xi|)(31)式中,s(x)為x點(diǎn)處的未知函數(shù)值;x為未知點(diǎn)坐標(biāo)x(x,y,z);xi為第i已知數(shù)據(jù)點(diǎn)坐標(biāo)xi(x,y,z);?為選用的徑向基函數(shù);αi為相應(yīng)于第i個數(shù)據(jù)點(diǎn)的待求系數(shù);函數(shù)|x-xi|為歐幾里得距離,對于三維空間可用r直接表示為:r=|x-xi|=√(x-xi)2+(y-yi)2+(z-zi)2(32)根據(jù)RBF插值算法,流體表面x、y、z方向的節(jié)點(diǎn)位移dfx、,dfy、dfz也由固體表面x、y、z方向的節(jié)點(diǎn)位移dsx、dsy、dsz左乘以傳遞矩陣H求得.得到傳遞矩陣H后,根據(jù)能量守恒原理即可編制相應(yīng)的計(jì)算程序,對流固耦合界面位移、速度、壓力等物理量的進(jìn)行傳遞.值得注意的是對自變量空間,可以找到一個截?cái)喽囗?xiàng)式,使得由其產(chǎn)生的徑向函數(shù)在給定的維數(shù)空間正定且滿足給定的連續(xù)性條件,即采用緊支徑向基函數(shù)(compactlysupportedradialbasisfunction,簡稱CSRBF).這類函數(shù)適用于處理大規(guī)模計(jì)算,因?yàn)橛纱藢?dǎo)出的線性方程組是一個稀疏矩陣.RBF插值方法允許CSD、CFD采用任意網(wǎng)格形式,因此在CSD和CFD計(jì)算程序之間可以開發(fā)獨(dú)立的界面信息傳遞的接口程序,很容易實(shí)現(xiàn)弱耦合分析,