李迎華，吳寶山，張 華

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

1 引 言

對(duì)于水下滑翔機(jī)[1-2]這一類以剩余浮力為主要驅(qū)動(dòng)方式的水下低速航行體而言，其在進(jìn)行縱傾與浮力調(diào)節(jié)時(shí)，所作的運(yùn)動(dòng)是非定常的，因此對(duì)其非定常操縱運(yùn)動(dòng)的預(yù)報(bào)研究，是指導(dǎo)其設(shè)計(jì)開(kāi)發(fā)的基礎(chǔ)。

由于水下低速航行體非定常操縱運(yùn)動(dòng)中其雷諾數(shù)低于臨界雷諾數(shù)，水動(dòng)力系數(shù)會(huì)隨時(shí)間發(fā)生變化，因而采用傳統(tǒng)的拘束模型試驗(yàn)[3-5]，以及基于準(zhǔn)定常假設(shè)的數(shù)值方法[6]對(duì)其進(jìn)行運(yùn)動(dòng)預(yù)報(bào)時(shí)，則需要對(duì)各種雷諾數(shù)下的水動(dòng)力分別進(jìn)行測(cè)試和計(jì)算，這顯然需要大量的工作，而且對(duì)于拘束模型試驗(yàn)來(lái)說(shuō)，其代價(jià)也是極其昂貴的。近年來(lái)，隨著商業(yè)軟件Fluent的不斷發(fā)展，基于RANS求解器的CFD數(shù)值預(yù)報(bào)方法已成為一種經(jīng)濟(jì)、有效的手段。因此，研究非定常操縱運(yùn)動(dòng)的CFD數(shù)值模擬技術(shù)對(duì)水下低速航行體設(shè)計(jì)研究具有重要的意義。

對(duì)于數(shù)值模擬時(shí)水動(dòng)力及運(yùn)動(dòng)參數(shù)需要實(shí)時(shí)求解的這類操縱運(yùn)動(dòng)，F(xiàn)luent專門提供了一種研究手段—?jiǎng)泳W(wǎng)格技術(shù)。在模擬時(shí)將運(yùn)動(dòng)方程通過(guò)用戶自定義函數(shù)（UDF）導(dǎo)入至Fluent中，即可實(shí)現(xiàn)水動(dòng)力和運(yùn)動(dòng)參數(shù)實(shí)時(shí)自主求解。因此，動(dòng)網(wǎng)格技術(shù)適合于水下低速航行體非定常操縱運(yùn)動(dòng)的直接模擬。

目前動(dòng)網(wǎng)格技術(shù)大多應(yīng)用于周期性操縱運(yùn)動(dòng)的數(shù)值計(jì)算，如垂蕩、升沉等，所采用的動(dòng)網(wǎng)格劃分形式一般是全非結(jié)構(gòu)網(wǎng)格[7]，由于這種網(wǎng)格劃分形式在近壁面粘性邊界層內(nèi)難以生成大展弦比的非結(jié)構(gòu)網(wǎng)格，且網(wǎng)格數(shù)量較結(jié)構(gòu)網(wǎng)格多出許多，因而將其推廣到非周期性操縱運(yùn)動(dòng)預(yù)報(bào)上有一定的局限性，因此，近年來(lái)，結(jié)合結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格優(yōu)勢(shì)的混合網(wǎng)格技術(shù)得到了重視和發(fā)展，張來(lái)平等[8]提出了三棱柱/四面體/矩形三維動(dòng)態(tài)混合網(wǎng)格劃分形式，并用無(wú)粘流的圓球超聲速繞流算例進(jìn)行了驗(yàn)證，得到了與定常計(jì)算相一致的結(jié)果。

本文在此基礎(chǔ)上對(duì)三維動(dòng)態(tài)混合網(wǎng)格模型進(jìn)行了發(fā)展，提出了六面體/四面體三維動(dòng)態(tài)混合網(wǎng)格模型，并以此建立了基于動(dòng)網(wǎng)格技術(shù)的CFD數(shù)值預(yù)報(bào)方法。本文首先以帶攻角均勻直航的滑行狀態(tài)為研究對(duì)象，應(yīng)用動(dòng)網(wǎng)格技術(shù)對(duì)一典型水下低速航行體縮比模型主體及帶附體的水動(dòng)力進(jìn)行了計(jì)算以驗(yàn)證其水動(dòng)力的計(jì)算精度，為便于比較，同時(shí)給出了靜態(tài)網(wǎng)格計(jì)算結(jié)果，進(jìn)而將其應(yīng)用于一水下滑翔機(jī)自航模非定常滑翔運(yùn)動(dòng)的預(yù)報(bào)，以驗(yàn)證動(dòng)網(wǎng)格技術(shù)的應(yīng)用能力。

