多極子面元法近水面橢球體興波時域研究

沈王剛，鄭堯坤，林志良

(上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室，上海 200240)

0 引 言

研究潛體近水面航行時的興波現象及其力學特性，一直是一項具有重要工程意義的課題。尤其近幾年各類附加水翼新概念船型的大量涌現，使得潛體興波研究顯得更為重要。由于計算效率的限制，以往對潛體興波研究主要采用定常方法。如Farell改進了Havelock源格林函數，推導得到潛航回轉橢球體在Neumann-Kelvin問題下的半解析理論解[1]；Doctors和Beck研究了潛體在Neumann-Kelvin問題下，采用Havelock源計算的收斂性[2]。隨著時代的發展，時域方法以其可以計算非定常興波，得到運動物體實時興波信息，便于仿真模擬和后處理等優點，吸引了越來越多學者的注意。

自從Hess和Smith[3]第一次應用面元法求解三維物體繞流問題以來，面元法在船舶與海洋工程領域得到了廣泛的應用。其中Rankine源形式簡單，具有較高靈活性，且易于推廣到非線性計算，已經有不少學者嘗試使用Rankine源處理興波問題[4–5]。

然而，傳統的面元法在求解問題時，需要計算每個面元與配置點之間的相互作用，生成一個N階方陣（N為離散面元個數）。一方面儲存這個N階稠密矩陣需要大量的內存空間，另一方面求解該矩陣也將耗去大量時間。例如采用高斯消去法直接求解，時間復雜度將達到O（N3），即使使用迭代法求解，時間復雜度也有O（N2）。計算效率上的缺陷將面元法的應用限制在了小規模流場，更遑論時域計算時需要進行多次時間步進迭代，計算規模將被進一步壓縮。為克服傳統面元法這一局限性，引入快速多極子法，可同時減小存儲量與時間復雜度。

快速多極子法最早由Rokhlin[6]于1985年提出，用來加速二維勢流問題的求解。此后Liu[7]，Nishimura[8]及Yoshida[9]等將快速多極子法和邊界元法結合，并對其進行深入研究。多極子法與邊界元法的結合已被證實可以大幅度提高其計算效率，適用于求解大規模勢流問題[10–11]。

本文結合多極子法與面元法，對三維潛體的興波特性進行時域研究，將所得結果與傳統方法結果進行比較，證明多極子面元法高精度與高效率的特性。

1 潛體興波問題描述

考慮如圖1所示的計算流場，O-xyz為三維笛卡爾直角坐標系，xy平面與自由面平行，潛體沿x軸方向以速度U運動，z軸垂直向上，整個流場邊界由潛體物面SB和自由面SF組成。潛體的特征長度為l，特征寬度為d，潛體中心距自由表面深度為h。對于時域興波問題，自由面條件中對流項起到主要作用，足夠遠處將自由水面截斷，截斷邊界上的反射對勢流解的影響并不敏感[12]。因此本文中采用截斷有限自由表面代替無限自由表面，其長為L，寬為D。假設流體無壓、無粘并且無旋，流場中存在一速度勢，滿足拉普拉斯方程：

在潛體表面SB上滿足不可穿透條件，即

在本文時域興波問題中，假設潛體由靜止快速啟動至勻速運動狀態，且在自由面SF上初始條件（t=0）為：

圖1 計算流場示意圖Fig. 1 The schematic of computational flow field

2 傳統面元法求解

首先將整個流場的邊界劃分成N個離散面元，其中前NB個為三角形常數面元，用于表征潛體表面；后NF個為四邊形常數面元，用于表征自由水面，則有N=NB+NF。第i個面元的的面積為Si，其上均勻分布強度為σi的面源。此時，可得到流場中任一點的誘導速度勢

誘導速度

此外，選取每一個常數面元的中點為配置點，在每一個配置點上滿足已知的邊界條件，例如在t=0初始時刻，潛體物面上的配置點應滿足式（2），而自由表面上的配置點應滿足式（5）中的據此可得到如下線性方程組，即

常數面元法的優點在于公式推導簡單，易于應用；但為了提高精度，需要增加離散面元數目，使得計算規模增大、計算效率降低。此外，不少學者嘗試使用高階面元法求解勢流問題，如Liu Y H[15]，Romate JE[16]等。盡管計算精度和效率有所提升，高階面元法仍無法有效克服傳統面元法求解大規模流場問題的缺陷。因此，本文采用快速多極子方法加速傳統常數面元法計算，求解大規模流場時域興波問題。

