直接邊界元法解勢流速度場問題

李井煜，盧曉平，趙鵬偉

直接邊界元法解勢流速度場問題

李井煜，盧曉平，趙鵬偉

海軍工程大學(xué)艦船工程系，湖北武漢430033

在船舶水動(dòng)力學(xué)中，大多采用以Hess-Smith方法為基礎(chǔ)的間接邊界元法求解勢流繞流問題，但Hess-Smith方法本質(zhì)上是基于物理直觀提出，在理論和數(shù)值計(jì)算上都存在著缺點(diǎn)。直接邊界元法雖然在船舶水動(dòng)力領(lǐng)域有著非常廣闊的應(yīng)用前景，但至今應(yīng)用較少。為推廣直接邊界元法在船舶水動(dòng)力學(xué)中的應(yīng)用，根據(jù)邊界積分法建立積分方程，采用直接邊界元法對無界勢流繞流問題予以求解，得出流場速度勢和物面上的速度分布，并通過與解析解的比較進(jìn)行誤差分析。對二維、三維問題的算例進(jìn)行數(shù)值計(jì)算。數(shù)值計(jì)算過程用Matlab編程實(shí)現(xiàn)。結(jié)果表明：直接邊界元法在求解船舶勢流繞流問題中具有足夠的精度和較高的效率，且數(shù)值計(jì)算實(shí)現(xiàn)過程更簡潔，可發(fā)展成為求解船舶興波等船舶水動(dòng)力學(xué)問題的通用方法。

直接邊界元法；勢流理論；數(shù)值積分；船舶水動(dòng)力學(xué)

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/42.1755.TJ.20150128.1201.005.html

期刊網(wǎng)址：www.ship-research.com

引用格式：李井煜，盧曉平，趙鵬偉.直接邊界元法解勢流速度場問題［J］.中國艦船研究，2015，10（1）：68-75. LI Jingyu，LU Xiaoping，ZHAO Pengwei.Direct boundary element method for the problem of potential flow velocity field［J］.Chinese Journal of Ship Research，2015，10（1）：68-75.

0　引言

船舶水動(dòng)力學(xué)中的許多問題都可以歸結(jié)為按初始條件和邊界條件求解偏微分方程的初值和邊值問題，但能求出解析解的僅限于極少數(shù)情況，故一般只用數(shù)值方法求出近似解。邊界元法是一種很好的求解船舶水動(dòng)力學(xué)問題的數(shù)值計(jì)算方法，經(jīng)過幾十年的發(fā)展，其已經(jīng)在諸多領(lǐng)域得到廣泛應(yīng)用。相比有限元法等其他數(shù)值計(jì)算方法，邊界元法的優(yōu)勢在于利用了“降維”功能，從而可以大大減少對計(jì)算機(jī)存儲的需求［1-2］。采用邊界元法僅需在計(jì)算域的邊界上進(jìn)行求解即可，其可將三維問題轉(zhuǎn)化為二維問題，或者將二維問題轉(zhuǎn)化為一維問題，且能方便處理無界區(qū)域問題，這也是邊界元法在求解船舶水動(dòng)力學(xué)勢流繞流問題中得到廣泛應(yīng)用的重要原因之一。

邊界元法分直接法和間接法2類。直接法是用物理意義明確的變量來建立積分方程，積分方程中的未知函數(shù)就是物理量在邊界上的值；間接法是用物理意義不一定明確的變量來建立積分方程，例如源分布或偶極子分布函數(shù)的概念，都是基于物理直觀提出，其在理論和數(shù)值計(jì)算上都存在著缺點(diǎn)，不便于在工程應(yīng)用中推廣。邊界元法在船舶水動(dòng)力問題上的應(yīng)用可追溯到Hess-Smith方法。該方法是一種以源分布函數(shù)為未知函數(shù)的間接方法，或許是基于這個(gè)原因，后續(xù)在船舶興波、耐波性和操縱性問題的勢流理論求解中大多采用以Hess-Smith法為基礎(chǔ)的間接法，長期以來，該方法在求解船舶水動(dòng)力學(xué)勢流問題中發(fā)揮著重要作用［3-5］。

