無(wú)限區(qū)域二維勢(shì)流直接邊界元法精度分析

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

無(wú)限區(qū)域二維勢(shì)流直接邊界元法精度分析

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

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

直接邊界元法；勢(shì)流理論；邊界離散；數(shù)值積分

0　引言

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

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

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

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

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

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

式中：φ為場(chǎng)點(diǎn)Q的待定函數(shù)，即φ(Q)；qˉ為邊界LQ上勢(shì)函數(shù)u的法向?qū)?shù)的已知值；V∞為無(wú)窮遠(yuǎn)處均勻來(lái)流速度；nQ為邊界上的單位內(nèi)法向量。

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

二維圓柱無(wú)界繞流場(chǎng)中任意位置的速度勢(shì)函數(shù)的解析解為

根據(jù)基本解的定義，二維拉普拉斯方程的基本解u*應(yīng)滿足

對(duì)于二維圓柱繞流，當(dāng)場(chǎng)點(diǎn)P位于邊界面上時(shí)，由物面條件和格林定理，可得場(chǎng)點(diǎn)P的速度勢(shì)為

式中，r為源點(diǎn)Q與場(chǎng)點(diǎn)P之間的距離。

2　邊界單元的劃分

在邊界離散方面，針對(duì)不同的具體問(wèn)題，對(duì)邊界量可以采用不同類型的插值函數(shù)，例如零次插值函數(shù)、線性插值函數(shù)和二次插值函數(shù)等，其相應(yīng)的邊界單元稱為常數(shù)單元、線性單元和二次單元。

2.1常數(shù)單元

對(duì)于常數(shù)單元，如圖2所示，節(jié)點(diǎn)取在邊界單元的中點(diǎn)，邊界上的函數(shù)值u和函數(shù)的法向?qū)?shù)值q在邊界單元上為常量，并等于節(jié)點(diǎn)處的值。

圖2　常數(shù)單元節(jié)點(diǎn)示意圖Fig.2 Nodes of constantelement

在每個(gè)單元上對(duì)邊界積分方程進(jìn)行離散，可得

式中：uj和qj分別為邊界單元節(jié)點(diǎn)處的函數(shù)值和函數(shù)的法向?qū)?shù)值；hij和gij分別為系數(shù)矩陣H 和G中的元素，其表達(dá)式為：當(dāng)i=j時(shí)，

式中：li為邊界單元的長(zhǎng)度；dij為節(jié)點(diǎn)Pi到邊界單元Li的垂直距離；rij為節(jié)點(diǎn)到積分單元上點(diǎn)的距離。從節(jié)點(diǎn)P1開(kāi)始，依次對(duì)所有節(jié)點(diǎn)的hij和gij進(jìn)行求解，利用Matlab對(duì)公式進(jìn)行編程并計(jì)算出結(jié)果。

2.2線性單元

對(duì)于線性單元，如圖3所示，節(jié)點(diǎn)取在邊界單元的兩端，邊界上的函數(shù)值u和函數(shù)的法向?qū)?shù)值q在邊界單元上服從線性分布，其離散化的邊界積分方程為

式中，uj，uj+1，qj，qj+1分別為邊界單元上2個(gè)端點(diǎn)處的函數(shù)值和函數(shù)的法向?qū)?shù)值。

圖3　線性單元節(jié)點(diǎn)示意圖Fig.3Nodesoflinearelement

每個(gè)線性單元可用無(wú)因次坐標(biāo)ξ表示為

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

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

式中：(xp，yp)為節(jié)點(diǎn)Pi的空間坐標(biāo)；為雅克比行列式，對(duì)于二維線性單元，有J=lj。從節(jié)點(diǎn)P1開(kāi)始，依次對(duì)各個(gè)節(jié)點(diǎn)的hij和gij進(jìn)行求解，求解過(guò)程與常數(shù)單元相同。

當(dāng)i=j時(shí)，邊界積分方程的形函數(shù)C由下式確定：