3 多極子面元法

快速多極子法與傳統面元法的結合過程和其與傳統邊界元法結合過程相同，詳細內容可參考文獻[7 – 9]。本節將簡要介紹快速多極子方法的主要思想、基本公式和與面元法結合后的計算步驟。

3.1 多極子法主要思想

以粒子在粒子群中受力為例，兩兩直接計算再求和策略的計算量級是O（N2），其中N表示粒子數量?？焖俣鄻O子法將所有粒子按其空間位置分成不同的集合（如圖2中的A和B），首先將粒子對外的作用力匯聚到集合內一點（如圖2中的①過程），其次計算集合與集合之間的相互作用力（如圖2中的②過程），最后將集合所受到的外界作用力分配到集合內每一個粒子（如圖2中的③過程）。對于相距很近的粒子，相互間作用力直接計算得到；而對于相距較遠的粒子，相互間的作用力則采用上述多極子方法展開得到，從而將計算量級減少至O（N·LogN)[7–9]。其中區分粒子間距關系，可采用樹狀結構劃分策略[17]加以實現。在傳統面元法中，主要過程在于計算邊界離散面元間的相互作用，因此可以采用快速多極子方法來提高計算效率。

圖2 多極子計算示意圖Fig. 2 The schematic of multipole calculation

圖3 多極子展開關鍵點Fig. 3 Related points in multipole expansion

3.2 多極子法基本公式

關于多極子法公式具體推導請參見文獻[7-9]，這里以面元法中三維Rankine源誘導速度勢核函數多極子展開為例，其展開形式如下：

得到核函數的展開形式后，即可得到積分方程的展開形式如下：

對于核函數的局部展開系數則定義如下：

3.3 多極子面元法計算步驟

在此只簡要介紹多極子面元法的計算步驟，具體請參見文獻[18]。

步驟1與傳統面元法一樣對邊界進行單元離散，布置面源，選取配置點。

步驟2劃分樹狀結構。以二維問題為例，利用逐層嵌套的正方形來劃分空間樹狀結構。首先利用一個足夠大的正方形包圍整個邊界S，并定義其層級數為0級?？疾斓贚層級的正方形單元C，將其均分為4個小正方形單元，層級為L+1，正方形單元C是它們的父單元。若某離散邊界單元的中點位于小正方形內，則認為其屬于該小正方形單元。對于4個小正方形單元，分成3種情況。若該小正方形單元內離散邊界單元數為0，則將它舍棄，不計入樹狀結構內；若該小正方形單元內離散邊界單元數大于0，且不大于事先給定的常數p時，將其定義為葉單元，不再往下劃分；若該小正方形單元內離散邊界單元數大于常數p時，則繼續將其往下劃分。如此嵌套割劃，直到最后一層均為葉單元，或達到規定最大層數則停止。對于正方形單元之間的空間關系，定義如下（見圖4）：若2個正方形單元共用一個角點，則定義它們為相鄰單元；若2個正方形單元不相鄰，但其父單元相鄰，則定義它們為相隔單元；若2個正方形單元的父單元也不相鄰，則定義它們為遠距單元。

圖4 正方形單元相互關系Fig. 4 Relations between square cells

步驟3向上傳遞。首先通過式（24）或式（25）計算各個葉單元的多極子動量矩，其次利用式（26）將葉單元的多極子動量矩轉移到其父單元中。以此類推，多極子動量矩不斷上聚，直至第2層。

步驟4向下傳遞。對于相隔單元，通過式（28）計算其作用；對于遠距單元，則利用式（29）將父單元通過得到的局部展開系數繼承到子單元，以此來計算其作用。如此類推，直到得到每個葉單元的局部展開系數。向上傳遞和向下傳遞過程參見圖5。

圖5 向上和向下傳遞過程Fig. 5 Upward and downward passes

步驟5近場遠場作用疊加。對于某葉單元C中的配置點，利用傳統面元法的直接作用計算該葉單元及相鄰單元中其他離散面元對其的作用，利用式（27）將局部展開系數從單元中心轉移到配置點以計算相隔或遠距單元中離散面元對其的作用。將近場和遠場的作用疊加，就得到了整個邊界所有離散面元對配置點的作用。