針對船舶水動(dòng)力學(xué)數(shù)值求解中的這種趨勢，提出采用直接邊界元法求解船舶水動(dòng)力的勢流繞流問題，直接將變量取為流場的速度勢，然后求解流動(dòng)區(qū)域或邊界上的速度場，并采用簡單的算例進(jìn)行數(shù)值計(jì)算。通過比較數(shù)值解和解析解［6］，分析計(jì)算精度。對于二維流動(dòng)，選取無限域的均勻來流為計(jì)算模型，對于三維流場，則選取無限域圓球繞流為計(jì)算模型。所做研究可為進(jìn)一步推廣使用直接邊界元法求解帶有興波面、附體或多片體等復(fù)雜邊界條件的船舶勢流繞流問題打下基礎(chǔ)。

1　二維勢流問題算例

1.1數(shù)學(xué)模型

考慮無窮遠(yuǎn)邊界的流場，有沿y軸正向的速度為1的均勻來流，流體為無旋、無粘性、不可壓縮的理想流體。研究邊長為1的正方形區(qū)域流場，定解問題如下：

根據(jù)第三格林公式，可得到相應(yīng)的邊界積分方程為

式中：Γ為邊界線；φ為速度勢；φ*為拉普拉斯算子的基本解（或稱格林函數(shù)），在二維問題中

其中，r（P，Q）為場點(diǎn)P與源點(diǎn)Q之間的距離；C為

其中，θ在二維中為邊界點(diǎn)P處邊界切線的夾角，在三維中為邊界點(diǎn)P的立體角。當(dāng)邊界Γ光滑時(shí)，C（P）=1/2。當(dāng)場點(diǎn)P取在邊界上時(shí)，式（2）即為求解邊界上速度勢φ的邊界積分方程，根據(jù)邊界上的速度勢，可進(jìn)一步解出流場中任意點(diǎn)的速度勢。

1.2離散與數(shù)值求解

要采用直接邊界元法求解式（1）所示的邊界積分方程，需對邊界劃分單元。如圖1所示，將正方形邊界劃分為20個(gè)單元，在每個(gè)單元上選取插值函數(shù)模式，將邊界積分方程離散后，即轉(zhuǎn)化為線性方程組，據(jù)此，便可對邊界上的速度勢φ進(jìn)行求解。

圖1　正方形網(wǎng)格區(qū)域Fig.1　A square flow field area

1.2.1常數(shù)單元

對于常數(shù)單元，單元節(jié)點(diǎn)位于單元中點(diǎn)處。如圖2所示，將每個(gè)單元上的物理量等）都設(shè)為常量，且取為節(jié)點(diǎn)上的值。

圖2　常數(shù)單元示意圖Fig.2　Constant element

此時(shí)，邊界積分方程（2）便轉(zhuǎn)為以下線性方程組：

或

對于節(jié)點(diǎn)（本例中為20個(gè)），線性方程組（3）可表示為如下矩陣形式：

式（5）中的系數(shù)矩陣H和G的計(jì)算方式如下：

1）當(dāng)i=j時(shí)，場點(diǎn)P分布在單元內(nèi)，位于單元節(jié)點(diǎn)的位置。因積分計(jì)算帶有奇異點(diǎn)，按柯西主值積分，可以推導(dǎo)得到

式中，l為單元長度。

2）當(dāng)i≠j時(shí)，

式中：hij和gij為場點(diǎn)在第i個(gè)單元時(shí)第j個(gè)單元產(chǎn)生的影響系數(shù)；R為場點(diǎn)P到源點(diǎn)Q的向量矢徑；n為單元的單位法向量，指向方形區(qū)域外部。積分式（9）和式（10）采用4點(diǎn)高斯積分進(jìn)行數(shù)值計(jì)算。本算例的計(jì)算結(jié)果如表1和表2所示。

在劃分20個(gè)網(wǎng)格單元，采用常數(shù)單元離散，積分方法采用4點(diǎn)高斯積分公式的情況下：速度勢的最大相對誤差為5.225 4%，速度勢的平均相對誤差為1.365 3%，速度的最大相對誤差為5.415 3%，速度的平均相對誤差為2.491 2%。

表1　常數(shù)單元速度勢計(jì)算結(jié)果Tab.1The calculation results of constant element velocity potential

