無限區域二維勢流直接邊界元法精度分析

韓玉超，盧曉平，王　中海軍工程大學艦船工程系，湖北武漢430033

無限區域二維勢流直接邊界元法精度分析

韓玉超，盧曉平，王中
海軍工程大學艦船工程系，湖北武漢430033

邊界元法作為一種重要的數值方法已在許多領域得到廣泛應用，但在船舶水動力勢流理論數值計算方面，有關直接邊界元法的研究并不充分，尤其是在船舶興波阻力勢流理論求解方面，以往的“面元法”通常是基于Hess-Smith法的間接法，這類方法在理論和數值計算上都存在著缺陷。針對船舶水動力勢流理論計算，采用直接邊界元法，對二維勢流無界繞流算例進行系統的數值計算，并根據二維勢流問題的計算結果詳細探討邊界單元離散形式和單元上的數值積分方法對計算精度的影響，各項數值計算均以Matlab軟件編程實現。結果表明，采用常數單元和龍貝格積分法能夠得到較準確的結果，且計算速度較快。

直接邊界元法；勢流理論；邊界離散；數值積分

0　引言

根據邊界積分方程中未知函數的不同，邊界元法可分為直接法和間接法2種。在以往的船舶興波阻力勢流理論問題求解方面，通常是采用基于Hess-Smith法的間接法［1］。在間接法中，邊界積分方程中的未知函數并不是待定函數本身，而是作用于所考察區域邊界相應線（面）上的某種虛擬的分布量，由于待定函數是通過求出虛擬分布量后間接求出，且沒有嚴格的數學定義，故間接邊界元法在理論和數值計算上都存在著缺陷。在直接法中，邊界積分方程建立了待定函數在區域內的值與邊界上的值之間的關系，具有直觀的數學意義，并且一旦待定函數的邊界值全部確定，其區域內的值也隨之確定［2-3］，這樣既減少了計算工作量，又便于結果的整理與分析。

通常，邊界元法誤差的主要來源是邊界的離散誤差和積分的計算誤差［4］。在邊界離散方面，單元的剖分是重要的一環。而在數值積分方面，誤差主要包括數值積分的誤差和奇異積分的處理，在奇異積分方面，許多學者已做了大量工作［5-7］。因此，找到一種快速且準確的邊界離散形式和數值積分方法對進一步提高邊界元法的計算精度將有很大的幫助。

本文將針對船舶水動力勢流理論計算［8-9］，采用直接邊界元法，先給出二維拉普拉斯方程邊值問題的邊界積分方程，然后分別采用常數單元、線性單元和二次單元對積分邊界進行離散。在數值計算方面，將分別采用多種常用的數值積分方法對問題進行求解分析［10］，然后將所得數值解與解析解進行比較，著重討論不同類型邊界單元和不同數值積分方法對直接邊界元法計算精度的影響，并選出一種適合于求解一般問題的邊界離散形式和數值積分方法，進而推廣至三維無限域問題，這對三維問題離散形式和積分方法的選取具有一定的指導意義。

1　數學模型

考慮在無限水深條件下無旋、非粘性、不可壓縮來流中的任意物體，取一外部控制面將其封閉在內。流域的邊界面由物面L組成，在邊界上存在擾動速度勢φ，在整個流場中滿足流體質量守恒定理的二維拉普拉斯方程：

對于二維圓柱無界繞流，如圖1所示，計算邊界僅為物體表面LQ，需滿足物面邊界條件：

式中：φ為場點Q的待定函數，即φ(Q)；qˉ為邊界LQ上勢函數u的法向導數的已知值；V∞為無窮遠處均勻來流速度；nQ為邊界上的單位內法向量。

圖1　二維圓柱繞流示意圖Fig.1 Two dimensional flow around cylinder

二維圓柱無界繞流場中任意位置的速度勢函數的解析解為

根據基本解的定義，二維拉普拉斯方程的基本解u*應滿足

對于二維圓柱繞流，當場點P位于邊界面上時，由物面條件和格林定理，可得場點P的速度勢為

式中，r為源點Q與場點P之間的距離。

2　邊界單元的劃分

在邊界離散方面，針對不同的具體問題，對邊界量可以采用不同類型的插值函數，例如零次插值函數、線性插值函數和二次插值函數等，其相應的邊界單元稱為常數單元、線性單元和二次單元。