2 基于動(dòng)網(wǎng)格技術(shù)的數(shù)值計(jì)算方法

2.1 控制方程

時(shí)均的不可壓縮連續(xù)性方程和N-S方程（RANS方程）如（1）、（2）式所示：

2.2 湍流模型

對(duì)于本文研究的水下低速航行體而言，其運(yùn)動(dòng)過(guò)程中縱向雷諾數(shù)較低。RNG k-ε湍流模型提供了針對(duì)低雷諾數(shù)的有效粘性的微分解析式，且具備數(shù)值穩(wěn)定性好和處理流線彎曲程度較大的流動(dòng)更優(yōu)等優(yōu)點(diǎn)[9]，本文數(shù)值計(jì)算即采用RNG k-ε湍流模型。

2.3 幾何守恒律

動(dòng)網(wǎng)格需滿足幾何守恒律，控制體積的時(shí)間導(dǎo)數(shù)由（3）式計(jì)算：

式中：nf是控制面的數(shù)目，是j面的面積矢量。

式中：δVj是控制面j在一個(gè)時(shí)間步長(zhǎng)內(nèi)的體積更新量。

2.4 邊界條件

（1）遠(yuǎn)場(chǎng)邊界條件：采用零剪應(yīng)力壁面邊界條件處理。

（2）物面條件：滿足壁面黏附條件，壁面處流體速度與運(yùn)動(dòng)邊界速度相同。

2.5 動(dòng)態(tài)混合網(wǎng)格生成方法

本文的六面體/四面體動(dòng)態(tài)混合網(wǎng)格模型如下建立：在近壁面處依次生成貼體六面體網(wǎng)格和四面體網(wǎng)格共同構(gòu)成運(yùn)動(dòng)區(qū)域隨運(yùn)動(dòng)邊界一起運(yùn)動(dòng) （四面體網(wǎng)格與六面體網(wǎng)格之間采用交界面過(guò)渡），外場(chǎng)生成四面體網(wǎng)格并時(shí)刻保持靜止，在運(yùn)動(dòng)區(qū)域和靜止區(qū)域之間采用四面體網(wǎng)格過(guò)渡構(gòu)成變形區(qū)域。當(dāng)物體位移較小時(shí)，采用“類彈簧原理”對(duì)變形區(qū)域網(wǎng)格節(jié)點(diǎn)進(jìn)行松弛，當(dāng)物體位移較大導(dǎo)致變形區(qū)域網(wǎng)格質(zhì)量較差時(shí)，則在局部重新生成網(wǎng)格。

本文采用上述改進(jìn)的動(dòng)態(tài)混合網(wǎng)格模型對(duì)一典型水下低速航行體縮比模型表面及其周圍流動(dòng)區(qū)域進(jìn)行了網(wǎng)格劃分，此模型主體為細(xì)長(zhǎng)回轉(zhuǎn)體，長(zhǎng)2 150mm，直徑300mm，在主體中后部布置有一對(duì)后掠的水平翼，翼兩端間距1 200mm。網(wǎng)格劃分的具體形式為：主體C-O型（縱向-周向），附體C-O型（展向—弦向）；計(jì)算域大小為：沿航行體首前端、航行體徑向分別延伸3L距離，沿航行體尾后端延伸4L距離（L為航行體縱向長(zhǎng)度）。單獨(dú)主體網(wǎng)格數(shù)65萬(wàn)，帶附體網(wǎng)格數(shù)120萬(wàn)。圖1、2為帶附體表面及近壁面流場(chǎng)區(qū)域網(wǎng)格劃分示意圖。為便于對(duì)比分析，也給出了靜態(tài)網(wǎng)格劃分示意圖，如圖3所示。單獨(dú)主體的靜態(tài)網(wǎng)格數(shù)量32萬(wàn)，帶附體的靜態(tài)網(wǎng)格數(shù)量57萬(wàn)。