表2　常數(shù)單元速度計(jì)算結(jié)果Tab.2The calculation results of constant element velocity

1.2.2線性單元

對于線性單元，其節(jié)點(diǎn)設(shè)在單元端點(diǎn)處。如圖3所示，將物理量在單元上進(jìn)行線性插值處理，建立單元上任意點(diǎn)的函數(shù)值與單元節(jié)點(diǎn)的關(guān)系。

式中：ξ為邊界元的局部坐標(biāo)，首尾兩節(jié)點(diǎn)的局部坐標(biāo)值分別為0和1；φ1，φ2和分別為首尾兩節(jié)點(diǎn)的值與值；Φ為插值基函數(shù)。

圖3　計(jì)算影響系數(shù)時(shí)需代入邊界條件的單元Fig.3　The element needing to concern boundary conditions when calculating the influence coefficients

將式（11）代入邊界積分方程（2）并進(jìn)行簡化，得

其中：

與常數(shù)單元類似，將式（13）簡化成如式（5）一樣的矩陣形式HU=GN后，代入已知邊界條件，即可求解得出邊界上未知的φ和

另需注意，對于線性單元，有些單元需進(jìn)行特殊處理。如圖3所示，對于本算例的角點(diǎn)單元：一個(gè)節(jié)點(diǎn)的已知條件屬于本質(zhì)邊界條件，另一節(jié)點(diǎn)屬于自然邊界條件。對這樣的單元做插值處理將會(huì)與已知條件沖突，因此在進(jìn)行程序設(shè)計(jì)時(shí)需特別注意，否則，將導(dǎo)致這類單元節(jié)點(diǎn)的計(jì)算結(jié)果出現(xiàn)較大誤差。

為便于與常數(shù)單元比較，積分方法依然選用高斯積分法，結(jié)果如表3和表4所示。

表3　線性單元速度勢計(jì)算結(jié)果Tab.3The calculation results of linear element velocity potential

表4　線性單元速度計(jì)算結(jié)果Tab.4The calculation results of linear element velocity

在劃分20個(gè)網(wǎng)格單元，采用線性單元離散，積分方法采用4點(diǎn)高斯積分公式的情況下：速度勢的最大相對誤差為0.543 9%，速度勢的平均相對誤差為0.137 3%，速度的最大相對誤差為5.415 3%，速度的平均相對誤差為2.491 2%。

1.3常數(shù)單元與線性單元比較

從以上結(jié)果可以看出，采用常數(shù)或線性單元的直接邊界元法解勢流速度場可以得出有效的計(jì)算結(jié)果。其誤差主要是由高斯積分法近似引起，用Matlab工具箱中的int符號積分函數(shù)［7］可以提高精度，但運(yùn)算效率較低。也可以通過使用更合適的數(shù)值積分方法來提高計(jì)算精度。

常數(shù)單元在數(shù)學(xué)處理和編程實(shí)現(xiàn)上都較為簡便，在單元數(shù)相同的條件下，常數(shù)單元較線性單元精度低，但是可以通過增加單元數(shù)來提高計(jì)算精度。線性單元的節(jié)點(diǎn)設(shè)在單元邊界上，需進(jìn)行角點(diǎn)處理，因而對復(fù)雜邊界的適應(yīng)性不強(qiáng)。經(jīng)權(quán)衡，初步認(rèn)為總體上采用常數(shù)單元效率較高。

2　三維問題計(jì)算

2.1數(shù)學(xué)模型

考慮無窮遠(yuǎn)邊界的流場［8-9］，有沿x軸負(fù)向的、速度為1的均勻來流。流體無旋、無粘性、不可壓縮，流場中有一個(gè)半徑為1的圓球。流場速度勢可表示為

式中：Φ0為總速度勢；V∞為無窮遠(yuǎn)處的來流速度；φ為圓球的擾動(dòng)速度勢。若直接求解Φ0，由于物面不可穿透條件，離散化之后的代數(shù)方程組為齊次線性方程組，不能得出唯一的非零解。因此，將求解的直接變量選擇為擾動(dòng)速度勢φ，擾動(dòng)速度勢φ滿足定解方程

式中：Ω為流場區(qū)域；S為物面邊界；nQ為Q點(diǎn)的單位法向量；nx為法向量在x軸上的分量。

