利用能量守恒和徑向基函數(shù)插值的流固耦合界面數(shù)據(jù)傳遞方法.pdf

1、第 43 卷 第 9 期西 安交 通 大 學(xué) 學(xué) 報(bào)Vol. 4392009 年 9 月JOU RNAL OFXI AN JIA OTONG U NIVERSITYSep. 2009利用能量守恒和徑向基函數(shù)插值的流固耦合界面數(shù)據(jù)傳遞方法蘇波1, 2 , 錢若軍1 , 袁行飛3 , 于猛1( 1. 同濟(jì)大學(xué)土木工程學(xué)院, 200092, 上海; 2. 江蘇大學(xué)理學(xué)院, 212013, 江蘇鎮(zhèn)江;3. 浙江大學(xué)土木工程學(xué)院, 310027, 杭州)摘要: 根據(jù)界面能量守恒原理, 將徑向基函數(shù)( RBF) 引入流固耦合分析領(lǐng)域, 提出了基于 RBF 插 值的流固耦合界面數(shù)據(jù)傳遞方法( RBF/ FS

2、I) . 推導(dǎo)了界面位移傳遞矩陣, 并根據(jù) RBF/ FSI 算法編 制了相應(yīng)的界面信息傳遞的計(jì)算程序. 以三維耦合界面的位移信息傳遞為例, 將數(shù)值計(jì)算結(jié)果和解 析解進(jìn)行了對(duì)比分析, 結(jié)果表明, RBF/ FSI 算法計(jì)算效率高, 計(jì)算結(jié)果準(zhǔn)確, 適合處理復(fù)雜耦合界 面的流固耦合信息傳遞. 此外, RBF/ FSI 算法允許 CSD( 計(jì)算結(jié)構(gòu)力學(xué)) 、CFD( 計(jì)算流體力學(xué)) 采用 任意網(wǎng)格形式, 因此在 CSD 和 CFD 計(jì)算程序之間可以開發(fā)獨(dú)立的界面信息傳遞接口程序, 很容易 實(shí)現(xiàn)弱耦合分析.關(guān)鍵詞: 流固耦合; 數(shù)據(jù)傳遞; 徑向基函數(shù); 能量守恒中圖分類號(hào): TU311 3; O24

3、2 2 文獻(xiàn)標(biāo)志碼: A 文章編號(hào): 0253 987X( 2009) 09 0114 06Data Exchange Method for Fluid-Structure Interaction Based on Energy Conservation and Interpolation Algorithm AdoptingRadial Basis FunctionSU Bo1, 2 , Q IAN Ruojun1 , YUAN Xingfei3 , YU Meng1( 1. College of Civil Engineering, T on gji U niver sity,S hang

4、 hai 200092,C hina; 2. Facu lty of S cience, Jiangsu U nivers ity,Zhenjiang, Jian gsu 212013, Chin a; 3. C ollege ofCivil E ngineering,Zh ejiang University, H angzhou 310027, C hina)Abstract: Fo llow ing the principle o f energ y conservation, a method adopting radial basis function ( RBF) is pro po

5、sed to per for m fluid- structure interaction ( FSI) analysis. T he displacem ent trans-fer m atr ix H is deduced and a corresponding prog ram by RBF/ FSI alg orithm is developed. Taking displacement data exchange w ith 3D interface as an example, the numerical solution by RBF/ FSI method is compar

6、ed w ith the analytical one. T he result show s that RBF/ FSI method is efficientand accurate, and is suitable fo r FSI data exchange oncomplicated interface. M oreover, themethod allow s CSD( Co mputational Structural Dy nam ics)and CFD( Com putational Fluid Dynam-ics) to choose ar bitrary g rid fo

7、rm to easily realize loose coupling method by developing an inde-pendent interface prog ram betw een CSD and CFD solv ers.Keywords: fluid- structur e interaction; data ex chang e; radial basis function; energy conservatio n通常情況下 CFD( 計(jì)算流體力學(xué)) 計(jì)算所要求傳遞問題, 因此尋求一種合理有效的數(shù)據(jù)交換方法的網(wǎng)格密度要比 CSD( 計(jì)算結(jié)構(gòu)力學(xué)) 的大得多, 由是實(shí)