2.6 動(dòng)網(wǎng)格更新算法

動(dòng)網(wǎng)格更新算法有彈簧近似光滑算法（Spring Based Smoothing）、動(dòng)態(tài)分層算法（Dynamic Layering）和局部重構(gòu)算法（Local Remeshing）。對(duì)于本文研究的四面體/六面體混合網(wǎng)格，選用彈簧近似光滑算法和局部重構(gòu)算法作為動(dòng)網(wǎng)格更新算法。

2.7 數(shù)值離散和求解

（1）時(shí)間項(xiàng)的離散：采用直接一階隱式離散。

（2）空間項(xiàng)的離散：其中擴(kuò)散項(xiàng)以中心差分格式差分，對(duì)流項(xiàng)采用二階迎風(fēng)格式。

應(yīng)用SIMPLE法處理壓力速度耦合問(wèn)題，離散方程以Gauss-Seidel迭代法求解。

2.8 時(shí)間步長(zhǎng)的選取

由于時(shí)間項(xiàng)采用一階隱式離散，是無(wú)條件穩(wěn)定的，因此時(shí)間步長(zhǎng)的選取可以不用考慮離散格式的穩(wěn)定性，理論上可以取較大的時(shí)間步長(zhǎng)，但過(guò)大的時(shí)間步長(zhǎng)將導(dǎo)致網(wǎng)格更新過(guò)快，進(jìn)而出現(xiàn)負(fù)體積。本文從網(wǎng)格更新角度推導(dǎo)出時(shí)間步長(zhǎng)計(jì)算公式如下：

不考慮控制面上的面積更新量則由（4）式可得對(duì)于單元j上的當(dāng)?shù)貢r(shí)間步長(zhǎng)：

式中：αj=hj/Δhj，Δhj為單元j上的高度更新量，hj為單元j上更新前的高度。

全場(chǎng)時(shí)間步長(zhǎng)應(yīng)該取全動(dòng)區(qū)域最小的即

式中：h1表示動(dòng)區(qū)域最小的第一層網(wǎng)格高度；vel為邊界運(yùn)動(dòng)的合速度；α為一系數(shù)。

2.9 雙振蕩圓板驗(yàn)證算例

為對(duì)動(dòng)網(wǎng)格技術(shù)的可行性進(jìn)行驗(yàn)證，本文選取文獻(xiàn)[10]中的雙強(qiáng)迫振蕩圓板算例進(jìn)行了計(jì)算，并與文獻(xiàn)給出的結(jié)果進(jìn)行了對(duì)比。采用全非結(jié)構(gòu)網(wǎng)格劃分形式對(duì)其進(jìn)行網(wǎng)格劃分，網(wǎng)格數(shù)量12萬(wàn)，如圖4所示。

采用彈簧近似光滑算法和局部重構(gòu)算法作為動(dòng)網(wǎng)格更新算法，將雙板振蕩速度以用戶自定義函數(shù)（UDF）形式導(dǎo)入到Fluent中，時(shí)間步長(zhǎng) Δt=0.1s。

本文計(jì)算得到的阻力系數(shù)變化周期T=25s，阻力系數(shù)幅值A(chǔ)=13.7；而文獻(xiàn)給出的阻力系數(shù)變化周期T=25s，阻力系數(shù)幅值A(chǔ)=14.2，可見(jiàn)本文采用的動(dòng)網(wǎng)格技術(shù)可以得到相對(duì)比較精確的計(jì)算結(jié)果，因此可以應(yīng)用于水下滑翔機(jī)非定常運(yùn)動(dòng)的預(yù)報(bào)。

3 水下滑翔機(jī)非定常運(yùn)動(dòng)的CFD計(jì)算研究

本文以中國(guó)船舶科學(xué)研究中心自主開(kāi)發(fā)的“前哨”號(hào)水下滑翔機(jī)自航模[2]為對(duì)象，應(yīng)用動(dòng)網(wǎng)格技術(shù)對(duì)其非定常操縱運(yùn)動(dòng)進(jìn)行了預(yù)報(bào)，以進(jìn)一步驗(yàn)證動(dòng)網(wǎng)格技術(shù)的應(yīng)用。