當(dāng)邊界光滑時(shí)，C（P）=1/2。

2.2離散與數(shù)值求解

（3）云計(jì)算技術(shù)的使用可以有效的提升各個(gè)采油廠的信息化水平，采油廠的信息化減少的投入和油田公司相差甚遠(yuǎn)，對此，就可以使用SaaS模式進(jìn)行處理，構(gòu)建采油廠的信息化平臺，讓油田公司的數(shù)據(jù)中心扮演采油廠的云服務(wù)商一角色，以各個(gè)采油廠的實(shí)際業(yè)務(wù)需求為基準(zhǔn)，實(shí)行基礎(chǔ)設(shè)施的托管服務(wù)以及SaaS系統(tǒng)服務(wù)。只有這樣才可以更為高效且快速的去打造各類采油廠信息化的平臺，從根源上滿足采油廠的實(shí)際應(yīng)用需求。

因三角形單元對復(fù)雜邊界的適應(yīng)性強(qiáng)，無需對物面的點(diǎn)進(jìn)行坐標(biāo)變換處理，所以選擇將物面劃分為三角形單元，如圖4所示。由于三維問題的計(jì)算量比二維問題的大，為提高編程和運(yùn)算效率，采用常數(shù)單元求解。

圖4　圓球劃分三角網(wǎng)格效果Fig.4　The effect of sphere divided into triangular meshes

由邊界積分方程（16）離散得出線性方程組的過程與二維問題類似，見式（3）～式（5）。最后，簡化為與式（5）一樣的矩陣形式HU=GN，然后把邊界條件·nx代入N中，即可求解得出圓球物面各個(gè)單元節(jié)點(diǎn)的速度勢，進(jìn)而可進(jìn)一步求解出物面的速度場。

2.3系數(shù)矩陣與速度求解

影響系數(shù)矩陣的表達(dá)式和計(jì)算方法的說明如下：

1）當(dāng)i=j時(shí)，場點(diǎn)P分布在單元內(nèi)，涉及到高維奇異積分的處理。由于在單元內(nèi)有，所以

方法1：采用3點(diǎn)高斯積分式。

因?yàn)槠纥c(diǎn)只在三角區(qū)域的中心處存在，而3點(diǎn)高斯積分法需要的積分點(diǎn)都是在三角形的邊界上，所以避開了奇點(diǎn)。如圖5所示。

式中：f為被積函數(shù)；A為三角形區(qū)域的面積；w為權(quán)系數(shù)；ξ1，ξ2，ξ3為三角形面積坐標(biāo)［1，10］。

圖5　高斯積分點(diǎn)示意圖Fig.5　Sketch of Gauss integral points

方法2：采用在四邊形上均勻分布奇點(diǎn)的誘導(dǎo)速度公式。

引入坐標(biāo)系oξηζ，原點(diǎn)一般取在四邊形ΔS的形心處，坐標(biāo)平面ξηζ即為四邊形ΔS所在平面，ζ軸的正向?yàn)棣的法線方向。ΔS的四個(gè)頂點(diǎn)p1，p2，p3，p4按逆時(shí)針方向排列，頂點(diǎn)pi的坐標(biāo)為（ξi，ηi，0）（i=1，2，3，4）。奇點(diǎn)的坐標(biāo)p的坐標(biāo)為（x，y，z）。給出單層位勢計(jì)算公式［5］：

式中：li，j為四邊形邊長，其中i，j為頂點(diǎn)編號；ri為各頂點(diǎn)到奇點(diǎn)p的距離。

把三角形視為某兩個(gè)頂點(diǎn)重合的四邊形進(jìn)行計(jì)算，發(fā)現(xiàn)與高斯積分法計(jì)算結(jié)果相比，其節(jié)點(diǎn)速度勢最大誤差提高了，平均誤差降低了。

方法3：采用Matlab函數(shù)quad2d（）處理匿名函數(shù)積分［7］。

2009版以后的Matlab軟件帶有二重積分函數(shù)quad2d（），可以近似計(jì)算奇異積分，運(yùn)算時(shí)的效率也挺高，但和高斯積分相比要慢。其精度與方法2的處理結(jié)果相近。

