基于rbf插值的流固耦合界面數(shù)據(jù)傳遞方法

基于rbf插值的流固耦合界面數(shù)據(jù)傳遞方法

通常，計算cfd（計算力學(xué)）所需的網(wǎng)格密度大于csd（計算結(jié)構(gòu)力學(xué)），因此在流固耦合界面中生成了兩種非組合網(wǎng)格之間的數(shù)據(jù)傳輸問題。因此，尋求合理有效的數(shù)據(jù)處理方法是流固耦合分析的關(guān)鍵。近年來，許多科學(xué)家對基于rbf的近似方法進(jìn)行了大量研究。其中，buhaman和吳宗敏在rbf的基礎(chǔ)理論和應(yīng)用研究方面發(fā)揮了重要作用。rbf的應(yīng)用于工程實踐的計算非常重要。目前，主要應(yīng)用如下：rbf搜索方法和大型離散數(shù)據(jù)處理。第一個可以用作固體或液體控制方程，第二個可以應(yīng)用于地形學(xué)、醫(yī)學(xué)和計算機圖形學(xué)領(lǐng)域。第三個可以使用于計算體積和流量控制方程。第二個可以應(yīng)用于地形學(xué)、醫(yī)學(xué)和計算機工程。f。比較了散列數(shù)據(jù)的不同插值方法。結(jié)論是，流固函數(shù)的插值結(jié)果最令人滿意。本文將RBF引入流固耦合分析領(lǐng)域,根據(jù)能量守恒原理,闡述如何采用RBF插值進(jìn)行界面信息傳遞,并利用本文提出的RBF/FSI算法編制了相應(yīng)的計算程序,對三維耦合界面上的位移信息傳遞進(jìn)行了計算分析.1耦合界面上固體、流體的虛位移耦合界面的能量守恒原理是指在耦合作用過程中,耦合界面上流體荷載(外力)、固體力(內(nèi)力)在界面位移上所做的虛功相等,即δW=δuΤs?fs=δuΤf?ff(1)δW=δuTs?fs=δuTf?ff(1)式中:δus、δuf分別為耦合界面上固體、流體的虛位移;fs、ff分別為耦合界面上固體、流體的表面力.耦合界面上流體、固體虛位移之間的關(guān)系可表示為δuf=Ηδus(2)δuf=Hδus(2)式中:H為界面位移傳遞矩陣,可由不同計算方法得到,本文采用基于RBF插值的計算方法.將式(2)代入式(1),可得壓力傳遞計算公式fs=ΗΤff(3)根據(jù)式(1)～式(3),可知H的求解至關(guān)重要.2rbf插值函數(shù)及其定解方法2.1rbf插值函數(shù)在d維歐幾里德空間給定一組位置不同的點X={x1…xN}?Rd,稱之為“中心點”.進(jìn)一步,如果知道在這些中心點上的標(biāo)量值g1,g2,…,gN,可以由此確定一個連續(xù)函數(shù),使函數(shù)通過這些中心點,當(dāng)連續(xù)函數(shù)采用式(4)表示的形式時,稱之為RBF插值函數(shù)s(x)=Ν∑i=1αi?(∥x-xi∥)+p(x)(4)式中:s(x)為x點處的未知函數(shù)值;x為未知點坐標(biāo)(x,y,z);xi為第i已知數(shù)據(jù)點坐標(biāo)(xi,yi,zi);?為選用的徑向基函數(shù);αi為相應(yīng)于第i個數(shù)據(jù)點的待求系數(shù);‖x-xi‖為歐幾里得距離,對于三維空間可用r直接表示為r=∥x-xi∥=[(x-xi)2+(y-yi)2+(z-zi)2]1/2(5)p(x)為低階d維多項式,對于小于等于2階的有條件正定基函數(shù),p(x)可以采用線性多項式p(x)=γ0+γxx+γyy+γzz(6)2.2型degx及廣義解析條件式(4)中的系數(shù)αi和p(x)可以根據(jù)式(7)和式(8)的定解條件求得s(xi)=gii=1?2???Ν(7)Ν∑i=1αiq(x)=0(8)式中:多項式q(x)為滿足deg(q(x))≤deg(p(x))的所有多項式.3流體表面節(jié)點位移的確定在流固耦合界面上,固體域上存在Ns個固體節(jié)點xsi=(xsi?ysi?zsi)∈R3i=1?2???Νs(9)流體域上存在Nf個流體節(jié)點xfi=(xfi?yfi?zfi)∈R3i=1?2???Νf(10)ds、df分別表示固體、流體邊界節(jié)點的位移矢量,位移傳遞的目標(biāo)旨在用ds(xsi)表示df(xfi).將式(9)代入式(7),并考慮式(8),經(jīng)過簡單的計算處理,Ns個固體節(jié)點位移ds在x方向的位移分量dsx可用下式表示Dsx=Cssαx(11)式中基函數(shù)為節(jié)點1和2間的距離函數(shù),余同.根據(jù)式(11),即可求得系數(shù)αxαx=C-1ssDsx(15)為確定流體表面節(jié)點的位移,首先建立如下點位矩陣Afs=[1xf1yf1zf1?f1,s1?f1,s2??f1,sΝs???????1xfΝfyfΝfzfΝf?fΝf?s1?fΝf,s2??fΝf,sΝs](16)流體表面節(jié)點的位移dfx可以表示為dfx=Afsαx=AfsC-1ssDsx(17)式(17)可進(jìn)一步處理,因為矩陣Css有零塊,將之進(jìn)行分塊處理Css=[0ΡΡΤΜ](18)矩陣αx可以分解成2部分,,其中αΡLYx、αRBFx分別表示為αΡLYx=[γx0γxxγxyγxz]；αRBFx=[αxs1αxs2?αxsΝs](19)上標(biāo)PLY、RBF分別表示屬于多項式p(x)和RBF插值函數(shù)中的系數(shù).將式(17)寫成分塊矩陣的形式[αΡLYxαRBFx]=C-1ss[0dsx](20)經(jīng)過矩陣運算,可求得αΡLYx=ΜpΡΜ-1dsx(21)αRBFx=(Μ-1-Μ-1ΡΤΜpΡΜ-1)dsx(22)式中:Mp=(PM-1PT)-1.由式(17)可得dfx=Afsαx=Afs[αΡLYxαRBFx]=Afs[ΜpΡΜ-1Μ-1-Μ-1ΡΤΜpΡΜ-1]dsx(23)式(23)建立了流體節(jié)點位移和固體節(jié)點位移間的關(guān)系,因此傳遞矩陣H可直接寫為Η=Afs[ΜpΡΜ-1Μ-1-Μ-1ΡΤΜpΡΜ-1](24)類似地,流體表面y、z方向的節(jié)點位移dfy、dfz也可由固體表面y、z方向的節(jié)點位移dsy、dsz左乘以傳遞矩陣H求得,寫成矩陣的形式為[dfxdfydfz]=[ΗΗΗ][dsxdsydsz](25)得到傳遞矩陣H后,根據(jù)能量守恒原理即可編制相應(yīng)的計算程序,對流固耦合界面的位移、速度、壓力等物理量進(jìn)行傳遞.4cfd和csd的模型構(gòu)建RBF插值函數(shù)應(yīng)用于流固耦合信息傳遞是一種整體插值的方法,對數(shù)據(jù)點來源并無特殊要求,可以是規(guī)則網(wǎng)格數(shù)據(jù)點(如有限元網(wǎng)格節(jié)點),也可以是無規(guī)則的散亂點.因此,在使用網(wǎng)格化方法時,CFD和CSD可以保持相對的獨立性,這有利于獨立耦合界面信息程序的開發(fā).如式(24)所示,采用RBF插值函數(shù)進(jìn)行數(shù)據(jù)傳遞時涉及傳遞矩陣H的計算,其中包括矩陣M的求逆運算,只要M是非奇異矩陣,就可采用RBF插值函數(shù)進(jìn)行信息傳遞.根據(jù)文獻(xiàn)中關(guān)于RBF插值函數(shù)存在性的相關(guān)定理,當(dāng)RBF函數(shù)?為正定函數(shù)或條件正定函數(shù)時,對于兩兩不同的數(shù)據(jù)點,M為非奇異矩陣.5x、y、z-dxa、dfza如圖1所示,耦合界面為一半圓柱殼體曲面,計算域取0≤x≤1,0.5≤y≤1,0≤z≤2.假定界面在x方向的位移場可分別由以下的解析式表示dx(x?y)=3[1+cos((x-0.5)π)][1+cos((y-1)π)][1+cos((z-1)π)]+0.5(26)dy(x?y)=2[1+cos((x-0.5)π)][1+cos((y-1)π)][1+cos((z-1)π)]+1(27)dz(x?y)=[1+cos((x-0.5)π)][1+cos((y-1)π)][1+cos((z-1)π)]+1(28)CSD、CFD界面均為三角形單元,考慮4種不同密度的網(wǎng)格形式,見圖2.根據(jù)式(26)～式(28),可求得流體界面網(wǎng)格點在x、y、z方向的位移解析解dfx,a、dfy,a、dfz,a;根據(jù)RBF/FSI算法,可求得流體界面網(wǎng)格點在x、y、z方向的數(shù)值解dfx,n、dfy,n、dfz,n.流體界面每個網(wǎng)格點解析解、數(shù)值解的相對誤差為Rx=|dfx?a-dfx?n|dfx?a；Ry=|dfy?a-dfy?n|dfy?aRz=|dfz?a-dfz?n|dfz?a考慮一般情況下CFD的網(wǎng)格密度大于CSD的網(wǎng)格密度,做3組界面網(wǎng)格計算組合,見表1,表中同時列出了每個組合在x、y、z方向的最大計算相對誤差.在x、y、z方向全部網(wǎng)格點對應(yīng)的相對誤差(取對數(shù))分別繪于圖3～圖5.先分析x方向位移的計算結(jié)果:從表1和圖3可以看出,對于計算分組1、2,隨著計算網(wǎng)格密度的增大,網(wǎng)格點的計算相對誤差和最大相對誤差在不斷減小的同時,網(wǎng)格點相對誤差的分布范圍變窄,這說明計算精度隨網(wǎng)格密度的增大而提高.但是,對比計算分組2、3可以看出,對于高密度網(wǎng)格,計算精度并沒有隨著網(wǎng)格密度的增加而提高,相反卻略有降低,但網(wǎng)格點相對誤差的分布范圍繼續(xù)變窄.從表1和圖4、圖5可以看出,y、z方向位移的計算結(jié)果與x方向位移的計算結(jié)果具有類似的性質(zhì),但對于z方向位移,計算精度并沒有隨著網(wǎng)格密度的增加而略有降低的現(xiàn)象.隨著網(wǎng)格密度的增加,網(wǎng)格面漸趨逼近殼體幾何面,因此計算精度隨之增加是合理的,但如表1所示,在計算分組2、3時,卻出現(xiàn)了“反?！爆F(xiàn)象,即在采用高密度網(wǎng)格時,計算精度沒有隨網(wǎng)格密度的增加而提高,反而在x、y方向出現(xiàn)了略有降低的現(xiàn)象.這是因為采用能量原理進(jìn)行數(shù)據(jù)傳遞時,需根據(jù)式(24)求得顯式的傳遞矩陣H,因此存在矩陣M的求逆運算,當(dāng)采用高密度網(wǎng)格時,大型矩陣的求逆運算導(dǎo)致了計算精度的不穩(wěn)定.如果不采用能量守恒原理,位移、壓力的傳遞就可以獨立進(jìn)行插值計算,位移的求解可不進(jìn)行矩陣求逆運算,直接求解式(11)的線性代數(shù)方程組.此時,對于計算分組2、3,采用高斯消除法求得的x方向的最大相對誤差分別為-2.15517和-2.32350,y方向的最大相對誤差分別為-2.24577和-2.41818,計算結(jié)果沒有出現(xiàn)上述的“反?！爆F(xiàn)象.根據(jù)各分組的整體計算結(jié)果來看,x、y、z方向絕大多數(shù)節(jié)點的位移相對誤差在10-2～10-6之間,計算精度較高.更多的計算結(jié)果表明,在進(jìn)行網(wǎng)格剖分時,固體、流體的界面網(wǎng)格邊長宜分別小于耦合界面最大邊長的1/20和1/40.6界面信息傳遞根據(jù)能量守恒原理,位移傳遞矩陣的計算在流固耦