步驟6利用GMRES求解式（10），若迭代誤差小于設定值，則終止運算，否則返回第3步進行下一步迭代。本文中迭代誤差設定值均取為10–4。

4 計算結果與分析

4.1 近水面回轉橢球體勻速直航興波問題

本節將以近水面回轉橢球體勻速直航產生的興波波形和興波升力、阻力為切入點，闡述如何選取合適的自由面計算域和自由面及物面網格密度，驗證使用多極子面元法計算潛體興波問題的可行性，并將多極子面元法與傳統面元法進行計算效率對比。

計算模型為回轉橢球體，長軸長l=1.0 m、短軸長d=0.2 m的回轉橢球在自由水面下（球心距靜水面距離為h=0.16 m）從靜止狀態突然啟動，達到常速U=1.566 m/s并沿軸正向勻速直航，Froude數為浸深比為橢球表面網格劃分如圖6所示，其中圖6（a）橢球面元數為528，圖6（b）橢球面元數為1 304，網格劃分時均為兩端密，中間疏；圖6（c）橢球面元數為2 640，網格劃分時兩端與中間均加密。

圖6 橢球網格劃分Fig. 6 The mesh generation of elipsoid

自由面的網格密度和時間步長會直接影響到計算量、波面的光滑程度以及興波水動力學參數的準確性，網格和時間步長的選取參照的是Dommermuth等人通過大量數值試驗得到的經驗公式[19]：

4.1.1 自由面靜網格

自由面為一足夠大的長L=20 m，寬D=6 m的矩形平面，對稱于坐標系原點。橢球球心初始坐標為x=–7.5，y=0，沿軸正向勻速直航。

橢球網格劃分如圖6（b）所示。自由面網格密度有3套，其中粗網格選取為100×30，網格間距為中網格選取為150×45，網格間距為細網格選取為200×60，網格間距時間步長間隔為0.04 s，步進次數200次，總模擬時間8 s。

圖7為t=8 s時波面圖，經過對比可發現，粗網格與細網格所描繪出的波面波形圖及等高線圖基本一致，但粗網格仍有一些波面鋸齒，而細網格的波面則光滑許多。對比圖8亦可發現，3種網格描繪出的自由面中剖線二維波形圖在第一個波谷處有比較明顯的差別，除此之外中網格和細網格的結果完全貼合，而粗網格則與其他2種網格差別較大。

圖7 t=8 s時波面圖Fig. 7 The wave surface at t=8 s

其次考察興波導致的橢球體升力、阻力。圖9（a）為升力系數隨時間變化情況，圖9（b）為興波阻力系數隨時間變化情況。圖中圓點直線為Farell解[1]，3種網格所計算得到的力學系數都能收斂到Farell解附近。然而粗網格結果，特別是其升力系數的波動很大，毛刺明顯；中網格和細網格的毛刺較粗網格有所減少，但仍然比較明顯。此外，粗網格計算得到的升力系數收斂值和其他2種網格比較吻合，但興波阻力收斂值有比較明顯的差別，顯示粗網格密度并不夠，計算精度低；同時隨著自由面網格不斷加密，其計算結果也將逐漸收斂。

圖8 t=8 s時自由面二維波形圖（y=0）Fig. 8 The 2D waveform of free surface at t=8 s (y=0)

圖9 興波水動力參數隨時間變化情況Fig. 9 The change of wave making hydrodynamic parameters over time

使用靜網格固然能描述潛體的整個尾后興波波面，但其取的自由面計算域過大，網格數量過多，時間效率低下；此外有大片自由面區域在長時間內興波影響很小，是一種計算資源浪費。為了能提高計算效率，減少計算資源浪費，使用能隨潛體一起前進的動網格來進行興波時域計算。

4.1.2 自由面動網格與網格密度

自由面為一個長L=8 m，寬D=5 m的矩形平面，初始時刻對稱于坐標系原點。橢球球心初始坐標為x=–7.5，y=0，沿軸正向勻速直航。

橢球網格劃分有3套（見圖6）。自由面網格密度選取為80×50，網格間距時間步長間隔為0.04 s，步進次數200次，總模擬時間8 s。在每一次步進后，若橢球中心與網格中心的距離大于，則將網格向軸正向移動一個網格間距，相當于消去最后一排網格，并在第一排網格前方再加上一排網格，其初始條件如式（5）所示。