8、現(xiàn)流固耦合分析的關(guān)鍵. 近年來, 數(shù)學(xué)界對(duì)以徑此產(chǎn)生流固耦合界面上 2 套非匹配網(wǎng)格之間的數(shù)據(jù)向基函數(shù)( RBF) 為基礎(chǔ)的近似方法進(jìn)行了大量的研收稿日期: 2009 02 01. 作者簡介: 蘇波( 1977- ) , 男, 博士生; 錢若軍( 聯(lián)系人) , 男, 教授, 博士生導(dǎo)師. 基金項(xiàng)目: 國 家自然科學(xué)基金資助項(xiàng)目( 50708092) ; 浙江省自然科學(xué)基金資助項(xiàng)目( R1080363) .第 9 期蘇波, 等: 利用能量守恒和徑向基函數(shù)插值的流固耦合界面數(shù)據(jù)傳遞方法115究, 其中 Buhm ann 1 和吳宗敏 2 3 對(duì) RBF 的基礎(chǔ)理 論和應(yīng)用研究起了重要的推動(dòng)作用.

9、將 RBF 應(yīng)用于 工程實(shí)踐計(jì)算具有重要意義, 目前主要有 2 個(gè)應(yīng)用 領(lǐng)域: 無網(wǎng)格法求解偏微分方程和大型散亂數(shù)據(jù)的 處理. 前者可以用來求解固體或流體控制方程, 后者 可應(yīng)用于地形學(xué)、醫(yī)學(xué)、計(jì)算機(jī)圖形學(xué)等領(lǐng)域. Frank 4 曾對(duì)散亂數(shù)據(jù)的各種插值方法做了實(shí)例比 較, 得到的結(jié)論是, 徑向基函數(shù)插值的結(jié)果最能使人 滿意.本文將 RBF 引入流固耦合分析領(lǐng)域, 根據(jù)能量 守恒原理, 闡述如何采用 RBF 插值進(jìn)行界面信息傳 遞, 并利用本文提出的 RBF/ FSI 算法編制了相應(yīng)的 計(jì)算程序, 對(duì)三維耦合界面上的位移信息傳遞進(jìn)行 了計(jì)算分析.1 界面信息傳遞能量守恒原理耦合界面的能量守恒原

10、理是指在耦合作用過程 中, 耦合界面上流體荷載( 外力) 、固體力( 內(nèi)力) 在界面位移上所做的虛功相等 5 7 , 即( 1)W = usf s = uf f fTT式中: us 、uf 分別為耦合界面上固體、流體的虛位 移; f s 、f f 分別為耦合界面上固體、流體的表面力. 耦合界面上流體、固體虛位移之間的關(guān)系可表示為u = H u( 2)fs式中: H 為界面位移傳遞矩陣, 可由不同計(jì)算方法得到, 本文采用基于 RBF 插值的計(jì)算方法.將式( 2) 代入式( 1) , 可得壓力傳遞計(jì)算公式f s= HT f f( 3)根據(jù)式( 1) 式( 3) , 可知 H 的求解至關(guān)重要.2 R

11、BF 插值函數(shù)及其定解條件2 1 RBF 插值函數(shù)在 d 維歐幾里德空間給定一組位置不同的點(diǎn)X= x1 xN Rd , 稱之為 中心點(diǎn) . 進(jìn)一步, 如 果知道在這些中心點(diǎn)上的標(biāo)量值 g1 , g2 , , gN , 可 以由此確定一個(gè)連續(xù)函數(shù), 使函數(shù)通過這些中心點(diǎn), 當(dāng)連續(xù)函數(shù)采用式( 4) 表示的形式時(shí), 稱之為 RBF 插值函數(shù)Ns( x) =i (x - xi ) + p ( x)i =1( 4)式中: s( x) 為 x 點(diǎn)處的未知函數(shù)值; x 為未知點(diǎn)坐標(biāo)( x , y , z ) ; xi 為第 i 已知數(shù)據(jù)點(diǎn)坐標(biāo)( x i , y i , z i ) ; 為選用的徑向基函數(shù);

12、 i 為相應(yīng)于第 i 個(gè)數(shù)據(jù)點(diǎn)的待求系數(shù); x- xi為歐幾里得距離, 對(duì)于三維空間可用 r 直接表示為r =x - xi = ( x - x i ) 2 +( y - y i ) 2 + ( z - z i ) 2 1/ 2( 5)p ( x) 為低階 d 維多項(xiàng)式, 對(duì)于小于等于 2 階的有條 件正定基函數(shù), p ( x) 可以采用線性多項(xiàng)式 1p ( x) = 0 + x x + y y + z z( 6)2 2 RBF 插值函數(shù)的定解條件式( 4) 中的系數(shù) i和 p ( x) 可以根據(jù)式( 7) 和式( 8) 的定解條件求得s( xi ) = g ii = 1, 2, , N( 7)

