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

沈王剛，鄭堯坤，林志良

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

0 引 言

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

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

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

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

本文結(jié)合多極子法與面元法，對三維潛體的興波特性進(jìn)行時(shí)域研究，將所得結(jié)果與傳統(tǒng)方法結(jié)果進(jìn)行比較，證明多極子面元法高精度與高效率的特性。

1 潛體興波問題描述

考慮如圖1所示的計(jì)算流場，O-xyz為三維笛卡爾直角坐標(biāo)系，xy平面與自由面平行，潛體沿x軸方向以速度U運(yùn)動，z軸垂直向上，整個流場邊界由潛體物面SB和自由面SF組成。潛體的特征長度為l，特征寬度為d，潛體中心距自由表面深度為h。對于時(shí)域興波問題，自由面條件中對流項(xiàng)起到主要作用，足夠遠(yuǎn)處將自由水面截?cái)?，截?cái)噙吔缟系姆瓷鋵萘鹘獾挠绊懖⒉幻舾衃12]。因此本文中采用截?cái)嘤邢拮杂杀砻娲鏌o限自由表面，其長為L，寬為D。假設(shè)流體無壓、無粘并且無旋，流場中存在一速度勢，滿足拉普拉斯方程：

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

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

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

2 傳統(tǒng)面元法求解

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

誘導(dǎo)速度

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

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

3 多極子面元法

快速多極子法與傳統(tǒng)面元法的結(jié)合過程和其與傳統(tǒng)邊界元法結(jié)合過程相同，詳細(xì)內(nèi)容可參考文獻(xiàn)[7 – 9]。本節(jié)將簡要介紹快速多極子方法的主要思想、基本公式和與面元法結(jié)合后的計(jì)算步驟。

3.1 多極子法主要思想

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

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

圖3 多極子展開關(guān)鍵點(diǎn)Fig. 3 Related points in multipole expansion

3.2 多極子法基本公式

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

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

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

3.3 多極子面元法計(jì)算步驟

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

步驟1與傳統(tǒng)面元法一樣對邊界進(jìn)行單元離散，布置面源，選取配置點(diǎn)。

步驟2劃分樹狀結(jié)構(gòu)。以二維問題為例，利用逐層嵌套的正方形來劃分空間樹狀結(jié)構(gòu)。首先利用一個足夠大的正方形包圍整個邊界S，并定義其層級數(shù)為0級?？疾斓贚層級的正方形單元C，將其均分為4個小正方形單元，層級為L+1，正方形單元C是它們的父單元。若某離散邊界單元的中點(diǎn)位于小正方形內(nèi)，則認(rèn)為其屬于該小正方形單元。對于4個小正方形單元，分成3種情況。若該小正方形單元內(nèi)離散邊界單元數(shù)為0，則將它舍棄，不計(jì)入樹狀結(jié)構(gòu)內(nèi)；若該小正方形單元內(nèi)離散邊界單元數(shù)大于0，且不大于事先給定的常數(shù)p時(shí)，將其定義為葉單元，不再往下劃分；若該小正方形單元內(nèi)離散邊界單元數(shù)大于常數(shù)p時(shí)，則繼續(xù)將其往下劃分。如此嵌套割劃，直到最后一層均為葉單元，或達(dá)到規(guī)定最大層數(shù)則停止。對于正方形單元之間的空間關(guān)系，定義如下（見圖4）：若2個正方形單元共用一個角點(diǎn)，則定義它們?yōu)橄噜弳卧?；?個正方形單元不相鄰，但其父單元相鄰，則定義它們?yōu)橄喔魡卧蝗?個正方形單元的父單元也不相鄰，則定義它們?yōu)檫h(yuǎn)距單元。

圖4 正方形單元相互關(guān)系Fig. 4 Relations between square cells

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

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

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

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

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

4 計(jì)算結(jié)果與分析

4.1 近水面回轉(zhuǎn)橢球體勻速直航興波問題

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

計(jì)算模型為回轉(zhuǎn)橢球體，長軸長l=1.0 m、短軸長d=0.2 m的回轉(zhuǎn)橢球在自由水面下（球心距靜水面距離為h=0.16 m）從靜止?fàn)顟B(tài)突然啟動，達(dá)到常速U=1.566 m/s并沿軸正向勻速直航，F(xiàn)roude數(shù)為浸深比為橢球表面網(wǎng)格劃分如圖6所示，其中圖6（a）橢球面元數(shù)為528，圖6（b）橢球面元數(shù)為1 304，網(wǎng)格劃分時(shí)均為兩端密，中間疏；圖6（c）橢球面元數(shù)為2 640，網(wǎng)格劃分時(shí)兩端與中間均加密。

圖6 橢球網(wǎng)格劃分Fig. 6 The mesh generation of elipsoid

自由面的網(wǎng)格密度和時(shí)間步長會直接影響到計(jì)算量、波面的光滑程度以及興波水動力學(xué)參數(shù)的準(zhǔn)確性，網(wǎng)格和時(shí)間步長的選取參照的是Dommermuth等人通過大量數(shù)值試驗(yàn)得到的經(jīng)驗(yàn)公式[19]：

4.1.1 自由面靜網(wǎng)格

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

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

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

圖7 t=8 s時(shí)波面圖Fig. 7 The wave surface at t=8 s

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

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

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

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

4.1.2 自由面動網(wǎng)格與網(wǎng)格密度

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

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

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

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

4.1.3 計(jì)算效率對比

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

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

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

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

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

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

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

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

5 結(jié) 語

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

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