從圖10中可以看到，雖然動網格的網格密度與靜網格中細網格密度相同，但其計算結果中的毛刺現象卻有大幅度改善，究其原因可能是靜網格所取的自由面計算域過大，潛體相對網格的位置變化明顯，在網格密度不夠絕對密的情況下毛刺就會比較多，從圖9中也可看出，當自由面靜網格加密后，毛刺現象有所改善。此外，1304面元的橢球網格2和2640面元的橢球網格3所計算得到的升力曲線完全一致，興波阻力曲線有略微差別，但只在1.8%左右；而528面元的橢球網格1則與其他2種網格計算結果差別較大，因此可認為隨著橢球網格數量的增加，計算結果將逐漸收斂，同時1304面元的橢球網格2精度已經足夠。

此外，考察自由面橫向網格密度變化率q的影響，q=1代表橫向網格均勻分布，q＞1代表中縱線附近網格密，橫向兩端網格疏。表1顯示了不同橫向網格密度變化率下計算所得結果，其中2種網格布置的橫向網格數相同，均為50個。由表中可看出橫向網格變化率對計算結果的影響很小，幾乎可以忽略不計。

4.1.3 計算效率對比

圖10 不同物面網格下興波水動力參數隨時間變化情況Fig. 10 The change of wave making hydrodynamic parameters over time in different object surface mesh

表1 不同橫向網格密度變化率下興波水動力學系數Tab. 1 The wave making hydrodynamic parameters in different changing rate of transverse mesh density

表2 多極子面元法和傳統面元法單步計算效率比較Tab. 2 Comparison of single-step computational efficiency between FMBEM&BEM

多極子法對傳統面元法的改進，首先體現在內存使用量上，傳統面元法需要占用O（N2）的內存，普通臺式計算機無法駕馭面元數量過萬的問題，但多極子面元法只需要占用O（N）的內存，使得普通臺式計算機也可以計算面元數量達到數百萬甚至數千萬的大規模問題。表2和圖11是分別采用多極子面元法和傳統面元法進行時域單步計算的效率對比，計算環境為1.9 GB內存和2.50 GHz*2處理器。從圖表中可以看出，當面元數較小時，使用多極子面元法和傳統面元法所需要的計算量相近，甚至當面元數只有數百個時，傳統面元法還要略優于多極子面元法；然而隨著面元數的增長，多極子面元法所需的計算量近似線性增長，而傳統面元法所需的計算量卻是以三次方的速度增長，當面元數上萬時，使用傳統面元法求解不僅在內存上受到限制，在計算時間上也顯得不現實。此外，由于時域問題較定常問題要多一個時間維度，因此利用多極子面元法來計算時域問題的可行性和高效性不言而喻，具有巨大的應用潛力。

圖11 多極子面元法和傳統面元法單步計算效率比較Fig. 11 Comparison of single-step computational efficiency between FMBEM&BEM

4.2 不同航速及潛深下橢球定常興波

本小節將以近水面橢球勻速直航達到穩定狀態后所產生的升力和興波阻力為切入點，說明不同航速和下潛深度對其定常興波的影響。計算模型與上節一樣為回轉橢球體，長軸長l=1.0 m、短軸長d=0.2 m，改變航速和潛深以研究其定常興波力學系數的變化情況。

興波水動力系數如圖12所示，其中圖12（a）是升力系數隨Froude數變化情況，圖12（b）是興波阻力系數隨Froude數變化情況。與Doctor解[2]相比，本文方法所得結果在時吻合較好；在時，高Froude數情況下的升力系數和阻力系數，以及Fr=0.45～0.5段的阻力系數會有較為明顯的差別，其余部分吻合較好。由圖中可以看出，浸深比下升力系數和阻力系數的值總小于浸深比下的值，升力系數在左右達到最大，在Fr=0.55～0.6后變為負升力；阻力系數在左右達到最大，在低Froude數時接近零。此外，下升力系數和阻力系數最大值所對應的Froude數較略微后移。

5 結 語

從計算結果來看，本文所開發的快速多極子面元法可以較準確地預報近水面潛體時域興波特性，定常興波計算結果與已有研究結果吻合良好，證明了其合理性與可靠性；且相較于傳統面元法，有效地克服了計算規模和效率的局限性，引進了運動網格提高計算精度和效率。本文方法將進一步改進和發展，以計算完全非線性興波問題、船體興波問題等，具有較大的應用前景。

圖12 興波水動力參數Fig. 12 The wave making hydrodynamic parameters