由于系數(shù)C只與節(jié)點(diǎn)Pi處邊界的幾何形狀有關(guān)，故兩線性單元間的夾角為非定值，如圖4所示。又因?qū)蔷€元素Hii的計(jì)算包含有系數(shù)C和Cauchy型積分的計(jì)算，因此，采用一種簡(jiǎn)便的均勻場(chǎng)的辦法便可間接計(jì)算C值。

圖4　線性單元節(jié)點(diǎn)處夾角示意圖Fig.4Anglearoundnodeoflinearelement

當(dāng)一個(gè)均勻位勢(shì)場(chǎng)作用于整個(gè)邊界時(shí)，函數(shù)的法向?qū)?shù)值必定為零。因此，系數(shù)矩陣H的任意一行中所有元素的總和均為零，即

對(duì)于無(wú)限域問(wèn)題，可設(shè)想所有內(nèi)部邊界被一個(gè)無(wú)限大的圓LR所包圍，則除了出現(xiàn)在邊界積分方程中的所有內(nèi)部邊界的積分外，還要加上在此大圓上的積分。

令R→∞，因

對(duì)邊界積分方程離散化后，可得

而Gii可推導(dǎo)出其積分的解析表達(dá)式為

2.3二次單元

對(duì)于二次單元，如圖5所示，節(jié)點(diǎn)分別取在邊界單元的兩端和中點(diǎn)，邊界上的u和q在邊界單元上服從二次分布。

圖5　二次單元節(jié)點(diǎn)示意圖Fig.5 Nodes of quadratic element

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

式中：uj，uj+1，uj+2分別為邊界單元上3個(gè)節(jié)點(diǎn)處的函數(shù)值；qj，qj+1，qj+2為其節(jié)點(diǎn)處函數(shù)的法向?qū)?shù)值。

每個(gè)二次單元可用無(wú)因次坐標(biāo)ξ表示為

式中，(x1，y1)，(x2，y2)和(x3，y3)分別表示二次單元上3個(gè)節(jié)點(diǎn)的坐標(biāo)值。

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

當(dāng)i=j時(shí)，與線性單元相同，Hii可由式（17）和式（19）求得。而Gii由于公式過(guò)于復(fù)雜，無(wú)法推導(dǎo)出解析解表達(dá)式，而且在單元端點(diǎn)處的積分存在奇異性，所以只能采取積分上的近似解法，即取一趨近于0的極小值ε。在對(duì)g1ii和gi3i積分時(shí)，積分區(qū)間取為[1-ε，1+ε]；對(duì)gi2i積分時(shí)，積分區(qū)間取為[0，0.5-ε]和[0.5+ε，1]，按上述區(qū)間對(duì)gii的積分項(xiàng)進(jìn)行積分，即可得到近似解。

值得說(shuō)明的是，由于高斯積分法所取的高斯積分點(diǎn)一般在積分區(qū)間內(nèi)，所以在積分過(guò)程中不需要求出帶奇異性的端點(diǎn)值，在不改變積分區(qū)間的前提下可直接避開(kāi)對(duì)積分奇異性的處理，并且多點(diǎn)高斯積分在計(jì)算精度方面也能滿足一定的要求。所以，高斯積分法在處理二次單元的問(wèn)題上存在著天然優(yōu)勢(shì)，該問(wèn)題在后文中將予以詳細(xì)討論。

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

3.1數(shù)值積分方法與計(jì)算精度

用前文介紹的邊界單元離散形式和數(shù)值積分方法對(duì)二維圓柱無(wú)界繞流問(wèn)題進(jìn)行求解。采用不同的節(jié)點(diǎn)數(shù)和數(shù)值積分方法對(duì)每種單元進(jìn)行計(jì)算，邊界單元離散的節(jié)點(diǎn)數(shù)依次從8取至200。之后，對(duì)每個(gè)節(jié)點(diǎn)處的計(jì)算值與解析值進(jìn)行比較，求出相對(duì)誤差并統(tǒng)計(jì)每次計(jì)算的時(shí)間，其結(jié)果如圖6～圖8所示。

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

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

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

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

