王 星，吳乘勝，趙 峰，王建春，金奕星

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

0 引 言

阻力性能是船舶最重要的水動力性能之一，阻力預(yù)報(bào)則是船舶設(shè)計(jì)過程中進(jìn)行航速預(yù)報(bào)時(shí)必不可少的工作。基于模型尺度的水池試驗(yàn)或模型試驗(yàn)結(jié)果預(yù)報(bào)實(shí)船阻力，是目前主要的技術(shù)手段，而將模型尺度的船模阻力換算成實(shí)船阻力的具體方法有多種，目前常用的主要包括二因次法和三因次法[1]。

最早的體系性的實(shí)船阻力換算方法是傅汝德提出的二因次法，該方法將船體總阻力分成摩擦阻力和剩余阻力，并認(rèn)為摩擦阻力僅與雷諾數(shù)（Re）有關(guān)，剩余阻力僅與傅汝德數(shù)（Fr）有關(guān)，且二者互不干擾。

由于二因次法無法考慮船體的三維形狀與平板的差異所產(chǎn)生的對粘性阻力的影響，在其基礎(chǔ)上，休斯提出將船體阻力分成與Re相關(guān)的粘性阻力和與Fr相關(guān)的興波阻力，并且認(rèn)為粘性阻力中的粘壓阻力系數(shù)與摩擦阻力系數(shù)的比值是常數(shù)k，1+k就是形狀因子，與船體形狀有關(guān)。這種實(shí)船阻力換算方法中引入了形狀因子，即為三因次法。由于三因次法考慮了船體三維形狀的影響，理論上更為合理，目前廣泛應(yīng)用于各大水池的船模阻力試驗(yàn)預(yù)報(bào)中。

使用三因次法換算實(shí)船阻力，關(guān)鍵問題是如何確定形狀因子1+k。傳統(tǒng)的形狀因子確定方法主要基于船模低速拖曳試驗(yàn)，采用Prohaska 方法或ITTC’1978 方法獲得[1]（這兩種方法分別將興波阻力系數(shù)表示為正比于Fr的4次和n次方的形式）。然而，通過水池模型試驗(yàn)獲得形狀因子并非易事，因?yàn)榈退偻弦窢顟B(tài)下，船模阻力較小，易受各種因素的干擾，導(dǎo)致測量結(jié)果相對誤差較大，在給模型試驗(yàn)帶來一定難度的同時(shí)，也會影響形狀因子計(jì)算的準(zhǔn)確度。

近年來，CFD應(yīng)用技術(shù)正高密度融入預(yù)報(bào)-評估-設(shè)計(jì)，為船舶水動力學(xué)性能研究和船型設(shè)計(jì)的創(chuàng)新發(fā)展提供了前所未有的利器。而在形狀因子的計(jì)算和研究中，CFD技術(shù)有著其獨(dú)特的優(yōu)點(diǎn)：一是可以采用疊模計(jì)算，可以完全排除興波的影響；二是可以在較大的航速范圍內(nèi)開展計(jì)算，研究Re數(shù)等對形狀因子的影響（模型試驗(yàn)由于必須在低速工況下開展，因而難以實(shí)現(xiàn)）。

國內(nèi)外不少研究人員都利用CFD 技術(shù)開展了形狀因子研究。Kouh 等[2]應(yīng)用CFD 技術(shù)，針對Wig?ley船、KVLCC2、DTMB5415、DTRC4621、一艘漁船和一條魚雷，開展了船舶形狀因子的尺度效應(yīng)研究，結(jié)果表明形狀因子隨Re大致呈線性增大的趨勢；Shen 等[3]針對細(xì)長型回轉(zhuǎn)體（長徑比11～13），運(yùn)用數(shù)值方法，研究并推導(dǎo)了由模型剩余阻力換算實(shí)船/艇阻力的外推公式，并提出模型尺度下的形狀因子數(shù)值與實(shí)船尺度下的呈現(xiàn)明顯差異；王金寶等[4]研究了肥大型船舶形狀因子的CFD 計(jì)算方法，通過對多艘肥大型船舶形狀因子的CFD 計(jì)算并與模型試驗(yàn)結(jié)果的比較，以及典型船舶航速預(yù)報(bào)與實(shí)船測試結(jié)果的比較，表明采用CFD 計(jì)算獲得的形狀因子（1+k）可以用于肥大型船舶實(shí)船阻力的三因次換算；江杰等[5]則應(yīng)用CFD 技術(shù)，針對DTMB5145和KCS，研究了航行姿態(tài)及尺度效應(yīng)對形狀因子的影響，結(jié)果表明二者均會導(dǎo)致形狀因子增大。

