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

，2

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

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

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

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

1 數學公式

1.1 數學模型

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

(1)

圖1 坐標

(2)

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

S——物面。

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

(3)

1.2 有關物理量的表達

有

(4)

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

dij——物面控制點網格；

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

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

wij——權因子；

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

(5)

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

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

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

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

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

(6)

(7)

2 數值方法

圖2 場點與面元的關系

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

3 數值計算與結果比較

3.1 平面積分

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

(8)

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

表1 兩種奇異積分計算結果對比

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

3.2 圓球

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

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

表2 球面速度計算結果的對比

3.3 橢球

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

圖3 橢球體附加質量系數 K22

4 結束語

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

[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] 趙成璧. 船舶曲面設計現代方法與軟件系統的研究[D].武漢： 武漢交通科技大學, 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高階面元法的有航速船舶輻射問題數值計算研究[D].武漢：武漢理工大學, 2005.