2）當(dāng)i≠j時(shí)，

式中，n為單元的單位法向量，指向流場外部（物面內(nèi)部）。該積分式因不含奇異性，故結(jié)果易收斂。

求得物面控制點(diǎn)處的勢函數(shù)值后，文獻(xiàn)［11］提供了一種簡便的用數(shù)值微分方法計(jì)算物面上速度分布的方法。設(shè)物面方程為y=y（x，z），勢函數(shù)在物面上的增量可以寫成

式中：。在控制點(diǎn)附近任意找兩點(diǎn)，兩次使用式（23），便可求解出u*，w*。用下式求解速度矢量：

式中：u，v，w為控制點(diǎn)速度的3個(gè)分量；l，m，k為控制點(diǎn)單位法向量的3個(gè)分量。

由于單元數(shù)對計(jì)算結(jié)果的影響很大，故本文選取不同的單元數(shù)分別采用上述方法進(jìn)行了試驗(yàn)。以7點(diǎn)的高斯積分法計(jì)算系數(shù)矩陣的結(jié)果如圖6～圖11所示。

當(dāng)單元數(shù)為360時(shí)，速度勢的最大誤差為14.683 7%，平均誤差為3.002 5%；速度的最大誤差為12.750 0%，平均誤差為2.080 0%。

基于網(wǎng)格的劃分方法，球的兩個(gè)極點(diǎn)附近的單元數(shù)較少，因而其誤差較其它位置節(jié)點(diǎn)處的突出。由于所使用二維積分式的數(shù)值積分精度不高，因此三維問題的計(jì)算結(jié)果不如二維問題的好，可以考慮結(jié)合采用改進(jìn)網(wǎng)格劃分和數(shù)值積分的方法提高計(jì)算精度。

圖6　360個(gè)網(wǎng)格速度勢計(jì)算結(jié)果Fig.6　The velocity potential calculation results of 360 grids

圖7　360個(gè)網(wǎng)格速度計(jì)算結(jié)果Fig.7　The velocity calculation results of 360 grids

圖9　 1230個(gè)網(wǎng)格速度勢計(jì)算結(jié)果Fig.9　The velocity potential calculation results of 3 120 grids

圖11　 3120個(gè)網(wǎng)格速度矢量Fig.11　The velocity vector results of 3 120 grids

圖10　 3120個(gè)網(wǎng)格速度計(jì)算結(jié)果Fig.10　The velocity calculation results of 3 120 grids

圖8　360個(gè)網(wǎng)格速度矢量Fig.8　The velocity vector results of 360 grids

速度幅值計(jì)算結(jié)果的精度是可以接受的，但該計(jì)算結(jié)果未能精確反映解析解給出的速度變化趨勢，例如，駐點(diǎn)和速度最大值的位置未能精確反映出來。這主要是由數(shù)值計(jì)算誤差所致，不僅物面速度計(jì)算存在誤差，而且速度勢函數(shù)計(jì)算本身也存在誤差。

當(dāng)單元數(shù)為3 120個(gè)時(shí)，速度勢最大誤差為5.807 8%，平均誤差為0.621 8%；速度最大誤差為7.090 0%，平均誤差為0.120 0%。隨著單元數(shù)的增加，計(jì)算精度顯著提高。

3　結(jié)語

對應(yīng)用于船舶水動(dòng)力學(xué)數(shù)值計(jì)算中的邊界元法，多采用以Hess-smith方法為基礎(chǔ)的間接邊界元法。Hess-smith方法是先求解各單元的源強(qiáng)密度，使其與系數(shù)矩陣相乘而得到速度勢，由于其借助了作為中間變量的源強(qiáng)密度，從而使得這類間接邊界元法在編程和計(jì)算效率上不如直接邊界元法。另外，相對于間接邊界元法，直接邊界元法的數(shù)學(xué)意義更嚴(yán)謹(jǐn)，物理意義更明確。

邊界元法的計(jì)算量主要集中在影響系數(shù)矩陣的計(jì)算上，誤差也大多源于此。本文對多種數(shù)值積分方法進(jìn)行了數(shù)值試驗(yàn)，分析了各種方法對最終結(jié)果誤差的影響。改進(jìn)本文方法或使用精度更高的數(shù)值積分算法，是提高計(jì)算結(jié)果精度的重要途徑。

