船舶操縱水動(dòng)力導(dǎo)數(shù)的數(shù)值求解及敏感度分析

孫寒冰，肖佳峰，王偉，劉偉杰，鄭興*

1 哈爾濱工程大學(xué) 船舶工程學(xué)院, 黑龍江 哈爾濱 150001 2 中國船舶工業(yè)系統(tǒng)工程研究院, 北京 100036

0 引 言

船舶操縱性是船舶重要的航行性能之一，它是船舶在控制裝置的作用下，按照駕駛者的意圖保持或改變船舶運(yùn)動(dòng)狀態(tài)（包括航向、航速和位置等）的能力[1]。隨著國際海事組織（IMO）關(guān)于船舶操縱性標(biāo)準(zhǔn)的實(shí)行，在船舶設(shè)計(jì)之初評(píng)估船舶的操縱性變得越來越重要[2-3]。

第25屆國際船模拖曳水池會(huì)議（ITTC）專業(yè)技術(shù)（操縱性）委員會(huì)將目前的幾類船舶操縱性預(yù)報(bào)方法進(jìn)行總結(jié)，分為了3類[4]：直接預(yù)報(bào)法、基于系統(tǒng)的操縱模擬方法和計(jì)算流體力學(xué)（CFD）數(shù)值模擬方法。直接預(yù)報(bào)法和基于系統(tǒng)的操縱模擬方法是通過經(jīng)驗(yàn)估算公式、船模以及實(shí)船試驗(yàn)來預(yù)報(bào)船舶的操縱性，存在船舶依賴性強(qiáng)、花費(fèi)的時(shí)間多、金錢多等問題，不適用于船舶初步設(shè)計(jì)階段。隨著計(jì)算流體力學(xué)的飛速發(fā)展，CFD數(shù)值模擬方法為在船舶設(shè)計(jì)階段快速預(yù)報(bào)操縱性提供了一種新的手段，大量學(xué)者針對(duì)CFD數(shù)值模擬方法進(jìn)行了研究?；贑FD數(shù)值模擬方法預(yù)報(bào)船舶操縱性的方法可以分為2大類：預(yù)報(bào)操縱運(yùn)動(dòng)參數(shù)的自航船模操縱模擬和預(yù)報(bào)水動(dòng)力導(dǎo)數(shù)的約束船模操縱模擬，其中自航運(yùn)動(dòng)包括直航、回轉(zhuǎn)、Z形等典型的操縱運(yùn)動(dòng)。

在CFD自航船模操縱模擬方面，Carrica等[5]采用重疊網(wǎng)格方法，通過求解URANS方程，對(duì)DTMB 5415船模進(jìn)行了回轉(zhuǎn)運(yùn)動(dòng)和Z形操縱運(yùn)動(dòng)的自航船模數(shù)值模擬。Wang等[6]采用重疊網(wǎng)格方法，對(duì)全附體雙槳雙舵ONRT船模進(jìn)行了35°舵角工況下的回轉(zhuǎn)運(yùn)動(dòng)數(shù)值模擬，預(yù)報(bào)的回轉(zhuǎn)操縱性參數(shù)同試驗(yàn)結(jié)果吻合良好，驗(yàn)證了naoe-FOAM-SJTU求解器對(duì)全附體船舶自航操縱問題的適用性。Mofidi等[7]通過將基于渦格法的升力面螺旋槳體積力求解模型PUF-14與CFD求解器RXE相結(jié)合，對(duì)KCS自航船舶進(jìn)行了操縱模擬，較為精確地捕捉了操縱運(yùn)動(dòng)過程中船?槳?舵之間的干擾。Jin等[8]通過求解非定常RANS方程，使用BFM和DPM這2種不同的推力模擬技術(shù)對(duì)DTMB 5415船模進(jìn)行了回轉(zhuǎn)運(yùn)動(dòng)和Z形操縱試驗(yàn)的自航數(shù)值模擬，并進(jìn)行了比較，從中可以看出，自航模擬能夠直接預(yù)報(bào)操縱運(yùn)動(dòng)參數(shù)，但需同時(shí)對(duì)船?槳?舵系統(tǒng)的復(fù)雜運(yùn)動(dòng)和流場(chǎng)進(jìn)行求解，且大部分要求使用二次開發(fā)的求解器才能較好地模擬船?槳?舵之間的干擾。然而，這種方法對(duì)計(jì)算資源的需求量很大，若對(duì)各工況逐個(gè)計(jì)算將導(dǎo)致時(shí)間成本非常高。

