一種處理高階面元法中奇異積分的方法*

，2

1.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 上海 200030 2.上海交通大學(xué) 海洋工程國家重點實驗室 上海 200030

自從Hess和Smith[1]用一階面元法求解三維勢流問題以來，面元法作為一種有效的數(shù)值方法已被廣泛地應(yīng)用于求解空氣動力學(xué)和水動力學(xué)中的三維勢流問題。最初的面元法中存在幾個缺點，如物面用平面四邊形面元離散，相鄰面元間存在間隙；在離散大曲率的曲面時必須要有大量的面元，為此需要大量的計算時間；曲面上的源(匯)強度在面元邊界上不連續(xù)。

為了克服低階面元法的不足，人們引入線性或二次的函數(shù)表達面元，并在面元上分布一階或二階多項式源項。近年來，由于三次樣條、B樣條以及非均勻有理B樣條(non-uniform rational B-splines, NURBS)在工業(yè)產(chǎn)品外形設(shè)計中的廣泛應(yīng)用，人們已經(jīng)把他們應(yīng)用到高階面元法中以表達物面幾何和速度勢[2-5]。對于高階面元法中的Rankine源項，Maniar和Danmeier以及Lee等[2-4]均運用三角形細分和級數(shù)展開的方式來近似求解。該方法在計算中較為復(fù)雜。趙成璧[5]的方法是離散第二類Fredholm積分方程，再直接高斯積分。該方法在場點離面元較遠時可獲得較好的結(jié)果，而當(dāng)場點離面元較近時，數(shù)值積分的精度將受到很大的影響。

在處理基于B樣條的高階面元法中的奇異積分時，根據(jù)場點到面元的距離把曲面積分分為遠場積分和近場積分。遠場積分中，直接應(yīng)用高斯-勒讓德公式；近場積分中，采用Telles等[6-7]介紹的方法，先對積分函數(shù)進行三次多項式變換，然后應(yīng)用高斯-勒讓德公式。

1 數(shù)學(xué)公式

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

考慮物體在無限靜止流體中以勻速U向前運動，定義一個固定在物體上的右手直角坐標(biāo)系o-xyz(圖1)。假設(shè)流體是不可壓縮理想流體，流動無旋，則存在定常擾動勢φ滿足下列邊值問題：

(1)

圖1 坐標(biāo)

(2)

式中：φ(q)——物面上源點q(ξ,η,ζ)處的源(匯)強度；

S——物面。

把式(2)代入式(1) 中得到確定源(匯)強度的第二類Fredholm積分方程。

(3)

1.2 有關(guān)物理量的表達

有

(4)

式中：u0、v0——代表NURBS曲面中的兩個參數(shù)方向；

dij——物面控制點網(wǎng)格；

m、n——u0、v0方向上物面控制點的數(shù)目；

Ni、g(u0)、Nj、h(v0)——定義在u0、v0上的B樣條基函數(shù)；

wij——權(quán)因子；

g、h——u0、v0上B樣條基函數(shù)的階數(shù)。

(5)

式中：u、v——B樣條曲面中的兩個參數(shù)方向；

m′、n′——u、v方向上控制點的數(shù)目；

Ni,g(u)、Nj,h(v)——B樣條基函數(shù)；

g、h——u,v上B樣條基函數(shù)的階數(shù)。

把式(4)、(5)代入式(3)中，可以得到用B樣條參數(shù)表達的積分方程：

(6)

(7)

2 數(shù)值方法

圖2 場點與面元的關(guān)系

在結(jié)合高斯-勒讓德公式積分時，三次多項式變換在不對面元進行任何細分的情況下，可以使高斯點自動向奇異點靠近，這個特性將有助于提高準奇異積分的精確性。因此，對于近場積分，先把式(6)中的積分函數(shù)進行三次多項式變換，然后再結(jié)合高斯-勒讓德公式積分。

3 數(shù)值計算與結(jié)果比較

3.1 平面積分

為了驗證三次多項式變換法的效果，考慮以下準奇異的二重積分

(8)

其中a,b為常數(shù)。在此使用8×8的高斯-勒讓德公式，分別采用直接積分和經(jīng)三次多項式變換后的積分。表1給出了計算結(jié)果和解析解的比較。

表1 兩種奇異積分計算結(jié)果對比

從表中可以看出，當(dāng)場點在(1.40，0)處時，由于它離面元相對較遠，直接積分法和多項式變換法的計算誤差都較小，但隨著場點離面元距離的減少，當(dāng)減少到(1.05, 0)處時直接積分法的誤差迅速增大到50.80%，而多項式變換法的誤差僅為8.66%，說明后者有著比前者更好的計算精度。

3.2 圓球

考慮一單位半徑的球體以單位速度沿x軸正向運動，分別采用直接積分法和文中所介紹的方法對定常擾動勢進行了計算。

表2給出了計算結(jié)果。從表2中可以看出，兩種方法的計算結(jié)果均和解析解比較吻合，但本文方法的計算誤差更小，尤其在靠近x=-1.0處。其中一個重要原因是在球的NURBS表達中，靠近x軸兩端的面元較小，場點離面元相對較近，積分的奇異性表現(xiàn)得非常明顯，直接積分將忽略這種特性，導(dǎo)致計算誤差，而本文的計算方法通過對積分函數(shù)進行三次多項式變換，提高了計算精度。

表2 球面速度計算結(jié)果的對比

3.3 橢球

采用文中介紹的方法，對一系列以單位速度平移的橢球體在第二模態(tài)下的附加質(zhì)量系數(shù)K22=3λ22(4πρabc)進行了計算，其中λ22是附加質(zhì)量，ρ是流體密度，a、b和c分別是橢球體的三個半軸長；在此a∶b∶c的值分別取7∶1∶1, 5∶1∶1, 3∶1∶1, 2∶1∶1和1∶1∶1。根據(jù)三半軸的比值，布置不同的面元數(shù)。計算結(jié)果與Lamb給出的解析解和文獻[7]中的計算結(jié)果(直接高斯積分法)進行了比較，如圖3所示。可以看出，本文的結(jié)果比文獻[7]的計算值更接近于解析解。

圖3 橢球體附加質(zhì)量系數(shù) K22

4 結(jié)束語

本方法為解決高階面元法中的奇異積分問題提供了一種有效的手段，在船舶運動性能預(yù)報相關(guān)的m項計算和水動力系數(shù)計算中具有較大的應(yīng)用價值。

[1] Hess J L， Smith A M. Calculation of nonlifting potential flow about arbitrary three-dimensional smooth bodies[J]. Journal of Ship Research, 1964(7):22-44.

[2] Maniar H D. A three dimensional higher order panel method based on B-splines[D], MIT, Massachusetts, 1995.

[3] Danmeier D G. A Higher-Order Panel Method for Larger-Amplitude Simulation of Bodies in Waves[D]. MIT, Massachusetts, 1999.

[4] Lee C H，Newman J N, Solution of Radiation Problems with Exact Geometry, Proc. of 16th International Workshop on Water Wavers and Floating Bodies[C], Hiroshima, Japan, 2001：93-96.

[5] 趙成璧. 船舶曲面設(shè)計現(xiàn)代方法與軟件系統(tǒng)的研究[D].武漢： 武漢交通科技大學(xué), 1999.

[6] Telles J C F，Oliveira R F. Third degree polynomial transformation for boundary element integrals: Further improvement[J]. Engineering Analysis with Boundary Elements, 1994(13)：135-141.

[7] 王化明. 基于NURBS高階面元法的有航速船舶輻射問題數(shù)值計算研究[D].武漢：武漢理工大學(xué), 2005.