當(dāng)前，CFD與EFD結(jié)合用于船舶航行性能研究和預(yù)報(bào)，已成為趨勢和共識的技術(shù)手段。通過CFD計(jì)算水面船形狀因子并用于實(shí)船阻力預(yù)報(bào)，便是一種典型的應(yīng)用。由于形狀因子中的k是個(gè)小量（通常在0.1～0.2左右），對CFD求解器和計(jì)算方法的精度、可靠性和分辨力等都有很高的要求。

本文介紹了針對水面船形狀因子所開展的CFD計(jì)算研究。由于計(jì)算域中船舶艏部和艉部附近不可避免地會存在低正交性、高扭曲率的低質(zhì)量網(wǎng)格，對計(jì)算結(jié)果影響很大，如處理不當(dāng)，會對計(jì)算精度和穩(wěn)定性產(chǎn)生不利影響。為解決這一問題，在CFD 求解器中引入了非正交修正和梯度限制器。本文針對KVLCC2、KCS 和“育鵬”輪三種典型船型，開展了水面船形狀因子CFD 計(jì)算研究，并對計(jì)算結(jié)果進(jìn)行了不確定度分析。結(jié)果表明，形狀因子CFD 計(jì)算結(jié)果以較高水平通過了不確定度分析的驗(yàn)證和確認(rèn)，說明CFD求解器較好地解決了關(guān)鍵局部區(qū)域低質(zhì)量網(wǎng)格帶來的問題，能夠以較高的精度計(jì)算水面船形狀因子。

1 數(shù)值計(jì)算方法

1.1 CFD求解器簡介

船舶水動力學(xué)CFD 求解器NaViiX（Naval Hydrodynamics Oriented CFD Solvers），由中國船舶科學(xué)研究中心獨(dú)立自主研發(fā)，具有完全自主知識產(chǎn)權(quán)，目前具備以下功能：（1）能夠?qū)崿F(xiàn)三維航行體單相、兩相（帶自由面）湍流繞流CFD 模擬；（2）支持結(jié)構(gòu)化網(wǎng)格、非結(jié)構(gòu)化網(wǎng)格、混合網(wǎng)格、交界面網(wǎng)格和滑移網(wǎng)格；（3）支持標(biāo)準(zhǔn)k-ε、RNGk-ε、k-ω和SSTk-ω湍流模型；（4）支持慣性坐標(biāo)系、非慣性坐標(biāo)系和多參考坐標(biāo)系求解；（5）支持MPI并行計(jì)算。

通過大量的CFD 計(jì)算研究[6]并結(jié)合以前的實(shí)踐經(jīng)驗(yàn)[7]，本文的數(shù)值計(jì)算與分析，采用基于RANS方程（Reynolds Averaged Navier-Stokes Equations）的CFD方法，其中湍流模型使用RNGk-ε模型。控制方程具體形式詳見文獻(xiàn)[6-7]。

CFD 求解器采用有限體積法（Finite Volume Method，F(xiàn)VM）離散控制方程，其中對流項(xiàng)采用二階迎風(fēng)差分格式，擴(kuò)散項(xiàng)采用中心差分格式，壓力速度耦合采用SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法解耦，代數(shù)方程組使用Gauss-Seidel 迭代求解，并使用代數(shù)多重網(wǎng)格（Multigrid）技術(shù)加速收斂。