在CFD約束船模操縱模擬方面，Kim等[9]采用自主開發(fā)的求解器SHIP_Motion，通過求解RANS方程，對(duì) KCS船模進(jìn)行了不同形式的平面運(yùn)動(dòng)機(jī)構(gòu)（planar motion mechanism，PMM）試驗(yàn)數(shù)值模擬，其靜態(tài)試驗(yàn)和橫蕩試驗(yàn)?zāi)M結(jié)果與試驗(yàn)結(jié)果吻合較好，但純艏搖試驗(yàn)的模擬結(jié)果與試驗(yàn)結(jié)果存在一定的誤差。Sakamoto等[10]采用自主開發(fā)的程序 CFD Ship-Iwoa對(duì)DTMB 5415船模進(jìn)行了靜態(tài)和動(dòng)態(tài)的PMM試驗(yàn)數(shù)值模擬，預(yù)報(bào)得到了相關(guān)的水動(dòng)力導(dǎo)數(shù)，并與試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。Roy-Choudhury等[11]基于商業(yè)軟件SHIPFLOW對(duì)KVLCC2船模的純艏搖和帶漂角艏搖等動(dòng)態(tài)約束船模試驗(yàn)進(jìn)行了數(shù)值模擬，并對(duì)不同水深下操縱運(yùn)動(dòng)的水動(dòng)力特性進(jìn)行了研究，結(jié)果證明在不同水深下采用CFD方法進(jìn)行操縱性模擬是可行的。程捷等[12]使用商業(yè)軟件Fluent對(duì)DTMB 5415船模進(jìn)行了動(dòng)態(tài)純橫蕩和純艏搖試驗(yàn)的數(shù)值模擬，計(jì)算得到的橫向力及回轉(zhuǎn)力矩的變化趨勢(shì)與實(shí)驗(yàn)結(jié)果一致，但忽略了自由面興波、船模的縱搖和垂蕩，精度還有待提高。李冬琴等[13]利用商用軟件Fine/Marine對(duì)DTMB 5415船模進(jìn)行了低頻純橫蕩及純艏搖運(yùn)動(dòng)的數(shù)值仿真，其數(shù)值結(jié)果與試驗(yàn)結(jié)果吻合較好，但忽略了自由面興波的影響。王建華等[14]采用重疊網(wǎng)格方法，通過基于OpenFOAM自主研發(fā)的水動(dòng)力學(xué)求解器naoe-FOAM-SJTU，對(duì)標(biāo)準(zhǔn)船模DTMB 5415進(jìn)行了純搖艏操縱運(yùn)動(dòng)模擬，結(jié)果顯示基于約束船模模擬預(yù)報(bào)的水動(dòng)力導(dǎo)數(shù)能夠代入操縱運(yùn)動(dòng)方程中用于預(yù)報(bào)操縱運(yùn)動(dòng)參數(shù)。Cura-Hochbaum[15]通過求解RANS方程，對(duì)KVLCC1船模的動(dòng)態(tài)約束試驗(yàn)進(jìn)行了模擬，并將回歸得到的水動(dòng)力導(dǎo)數(shù)代入操縱運(yùn)動(dòng)方程，得到了回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)的軌跡，其結(jié)果與試驗(yàn)結(jié)果大體一致，表明基于RANS方法數(shù)值預(yù)報(bào)船舶操縱性也是可行的。Guo等[16]以雙槳雙舵ONRT船模為研究對(duì)象，使用商業(yè)軟件STAR CCM+進(jìn)行了旋轉(zhuǎn)試驗(yàn)、斜拖試驗(yàn)和純橫傾試驗(yàn)的數(shù)值模擬，預(yù)報(bào)了相應(yīng)的操縱性水動(dòng)力導(dǎo)數(shù)，并結(jié)合四自由度的MMG（manoeuvring model group）模型對(duì)回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)進(jìn)行了仿真，其結(jié)果與自航試驗(yàn)結(jié)果吻合較好，證明采用CFD方法預(yù)報(bào)操縱性參數(shù)可靠。Hajivand等[17]對(duì)DTMB 5512船模進(jìn)行了不同形式的PMM試驗(yàn)數(shù)值模擬，并將獲得的水動(dòng)力導(dǎo)數(shù)代入操縱性模型模擬得到了自由回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)的軌跡。Sukas等[18]基于URANS方法，對(duì)DTMB 5415約束船模的操縱運(yùn)動(dòng)試驗(yàn)、自推進(jìn)和舵力試驗(yàn)進(jìn)行模擬，預(yù)報(bào)了水動(dòng)力導(dǎo)數(shù)、推力和舵力等相關(guān)參數(shù)，基于操縱運(yùn)動(dòng)方程模擬了自由回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)的軌跡，并對(duì)比了采用單次工況擬合和二次擬合這2種方法擬合水動(dòng)力導(dǎo)數(shù)的區(qū)別，從中可以看到提高水動(dòng)力導(dǎo)數(shù)預(yù)報(bào)的精度，進(jìn)而提高操縱預(yù)報(bào)精度是目前關(guān)注的要點(diǎn)。