操縱運(yùn)動(dòng)預(yù)報(bào)的基礎(chǔ)在于水動(dòng)力的預(yù)報(bào)，因此本文首先以帶攻角均勻直航的滑行狀態(tài)為研究對(duì)象，應(yīng)用動(dòng)網(wǎng)格技術(shù)對(duì)前述典型水下低速航行體縮比模型的主體及帶附體水動(dòng)力的預(yù)報(bào)精度進(jìn)行了驗(yàn)證，為便于比較分析，同時(shí)給出了靜態(tài)網(wǎng)格計(jì)算結(jié)果。為分析水動(dòng)力計(jì)算精度對(duì)運(yùn)動(dòng)預(yù)報(bào)的影響，本文還基于帶附體動(dòng)網(wǎng)格計(jì)算的水動(dòng)力對(duì)定?；柽\(yùn)動(dòng)參數(shù)進(jìn)行了分析。

3.1 基于動(dòng)網(wǎng)格技術(shù)的水動(dòng)力預(yù)報(bào)驗(yàn)證

通過(guò)用戶自定義函數(shù) （UDF）將水下低速航行體運(yùn)動(dòng)速度實(shí)時(shí)導(dǎo)入到Fluent中，時(shí)間步長(zhǎng)Δt=0.002s，應(yīng)用基于動(dòng)態(tài)混合網(wǎng)格模型的動(dòng)網(wǎng)格技術(shù)依次計(jì)算了雷諾數(shù)Re=2.4×106，攻角分別為6°、8°、10°、12°四種工況下主體及帶附體模型的水動(dòng)力，得到了無(wú)因次水動(dòng)力系數(shù)和水動(dòng)力中心作用點(diǎn)位置lα′（lα′=-M′/Z′），并與風(fēng)洞模型試驗(yàn)結(jié)果[12]進(jìn)行了對(duì)比。

（1）主體模型水動(dòng)力計(jì)算結(jié)果

主體模型靜態(tài)網(wǎng)格和動(dòng)網(wǎng)格計(jì)算結(jié)果如表1所示，水動(dòng)力隨攻角變化曲線如圖5-8所示。

表1 水下航行體主體模型水動(dòng)力計(jì)算結(jié)果（Re=2.4×106）Tab.1 The calculated results of the hydrodynamic forces for the main body of the underwater vehicle model（Re=2.4×106）

由計(jì)算結(jié)果可知，除個(gè)別工況外，主體模型靜態(tài)網(wǎng)格計(jì)算的軸向力系數(shù)相差16%以內(nèi)，垂向力系數(shù)相差8%以內(nèi)，俯仰力矩系數(shù)相差5%以內(nèi)，水動(dòng)力作用中心位置相差5%以內(nèi)；動(dòng)網(wǎng)格計(jì)算的軸向力系數(shù)相差5%以內(nèi)，垂向力系數(shù)、俯仰力矩系數(shù)均相差10%以內(nèi)，水動(dòng)力作用中心位置相差13%以內(nèi)。

（2）帶附體模型水動(dòng)力計(jì)算結(jié)果及對(duì)運(yùn)動(dòng)預(yù)報(bào)的影響

帶附體模型操縱性水動(dòng)力計(jì)算結(jié)果如表2所示，水動(dòng)力計(jì)算結(jié)果隨攻角變化曲線如圖9-12所示（由于靜態(tài)網(wǎng)格在12°時(shí)失速，因此只給出前三個(gè)攻角對(duì)應(yīng)的水動(dòng)力值）。

由計(jì)算結(jié)果可知，除個(gè)別工況外，帶附體模型靜態(tài)網(wǎng)格計(jì)算的軸向力系數(shù)相差16%以內(nèi)，垂向力系數(shù)和俯仰力矩系數(shù)相差14%以內(nèi)，水動(dòng)力作用中心位置相差5%以內(nèi)；動(dòng)網(wǎng)格計(jì)算的軸向力系數(shù)相差10%以內(nèi)，垂向力系數(shù)相差9%以內(nèi)，俯仰力矩系數(shù)相差13%以內(nèi)，水動(dòng)力作用中心位置lα′相差10%以內(nèi)。

表2 水下航行體帶附體模型水動(dòng)力計(jì)算結(jié)果（Re=2.4×106）Tab.2 The calculated results of the hydrodynamic forces for the main body with appendages of the underwater vehicle model（Re=2.4×106）