前面提到，由于船舶形狀因子中的k是個(gè)小量（通常在0.1～0.2左右），對CFD求解器和計(jì)算方法的精度、可靠性和分辨力等要求很高，因而常常需要使用正交性高、扭曲率低的高質(zhì)量結(jié)構(gòu)化網(wǎng)格。但水面船幾何外形是復(fù)雜的三維曲面，在生成結(jié)構(gòu)化網(wǎng)格的過程中，難免會存在高扭曲率、低正交性的低質(zhì)量網(wǎng)格，這些網(wǎng)格通常出現(xiàn)在船舶艏部和艉部附近。而船舶艏部和艉部通常曲率變化大，空間流動復(fù)雜，速度、壓力等物理量時(shí)空變化劇烈，且對計(jì)算結(jié)果影響巨大甚至是決定性的，如果處理不當(dāng)，不但會影響計(jì)算精度，還會影響計(jì)算穩(wěn)定性，甚至?xí)斐捎?jì)算發(fā)散。

為解決上述問題，在CFD求解器中引入了非正交修正和梯度限制器，以下分別進(jìn)行簡要介紹。

（1）非正交修正

非正交修正具體是針對控制方程擴(kuò)散項(xiàng)的離散。輸運(yùn)方程中，擴(kuò)散項(xiàng)表達(dá)式如下：

式中，γ是控制體單元界面的擴(kuò)散率，S則為界面的面外法向矢量。擴(kuò)散項(xiàng)離散后可分為正交部分和非正交部分[8]，見式（2）。對于網(wǎng)格正交性較好的部分，第二部分可直接忽略，但當(dāng)網(wǎng)格單元正交性低、扭曲率高時(shí)，加入非正交部分可有效保證計(jì)算精度和穩(wěn)定性。

式中，連接向量dNP=xN-xP，xP為控制體單元中心坐標(biāo)，xN則是控制體單元界面f所對應(yīng)的相鄰單元中心坐標(biāo)可利用控制體單元與相鄰單元的變量梯度值通過線性插值得到。

（2）梯度限制器

在同位網(wǎng)格中，流場變量值存儲在控制體單元中心，而單元界面的變量值φf,P則需要通過控制體單元φP重構(gòu)獲得，表達(dá)式為

式中，連接向量d=xf-xP，xP和xf分別是控制體單元中心和界面中心坐標(biāo)，( )?φP為單元中心處變量梯度。

在重構(gòu)過程中，為了避免出現(xiàn)新的極大值和極小值，需要引入梯度限制器。常見的梯度限制器有Barth and Jespersen限制器和Venkatakrishnan限制器[9]。本求解器使用改進(jìn)的Venkatakrishnan限制器，該限制器在精細(xì)流動以及非定常流動模擬中具有更小的數(shù)值耗散。具體表達(dá)式為

式中，

圖1 給出了針對一典型算例——CFD 求解器使用非正交修正和梯度限制器前后阻力收斂歷程的對比。從圖中可以看出：在使用前（左圖），整個(gè)計(jì)算過程中壓差阻力收斂歷程呈不規(guī)則振蕩，且振蕩幅度很大，摩擦阻力收斂歷程也存在一定振蕩，不過幅度較小而已；在使用后（右圖），壓差阻力和摩擦阻力都是在經(jīng)歷計(jì)算初始階段的振蕩后，很快趨于平穩(wěn)收斂的狀態(tài)。

圖1 阻力收斂歷程比較Fig.1 Comparison of convergence histories of resistance

可見，CFD求解器在引入非正交修正和梯度限制器之后，計(jì)算更加穩(wěn)定，同時(shí)收斂速度也有所提高。

1.2 計(jì)算區(qū)域、網(wǎng)格劃分與邊界條件