本文將主要考慮如何在船舶設(shè)計(jì)初期開展操縱性能預(yù)報(bào)，在不降低操縱運(yùn)動(dòng)預(yù)報(bào)精度的前提下盡可能減少數(shù)值計(jì)算工作，提高船舶操縱運(yùn)動(dòng)預(yù)報(bào)效率。首先 ，以DTMB 5415船模為研究對(duì)象，開展拘束船模數(shù)值模擬計(jì)算，通過直航、純橫蕩以及純艏搖運(yùn)動(dòng)獲得線性水動(dòng)力導(dǎo)數(shù)，并通過經(jīng)驗(yàn)公式得到非線性水動(dòng)力導(dǎo)數(shù)，進(jìn)而求解船舶操縱運(yùn)動(dòng)方程，進(jìn)行回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)仿真，并將仿真結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證所提預(yù)報(bào)方法的可行性；然后，開展水動(dòng)力導(dǎo)數(shù)對(duì)于操縱性參數(shù)的敏感度分析，進(jìn)一步驗(yàn)證預(yù)報(bào)模型的合理性；最后，采用本文方法并結(jié)合數(shù)值計(jì)算方法與計(jì)算機(jī)仿真，針對(duì)特定船型預(yù)報(bào)其水動(dòng)力導(dǎo)數(shù)，以提高水動(dòng)力導(dǎo)數(shù)預(yù)報(bào)的精度，同時(shí)使用經(jīng)驗(yàn)公式計(jì)算敏感度較低的水動(dòng)力導(dǎo)數(shù)，用以降低計(jì)算成本。

1 數(shù)學(xué)模型的建立

1.1 三自由度操縱運(yùn)動(dòng)數(shù)學(xué)模型

本文采用日本拖曳水池委員會(huì)提出的MMG分離型數(shù)學(xué)模型，將作用于船舶上的流體動(dòng)力和力矩按照其物理意義，分離成船體、槳、舵的單獨(dú)受力及其之間的相互干涉。船舶在水平面上的操縱運(yùn)動(dòng)方程如下：

式中：m為船體質(zhì)量；mx，my分別為船體在x軸和y軸方向的附加質(zhì)量；X(u)為船體直航時(shí)的阻力；u為船體的縱向速度；v為船體的橫向速度；r為水平面內(nèi)的轉(zhuǎn)艏角速度；Izz為船體質(zhì)量繞z軸的轉(zhuǎn)動(dòng)慣量；Jzz為船體繞z軸轉(zhuǎn)動(dòng)的附加慣性質(zhì)量；Xvv，Xrr，Xvr，Yv，Yr，Y|v|v，Y|r|r，Yvrr，Yvvr，Nv，Nr，N|v|v，N|r|r，Nvrr，Nvvr為采用貴島數(shù)學(xué)模型的黏性流體動(dòng)力和力矩導(dǎo)數(shù)；XP為 螺旋槳產(chǎn)生的推力；XR，YR，NR為舵產(chǎn)生的水動(dòng)力和力矩。

1.1.1 螺旋槳模型

縱向螺旋槳推力XP的表達(dá)式如下：

1.1.2 舵力模型

舵力的表達(dá)式如下：

式中：tR為 推力減額系數(shù)；αH為船體橫向力受到操舵影響的系數(shù)；xH為操舵誘導(dǎo)的船體橫向力作用中心至船舶重心的距離；xR為舵中心的縱向坐標(biāo)；，分別為左、右舵中心的橫向坐標(biāo)； δP，δS分別為左、右舵的舵角，其中右舵為正；，分別為左、右舵的法向力，具體的表達(dá)式如下：

式中：AR為舵面積；fα為 舵法向力系數(shù)，UR，αR分別為計(jì)入船體、螺旋槳影響后舵前的有效流速和有效舵角。

1.2 純橫蕩運(yùn)動(dòng)數(shù)學(xué)模型

船模做純橫蕩運(yùn)動(dòng)時(shí)，艏向角在運(yùn)動(dòng)過程中不發(fā)生變化，而是在沿縱向勻速運(yùn)動(dòng)的同時(shí)船模疊加了一個(gè)正弦變化的橫向位移。式（5）為純橫蕩運(yùn)動(dòng)方程。

式中：ψ為艏向角；a為橫向簡(jiǎn)諧振蕩運(yùn)動(dòng)幅值；ymax，vmax，v˙max分別為橫向位移、橫向速度和橫向加速度最大幅值；ω為簡(jiǎn)諧運(yùn)動(dòng)頻率， ω=2πf，其中f為運(yùn)動(dòng)頻率；t為時(shí)間。

1.3 純艏搖運(yùn)動(dòng)數(shù)學(xué)模型

船模做純艏搖運(yùn)動(dòng)時(shí)，漂角在運(yùn)動(dòng)過程中不發(fā)生變化，而在沿縱向勻速運(yùn)動(dòng)的同時(shí)疊加了一個(gè)余弦變化的艏向角，模型運(yùn)動(dòng)軌跡與其縱軸處相切。式（6）為純艏搖運(yùn)動(dòng)方程。

式中： ψmax，rmax，r˙max分別為艏向角、轉(zhuǎn)艏角速度和轉(zhuǎn)艏角加速度幅值；U為速度。

1.4 控制方程

本文使用STAR-CCM+軟件對(duì)船舶受到的水動(dòng)力進(jìn)行模擬，數(shù)值求解黏性流場(chǎng)中的Navier-Stokes方程。通過將黏性流場(chǎng)中的隨機(jī)脈動(dòng)項(xiàng)進(jìn)行時(shí)間平均，進(jìn)而將非定常問題轉(zhuǎn)化為對(duì)定常問題的求解。黏性流動(dòng)的連續(xù)方程如下所示：

