李井煜，盧曉平，趙鵬偉

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

李井煜，盧曉平，趙鵬偉

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

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

直接邊界元法；勢流理論；數(shù)值積分；船舶水動力學(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　引言

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

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

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

1　二維勢流問題算例

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

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

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

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

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

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

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

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

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

1.2.1常數(shù)單元

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

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

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

或

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

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

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

式中，l為單元長度。

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

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

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

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

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

1.2.2線性單元

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

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

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

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

其中：

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

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

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

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

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

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

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

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

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

2　三維問題計算

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

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

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

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

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

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

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

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

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

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

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

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

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

方法1：采用3點高斯積分式。

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

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

圖5　高斯積分點示意圖Fig.5　Sketch of Gauss integral points

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

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

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

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

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

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

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

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

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

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

式中：u，v，w為控制點速度的3個分量；l，m，k為控制點單位法向量的3個分量。

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

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

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

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

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

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

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

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

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

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

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

3　結(jié)語

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

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

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

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

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

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

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

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

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

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

［8］余靈，李干洛.球尾漁船表面流場的數(shù)值計算［J］.華南理工學(xué)報（自然科學(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］余靈，李干洛，羊少剛.船體表面流場的理論計算與數(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］王獻孚.計算船舶流體力學(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ò)出版時間：2015-1-28 12:01

國家部委基金資助項目

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

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