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

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

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

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

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

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

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

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

1.1 興波理論

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

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

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

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

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

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

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

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

1.2 數(shù)值離散

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

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

圖1 計(jì)算模型示意

1.3 貼體網(wǎng)格

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

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

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

式中J為雅可比矩陣。

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

2 數(shù)值計(jì)算

2.1 航行體模型

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

圖3 Wigley船型

圖4 KCS船型

表1 船型參數(shù)

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

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

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

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

圖5 自由面區(qū)域

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

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

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

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

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

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

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

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

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

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

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

表2 計(jì)算模型數(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）等硬件上進(jìn)行計(jì)算。C++AMP的使命是將GPU編程帶給每一位開發(fā)者，由于支持C++AMP的顯卡沒(méi)有限制，并且它包含了近乎全面的代碼庫(kù)，使用這些庫(kù)的代碼無(wú)需了解底層代碼就能獲得加速，這一點(diǎn)對(duì)于創(chuàng)建面向船體興波阻力預(yù)報(bào)的代碼庫(kù)而言很重要。

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

圖12 Dawson方法求解過(guò)程

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

圖13 Dawson方法加速計(jì)算分析

4 數(shù)值方法驗(yàn)證

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

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

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

圖15 Fr=0.316時(shí)Wigley船型波高分布

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

圖16 Fr=0.316時(shí)船體表面波切

圖17 阻力計(jì)算

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

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

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

圖19 KCS船型波高分布

圖20 KCS船型阻力計(jì)算

5 結(jié)論

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

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

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

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

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