雷諾平均的納維?斯托克斯（RANS）方程為

式中： ρuiuj為雷諾應(yīng)力項(xiàng)；上劃線“—”代表對(duì)物理量進(jìn)行時(shí)間平均；μ為動(dòng)力黏性系數(shù)；p為壓力時(shí)均值；ui,uj為速度分量時(shí)均值，xi,xj為位移分量時(shí)均值；si為 單位體積力；i，j=1，2，3。

2 CFD數(shù)值模擬

2.1 研究對(duì)象

本文的計(jì)算對(duì)象為DTMB 5415標(biāo)準(zhǔn)船模，為美國海軍艦艇模型，也是ITTC推薦的標(biāo)準(zhǔn)船模之一。DTMB 5415標(biāo)準(zhǔn)船模的模型主要參數(shù)如表1所示，幾何模型如圖1所示。表中，“全尺度模型”指DTMB 5415的實(shí)船尺度，“MARIN模型”為荷蘭海事研究所實(shí)驗(yàn)用船模。

表1 DTMB 5 415船模主要參數(shù)Table 1 Main parameters of DTMB 5 415 ship model

圖1 DTMB 5 415模型示意圖Fig.1 Schematic diagram of DTMB 5 415 ship model

2.2 數(shù)值方案與網(wǎng)格收斂性分析

本文采用MARIN模型進(jìn)行數(shù)值模擬。由于船模的對(duì)稱性，在直航模擬時(shí)為了節(jié)省計(jì)算量并提高計(jì)算效率，只對(duì)y>0的半個(gè)計(jì)算域進(jìn)行了計(jì)算。在進(jìn)行PMM模擬時(shí)，則對(duì)整個(gè)計(jì)算域進(jìn)行了計(jì)算。

進(jìn)行數(shù)值計(jì)算時(shí)，流場(chǎng)采用非定常RANS方程模擬，參考文獻(xiàn)[14]，選擇SSTk-ω湍流模型并結(jié)合Ally+Wall Treatment壁面處理方式，湍流動(dòng)能方程和湍流耗散比率則均采用二階迎風(fēng)格式求解。利用流體體積（VOF）法捕捉自由液面，采用剛體運(yùn)動(dòng)DFBI （dynamic fluid body interaction）模塊對(duì)船模約束運(yùn)動(dòng)進(jìn)行數(shù)值模擬。

邊界條件：船體使用無滑移的壁面邊界條件；計(jì)算流域入口選用速度進(jìn)口邊界；流域出口選用壓力出口邊界；流域頂部和底部選用速度進(jìn)口邊界；流域側(cè)面在進(jìn)行直航模模擬時(shí)選用對(duì)稱平面邊界，進(jìn)行約束模模擬時(shí)選用速度進(jìn)口邊界。進(jìn)行直航運(yùn)動(dòng)模擬時(shí)，在除頂部、底部以及右側(cè)面之外的流域表面設(shè)置VOF法波阻尼；進(jìn)行約束運(yùn)動(dòng)模擬時(shí)，在除頂部和底部之外的流域表面設(shè)置VOF法波阻尼，用以減小波的振蕩。

計(jì)算流域選為長方體，邊界入口距船艏1.5L（L為船長），邊界出口距船艉3.5L。周向邊界中的上、下邊界分別在距龍骨2/3L和4/3L處，側(cè)面邊界在距船中縱剖面1.5L處，船模在流域中的相對(duì)位置如圖2所示。

圖2 DTMB 5 415船模在流域中的位置及邊界條件示意圖Fig.2 Computational domain and boundary conditions of DTMB 5415 ship model

計(jì)算流域采用切割體進(jìn)行網(wǎng)格劃分，其基本原則是在兼顧計(jì)算精度的情況下提高計(jì)算效率。因此，對(duì)計(jì)算模型周圍區(qū)域和自由液面進(jìn)行網(wǎng)格加密，以精確捕捉船舶運(yùn)動(dòng)過程中周圍流體的流動(dòng)特征；在離船舶較遠(yuǎn)處，則采用較為稀疏的網(wǎng)格。另外，為了細(xì)化船體的銳角或銳邊，對(duì)船艏和船艉進(jìn)行網(wǎng)格加密。在沿船體壁面的法線方向上，采用棱柱層網(wǎng)格劃分邊界層，用以捕捉近壁面的流動(dòng)細(xì)節(jié)。棱柱層網(wǎng)格設(shè)為8層，船體表面第1層網(wǎng)格與壁面的距離y+取為100。模型的船體表面網(wǎng)格劃分及計(jì)算域體網(wǎng)格劃分如圖3所示。

圖3 網(wǎng)格劃分方案Fig.3 Mesh division schemes

為了確定網(wǎng)格尺寸對(duì)計(jì)算結(jié)果的影響并驗(yàn)證網(wǎng)格方案的有效性，以Fr= 0.28時(shí)的直航運(yùn)動(dòng)為研究對(duì)象，進(jìn)行網(wǎng)格收斂性分析。將網(wǎng)格尺寸以為比例系數(shù)進(jìn)行加密和加粗。