由主體及帶附體模型水動(dòng)力計(jì)算結(jié)果可知，基于動(dòng)態(tài)混合網(wǎng)格模型的動(dòng)網(wǎng)格技術(shù)水動(dòng)力預(yù)報(bào)結(jié)果與試驗(yàn)吻合較好，且可以達(dá)到與靜態(tài)網(wǎng)格相類似的精度。

同時(shí)從計(jì)算結(jié)果也可看出：數(shù)值計(jì)算結(jié)果與試驗(yàn)仍存在一定差距，這在很大程度上是因?yàn)閷?shí)際流動(dòng)處于層流和湍流的混合區(qū)，而本文采用適用于充分發(fā)展湍流（Re>1×107）的湍流模型去解算，因而會(huì)帶來(lái)一定的誤差。

為分析水動(dòng)力計(jì)算精度對(duì)運(yùn)動(dòng)預(yù)報(bào)的影響，本文基于上述動(dòng)網(wǎng)格計(jì)算的水動(dòng)力和試驗(yàn)測(cè)得的水動(dòng)力，分別計(jì)算了給定雷諾數(shù)Re=2.4×106下帶附體模型穩(wěn)定滑翔時(shí)不同剩余浮力系數(shù)k（k=ΔB/B，ΔB為剩余浮力，B為水下靜均衡后的浮力） 對(duì)應(yīng)的潛浮角 χ （χ=θ-α，θ為縱傾角，α 為攻角）、水平速度分量uξ和垂直速度分量uζ，計(jì)算結(jié)果曲線如圖13-15所示。

由計(jì)算結(jié)果來(lái)看：動(dòng)網(wǎng)格預(yù)報(bào)的結(jié)果與試驗(yàn)預(yù)報(bào)結(jié)果相比，潛浮角相差0.4°以內(nèi)，水平速度和垂直速度相差0.007m/s以內(nèi)，表明了上述動(dòng)網(wǎng)格水動(dòng)力計(jì)算誤差對(duì)運(yùn)動(dòng)預(yù)報(bào)的影響比較小。

3.2 基于動(dòng)網(wǎng)格技術(shù)的非定常運(yùn)動(dòng)的CFD預(yù)報(bào)

對(duì)“前哨”號(hào)水下滑翔機(jī)自航模垂直面內(nèi)滑翔運(yùn)動(dòng)的預(yù)報(bào)可采用如圖16所示的兩類坐標(biāo)系：地面坐標(biāo)系E-ξζ和運(yùn)動(dòng)坐標(biāo)系o-xz，其中原點(diǎn)o位于水下航行體浮心處，ox軸指向首部為正，oz軸垂直ox軸指向下方為正。圖中X，Z，M分別為此水下航行體所受的軸向、垂向水動(dòng)力和縱傾力矩。

此水下滑翔機(jī)垂直面內(nèi)非定?；柽\(yùn)動(dòng)方程可簡(jiǎn)化為如（7）～（12）式所示的形式：

式中：W（t）為水下滑翔機(jī)重量；xg（t）為水下滑翔機(jī)重心縱向坐標(biāo)；zg（t）為水下滑翔機(jī)重心垂向坐標(biāo)；B為水下滑翔機(jī)所受浮力。

將滑翔運(yùn)動(dòng)方程（7）～（12）進(jìn)行時(shí)間一階向前差分離散，并把離散后的方程以UDF形式導(dǎo)入到Fluent中，即可實(shí)現(xiàn)Fluent中每一時(shí)間步的水動(dòng)力和運(yùn)動(dòng)參數(shù)實(shí)時(shí)自主求解。其求解流程如圖17所示。

本文采用基于動(dòng)態(tài)混合網(wǎng)格模型的動(dòng)網(wǎng)格技術(shù)對(duì)“前哨”號(hào)水下滑翔機(jī)非定常下潛運(yùn)動(dòng)中的一段進(jìn)行了數(shù)值模擬，結(jié)果如圖18、19所示。

從數(shù)值計(jì)算結(jié)果與試驗(yàn)的比較來(lái)看：基于動(dòng)網(wǎng)格技術(shù)數(shù)值方法的預(yù)報(bào)結(jié)果與試驗(yàn)符合較好，且變化趨勢(shì)一致。