為排除興波的干擾，CFD計(jì)算采用重疊模，即靜水面作為對稱面處理；同時(shí)由于船舶左右對稱，且流動是定?；驕?zhǔn)定常的，同樣具有對稱性，CFD 模擬時(shí)只需計(jì)算一半?yún)^(qū)域（如圖2 所示）。計(jì)算區(qū)域范圍參考了相關(guān)文獻(xiàn)[7]的研究成果，其邊界如下：（a）前端——模型首部前約1.2 倍船長處，設(shè)置為入口邊界條件；（b）后端——模型尾部后約2.4 倍船長處，設(shè)置為出口邊界條件；（c）外邊界——模型外側(cè)約1.2 倍船長處，設(shè)置為入口邊界條件；（d）對稱面——模型中縱剖面的延展面和靜水面，設(shè)置為對稱面邊界條件；（e）船舶模型表面——設(shè)置為無滑移壁面邊界條件。

計(jì)算網(wǎng)格為多塊分區(qū)的H-O 型結(jié)構(gòu)化網(wǎng)格（縱向H 型、橫向O 型，如圖2 所示），網(wǎng)格劃分的基本原則為：船模艏部和艉部等曲率變化較大部分網(wǎng)格加密、舯部網(wǎng)格相對較為稀疏；模型近壁面附近網(wǎng)格加密，其中近壁面第一層網(wǎng)格法向高度根據(jù)y+確定（y+平均約為45）。

圖2 計(jì)算區(qū)域及網(wǎng)格示意圖Fig.2 Computation domain and meshes

CFD 計(jì)算中，邊界條件的具體設(shè)置如下：（a）在入口邊界上，根據(jù)船模前進(jìn)速度，給定入口流動速度；（b）出口邊界距離船模足夠遠(yuǎn)，其壓力分布設(shè)置為靜水壓力；（c）在船模表面，引入標(biāo)準(zhǔn)壁面函數(shù)；（d）在對稱面上，滿足對稱條件。

2 研究對象與計(jì)算工況

2.1 研究對象

CFD計(jì)算研究針對三型船模開展，分別為KVLCC2、KCS和“育鵬”輪。這三型船模很具有典型性：KVLCC2 是由韓國船舶與海洋工程研究所（Korean Research Institute for Ship and Ocean Engineering，KRISO）設(shè)計(jì)的30 萬噸VLCC 船型（實(shí)船并未建造），共有兩種線型，目前國際通常采用的是第二種線型，即KVLCC2，是典型的低速大方形系數(shù)船型；KCS 是一艘?guī)虮囚嫉募b箱船（實(shí)船并未建造），也由KRISO 設(shè)計(jì)，是典型的中高速中等方形系數(shù)船型；“育鵬”輪是一艘多用途教學(xué)實(shí)習(xí)船，航速和方形系數(shù)介于KVLCC2和KCS之間，也代表了近年來集裝箱船低速化的一種趨勢。

CFD 計(jì)算研究的三型船模主尺度參數(shù)見表1，圖3 給出了KVLCC2、KCS 和“育鵬”輪船模幾何外形。

表1 KVLCC2、KCS和“育鵬”輪主尺度參數(shù)Tab.1 Main particulars of KVLCC2,KCS and YUPENG

圖3 船模幾何外形Fig.3 Geometry of ship models

2.2 計(jì)算工況

CFD計(jì)算分別針對三型船模設(shè)計(jì)吃水狀態(tài)下的兩個(gè)航速（設(shè)計(jì)航速和較低航速）開展，計(jì)算工況列于表2中。

表2 船模形狀因子CFD計(jì)算工況Tab.2 Computational cases for form factor of surface ships

對三型船模，都采用三套網(wǎng)格來進(jìn)行數(shù)值預(yù)報(bào)，網(wǎng)格在軸向、周向和徑向三個(gè)方向上按 2 的加細(xì)比細(xì)化：對KVLCC2船模，三套網(wǎng)格的網(wǎng)格數(shù)分別為35萬、100萬和260萬；對KCS船模，三套網(wǎng)格的網(wǎng)格數(shù)分別為31萬、90萬和240萬；對“育鵬”輪船模，三套網(wǎng)格的網(wǎng)格數(shù)分別為33萬、91萬和240萬；網(wǎng)格細(xì)化過程中，近船模表面第一層網(wǎng)格法向高度保持不變（y+平均約為45）。