在進(jìn)行網(wǎng)格收斂性分析（表2）時(shí)，通過收斂性參數(shù)RG來判斷計(jì)算結(jié)果的穩(wěn)定性，定義如下：

表2 網(wǎng)格收斂性分析Table 2 Grid convergence analysis

式中：S1，S2，S3分別為細(xì)網(wǎng)格、中等網(wǎng)格以及粗網(wǎng)格的計(jì)算結(jié)果。當(dāng)0

由表2可見，計(jì)算結(jié)果一致收斂，其隨網(wǎng)格數(shù)量的增大逐漸穩(wěn)定，相鄰網(wǎng)格計(jì)算結(jié)果之間的差距也在逐漸減小，且不同網(wǎng)格的計(jì)算結(jié)果與試驗(yàn)結(jié)果間吻合較好，誤差均在3%以內(nèi)，計(jì)算精度較高。為了在兼顧計(jì)算精度的情況下提高計(jì)算效率，本文采用中等網(wǎng)格進(jìn)行了模擬。

2.3 直航阻力驗(yàn)證

表3 DTMB 5 415船模的總阻力CFD計(jì)算值與試驗(yàn)值比較Table 3 Comparison of the CFD calculated results and experimental data of total resistance for DTMB 5 415 ship model

將表3中的CFD計(jì)算值繪制成總阻力曲線，如圖4所示，并擬合出DTMB 5415船模的X(u)函數(shù)為：

圖4 CFD阻力計(jì)算值、試驗(yàn)值及擬合曲線Fig.4 The calculated results, experimental data and fitting curveof CFD resistance

2.4 約束運(yùn)動(dòng)模擬

對(duì)多個(gè)工況的船舶純橫蕩運(yùn)動(dòng)和純艏搖運(yùn)動(dòng)進(jìn)行模擬，計(jì)算船舶受到的水動(dòng)力，對(duì)計(jì)算結(jié)果進(jìn)行處理后即可獲得水動(dòng)力導(dǎo)數(shù)。約束運(yùn)動(dòng)模擬的工況如表4和表5所示。表中：為橫向位移的無因次化值，為轉(zhuǎn)艏角速度的無因次化值，

表4 純橫蕩計(jì)算工況Table 4 Calculation conditions of pure sway motion

表5 純艏搖計(jì)算工況Table 5 Calculation conditions of pure yaw motion

純橫蕩運(yùn)動(dòng)和純艏搖運(yùn)動(dòng)的橫向力Y及回轉(zhuǎn)力矩N的計(jì)算結(jié)果與SIMMAN 2008[3]中試驗(yàn)結(jié)果的對(duì)比如圖5和圖6所示。

由于船舶進(jìn)行約束運(yùn)動(dòng)時(shí)速度和角速度是正弦變化的，因此受到的力和力矩也是正弦變化的，且變化趨勢(shì)相似。從圖5和圖6的比較結(jié)果可以看出，計(jì)算出的橫向力及回轉(zhuǎn)力矩均是按照三角函數(shù)的規(guī)律變化的，變化趨勢(shì)相似，變化周期等于試驗(yàn)給定的船舶運(yùn)動(dòng)周期，最大值與最小值的絕對(duì)值也相等，幅值隨著運(yùn)動(dòng)幅值的增大而增大，且計(jì)算結(jié)果的相位與試驗(yàn)結(jié)果也吻合良好；幅值上橫向力與回轉(zhuǎn)力矩的誤差約為10%，其整體趨勢(shì)與試驗(yàn)結(jié)果吻合，計(jì)算效果較好。

圖5 純橫蕩運(yùn)動(dòng)橫向力、回轉(zhuǎn)力矩的計(jì)算值與試驗(yàn)值比較Fig.5 Comparison of the calculated results and experimental data of sway force and yaw moment in pure sway motion

圖6 純艏搖運(yùn)動(dòng)橫向力、回轉(zhuǎn)力矩的計(jì)算值與試驗(yàn)值比較Fig.6 Comparison of the calculated and experimental data of sway force and yaw moment of pure yaw motion

圖7所示為純橫蕩運(yùn)動(dòng)過程中一個(gè)運(yùn)動(dòng)周期T內(nèi)4個(gè)不同時(shí)刻的波形圖，圖右上角為對(duì)應(yīng)時(shí)刻下的船體表面水體積分?jǐn)?shù)正視圖。圖8所示為純艏搖運(yùn)動(dòng)過程中一個(gè)運(yùn)動(dòng)周期T內(nèi)不同時(shí)刻的波形圖，以及隨體坐標(biāo)系x方向的船體表面水體積分?jǐn)?shù)圖。

由圖7、圖8可以發(fā)現(xiàn)，采用本文中的數(shù)值模擬方法能夠準(zhǔn)確捕捉純橫蕩運(yùn)動(dòng)和純艏搖運(yùn)動(dòng)船體的興波，船后興波軌跡很好地反映了船的運(yùn)動(dòng)軌跡；波形圖里，興波的趨勢(shì)與船體表面水體積分?jǐn)?shù)圖中船體的水線變化一致，模擬效果較好。