4 結(jié)論和展望

本文基于傳統(tǒng)的動(dòng)態(tài)混合網(wǎng)格模型，發(fā)展了一種改進(jìn)的三維動(dòng)態(tài)混合網(wǎng)格模型，以此建立了非定常操縱運(yùn)動(dòng)數(shù)值預(yù)報(bào)方法，并以帶攻角均勻直航的滑行狀態(tài)為研究對(duì)象，應(yīng)用動(dòng)網(wǎng)格技術(shù)對(duì)一典型水下低速航行體縮比模型主體及帶附體的水動(dòng)力進(jìn)行了計(jì)算，結(jié)果與風(fēng)洞模型試驗(yàn)結(jié)果吻合較好，從而驗(yàn)證了基于動(dòng)態(tài)混合網(wǎng)格模型的動(dòng)網(wǎng)格技術(shù)對(duì)操縱運(yùn)動(dòng)水動(dòng)力的預(yù)報(bào)精度；本文進(jìn)而將其應(yīng)用于一水下滑翔機(jī)自航模非定常運(yùn)動(dòng)的預(yù)報(bào)，預(yù)報(bào)結(jié)果與試驗(yàn)符合較好，初步驗(yàn)證了動(dòng)網(wǎng)格技術(shù)的應(yīng)用能力和本文建立的方法的有效性。同時(shí)動(dòng)網(wǎng)格技術(shù)具有更廣闊的應(yīng)用前景，主要表現(xiàn)在：

（1）動(dòng)網(wǎng)格技術(shù)可以應(yīng)用于水下多體相對(duì)運(yùn)動(dòng)，潛艇的坐底、靠岸等的數(shù)值模擬。

（2）動(dòng)網(wǎng)格技術(shù)可以應(yīng)用于非定常操縱運(yùn)動(dòng)過(guò)程中流體動(dòng)力的精細(xì)研究。

[1]馬冬梅,馬 錚.水下滑翔機(jī)運(yùn)動(dòng)及水動(dòng)力特性研究[R].無(wú)錫:中國(guó)船舶科學(xué)研究中心科技報(bào)告,2007.

[2]李 龍,張 華,陳季軍,馬 錚,吳寶山.水下滑翔機(jī)自航模演示試驗(yàn)研究[C].中國(guó)造船工程學(xué)會(huì)學(xué)術(shù)論文集,2008,No.200804.

[3]Ma Ling,Cui Weicheng.Simulation of dive motion of a deep manned submersible[J].Journal of Ship Mechanics,2004,8(3):31-38.

[4]張 華,吳寶山.扁平型潛水器空間大機(jī)動(dòng)數(shù)學(xué)模型研究[C].中國(guó)造船工程學(xué)會(huì)學(xué)術(shù)論文集,2004,No.200405.

[5]謝俊元,馬 嶺,胡 震.載人深潛器運(yùn)動(dòng)操縱模擬[J].中國(guó)造船,2006,47(2):62-69.

[6]王長(zhǎng)濤,俞建成,吳利紅,封錫盛.水下滑翔機(jī)器人運(yùn)動(dòng)機(jī)理仿真與實(shí)驗(yàn)[J].海洋工程,2007,25(1):64-69.

[7]黃昆侖,龐永杰,蘇玉民,朱 軍.潛器線性水動(dòng)力系數(shù)計(jì)算方法研究[J].船舶力學(xué),2008,12(5):697-703.

[8]張來(lái)平,張涵信.三維動(dòng)態(tài)混合網(wǎng)格技術(shù)及非定常計(jì)算方法[C].全國(guó)流體力學(xué)青年研討會(huì)論文集,2005.

[9]Fluent 6.3.26 User’s Guide.Fluent Inc.[CP].

[10]顧 正.二維單圓柱、雙圓柱繞流問(wèn)題及三維振蕩板運(yùn)動(dòng)的數(shù)值模擬[D].上海:上海交通大學(xué),2007.

[11]施生達(dá).潛艇操縱性[M].北京:國(guó)防工業(yè)出版社,1995.

[12]倪歆韻,潘子英.低速水下航行體風(fēng)洞模型試驗(yàn)[R].無(wú)錫:中國(guó)船舶科學(xué)研究中心科技報(bào)告,2008.