陳帥，王超，周廣利，蔣彩霞

1. 中國船舶科學研究中心，江蘇 無錫 214082

2. 哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001

船體興波阻力快速預報是船型優(yōu)化設計方法研究的關鍵技術之一，它可以通過求解定常興波問題得到?，F在主流的船舶興波阻力預報的理論方法有基于勢流理論的Rankine源面元法（Dawson法）、Neumann-Michell （NM）理論[1-4]和考慮黏性的計算流體動力學(computational fluid dynamics,CFD)技術。CFD方法由于求解Navier-Stokes方程需要數小時或更長時間，其適用于評估設計，很難滿足船型優(yōu)化高效的要求。在上兩個世紀中，出現了一些經典的興波勢流理論，如Michell理論、Kelvin源（或Havelock源）、Neumann-Kelvin（NK）理論和新細長體理論等，這些理論為后續(xù)的興波問題研究奠定了基礎。

Dawson[5]在自由面和船體表面布置Rankine源1/r來求解定常興波問題，使得Dawson法在船舶興波阻力預報上取得了巨大的成功。Tarafder等[6]提出了一種改進的Dawson方法計算了Wigley和S60船型的阻力和波高分布，自由面采用貼體網格代替流線網格，其預報值與實驗值誤差較小。Rabaud等[7]基于Dawson法考慮了非線性的影響，預報了Wigley和KCS船型的阻力、縱傾、升沉和自由面波高，其預報值與實驗值誤差較小。

在國內，張寶吉等[8-10]利用Dawson法對船體興波阻力進行數值預報，并應用于某高速水面艦船進行優(yōu)化設計。程京普等[11]提出了一種改進的Dawson方法，采用靜水面貼體網格代替流線網格，對方尾船型的興波阻力進行預報并驗證。范井峰等[12]利用Dawson方法預報了小水線面雙體船在不同航速下的興波阻力、升沉和縱傾等，并與實驗結果進行了對比。何廣華等[13]建立了一種水下航行體興波尾跡和興波阻力的快速分析模型，其Rankine源項采用Newman解析法進行求解。

國內外更多針對Dawson方法的應用以及考慮非線性因素等，缺乏對求解過程的計算效率分析和研究，鑒于此，本文編寫了基于貼體網格的Dawson方法計算程序，引入加速大規(guī)模并行計算(accelerated massive Parallelism, AMP)技術，建立了船體興波阻力快速預報模型。本文開發(fā)的程序只需數十秒就可以計算一次興波，適用于船舶初步設計和船型優(yōu)化。

1 理論與數值離散模型

1.1 興波理論

本文采用C++語言自編程開發(fā)，基于Dawson法的興波阻力計算步驟[2]：1) 在船體表面上用面元法計算疊模擾動流；2) 在船體表面和自由面上利用物面邊界條件、自由邊界條件和輻射條件計算興波擾動流，把求解的疊模繞流和興波擾動流問題轉化為求解船體表面和自由面分布源點強度問題：

式中：U為航速；Φ為總速度勢，滿足拉普拉斯方程；φr為疊模擾動速度勢；φw為興波擾動速度勢；為疊模表面源強；σB和σF分別為船體表面和自由面源強；SB為船體表面；SF為自由面；rpqb、rpqw為場點p到船體表面上點qb和自由面上點qw的距離。

除了滿足拉普拉斯方程和遠前方無波的輻射條件，還在船體表面上滿足物面條件，在自由液面z=ζ(x,y)上滿足自由面邊界條件：

式中：下標x、y、z為該變量沿x、y、z方向的偏導數，Φx為總速度勢沿x方向的偏導數。

自由面邊界條件可推導[2]為

式中l(wèi)為該變量沿流線方向的導數。滿足線性化邊界條件的Bernoulli方程由船體周圍的速度表達，略去φw的二階項，船體表面無量綱化壓力系數為

濕表面積無量綱化的興波阻力系數為

式中：ΔSi為船體離散面元的面積，S為船體濕表面積，為面元單位法向量的x方向分量。水平面波高起伏為

1.2 數值離散

為獲得積分方程數值解，采用NURBS技術將船體表面和自由面離散成一系列小單元，該方法實現過程見文獻[14]。如圖1所示，考慮到船體具有對稱性，采用半自由面和半船體進行計算。根據式(1)～式(3)邊界條件，可得船體表面與自由面離散的線性方程組[14]：

式中：NB和NF分別為半船表面和半自由面的面元數，影響系數Aij、b(i)，j=1,2,···,NF，CBB、CFB、CBF、CFF和B(i)由文獻[5]提供計算方法。運用Gauss消去法可直接求解SB和SF上離散分布源的密度，進一步求解場點處興波引起的擾動速度?φw，最后計算出興波阻力系數Cws和波高起伏ζ(x,y)。