圖8 純艏搖運(yùn)動(dòng)興波Fig.8 The wave-making of pure yaw motion

觀察圖7還可見，由于選取的是計(jì)算穩(wěn)定后的周期，所以在t=0時(shí)刻船艏兩側(cè)波高不對(duì)稱；在t=0.25T時(shí)刻船艏兩側(cè)波高基本對(duì)稱；t=0.125T時(shí)刻，船艏兩側(cè)的水氣分布與t=0.375T時(shí)刻正好相反。另外，觀察兩側(cè)波高的變化，發(fā)現(xiàn)在0.5T之前是從左舷興波波高較大逐漸過渡到左右對(duì)稱，而后右舷興波波高更大，這與實(shí)際情況相符。

圖7 純橫蕩運(yùn)動(dòng)興波Fig.7 The wave-making of pure sway motion

觀察圖8還可見，在純艏搖運(yùn)動(dòng)下，其船艏興波與純橫蕩運(yùn)動(dòng)下的明顯不同，在0.5T之前始終是右舷興波波高更高，符合物理運(yùn)動(dòng)規(guī)律。

表6所示為通過最小二乘法二次擬合求出的水動(dòng)力導(dǎo)數(shù)的結(jié)果。從中可見，二次擬合后的水動(dòng)力導(dǎo)數(shù)與試驗(yàn)結(jié)果吻合較好，除Yν外，其他導(dǎo)數(shù)的誤差均在7%以內(nèi)，計(jì)算結(jié)果較精確，可見采用本文中數(shù)值計(jì)算方法能夠較好地預(yù)報(bào)船舶的線性水動(dòng)力導(dǎo)數(shù)。產(chǎn)生誤差的原因可能是，在運(yùn)動(dòng)過程中，船體伴流會(huì)產(chǎn)生嚴(yán)重的流動(dòng)分離現(xiàn)象，而本文數(shù)值模擬采用的雷諾時(shí)均方法無法準(zhǔn)確描述此現(xiàn)象，導(dǎo)致水動(dòng)力時(shí)歷曲線與試驗(yàn)有所差異，從而使得回歸得到的水動(dòng)力導(dǎo)數(shù)出現(xiàn)誤差。

表6 二次擬合的水動(dòng)力導(dǎo)數(shù)Table 6 Hydrodynamic derivatives obtained by quadratic fitting

2.5 非線性水動(dòng)力導(dǎo)數(shù)估算

采用數(shù)值模擬方法計(jì)算出線性水動(dòng)力導(dǎo)數(shù)后，為了預(yù)報(bào)船舶的操縱運(yùn)動(dòng)，在保證一定精度的情況下，在快速獲得水動(dòng)力導(dǎo)數(shù)后使用貴島數(shù)學(xué)模型的經(jīng)驗(yàn)公式估算無因次非線性水動(dòng)力導(dǎo)數(shù)，具體的經(jīng)驗(yàn)公式及無因次處理參考文獻(xiàn)[19]。采用經(jīng)驗(yàn)公式估算的無因次非線性水動(dòng)力導(dǎo)數(shù)如表7所示。

表7 經(jīng)驗(yàn)公式計(jì)算的水動(dòng)力導(dǎo)數(shù)Table 7 Hydrodynamic derivatives calculated by empirical formulas

3 船舶操縱性運(yùn)動(dòng)仿真及驗(yàn)證

將計(jì)算得到的水動(dòng)力導(dǎo)數(shù)代入運(yùn)動(dòng)方程后，采用四階定步長龍格?庫塔法進(jìn)行求解，可得到船舶任意時(shí)刻的加速度，然后對(duì)加速度進(jìn)行數(shù)值積分，即可獲得任意時(shí)刻的位置和速度，進(jìn)而預(yù)報(bào)船舶運(yùn)動(dòng)。至此，便完成了仿真模型的搭建。選取Fr=0.248進(jìn)行?35°舵角的回轉(zhuǎn)運(yùn)動(dòng)和?20°/?20°的Z形運(yùn)動(dòng)。為便于與SIMMAN 2008的試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，將結(jié)果換算成了實(shí)船尺度，對(duì)比結(jié)果如圖9～圖10所示。

圖10 Z形運(yùn)動(dòng)仿真值與試驗(yàn)值的對(duì)比Fig.10 Comparison of the simulated results and experimental data for zigzag motion

由圖9及表7可見，仿真與試驗(yàn)運(yùn)動(dòng)的軌跡及角速度吻合較好，仿真的回轉(zhuǎn)運(yùn)動(dòng)參數(shù)與試驗(yàn)值相比誤差大多在10%以內(nèi)，為可接受范圍，說明仿真結(jié)果能較好地反映試驗(yàn)結(jié)果。

圖9 回轉(zhuǎn)運(yùn)動(dòng)中仿真值與試驗(yàn)值的對(duì)比Fig.9 Comparison of the simulated results and experimental data for turning motion

