水面船形狀因子的CFD計算研究

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

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

0 引 言

阻力性能是船舶最重要的水動力性能之一，阻力預報則是船舶設計過程中進行航速預報時必不可少的工作?；谀Ｐ统叨鹊乃卦囼灮蚰Ｐ驮囼灲Y(jié)果預報實船阻力，是目前主要的技術手段，而將模型尺度的船模阻力換算成實船阻力的具體方法有多種，目前常用的主要包括二因次法和三因次法[1]。

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

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

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

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

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

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

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

1 數(shù)值計算方法

1.1 CFD求解器簡介

船舶水動力學CFD 求解器NaViiX（Naval Hydrodynamics Oriented CFD Solvers），由中國船舶科學研究中心獨立自主研發(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）支持標準k-ε、RNGk-ε、k-ω和SSTk-ω湍流模型；（4）支持慣性坐標系、非慣性坐標系和多參考坐標系求解；（5）支持MPI并行計算。

通過大量的CFD 計算研究[6]并結(jié)合以前的實踐經(jīng)驗[7]，本文的數(shù)值計算與分析，采用基于RANS方程（Reynolds Averaged Navier-Stokes Equations）的CFD方法，其中湍流模型使用RNGk-ε模型?？刂品匠叹唧w形式詳見文獻[6-7]。

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

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

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

（1）非正交修正

非正交修正具體是針對控制方程擴散項的離散。輸運方程中，擴散項表達式如下：

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

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

（2）梯度限制器

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

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

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

式中，

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

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

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

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

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

計算網(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 計算區(qū)域及網(wǎng)格示意圖Fig.2 Computation domain and meshes

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

2 研究對象與計算工況

2.1 研究對象

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

CFD 計算研究的三型船模主尺度參數(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 計算工況

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

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

對三型船模，都采用三套網(wǎng)格來進行數(shù)值預報，網(wǎng)格在軸向、周向和徑向三個方向上按 2 的加細比細化：對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)格細化過程中，近船模表面第一層網(wǎng)格法向高度保持不變（y+平均約為45）。

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

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

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

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

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

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

3.1 計算結(jié)果

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

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

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

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

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

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

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

（1）驗證流程

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

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

由于0

修正因子CG的計算如下，其中按照文獻推薦Pth=2：

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

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

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

（2）確認流程

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

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

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

同樣以KVLCC2設計航速計算工況為例，比較誤差和確認不確定度的計算為

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

由于 |E|

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

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

4 結(jié) 語

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

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