13、Ni q ( x) = 0( 8)i= 1式中: 多項(xiàng)式 q( x) 為滿足 deg ( q ( x) )deg ( p ( x) )的所有多項(xiàng)式.3 界面位移傳遞矩陣的計(jì)算在流固耦合界面上, 固體域上存在 N s 個(gè)固體節(jié)點(diǎn)xsi = ( x s i , y si , z si )R3i = 1, 2, , N s( 9)流體域上存在 N f 個(gè)流體節(jié)點(diǎn)xfi = ( x fi , y fi , z fi )R3i = 1, 2, , N f( 10)ds 、df 分別表示固體、流體邊界節(jié)點(diǎn)的位移矢量, 位 移傳遞的目標(biāo)旨在用 ds ( xs i ) 表示 df ( xfi ) .將式( 9

14、) 代入式( 7) , 并考慮式( 8) , 經(jīng)過簡單的 計(jì)算處理, N s 個(gè)固體節(jié)點(diǎn)位移 ds 在 x 方向的位移 分量ds x 可用下式表示Dsx=Css x( 11)式中x0x0xx0d s1D=;d=;=yzsx0sxxx0dsN sxs1dsxxsN s( 12)http:w w w. jdxb. cn116西 安 交 通 大 學(xué) 學(xué) 報(bào)第 43 卷Css =00001110000x s1x s2x sN s0000y s1y s2y s N s0000z s1z s2z sN s1x s1y s1z s1s11s12s1Ns1x sN sy sN sz sN ssN s 1sN

15、s 2sN s N s( 13)s12= (xs1 - xs2)( 14)基函數(shù)為節(jié)點(diǎn) 1 和 2 間的距離函數(shù),余同.根據(jù)式( 11) , 即可求得系數(shù) x( 15)x = CssDsx- 1為確定流體表面節(jié)點(diǎn)的位移, 首先建立如下點(diǎn)位矩陣Afs =1 x f1y f1z f1f1, s1f1, s2f1, sN s1 x fN fy fN fz fNff N f , s1f N f , s2f N f , sN s( 16)流體表面節(jié)點(diǎn)的位移 dfx 可以表示為( 17)dfx = Afsx = Afs Css Dsx- 1式( 17) 可進(jìn)一步處理, 因?yàn)榫仃?Css 有零塊, 將之進(jìn)行

16、分塊處理Css =0P( 18)PTMPLY矩陣 x 可以分解成 2 部分,x =x, 其中 x 、RBFPLYRBF分別表示為xxxx0s1xxPLY=x;RBF=s2( 19)xxxyxxzs N s上標(biāo) PLY、RBF 分別表示屬于多項(xiàng)式 p ( x) 和 RBF 插值函數(shù)中的系數(shù).將式( 17) 寫成分塊矩陣的形式PLY0x- 1( 20)RBF= Cssdsxx經(jīng)過矩陣運(yùn)算, 可求得ds x( 21)x= Mp PMPLY- 1x= ( M- MPMp PM) dsx( 22)RBF- 1- 1T- 1式中: Mp = ( PM- 1 PT ) - 1 .由式( 17) 可得PLYd

17、fx = Afsx = Afsx=RBFAfsMp PM- 1xdsx( 23)- 1- 1- 1 TM- M P Mp PM式( 23) 建立了流體節(jié)點(diǎn)位移和固體節(jié)點(diǎn)位移間的關(guān)系, 因此傳遞矩陣 H 可直接寫為Mp PM- 1M- 1 - M- 1 PT Mp PM- 1類似地, 流體表面 y 、z 方向的節(jié)點(diǎn)位移 dfy 、dfz也可由固體表面 y 、z方向的節(jié)點(diǎn)位移 d sy 、dsz 左乘以傳遞矩陣H 求得, 寫成矩陣的形式為dfxHdsxdfy=Hdsy( 25)dfzHdsz得到傳遞矩陣 H 后, 根據(jù)能量守恒原理即可編制相 應(yīng)的計(jì)算程序, 對(duì)流固耦合界面的位移、速度、壓力 等物理

18、量進(jìn)行傳遞.4 RBF 插值函數(shù)在流固耦合界面 信息傳遞中的應(yīng)用條件RBF 插值函數(shù)應(yīng)用于流固耦合信息傳遞是一 種整體插值的方法, 對(duì)數(shù)據(jù)點(diǎn)來源并無特殊要求, 可 以是規(guī)則網(wǎng)格數(shù)據(jù)點(diǎn)( 如有限元網(wǎng)格節(jié)點(diǎn)) , 也可以 是無規(guī)則的散亂點(diǎn). 因此, 在使用網(wǎng)格化方法時(shí), CFD 和 CSD 可以保持相對(duì)的獨(dú)立性, 這有利于獨(dú) 立耦合界面信息程序的開發(fā).如式( 24) 所示, 采用 RBF 插值函數(shù)進(jìn)行數(shù)據(jù)傳 遞時(shí)涉及傳遞矩陣 H 的計(jì)算, 其中包括矩陣 M 的 求逆運(yùn)算, 只要 M 是非奇異矩陣, 就可采用 RBF 插 值函數(shù)進(jìn)行信息傳遞. 根據(jù)文獻(xiàn) 8 中關(guān)于 RBF 插 值函數(shù)存在性的相關(guān)定理