2.1常數單元

對于常數單元，如圖2所示，節點取在邊界單元的中點，邊界上的函數值u和函數的法向導數值q在邊界單元上為常量，并等于節點處的值。

圖2　常數單元節點示意圖Fig.2 Nodes of constantelement

在每個單元上對邊界積分方程進行離散，可得

式中：uj和qj分別為邊界單元節點處的函數值和函數的法向導數值；hij和gij分別為系數矩陣H 和G中的元素，其表達式為：當i=j時，

式中：li為邊界單元的長度；dij為節點Pi到邊界單元Li的垂直距離；rij為節點到積分單元上點的距離。從節點P1開始，依次對所有節點的hij和gij進行求解，利用Matlab對公式進行編程并計算出結果。

2.2線性單元

對于線性單元，如圖3所示，節點取在邊界單元的兩端，邊界上的函數值u和函數的法向導數值q在邊界單元上服從線性分布，其離散化的邊界積分方程為

式中，uj，uj+1，qj，qj+1分別為邊界單元上2個端點處的函數值和函數的法向導數值。

圖3　線性單元節點示意圖Fig.3Nodesoflinearelement

每個線性單元可用無因次坐標ξ表示為

式中，(x1，y1)，(x2，y2)分別為線性單元上起點與終點的坐標值。

將邊界積分方程無因次化后可得hij和gij的表達式分別為：

式中：(xp，yp)為節點Pi的空間坐標；為雅克比行列式，對于二維線性單元，有J=lj。從節點P1開始，依次對各個節點的hij和gij進行求解，求解過程與常數單元相同。

當i=j時，邊界積分方程的形函數C由下式確定：

由于系數C只與節點Pi處邊界的幾何形狀有關，故兩線性單元間的夾角為非定值，如圖4所示。又因對角線元素Hii的計算包含有系數C和Cauchy型積分的計算，因此，采用一種簡便的均勻場的辦法便可間接計算C值。

圖4　線性單元節點處夾角示意圖Fig.4Anglearoundnodeoflinearelement

當一個均勻位勢場作用于整個邊界時，函數的法向導數值必定為零。因此，系數矩陣H的任意一行中所有元素的總和均為零，即

對于無限域問題，可設想所有內部邊界被一個無限大的圓LR所包圍，則除了出現在邊界積分方程中的所有內部邊界的積分外，還要加上在此大圓上的積分。

令R→∞，因

對邊界積分方程離散化后，可得

而Gii可推導出其積分的解析表達式為

2.3二次單元

對于二次單元，如圖5所示，節點分別取在邊界單元的兩端和中點，邊界上的u和q在邊界單元上服從二次分布。

圖5　二次單元節點示意圖Fig.5 Nodes of quadratic element

離散化后的邊界積分方程為

式中：uj，uj+1，uj+2分別為邊界單元上3個節點處的函數值；qj，qj+1，qj+2為其節點處函數的法向導數值。

每個二次單元可用無因次坐標ξ表示為

式中，(x1，y1)，(x2，y2)和(x3，y3)分別表示二次單元上3個節點的坐標值。

用二次單元離散后的邊界積分方程的求解步驟與前2種單元均相同，只是邊界單元中的積分函數有所不同，將邊界積分方程無因次化后，可得hij和gij的表達式分別為：

當i=j時，與線性單元相同，Hii可由式（17）和式（19）求得。而Gii由于公式過于復雜，無法推導出解析解表達式，而且在單元端點處的積分存在奇異性，所以只能采取積分上的近似解法，即取一趨近于0的極小值ε。在對g1ii和gi3i積分時，積分區間取為[1-ε，1+ε]；對gi2i積分時，積分區間取為[0，0.5-ε]和[0.5+ε，1]，按上述區間對gii的積分項進行積分，即可得到近似解。

值得說明的是，由于高斯積分法所取的高斯積分點一般在積分區間內，所以在積分過程中不需要求出帶奇異性的端點值，在不改變積分區間的前提下可直接避開對積分奇異性的處理，并且多點高斯積分在計算精度方面也能滿足一定的要求。所以，高斯積分法在處理二次單元的問題上存在著天然優勢，該問題在后文中將予以詳細討論。

3　圓柱繞流算例及其精度分析

3.1數值積分方法與計算精度

