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

李井煜，盧曉平，趙鵬偉

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

李井煜，盧曉平，趙鵬偉

海軍工程大學艦船工程系，湖北武漢430033

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

直接邊界元法；勢流理論；數值積分；船舶水動力學

網絡出版地址：http://www.cnki.net/kcms/detail/42.1755.TJ.20150128.1201.005.html

期刊網址：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　引言

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

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

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

1　二維勢流問題算例

1.1數學模型

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

根據第三格林公式，可得到相應的邊界積分方程為

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

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

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

1.2離散與數值求解

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

圖1　正方形網格區域Fig.1　A square flow field area

1.2.1常數單元

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

圖2　常數單元示意圖Fig.2　Constant element

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

或

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

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

1）當i=j時，場點P分布在單元內，位于單元節點的位置。因積分計算帶有奇異點，按柯西主值積分，可以推導得到

式中，l為單元長度。

2）當i≠j時，

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

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

表1　常數單元速度勢計算結果Tab.1The calculation results of constant element velocity potential

表2　常數單元速度計算結果Tab.2The calculation results of constant element velocity

1.2.2線性單元

對于線性單元，其節點設在單元端點處。如圖3所示，將物理量在單元上進行線性插值處理，建立單元上任意點的函數值與單元節點的關系。

式中：ξ為邊界元的局部坐標，首尾兩節點的局部坐標值分別為0和1；φ1，φ2和分別為首尾兩節點的值與值；Φ為插值基函數。

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

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

其中：

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

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

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

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

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

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

1.3常數單元與線性單元比較

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

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

2　三維問題計算

2.1數學模型

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

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

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

當邊界光滑時，C（P）=1/2。

2.2離散與數值求解

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

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

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

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

2.3系數矩陣與速度求解

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

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

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

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

式中：f為被積函數；A為三角形區域的面積；w為權系數；ξ1，ξ2，ξ3為三角形面積坐標［1，10］。

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

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

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

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

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

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

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

2）當i≠j時，

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

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

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

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

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

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

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

圖6　360個網格速度勢計算結果Fig.6　The velocity potential calculation results of 360 grids

圖7　360個網格速度計算結果Fig.7　The velocity calculation results of 360 grids

圖9　 1230個網格速度勢計算結果Fig.9　The velocity potential calculation results of 3 120 grids

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

圖10　 3120個網格速度計算結果Fig.10　The velocity calculation results of 3 120 grids

圖8　360個網格速度矢量Fig.8　The velocity vector results of 360 grids

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

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

3　結語

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

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

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

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

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

［3］鄧小敏.多體船計算及船型優化研究［D］.大連：大連理工大學，2012.

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

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

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

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

［8］余靈，李干洛.球尾漁船表面流場的數值計算［J］.華南理工學報（自然科學版），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］.廣東造船，1994（3）：1-7.

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

［11］王獻孚.計算船舶流體力學［M］.上海：上海交通大學，1992：244-245.

［責任編輯：盧圣芳］

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

網絡出版時間：2015-1-28 12:01

國家部委基金資助項目

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

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