圖1 計算模型示意

1.3 貼體網格

圖2 自由面網格圖

如圖2，靜水面上流線網格是通過上游點的速度追蹤求得，形成了一條條流線；貼圖網格是使用水線和計算域邊界采用NURBS曲面[14]求得。本文Dawson法的自由面離散采用貼體網格以代替?zhèn)鹘y(tǒng)的疊模流線網格，進而式(4)中的Φll數值離散的過程中，二階導數的計算轉換到大地坐標系下進行。有

式中：(ξ，η)為物理坐標系下的點，(x，y)為大地坐標系下的點?？v向和橫向的導數都采用單邊上游有限差分算子：

式中J為雅可比矩陣。

Φll求出之后可以解出式(5)～式(7)中的影響系數。同時，在使用貼體網格求解的過程中，減少了使用疊模流場求解自由面離散網格這一步驟，進一步減少了興波問題數值求解的計算時間。

2 數值計算

2.1 航行體模型

如圖3和圖4所示，采用2種航行體，船型參數如表1，一種是直接由數學公式來描述的Wigley船型，被廣泛地用于研究常規(guī)船型的興波問題；另一種為具有代表性的集裝箱船（KRISO container ship，KCS）。采用本文開發(fā)的程序計算其興波阻力、波高分布和船側波高，并與相關文獻的結果進行比較分析，驗證該計算方法的可靠性。

圖3 Wigley船型

圖4 KCS船型

表1 船型參數

2.2 數值模型收斂性驗證

為了研究數值模型對興波阻力預報精度的影響，本文模擬了Wigley船型以Fr=0.3均速航行時的興波運動，主要以航行體網格、自由面區(qū)域大小和自由面網格劃分3方面關鍵參數進行收斂性分析。如圖5～圖10所示，在其他參數值不作特別說明時，取值如下：自由液面下游le=2L，自由液面上游Lf=0.5L，自由液面寬度Lw=1L，數值計算離散模型均采用非均勻網格，網格沿著各方向以系數1.1的增長率劃分，船體模型網格數量NX=50，NY=20，自由液面網格數量NW=30，Ne=60，Nf=15。

2.2.1 自由面區(qū)域收斂性驗證

為了對自由液面區(qū)域大小進行收斂性驗證，如圖5和圖6所示，分別對計算域上游lf、下游le和寬度lw進行收斂性分析。從圖6可以看出，當自由液面寬度lw∈[1.5L,2.2L]、上游lf∈[0.5L,1.0L]和下游le∈[2.5L,3.3L]時，興波阻力計算結果十分接近。為了獲取更高的預報精度，在下面的算例中，計算域均大小均采用le=2.5L、lf=0.5L和lw=1.5L。

圖5 自由面區(qū)域

圖6 興波阻力隨自由面區(qū)域大小的變化

2.2.2 航行體網格劃分收斂性驗證

同上，改變Wigley船型X和Y方向網格數量進行收斂性驗證，船艏艉0.2L處區(qū)域網格進行加密，如圖7和圖8所示，分別對船體垂向網格數NY和縱向網格數NX進行收斂性分析。從圖8可以看出，NY對興波阻力預報的結果影響較小且具有較好的收斂性，縱向網格數NX∈[50,68]時，興波阻力計算結果具有較小波動。為了獲取更高的預報精度，在下面的算例中，船體網格劃分采用NX=50和NY=20。

圖7 航行體網格

圖8 興波阻力隨船體網格數量的變化

2.2.3 自由面網格劃分收斂性驗證

同理，對自由液面網格劃分進行收斂性分析，如圖9和圖10所示，加密增長率Ge、Gf和Gw均為1.1。從圖10可以看出，當自由液面寬度網格數較大時興波阻力預報精度較高，而自由液面縱向網格數對計算精度波動較大，但在Nf∈[10,20]和Ne＞40時，興波阻力計算結果十分接近。綜上，自由液面網格劃分采用Nw=35、Ne=40和Nf=15具 有較高的計算精度。

圖9 自由面網格劃分示意

圖10 興波阻力隨自由面區(qū)域大小的變化

2.3 數值模型網格劃分策略

同理，對KCS船型進行了收斂性分析，在下面的算例中，船體網格、自由面區(qū)域大小和自由面網格劃分如表2。由于KCS船體曲面較為復雜，為了提高計算精度，本文采用NURBS技術[14]對艏艉進行加密，圖11為船體和自由液面劃分3 830個面元的計算結果，從圖11中可以看出船艏艉沿縱向的網格數更為密集。

表2 計算模型數值離散網格值

圖11 KCS船體及自由面網格示意

3 并行技術