用前文介紹的邊界單元離散形式和數值積分方法對二維圓柱無界繞流問題進行求解。采用不同的節點數和數值積分方法對每種單元進行計算，邊界單元離散的節點數依次從8取至200。之后，對每個節點處的計算值與解析值進行比較，求出相對誤差并統計每次計算的時間，其結果如圖6～圖8所示。

由常數單元的計算結果可知，大部分數值積分方法從低節點數開始就能獲得較高的精度，且隨著節點數的增加，計算精度不斷提高，但相應的計算時間也呈幾何增長。其中，5點高斯積分公式法、龍貝格積分法等在精度和耗時方面都能獲得比較滿意的結果。在節點數為200時，其平均誤差均能達到10-3～10-5量級的精度，可以滿足工程應用需要。

圖6　采用常數單元時各種積分方法的相對誤差和計算時間Fig.6 Relative errors and calculation time of various integralmethodswhen using constantelement

圖7　采用線性單元時各種積分方法的相對誤差和計算時間Fig.7 Relative errors and calculation time of various integralmethodswhen using linear element

圖8　采用二次單元時各種積分方法的相對誤差和計算時間Fig.8 Relative errors and calculation time ofvarious integralmethodswhen using quadratic element

對于線性單元和二次單元，一部分數值積分方法已不適用于邊界積分方程的求解，有的計算結果會出現較大的波動和偏差，有的甚至會發散。尤其是在二次單元時，例如，逐次半分區間梯形公式法等已經無法保證計算精度，平均誤差達到了5%以上，且計算時間也呈指數增長，已超出邊界元法可接受的范圍。

綜上所述，本文認為龍貝格積分法比較適合邊界元法邊界積分方程的求解，其在計算精度和計算時間上均能獲得比較理想的結果。

3.2邊界單元離散形式與計算精度

邊界元法誤差的另一個重要來源是邊界的離散，單元的選擇和劃分對于計算結果的精度有著很大的影響，甚至超過了積分方法對精度的影響。在邊界離散方面，一般公認的觀點是邊界單元數劃分越多則計算精度越高，但邊界單元的離散形式對精度的影響并沒有一個定性的認識。本節將繼續采用二維勢流問題的計算結果，對使用相同數值積分方法不同離散形式下的計算精度和計算時間進行比較，比較結果如圖9所示。

圖9　不同離散方式下的相對誤差和計算時間比較Fig.9 Relative errors and calculation time of various boundary discretizations

由計算結果可知，在3種邊界單元離散形式中，常數單元的計算誤差一般都要小于其余2種，且隨著節點數的增多，常數單元的精度優勢也更加明顯。當劃分的節點數達到200時，常數單元的誤差基本上要比另2種單元的小一個數量級。此外，常數單元在計算時間上也遠遠優于其他單元，所耗時間通常只是其他單元的5%～10%。

從解析的角度分析，常數單元之所以優于線性單元和二次單元，很大一部分原因在于其對邊界函數的近似程度。常數單元雖然將單元上的函數值設為了定值，但在網格足夠多的前提下，常數單元反而能更好地表達離散物體表面的變化趨勢。而線性單元和二次單元只是直接把單元上的函數值規定為按線性或二次變化，并不一定能完全反映函數實際的變化趨勢。且線性單元和二次單元在具體的計算中會受到數值計算方法的限制，在含奇異積分的計算中無法推導出解析解，只能采用逼近奇點的近似解法。綜合考慮單元數和離散形式，常數單元可以得出更準確的計算結果。

因此，本文認為在網格單元劃分足夠多的情況下，常數單元可以取代線性單元和二次單元對邊界進行離散，并且計算誤差和計算時間均可以滿足邊界元法在工程實際中的應用。

4　結語

邊界元法由于采用了解析基本解，在計算精度上有限元法等數值方法難以與其相比。但許多學者由于對邊界元法的誤差來源了解不深，往往習慣性地采用高次單元，單純地認為離散單元階次越高，越能符合邊界條件的變化，但這樣反而不能把邊界元法在計算精度和速度上的優勢充分發揮出來。

