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

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

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

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

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

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

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

1 理論與數(shù)值離散模型

1.1 興波理論

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

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

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

式中：下標(biāo)x、y、z為該變量沿x、y、z方向的偏導(dǎo)數(shù)，Φx為總速度勢沿x方向的偏導(dǎo)數(shù)。

自由面邊界條件可推導(dǎo)[2]為

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

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

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

1.2 數(shù)值離散

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

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

圖1 計算模型示意

1.3 貼體網(wǎng)格

圖2 自由面網(wǎng)格圖

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

式中：(ξ，η)為物理坐標(biāo)系下的點，(x，y)為大地坐標(biāo)系下的點。縱向和橫向的導(dǎo)數(shù)都采用單邊上游有限差分算子：

式中J為雅可比矩陣。

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

2 數(shù)值計算

2.1 航行體模型

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

圖3 Wigley船型

圖4 KCS船型

表1 船型參數(shù)

2.2 數(shù)值模型收斂性驗證

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

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

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

圖5 自由面區(qū)域

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

2.2.2 航行體網(wǎng)格劃分收斂性驗證

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

圖7 航行體網(wǎng)格

圖8 興波阻力隨船體網(wǎng)格數(shù)量的變化

2.2.3 自由面網(wǎng)格劃分收斂性驗證

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

圖9 自由面網(wǎng)格劃分示意

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

2.3 數(shù)值模型網(wǎng)格劃分策略

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

表2 計算模型數(shù)值離散網(wǎng)格值

圖11 KCS船體及自由面網(wǎng)格示意

3 并行技術(shù)

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

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

圖12 Dawson方法求解過程

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

圖13 Dawson方法加速計算分析

4 數(shù)值方法驗證

4.1 Wigley船型的興波阻力預(yù)報

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

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

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

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

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

圖17 阻力計算

4.2 KCS船型的興波阻力預(yù)報

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

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

圖19 KCS船型波高分布

圖20 KCS船型阻力計算

5 結(jié)論

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

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

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

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

4）在Dawson方法中引入貼體網(wǎng)格和AMP并行技術(shù)，建立了船體興波阻力快速預(yù)報方法，在兩條標(biāo)模船模興波阻力的預(yù)報上取得了實效。