對(duì)于線性單元和二次單元，一部分?jǐn)?shù)值積分方法已不適用于邊界積分方程的求解，有的計(jì)算結(jié)果會(huì)出現(xiàn)較大的波動(dòng)和偏差，有的甚至?xí)l(fā)散。尤其是在二次單元時(shí)，例如，逐次半分區(qū)間梯形公式法等已經(jīng)無(wú)法保證計(jì)算精度，平均誤差達(dá)到了5%以上，且計(jì)算時(shí)間也呈指數(shù)增長(zhǎng)，已超出邊界元法可接受的范圍。

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

3.2邊界單元離散形式與計(jì)算精度

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

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

由計(jì)算結(jié)果可知，在3種邊界單元離散形式中，常數(shù)單元的計(jì)算誤差一般都要小于其余2種，且隨著節(jié)點(diǎn)數(shù)的增多，常數(shù)單元的精度優(yōu)勢(shì)也更加明顯。當(dāng)劃分的節(jié)點(diǎn)數(shù)達(dá)到200時(shí)，常數(shù)單元的誤差基本上要比另2種單元的小一個(gè)數(shù)量級(jí)。此外，常數(shù)單元在計(jì)算時(shí)間上也遠(yuǎn)遠(yuǎn)優(yōu)于其他單元，所耗時(shí)間通常只是其他單元的5%～10%。

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

因此，本文認(rèn)為在網(wǎng)格單元?jiǎng)澐肿銐蚨嗟那闆r下，常數(shù)單元可以取代線性單元和二次單元對(duì)邊界進(jìn)行離散，并且計(jì)算誤差和計(jì)算時(shí)間均可以滿足邊界元法在工程實(shí)際中的應(yīng)用。

4　結(jié)語(yǔ)

邊界元法由于采用了解析基本解，在計(jì)算精度上有限元法等數(shù)值方法難以與其相比。但許多學(xué)者由于對(duì)邊界元法的誤差來(lái)源了解不深，往往習(xí)慣性地采用高次單元，單純地認(rèn)為離散單元階次越高，越能符合邊界條件的變化，但這樣反而不能把邊界元法在計(jì)算精度和速度上的優(yōu)勢(shì)充分發(fā)揮出來(lái)。

本文首先通過(guò)算例證明了邊界元法的主要誤差是邊界離散誤差和數(shù)值積分誤差，可見(jiàn)，正確地選擇邊界單元離散形式和數(shù)值積分方法是保證邊界元法計(jì)算精度的第一要?jiǎng)?wù)。隨后，驗(yàn)證了不同數(shù)值積分方法在求解邊界積分方程中的適用性，通過(guò)與常用的高斯積分法進(jìn)行對(duì)比，最終選出了一種精度較高、速度較快的積分方法——龍貝格積分法作為后續(xù)計(jì)算采用的方法。最后，通過(guò)對(duì)3種邊界離散形式計(jì)算精度的分析，認(rèn)為在網(wǎng)格單元?jiǎng)澐肿銐蚨嗟那闆r下，常數(shù)單元可以取代線性單元和二次單元作為一般二維邊界元問(wèn)題的邊界離散形式。從理論上說(shuō)，該結(jié)論也可以推廣到三維問(wèn)題中，即常數(shù)單元在三維邊界元問(wèn)題中也能獲得最佳的精度，但該推論的正確性還有待后續(xù)的研究驗(yàn)證。

［1］陳材侃.計(jì)算流體力學(xué)［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，F(xiàn)ERNANDEZ 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］周煥林，牛忠榮，程長(zhǎng)征，等.各向異性位勢(shì)問(wèn)題邊界元法中幾乎奇異積分的解析算法［J］.計(jì)算力學(xué)學(xué)報(bào)，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］馬榮，劉繼朝，石建省.邊界元法中角點(diǎn)問(wèn)題的一種計(jì)算方法［J］.中國(guó)科學(xué)院研究生院學(xué)報(bào)，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，F(xiàn)INK K D.Numericalmethods using Matlab［M］.NJ：Prentice Hall，2004.

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

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網(wǎng)絡(luò)出版時(shí)間：2015-7-28 17:25:18

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

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