本文首先通過算例證明了邊界元法的主要誤差是邊界離散誤差和數值積分誤差，可見，正確地選擇邊界單元離散形式和數值積分方法是保證邊界元法計算精度的第一要務。隨后，驗證了不同數值積分方法在求解邊界積分方程中的適用性，通過與常用的高斯積分法進行對比，最終選出了一種精度較高、速度較快的積分方法——龍貝格積分法作為后續計算采用的方法。最后，通過對3種邊界離散形式計算精度的分析，認為在網格單元劃分足夠多的情況下，常數單元可以取代線性單元和二次單元作為一般二維邊界元問題的邊界離散形式。從理論上說，該結論也可以推廣到三維問題中，即常數單元在三維邊界元問題中也能獲得最佳的精度，但該推論的正確性還有待后續的研究驗證。

［1］陳材侃.計算流體力學［M］.重慶：重慶出版社，1992.

［2］BREBBIA C A.The boundary element method for engineers［M］.London：Pentech Press，1978.

［3］BREBBIA C A，TELLES JC F，WROBEL L.Boundary element techniques：theory and applications in engineering［M］.Berlin：Springer-Verlag，1984.

［4］姚振漢，王海濤.邊界元法［M］.北京：高等教育出版社，2010.

［5］KLASEBOER E，FERNANDEZ C R，KHOO B C.A note on true desingularisation of boundary integral methods for three-dimensional potential problems［J］. Engineering Analysis With Boundary Elements，2009，33（6）：796-801.

［6］周煥林，牛忠榮，程長征，等.各向異性位勢問題邊界元法中幾乎奇異積分的解析算法［J］.計算力學學報，2008，25（3）：333-338. ZHOU Huanlin，NIU Zhongrong，CHENG Changzheng，et al.Analytical algorithm of the nearly singular integrals in boundary elementmethod to anisotropic potential problems［J］.Chinese Journal of Computational Mechanics，2008，25（3）：333-338.

［7］馬榮，劉繼朝，石建省.邊界元法中角點問題的一種計算方法［J］.中國科學院研究生院學報，2011，28 （3）：315-321. MA Rong，LIU Jichao，SHI Jiansheng.A new algorithm to solve corner problem in the boundary element method［J］.Journal of the Graduate School of the Chinese Academy of Sciences，2011，28（3）：315-321.

［8］MIAO Q M，CHWANG A T.Ship waves by direct boundary elementmethod［J］.Journal of Ship Mechanics，2003，7（6）：27-36.

［9］BELIBASSAKISK A，GEROSTATHIS T P，KOSTAS K V，et al.A BEM-isogeometric method for the ship wave-resistanceproblem［J］.Ocean Engineering，2013，60（1）：53-67.

［10］MATHEWS JH，FINK K D.Numericalmethods using Matlab［M］.NJ：Prentice Hall，2004.

［責任編輯：盧圣芳］

Precision analysis of the two-dimensional potential flow problem in an in finite region with the direct boundary element method

HAN Yuchao，LU Xiaoping，WANGZhong Department of Naval Architecture Engineering，Naval University of Engineering，Wuhan 430033，China

The Boundary Element Method（BEM），as a key numerical method，has been widely applied in many fields.However，the research on the Direct Boundary Element Method（DBEM）for ship hydrodynamic numerical calculation problems is still insufficient，especially when it comes to the ship hydrodynamic potential flow theory.The general method-‘panel method'-is based on Hess-Smith method，which is an Indirect Boundary ElementMethod（IBEM）whosemajor flaws exist in both theory and numerical calculation.This paper，based on the ship hydrodynamic potential flow theory，adopts DBEM to calculate the example of two-dimensional unbounded potential flow around a cylinder，and analyzes the influence of the boundary element discrete forms and the numerical integralmethods on the calculation accuracy.The results carried out by Matlab clearly indicate that using the constant element and Romberg algorithm method could yield high calculation speed and accuracy.

Direct Boundary Element Method（DBEM）；potential flow theory；boundary discretization；numerical integration

U661.1

A

10.3969/j.issn.1673-3185.2015.04.006

2014-09-04網絡出版時間：2015-7-28 17:25:18

國家部委基金資助項目

韓玉超，男，1989年生，碩士生。研究方向：艦船水動力性能。E-mail：hanyuchaoa@163.com盧曉平（通信作者），男，1957年生，博士，教授。研究方向：艦船水動力性能。E-mail：Luxiaoping100@163.com王中，男，1981年生，博士，講師。研究方向：計算機輔助船舶設計制造。E-mail：wangzhonghj@sohu.com