另外，表8還對(duì)3篇文獻(xiàn)中的結(jié)果進(jìn)行了對(duì)比，其中SB指基于CFD進(jìn)行約束船模模擬，得到所有水動(dòng)力導(dǎo)數(shù)后預(yù)報(bào)得到的回轉(zhuǎn)運(yùn)動(dòng)參數(shù)；CFDB指基于CFD進(jìn)行自航船模模擬后得到的回轉(zhuǎn)運(yùn)動(dòng)參數(shù)；EMP指完全通過經(jīng)驗(yàn)公式預(yù)報(bào)得到的回轉(zhuǎn)運(yùn)動(dòng)參數(shù)。將每種方法的結(jié)果與試驗(yàn)值相比較，計(jì)算平均誤差，從中可以看到，預(yù)報(bào)誤差的平均誤差順序?yàn)镋MP >SB>CFDB>本文方法，其中經(jīng)驗(yàn)公式方法的誤差最大，其他方法的預(yù)報(bào)誤差均小于10%，精度可以滿足工程需要。表8中，Ad為 縱距，Tr為 橫 距，DT為戰(zhàn)術(shù)回轉(zhuǎn)直 徑，Do為穩(wěn)定回轉(zhuǎn)直徑，rs為穩(wěn)定回轉(zhuǎn)角速度，Us為穩(wěn)定回轉(zhuǎn)速度。

表8 回轉(zhuǎn)運(yùn)動(dòng)參數(shù)仿真值與試驗(yàn)值的對(duì)比Table 8 Comparison of the simulated results and experimental data for turning motion parameters

由圖10及表9可見，仿真與試驗(yàn)曲線的走向大體一致，誤差隨著時(shí)間的積累而增大；仿真得到的的Z形操縱運(yùn)動(dòng)參數(shù)誤差大部分在10%以內(nèi)，超越角的誤差較大，約達(dá)20%，這是因?yàn)槌浇堑慕^對(duì)值較小，約為5°，而1°的誤差就可以造成 20%的相對(duì)誤差，仿真得到的超越角與試驗(yàn)值間的誤差在1°左右，絕對(duì)誤差量級(jí)較小，可以接受。另外從表9中可以看到，采用3種方法預(yù)報(bào)的結(jié)果其平均誤差順序?yàn)镃FDB>本方法>SB，說明使用CFD方法對(duì)約束船模進(jìn)行模擬得到的所有水動(dòng)力導(dǎo)數(shù)能夠進(jìn)一步提高Z形運(yùn)動(dòng)的預(yù)報(bào)精度，但也意味著要進(jìn)行大量的數(shù)值計(jì)算。表9中， ψmax1為 第1最大艏向角， ψmax2為第2最大艏向角， ψOV1為 第1超越角， ψOV2為第2超越角，tα為初轉(zhuǎn)期，tl為轉(zhuǎn)艏滯后。

表9 Z形操縱運(yùn)動(dòng)參數(shù)仿真值與試驗(yàn)值的對(duì)比Table 9 Comparison of simulated results and experiment data for zigzag motion parameters

以上仿真結(jié)果表明，本文所提方法雖然犧牲了一定的預(yù)報(bào)精度，但大大節(jié)約了計(jì)算成本，且預(yù)報(bào)精度也能滿足工程應(yīng)用，能夠較為準(zhǔn)確地預(yù)報(bào)出DTMB 5415船模的自由回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)軌跡，更適用于船舶設(shè)計(jì)階段操縱性的驗(yàn)證及優(yōu)化。

4 水動(dòng)力導(dǎo)數(shù)敏感度分析

基于建立的仿真模型，對(duì)水動(dòng)力導(dǎo)數(shù)的敏感度進(jìn)行分析，研究線性水動(dòng)力導(dǎo)數(shù)和非線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)參數(shù)的影響，并對(duì)該預(yù)報(bào)方法的合理性進(jìn)行驗(yàn)證。

敏感度分析的公式如下：

式中：S為相對(duì)敏感度；R?，R分別為改變前、后的操縱運(yùn)動(dòng)參數(shù)；H?，H分別為改變前、后的水動(dòng)力導(dǎo)。取?100%～100%的變化范圍，每50%變化一次，相對(duì)敏感度取4次變化的平均值。水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)參數(shù)的敏感度如表10和表11所示。

由表10和表11可見，大部分非線性水動(dòng)力導(dǎo)數(shù)的平均相對(duì)敏感度在3%以下，線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響較大，幾乎都在3%以上。因此，從一定程度上可以說，非線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)模擬的影響比線性水動(dòng)力導(dǎo)數(shù)小，對(duì)非線性水動(dòng)力導(dǎo)數(shù)精度的要求比線性水動(dòng)力導(dǎo)數(shù)的更低，因而采用經(jīng)驗(yàn)公式估算非線性水動(dòng)力導(dǎo)數(shù)是合理的，這也是本文預(yù)報(bào)方法中非線性水動(dòng)力導(dǎo)數(shù)采用經(jīng)驗(yàn)公式能快速獲得的原因。

表10 水動(dòng)力導(dǎo)數(shù)對(duì)回轉(zhuǎn)運(yùn)動(dòng)參數(shù)的敏感度Table 10 Sensitivity of hydrodynamic derivatives to turning motion parameters