表3給出了三型船模三套網(wǎng)格中高質(zhì)量網(wǎng)格（扭曲率0.3以下）和低質(zhì)量網(wǎng)格（扭曲率0.7以上）的比例。網(wǎng)格扭曲率定義如下（以四邊形面單元為例）：

表3 計(jì)算域網(wǎng)格單元扭曲率Tab.3 Skewness coefficients of grid cells in computational domain

式中，θ為網(wǎng)格單元的最大或最小內(nèi)角。

從表3中可以看出：計(jì)算域中高質(zhì)量網(wǎng)格占大多數(shù)，比例在85%左右；低質(zhì)量網(wǎng)格比例很小，占比不超過3%，且大多在2%以下。高扭曲率的低質(zhì)量網(wǎng)格比例雖然不大，但絕大部分位于船模艏部和艉部附近這些關(guān)鍵區(qū)域，如圖4 所示（圖中紅色區(qū)域表示網(wǎng)格單元扭曲率高），對數(shù)值計(jì)算的影響很大，也給船模形狀因子的準(zhǔn)確、穩(wěn)定計(jì)算帶來了很大的挑戰(zhàn)。

圖4 船體附近網(wǎng)格單元扭曲率分布Fig.4 Distribution of skewness coefficients of grid cells in the vicinity of ship hull

3 計(jì)算結(jié)果與分析

3.1 計(jì)算結(jié)果

通過疊模CFD計(jì)算獲取水面船形狀因子，其定義為

表4給出了KVLCC2、KCS和“育鵬”輪船模在設(shè)計(jì)航速和較低航速下的阻力系數(shù)和形狀因子計(jì)算結(jié)果，表中同時(shí)給出了設(shè)計(jì)航速下的形狀因子模型試驗(yàn)結(jié)果[10-11]。圖5給出了三條船模在不同航速下形狀因子隨網(wǎng)格數(shù)量的變化曲線，圖中同樣也給出了設(shè)計(jì)航速下的形狀因子模型試驗(yàn)結(jié)果（圖中虛線所示）。

表4 船模阻力及形狀因子計(jì)算結(jié)果Tab.4 Computational results of resistance and form factor

圖5 船模形狀因子計(jì)算結(jié)果隨網(wǎng)格數(shù)量變化Fig.5 Computational results of form factor varying with grid cell number

由表4和圖5可見，三條船模的總阻力系數(shù)和形狀因子計(jì)算結(jié)果隨網(wǎng)格數(shù)的增加逐漸減小并趨于收斂——細(xì)網(wǎng)格和中網(wǎng)格計(jì)算結(jié)果之間的差異基本都在1%以內(nèi)，且形狀因子CFD計(jì)算結(jié)果收斂于模型試驗(yàn)結(jié)果。

3.2 計(jì)算結(jié)果不確定度分析

根據(jù)ITTC 推薦的CFD 不確定度分析規(guī)程[12]及相關(guān)指南[13-14]，對形狀因子計(jì)算結(jié)果進(jìn)行不確定度分析，包括驗(yàn)證流程和確認(rèn)流程。在分析過程中，定義下標(biāo)“1”代表細(xì)網(wǎng)格，下標(biāo)“2”代表中網(wǎng)格，下標(biāo)“3”代表粗網(wǎng)格。

（1）驗(yàn)證流程

以KVLCC2設(shè)計(jì)航速工況為例，相鄰兩套網(wǎng)格對應(yīng)的形狀因子計(jì)算結(jié)果之差為

網(wǎng)格收斂因子為

由于0

修正因子CG的計(jì)算如下，其中按照文獻(xiàn)推薦Pth=2：

CG= 4.800顯著大于1，則不確定度UG為