C++ AMP[15]是Microsoft Corporation開發(fā)的一套C++應用程序加速技術，使用AMP并行技術編寫的程序可在圖形處理器（graphics processing unit，GPU）等硬件上進行計算。C++AMP的使命是將GPU編程帶給每一位開發(fā)者，由于支持C++AMP的顯卡沒有限制，并且它包含了近乎全面的代碼庫，使用這些庫的代碼無需了解底層代碼就能獲得加速，這一點對于創(chuàng)建面向船體興波阻力預報的代碼庫而言很重要。

C++AMP的“入口”即parallel_for_each (array.extent, [=] (index＜1＞ idx) restrict (amp){kernel執(zhí)行體}，主要用到array容器，它用于儲存加速器上同時計算的一組數據，index＜N＞用于索引指定的某個位置，其整數值N為數組維度，extent代表執(zhí)行線程的數量，kernel執(zhí)行體為計算函數，也是AMP的核心，為了調用函數庫amp.h，必須使用restrict(amp)標記。計算環(huán)境運行在一臺配置為AMD 銳龍5 3500U CPU @2.1 GHz；操作系統(tǒng)為Windows 10；運行內存8 GB的筆記本上。如圖12，在Dawson方法的求解過程中可以在多個部分用到并行技術，主要是影響系數的計算。

圖12 Dawson方法求解過程

如圖13，在指定的面元數下，正常計算需要70～300 s，通過AMP方法加速計算后，基本在10～50 s左右就能算完。綠色虛線代表正常計算與AMP并行計算消耗時間的比例，總時間提速5～6倍。

圖13 Dawson方法加速計算分析

4 數值方法驗證

4.1 Wigley船型的興波阻力預報

圖14為以自由面橫向網格數為變量的興波阻力系數收斂圖，從圖14中可以看出，貼體網格的計算速度明顯快于流線網格，隨著面元數的增加興波阻力系數逐漸收斂。當Fr=0.316時，如圖15波高分布和圖16船體吃水處波切圖，從圖15和圖16中較為明顯地看出有首尾波興起，船體兩側出現開爾文波形，符合船體航行特征。

圖14 不同面元數下的興波阻力系數收斂圖

圖15 Fr=0.316時Wigley船型波高分布

如圖17為Wigley船型興波阻力系數隨不同Fr數下航行的變化曲線，對比值來自Kim[16]和陳紀康等[17]。由圖17可以看出，本文計算結果隨Fr數波動趨勢一致，與文獻值對比有一定的差距，考慮這種差距是由網格劃分的差異以及數值計算過程中本身的誤差所致。針對船體近水線網格劃分，采用NURBS曲面的離散精度較低，而流線網格是通過忽略黏性的疊模流場計算出來的，計算精度存在誤差疊加。從本文流線網格和貼體網格程序的計算結果對比來看，后者的計算速度提升30%左右，波高的幅值較前者略大，但興波阻力的預報結果較為接近，且貼體網格適用于方尾船型、多體船等復雜船型的興波阻力預報[11-12]，因此采用AMP并行技術和貼體網格的結合可實現船體興波阻力的快速預報。

圖16 Fr=0.316時船體表面波切

圖17 阻力計算

4.2 KCS船型的興波阻力預報

采用本文建立的快速預報方法對KCS船型進行計算，圖18為不同面元數下的興波阻力系數收斂圖。從圖18中可以看出，興波阻力系數具有較好的收斂值。當Fr為0.26和0.35時，如圖19為不同航速下的船體航行波高分布圖，從圖19中較為明顯地看出有首尾波興起，船體兩側出現開爾文波形，Fr數越大自由面興波起伏越大。如圖20為KCS船型興波阻力系數隨不同Fr數下航行的變化曲線，其中對比值來自NM理論[3]和高階Rankine源法[18]?？梢钥闯觯嬎憬Y果雖有些不同，但隨Fr數變化較為一致，說明本文的開發(fā)程序計算結果可靠。

圖18 不同面元數下的興波阻力系數收斂圖

圖19 KCS船型波高分布

圖20 KCS船型阻力計算

5 結論

本文編寫了基于Dawson方法的船體興波阻力計算程序，實現了典型船型的興波阻力預報，并對興波阻力的快速預報與精確求解展開研究與分析，可得出如下結論：

1）以航行體網格、自由面區(qū)域大小和自由面網格劃分3方面關鍵參數對興波阻力計算結果進行收斂性驗證，確定了Wigley和KCS船型數值離散的網格劃分策略；

2）在Dawson方法的求解過程中多處引入并行技術，采用AMP方法興波阻力計算提速5～6倍，數十秒完成一次預報，計算效率顯著提升；

3）自由面離散采用體貼網格代替流線網格，計算速度提升30%左右，波高的幅值較前者略大，興波阻力的預報結果收斂性較好；

4）在Dawson方法中引入貼體網格和AMP并行技術，建立了船體興波阻力快速預報方法，在兩條標模船模興波阻力的預報上取得了實效。