本文研究了無限域的勢流問題，通過求解經(jīng)典算例并與解析解的比較，證明了采用直接邊界元法求解船舶水動(dòng)力學(xué)勢流繞流問題的可行性。直接邊界元法的數(shù)值計(jì)算和程序設(shè)計(jì)過程具有通用性，適用于任何無限域三維勢流的速度勢和速度的求解。若進(jìn)一步引入興波條件和多片體干擾效應(yīng)等，還可求解更廣泛的船舶水動(dòng)力學(xué)勢流繞流數(shù)值計(jì)算問題。

［1］王元淳.邊界元法基礎(chǔ)［M］.上海：上海交通大學(xué)出版社，1988：6-33.

［2］姚壽廣.邊界元數(shù)值方法及其工程應(yīng)用［M］.北京：國防工業(yè)出版社，1995：1-24.

［3］鄧小敏.多體船計(jì)算及船型優(yōu)化研究［D］.大連：大連理工大學(xué)，2012.

［4］李子如.用面元法預(yù)報(bào)船體興波阻力［D］.武漢：武漢理工大學(xué)，2006.

［5］戴遺山，段文洋.船舶在波浪中運(yùn)動(dòng)的勢流理論［M］.北京：國防工業(yè)出版社，2008：24-28.

［6］張兆順，崔桂香.流體力學(xué)［M］.2版.北京：清華大學(xué)出版社，2006：130-131.

［7］陳澤，占海明.詳解MATLAB在科學(xué)計(jì)算中的應(yīng)用［M］.北京：電子工業(yè)出版社，2011.

［8］余靈，李干洛.球尾漁船表面流場的數(shù)值計(jì)算［J］.華南理工學(xué)報(bào)（自然科學(xué)版），1995，23（10）：155-163. YU Ling，LI Ganluo.Numerical calculation of flow fields over ship hulls of bulbous stern fishing boat［J］. Journal of South China University of Technology（Natu?ral Science），1995，23（10）：155-163.

［9］余靈，李干洛，羊少剛.船體表面流場的理論計(jì)算與數(shù)值分析［J］.廣東造船，1994（3）：1-7.

［10］HIRSCH C.Numerical computation of internal and external flows［M］.2nd ed.Butterworth-Heinemann，2007：221-223.

［11］王獻(xiàn)孚.計(jì)算船舶流體力學(xué)［M］.上海：上海交通大學(xué)，1992：244-245.

［責(zé)任編輯：盧圣芳］

Direct Boundary Element Method for the Problem of Potential Flow Velocity Field

LI Jingyu，LU Xiaoping，ZHAO Pengwei
Department of Naval Architecture Engineering，Naval University of Engineering，Wuhan 430033，China

In the field of ship hydrodynamics，scholars usually adopt Hess-Smith method to solve the po?tential flow problem.However，this method is essentially based on the physical intuition，and there are shortcomings concerning this theory in numerical calculation.Since the direct boundary element method is rarely used in ship hydrodynamic problems，the presented method in this paper has broad application pros?pects in the field of ship hydrodynamics and may promote the direct boundary element method.According to the boundary integral method，an integral equation is established to solve the unbounded potential flow problem，and the flow field velocity distribution on the surface of the velocity potential is then obtained.A comparison with the analytical solution is finally conducted for error analysis.In this paper，both 2D and 3D numerical examples are provided and programmed to realize the numerical calculation process in Mat?lab.The results show that the direct boundary element method has decent precision and efficiency，and the numerical implementation is even more concise.Moreover，this method can be developed into general forms to solve other ship dynamic problems.

direct boundary element method；potential flow theory；numerical integration；ship hydrody?namics

中國分類號：U661.1A

10.3969/j.issn.1673-3185.2015.01.010

2014-08-18

網(wǎng)絡(luò)出版時(shí)間：2015-1-28 12:01

國家部委基金資助項(xiàng)目

李井煜，男，1990年生，碩士生。研究方向：艦船流體動(dòng)力性能。E?mail：935228691@qq.com

盧曉平（通信作者），男，1957年生，博士，教授。研究方向：艦船流體動(dòng)力性能