類似地，可以進(jìn)行其他計(jì)算工況和其他船模計(jì)算結(jié)果的驗(yàn)證。KVLCC2、KCS 和“育鵬”輪船模形狀因子計(jì)算結(jié)果的驗(yàn)證列于表5 中。從表中可以看出，三條船模形狀因子計(jì)算結(jié)果的驗(yàn)證水平比較高，最大也不過在1+k的2%左右。

表5 形狀因子計(jì)算結(jié)果的驗(yàn)證Tab.5 Verification of computational results of form factor

（2）確認(rèn)流程

形狀因子CFD 計(jì)算結(jié)果的高水平確認(rèn)并非易事，原因不在于CFD 計(jì)算而在于模型試驗(yàn)。前面說過，通過水池模型試驗(yàn)獲得形狀因子并不容易，測量結(jié)果相對誤差較大，也就是試驗(yàn)結(jié)果的不確定度較大，一般情況下也很少給出形狀因子試驗(yàn)結(jié)果的不確定度。

為了進(jìn)行CFD 計(jì)算結(jié)果的確認(rèn)，必須要有模型試驗(yàn)結(jié)果的不確定度。為此，這里根據(jù)式（1）形狀因子的定義，由船模（疊模）總阻力模型試驗(yàn)結(jié)果的不確定度，推算形狀因子的不確定度。

假設(shè)船模（疊模）總阻力系數(shù)的不確定度為1%，則根據(jù)式（7），形狀因子模型試驗(yàn)結(jié)果的不確定度為（1+k）%。需要說明的是：對于船模（疊模）阻力模型試驗(yàn)而言，1%的不確定度是正常水平；但形狀因子通常是由船模低速拖曳試驗(yàn)確定的，由于此類試驗(yàn)測量結(jié)果相對誤差較大，總阻力1%的不確定度水平可能難以達(dá)到，從而導(dǎo)致形狀因子的不確定度較大。因此，對形狀因子CFD 計(jì)算結(jié)果而言，確認(rèn)的標(biāo)準(zhǔn)偏于嚴(yán)格。

同樣以KVLCC2設(shè)計(jì)航速計(jì)算工況為例，比較誤差和確認(rèn)不確定度的計(jì)算為

式中，D為模型試驗(yàn)結(jié)果，對于本算例，D=1.219。

由于 |E|

類似地，可以進(jìn)行其他船模計(jì)算結(jié)果的確認(rèn)。需要說明的是，由于沒有較低航速工況下的形狀因子模型試驗(yàn)結(jié)果，無法進(jìn)行對應(yīng)工況下CFD 計(jì)算結(jié)果的確認(rèn)，因此這里僅給出KVLCC2、KCS 和“育鵬”輪船模設(shè)計(jì)航速工況下計(jì)算結(jié)果的確認(rèn)，結(jié)果列于表6 中。從表中可以看出，三條船模形狀因子CFD計(jì)算結(jié)果的比較誤差均小于確認(rèn)不確定度，也就是都通過了UV水平的確認(rèn)。

表6 形狀因子計(jì)算結(jié)果的確認(rèn)Tab.6 Validation of computational results of form factor

4 結(jié) 語

通過對三種典型船型形狀因子的CFD 計(jì)算研究與分析可知：船模形狀因子CFD 計(jì)算結(jié)果隨網(wǎng)格數(shù)增加逐漸減小，并收斂于模型試驗(yàn)結(jié)果；船模形狀因子CFD計(jì)算結(jié)果以較高水平通過了不確定度分析的驗(yàn)證和確認(rèn)。

綜上，通過本文的研究工作可見，通過引入非正交修正和梯度限制器，自主研發(fā)的CFD 求解器較好地解決了船模艏部和艉部附近這些關(guān)鍵區(qū)域低質(zhì)量網(wǎng)格給數(shù)值計(jì)算帶來的問題，能夠以較高的精度計(jì)算水面船形狀因子，可以服務(wù)于船舶快速性性能的研究和預(yù)報(bào)。