表11 水動(dòng)力導(dǎo)數(shù)對(duì)Z形操縱運(yùn)動(dòng)參數(shù)的敏感度Table 11 Sensitivity of hydrodynamic derivatives to zigzag motion parameters

另外，由表10和表11中各個(gè)水動(dòng)力導(dǎo)數(shù)敏感度的平均值可知，Nv，Nr，Nrr對(duì)回轉(zhuǎn)運(yùn)動(dòng)具有較高的影響，相對(duì)敏感度在10%以上；Yvv，Nvvr對(duì)回轉(zhuǎn)運(yùn)動(dòng)具有中等的影響，相對(duì)敏感度在5%～10%之間，Nr˙，Xvv，Xrr，Yrr，Yvvr，Nvrr對(duì)回轉(zhuǎn)運(yùn)動(dòng)的影響極其微弱，相對(duì)敏感度在1%以下；Yv，Nv，Nr，Nr˙，Yvv對(duì)Z形運(yùn)動(dòng)具有較高的影響，Yν˙，Yr，Nrr對(duì)Z形運(yùn)動(dòng)具有中等的影響，Xvr，Xvv，Xrr，Yrr，Yvvr，Nvrr對(duì)Z形運(yùn)動(dòng)的影響極其微弱。尤其是Nv和Nr，其對(duì)操縱運(yùn)動(dòng)具有很高的敏感度，相對(duì)敏感度約為30%。當(dāng)然，非線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響也不可忽略，這也是采用本文方法時(shí)產(chǎn)生誤差的原因之一。

從總體上看，回轉(zhuǎn)力矩項(xiàng)的水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響最大，橫向力項(xiàng)次之，縱向力項(xiàng)的非線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響極其微弱；在Z形操縱運(yùn)動(dòng)參數(shù)方面，水動(dòng)導(dǎo)數(shù)對(duì)轉(zhuǎn)艏滯后的影響程度要遠(yuǎn)大于Z形操縱運(yùn)動(dòng)的其他參數(shù)。

5 結(jié) 論

本文將CFD技術(shù)與理論計(jì)算相結(jié)合，基于CFD技術(shù)對(duì)DTMB 5415船模的直航、純橫蕩以及純艏搖運(yùn)動(dòng)進(jìn)行數(shù)值模擬，回歸得到了線性水動(dòng)力導(dǎo)數(shù)，隨后又使用經(jīng)驗(yàn)公式計(jì)算出了非線性水動(dòng)力導(dǎo)數(shù)。然后將獲得的水動(dòng)力導(dǎo)數(shù)代入建立的操縱運(yùn)動(dòng)仿真模型，對(duì)船舶的回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)進(jìn)行了仿真，并對(duì)水動(dòng)力導(dǎo)數(shù)的敏感度進(jìn)行了分析，主要得到如下結(jié)論：

1） 基于RNAS的數(shù)值模擬方法可以較好地預(yù)報(bào)DTMB 5415船的線性水動(dòng)力導(dǎo)數(shù)，數(shù)值模擬得到的水動(dòng)力曲線與試驗(yàn)吻合較好，通過回歸方法得到的水動(dòng)力導(dǎo)數(shù)與試驗(yàn)相比誤差大部分小于10%。

2） 所提預(yù)報(bào)方法能夠快速且較準(zhǔn)確地預(yù)報(bào)DTMB 5415船的操縱性，回轉(zhuǎn)和Z形操縱運(yùn)動(dòng)的仿真軌跡與試驗(yàn)吻合較好，操縱性參數(shù)的誤差大部分小于10%。

3） 由水動(dòng)力導(dǎo)數(shù)敏感度分析發(fā)現(xiàn)，水動(dòng)力導(dǎo)數(shù)對(duì)不同操縱運(yùn)動(dòng)的敏感度各不相同，在大部分情況下，線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響較大，非線性水動(dòng)力導(dǎo)數(shù)對(duì)操縱運(yùn)動(dòng)的影響較小，驗(yàn)證了本文所提預(yù)報(bào)方法是有理論依據(jù)的。

本文基于CFD技術(shù)和經(jīng)驗(yàn)公式，結(jié)合MMG船舶操縱運(yùn)動(dòng)方程，對(duì)DTMB 5415船的操縱運(yùn)動(dòng)進(jìn)行了完整的預(yù)報(bào)，為在船舶設(shè)計(jì)階段快速預(yù)報(bào)操縱性提供了一種有效的方案，也為船舶操縱運(yùn)動(dòng)預(yù)報(bào)提供了新的思路。雖然采用經(jīng)驗(yàn)公式預(yù)報(bào)非線性水動(dòng)力導(dǎo)數(shù)能較快地獲得具有一定精度的非線性水動(dòng)力導(dǎo)數(shù)，但由于經(jīng)驗(yàn)公式對(duì)船型的依賴性強(qiáng)，有可能產(chǎn)生誤差，故需進(jìn)一步針對(duì)更多的船型開展研究，考慮不同操縱運(yùn)動(dòng)所必需的水動(dòng)力導(dǎo)數(shù)，進(jìn)而有針對(duì)性地開展CFD約束模數(shù)值模擬。