19、, 當(dāng) RBF 函數(shù) 為正定函 數(shù)或條件正定函數(shù)時(shí), 對(duì)于兩兩不同的數(shù)據(jù)點(diǎn), M 為非奇異矩陣.5 算例分析: 三維殼體耦合界面的 位移傳遞如圖 1 所示, 耦合界面為一半圓柱殼體曲面, 計(jì)算域取 0 x 1, 0 5 y 1, 0 z 2.假定界面在 x 方向的位移場(chǎng)可分別由以下的解析式表示 9dx ( x , y ) = 3 1+ co s( ( x - 0 5) ) 1 + cos( ( y - 1) ) 1+ cos( ( z - 1) ) + 0 5( 26)http:w w w. jdxb. cn第 9 期蘇波, 等: 利用能量守恒和徑向基函數(shù)插值的流固耦合界面數(shù)據(jù)傳遞方法117d

20、y ( x , y ) = 2 1 + co s( ( 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)圖 1半圓柱殼體曲面幾何CSD、CFD 界面均為三角形單元, 考慮 4 種不同密度的網(wǎng)格形式, 見圖 2.( a) 網(wǎng)格 1: 節(jié)點(diǎn)數(shù) 231, 網(wǎng)格數(shù) 400( b) 網(wǎng)格 2: 節(jié)點(diǎn)數(shù) 861, 網(wǎng)格數(shù) 1 600( c) 網(wǎng)格 3: 節(jié)點(diǎn)數(shù) 1 271, 網(wǎng)格數(shù) 2 4

21、00( d) 網(wǎng)格 4: 節(jié)點(diǎn)數(shù) 3 321, 網(wǎng)格數(shù) 6 400圖 2 4 種不同密度的界面網(wǎng)格根據(jù)式( 26) 式( 28) , 可求得流體界面網(wǎng)格點(diǎn) 在 x 、y、z 方向的位移解析解 d fx, a 、dfy, a 、dfz , a; 根據(jù) RBF/ FSI 算法, 可求得流體界面網(wǎng)格點(diǎn)在 x 、y、z 方 向的數(shù)值解d fx, n 、d fy , n 、dfz , n . 流體界面每個(gè)網(wǎng)格點(diǎn)解析解、數(shù)值解的相對(duì)誤差為R x = | dfx, a- dfx , n| ; Ry =| dfy , a- dfy, n|dfx, adfy, aRz =| dfz , a - dfz, n |

22、dfz, a考慮一般情況下 CFD 的網(wǎng)格密度大于 CSD 的 網(wǎng)格密度, 做 3 組界面網(wǎng)格計(jì)算組合, 見表 1, 表中 同時(shí)列出了每個(gè)組合在 x 、y 、z 方向的最大計(jì)算相 對(duì)誤差. 在 x 、y、z 方向全部網(wǎng)格點(diǎn)對(duì)應(yīng)的相對(duì)誤差 ( 取對(duì)數(shù)) 分別繪于圖 3 圖 5.表 1位移傳遞計(jì)算分組及最大相對(duì)誤差分界面網(wǎng)格形式最大相對(duì)誤差 max( lgR)組CSDCFDx 方向y 方向z 方向1網(wǎng)格 1網(wǎng)格 2- 1603 33 - 1706 81 - 1492 032網(wǎng)格 2網(wǎng)格 3- 2245 12 - 2355 38 - 2088 993 網(wǎng)格 3 網(wǎng)格 4 - 2 107 27 - 2 220 12 - 2 203 34先分析 x 方向位移的計(jì)算結(jié)果: 從表 1 和圖 3 可以看出, 對(duì)于計(jì)算分組 1、2, 隨著計(jì)算網(wǎng)格密度的 增大, 網(wǎng)格點(diǎn)的計(jì)算相對(duì)誤差和最大相對(duì)誤差在不 斷減小的同時(shí), 網(wǎng)格點(diǎn)相對(duì)誤差的分布范圍變窄, 這 說明計(jì)算精度隨網(wǎng)格密度的增大而提高. 但是, 對(duì)比 計(jì)