祁江濤，郭海鵬，張起敏，李廣年

(1.中國(guó)船舶科學(xué)研究中心 水動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇 無(wú)錫 214082；2.寧波大學(xué) 海運(yùn)學(xué)院，浙江 寧波 315211；3.寧波大學(xué) 東海戰(zhàn)略研究院，浙江 寧波 315211)

0 引言

針對(duì)實(shí)尺度船舶的水動(dòng)力性能研究，當(dāng)前廣泛采用的技術(shù)手段是利用模型尺度的水池試驗(yàn)或者數(shù)值模擬數(shù)據(jù)外推得到實(shí)船的水動(dòng)力性能。然而，船模尺度與實(shí)船尺度之間往往不滿(mǎn)足某些相似準(zhǔn)則，其中最主要的是表征流體慣性力與粘性力比值的雷諾數(shù)。對(duì)于不同尺度下的船模試驗(yàn)或數(shù)值計(jì)算，不能同時(shí)滿(mǎn)足傅汝德數(shù)與雷諾數(shù)相等，因此在水動(dòng)力外推預(yù)報(bào)時(shí)存在“尺度效應(yīng)”[1]。對(duì)于水下航行體而言，“尺度效應(yīng)”問(wèn)題是影響其水動(dòng)力性能預(yù)報(bào)精度的重要因素，給實(shí)尺度船舶航行性能評(píng)估帶來(lái)一定困擾。隨著計(jì)算流體動(dòng)力學(xué)（computational fluid dynamics,CFD）的發(fā)展，應(yīng)用CFD 方法對(duì)船舶及水下航行體水動(dòng)力的尺度效應(yīng)分析并預(yù)報(bào)實(shí)船水動(dòng)力性能是當(dāng)前船舶水動(dòng)力學(xué)研究領(lǐng)域的熱點(diǎn)之一。

操盛文等[2]對(duì)高雷諾數(shù)下SUBOFF 潛艇的全附體模型繞流進(jìn)行數(shù)值模擬，研究了網(wǎng)格數(shù)量及艇體尺度對(duì)高雷諾數(shù)條件下潛艇阻力計(jì)算結(jié)果的影響。吳方良等[3]對(duì)SUBOFF 潛艇的主艇體模型和全附體模型在不同雷諾數(shù)條件下的流場(chǎng)進(jìn)行了數(shù)值模擬，分析了潛艇粘壓阻力系數(shù)隨雷諾數(shù)變化的規(guī)律。司朝善等[4]以SUBOFF 潛艇全附體模型為對(duì)象，通過(guò)改變潛艇尺度實(shí)現(xiàn)不同雷諾數(shù)下的潛艇直航繞流數(shù)值模擬，分析了模型尺度對(duì)潛艇阻力性能和繞流場(chǎng)的影響規(guī)律。王展智等[5]以DTMB5415 裸船體和全附體模型為研究對(duì)象，計(jì)算了包含實(shí)尺度在內(nèi)的多種縮尺模型的流場(chǎng)，詳細(xì)分析了雙槳水面艦船附體阻力的尺度效應(yīng)。張恒等[6]以KCS 船型作為研究對(duì)象，計(jì)算了不同雷諾數(shù)下的阻力并對(duì)船體尾部流場(chǎng)進(jìn)行分析，探討了改變介質(zhì)的粘性系數(shù)以滿(mǎn)足全相似在實(shí)尺度預(yù)報(bào)方面的可行性。蔡博奧等[7]針對(duì)不同縮尺比的三體船模型進(jìn)行了各阻力成分的數(shù)值預(yù)報(bào)，探討了三體船阻力成分關(guān)于水動(dòng)力系數(shù)的變化規(guī)律并提出了三體船阻力換算的改進(jìn)方法。師超等[8]開(kāi)展了不同船體尺度下的裸船體斜航和定?；剞D(zhuǎn)數(shù)值計(jì)算，探討了不同尺度下船體操縱運(yùn)動(dòng)水動(dòng)力的變化規(guī)律，并結(jié)合船舶操縱運(yùn)動(dòng)數(shù)學(xué)模型分析了尺度效應(yīng)對(duì)操縱性的影響。Dogrul等[9]對(duì)不同尺度下的KCS 和KVLCC2 船的阻力試驗(yàn)進(jìn)行了數(shù)值模擬，探討了不同阻力成分的尺度效應(yīng)問(wèn)題，并據(jù)此分析了不同外推方法在船舶實(shí)尺度阻力預(yù)報(bào)方面的有效性。宋科委等[10]以DTMB5415 船型為研究對(duì)象，預(yù)報(bào)了多種尺度、多種工況及多種船體壁面條件下的阻力性能，對(duì)船體形狀因子、興波阻力及摩擦阻力的尺度效應(yīng)問(wèn)題進(jìn)行了分析，探討了多種實(shí)船阻力性能預(yù)報(bào)方法存在差異的原因。Sezen等[11]以SUBOFF 潛艇模型為對(duì)象，研究了尺度效應(yīng)對(duì)潛艇阻力等水動(dòng)力性能參數(shù)的影響，通過(guò)對(duì)實(shí)尺度潛艇阻力試驗(yàn)的數(shù)值模擬，驗(yàn)證了基于傅汝德外推法的實(shí)尺度阻力預(yù)報(bào)方面的有效性。

從以上研究可以看出，針對(duì)船舶及水下航行體水動(dòng)力性能尺度效應(yīng)的研究大多集中在阻力性能方面，而針對(duì)船舶及水下航行體操縱性能尺度效應(yīng)問(wèn)題的研究還十分有限。本文針對(duì)水下航行體的操縱性能評(píng)估需求，開(kāi)展水下航行體操舵水動(dòng)力性能的尺度效應(yīng)問(wèn)題數(shù)值研究，應(yīng)用CFD 方法開(kāi)展不同尺度的水下航行體直航操舵試驗(yàn)數(shù)值模擬，探討艇體繞流場(chǎng)及水動(dòng)力特性隨雷諾數(shù)的變化規(guī)律，為構(gòu)建實(shí)尺度水下航行體操舵水動(dòng)力性能的外推預(yù)報(bào)方法提供參考。

1 數(shù)值方法

考慮到水下航行體操舵工況下存在的大尺度流動(dòng)分離現(xiàn)象，應(yīng)用分離渦模擬（Detached Eddy Simulation,DES）求解器開(kāi)展數(shù)值計(jì)算。求解器所采用的DES 方法基于剪切應(yīng)力輸運(yùn)（Shear Stress Transfer,SST）k?ω湍流模型進(jìn)行構(gòu)造。SSTk?ω兩方程湍流模型如下：

式中：ρ為流體密度；uj為速度矢量；xj為位置矢量；k和 ω分別為湍流動(dòng)能和湍流比耗散率；μ和 μt分別為層流和湍流黏性系數(shù)。源項(xiàng)P和函數(shù)F1的具體形式以及系數(shù) β?，γ，σk，σω和σω2參考文獻(xiàn)[12]。

將RANS 湍流模型中的長(zhǎng)度尺度采用DES 的長(zhǎng)度尺度代替，即可獲得基于該湍流模型的DES 方法，DES 的長(zhǎng)度尺度表示如下：

其中，RANS 長(zhǎng)度尺度為最近壁面距離，表達(dá)如下：

此外，LES 長(zhǎng)度尺度僅與網(wǎng)格間距有關(guān)，表達(dá)如下：

式中，?為局部網(wǎng)格間距，一般取網(wǎng)格單元與相鄰網(wǎng)格單元中心連線長(zhǎng)度的最大值。

由DES 模型的構(gòu)造可知，當(dāng)RANS 長(zhǎng)度尺度激活時(shí)，以RANS 湍流模型的方式計(jì)算；而當(dāng)LES 長(zhǎng)度尺度激活時(shí)，湍流模型轉(zhuǎn)換為經(jīng)典的Smagorisnky 亞格子模型，從而在求解流場(chǎng)時(shí)實(shí)現(xiàn)RANS 方法與LES 方法的耦合。

2 研究對(duì)象及工況

本文以美國(guó)DARPA 潛艇模型SUBOFF-8 為研究對(duì)象，潛艇幾何形狀和主要尺度如圖1 和表1 所示。針對(duì)不同縮尺比的潛艇模型，開(kāi)展直航工況下的操舵試驗(yàn)數(shù)值模擬，工況見(jiàn)表1。數(shù)值模擬工況滿(mǎn)足傅汝德數(shù)相似，對(duì)應(yīng)傅汝德數(shù)Fr=0.512。

表1 研究對(duì)象及工況Tab.1 Study object and working conditions

圖1 DARPA 潛艇模型SUBOFF-8 模型示意圖Fig.1 Schematic diagram of DARPA submarine model SUBOFF-8 model

主要關(guān)注對(duì)水下航行體操縱性能有顯著影響的縱向力X、橫向力Y及轉(zhuǎn)首力矩N。為便于對(duì)比分析，水動(dòng)力及力矩均采用無(wú)因次形式描述，表達(dá)如下：

3 計(jì)算域與網(wǎng)格生成

計(jì)算域?yàn)殚L(zhǎng)方體，如圖2 所示。坐標(biāo)系原點(diǎn)位于潛艇尾端，計(jì)算域入口（Inlet）邊界距坐標(biāo)系原點(diǎn)2.5 倍艇長(zhǎng)LOA，計(jì)算域出口（Outlet）邊界距坐標(biāo)系原點(diǎn)4.5 倍艇長(zhǎng)LOA，計(jì)算域左右邊界（Left、Right）與上下邊界（Top、Bottom）距坐標(biāo)系原點(diǎn)2 倍艇長(zhǎng)LOA。邊界條件方面，出口邊界設(shè)為壓力出口，其他邊界均設(shè)置為速度入口，艇體表面設(shè)為無(wú)滑移壁面。

圖2 計(jì)算域劃分及邊界條件Fig.2 Division of the computational domain and boundary conditions

計(jì)算域采用非結(jié)構(gòu)六面體網(wǎng)格、棱柱層網(wǎng)格以及表面重構(gòu)的方法進(jìn)行網(wǎng)格劃分。為了保證潛艇附近的網(wǎng)格質(zhì)量，對(duì)其表面進(jìn)行了加密。通過(guò)使用切割體網(wǎng)格單元，在艇體附近及其尾流進(jìn)行網(wǎng)格密化，從而更好地捕捉流動(dòng)細(xì)節(jié)，獲得更準(zhǔn)確的流場(chǎng)信息。通過(guò)調(diào)整棱柱網(wǎng)格層的厚度、層數(shù)及增長(zhǎng)率來(lái)控制近壁面第一層網(wǎng)格高度，使無(wú)因次壁面距離y+值在30～300 之間。圖3 和圖4 給出了潛艇表面和計(jì)算域的網(wǎng)格分布情況。

圖3 潛艇表面網(wǎng)格分布情況Fig.3 Grid distribution on the submarine surface

圖4 計(jì)算域網(wǎng)格分布Fig.4 Grid distribution in the computational domain

4 網(wǎng)格收斂性分析

采用Stern等[13]提出的方法開(kāi)展網(wǎng)格收斂性分析：

式中：SG為該網(wǎng)格下的參數(shù)計(jì)算結(jié)果；εG為相鄰2 套網(wǎng)格對(duì)應(yīng)參數(shù)值的差；RG為網(wǎng)格收斂率。

收斂性分析結(jié)果存在3 種情況：1）01，網(wǎng)格發(fā)散。當(dāng)網(wǎng)格發(fā)散時(shí)，無(wú)法估計(jì)不確定性。網(wǎng)格振蕩收斂時(shí)，不確定性由下式估計(jì)：

式中：SU，SL分別為最密網(wǎng)格計(jì)算過(guò)程中的上、下限；UG為網(wǎng)格不確定性。

網(wǎng)格收斂時(shí)，采用廣義理查森外推法計(jì)算誤差的單項(xiàng)估算值和準(zhǔn)確度階數(shù)：

式中，rG為網(wǎng)格細(xì)化比。

式中，pGest為所采用數(shù)值方法的極限或理論精度，一般可以選擇離散精度階數(shù)。

如果修正因子遠(yuǎn)離1，缺少置信度，數(shù)值不確定度可以定義為：

在進(jìn)行網(wǎng)格收斂性分析時(shí)，采用粗（S3）、中（S2）、細(xì)（S1）三套網(wǎng)格，網(wǎng)格細(xì)化比不同縮尺比下的網(wǎng)格數(shù)量如表2 所示。不同縮尺比下X，Y和N的網(wǎng)格不確定性分析如表3～表7 所示。由表可知，當(dāng)L/LOA=1 時(shí)，X，Y和N的網(wǎng)格收斂率均滿(mǎn)足RG<0，即網(wǎng)格振蕩收斂，不確定度為2.83%～8.48%；當(dāng)L/LOA=2 時(shí)，X，Y和N的網(wǎng)格收斂率均滿(mǎn)足01，即網(wǎng)格發(fā)散；當(dāng)L/LOA=16 時(shí)，X的網(wǎng)格收斂率RG<0，即網(wǎng)格振蕩收斂，不確定度為4.09%，但是Y和N的網(wǎng)格收斂率RG>1，即網(wǎng)格發(fā)散。

表2 網(wǎng)格收斂性研究網(wǎng)格數(shù)量Tab.2 Grid number in the grid convergence study

表3 L/LOA=1 時(shí)準(zhǔn)確度階數(shù)和網(wǎng)格不確定度分析Tab.3 Estimated order of accuracy and grid uncertainty when L/LOA=1

表4 L/LOA=2 時(shí)準(zhǔn)確度階數(shù)和網(wǎng)格不確定度分析Tab.4 Estimated order of accuracy and grid uncertainty when L/LOA=2

表5 L/LOA=4 時(shí)準(zhǔn)確度階數(shù)和網(wǎng)格不確定度分析Tab.5 Estimated order of accuracy and grid uncertainty when L/LOA=4

表6 L/LOA=8 時(shí)準(zhǔn)確度階數(shù)和網(wǎng)格不確定度分析Tab.6 Estimated order of accuracy and grid uncertainty when L/LOA=8

表7 L/LOA=16 時(shí)準(zhǔn)確度階數(shù)和網(wǎng)格不確定度分析Tab.7 Estimated order of accuracy and grid uncertainty when L/LOA=16

對(duì)于網(wǎng)格收斂性研究中存在網(wǎng)格發(fā)散的情況，其原因是多方面的：潛艇大舵角操舵狀態(tài)下舵附近流場(chǎng)往往存在大尺度流動(dòng)分離現(xiàn)象，流場(chǎng)具有顯著的非線性和非定常性；所采用的DES 方法是一種混合算法，其湍流建模與網(wǎng)格尺度存在關(guān)聯(lián)性，不同密度網(wǎng)格其流場(chǎng)內(nèi)局部流動(dòng)的湍流求解方式也有所差異，導(dǎo)致難以呈現(xiàn)出較好的網(wǎng)格收斂性。

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

5.1 泄出渦系結(jié)構(gòu)

圖5～圖9 給出了不同縮尺比潛艇的艇體及舵附近的無(wú)量綱Q準(zhǔn)則瞬時(shí)等值面。等值面取Q=100s?2，并采用無(wú)因次縱向速度Ux/U0著色?？梢钥闯?，垂直舵附近產(chǎn)生了明顯的泄出渦結(jié)構(gòu)，這主要是由于舵背風(fēng)面處由于舵角所產(chǎn)生的流動(dòng)分離現(xiàn)象所致。通過(guò)觀察不同縮尺比的等值面圖可以發(fā)現(xiàn)，舵附近的泄出渦規(guī)模隨著縮尺比的改變而逐漸變化。在縮尺比較小（L/LOA=1）時(shí)，舵背風(fēng)面產(chǎn)生了大規(guī)模的流動(dòng)分離現(xiàn)象，泄出渦脫落起始位置接近舵導(dǎo)緣位置。同時(shí)，舵葉梢部泄出渦與背風(fēng)面泄出渦發(fā)生明顯的相互作用，導(dǎo)致梢部泄出渦結(jié)構(gòu)不明顯。隨著縮尺比增大（L/LOA=2，4，8，16），泄出渦脫落起始位置向尾緣移動(dòng)，背風(fēng)面泄出渦規(guī)模逐漸減小。同時(shí)，舵葉梢部泄出渦結(jié)構(gòu)逐漸清晰，受背風(fēng)面泄出渦的影響逐漸減小?？傮w而言，在舵角不變的條件下，舵背風(fēng)面流動(dòng)分離現(xiàn)象隨著雷諾數(shù)的增大逐漸減弱。

圖5 L/LOA=1 時(shí)艇體及舵附近的泄出渦結(jié)構(gòu)Fig.5 Vortex structure around the submarine body and rudder with L/LOA=1

圖6 L/LOA=2 時(shí)艇體及舵附近的泄出渦結(jié)構(gòu)Fig.6 Vortex structure around the submarine body and rudder with L/LOA=2

圖7 L/LOA=4 時(shí)艇體及舵附近的泄出渦結(jié)構(gòu)Fig.7 Vortex structure around the submarine body and rudder with L/LOA=4

圖8 L/LOA=8 時(shí)艇體及舵附近的泄出渦結(jié)構(gòu)Fig.8 Vortex structure around the submarine body and rudder with L/LOA=8

圖9 L/LOA=16 時(shí)艇體及舵附近的泄出渦結(jié)構(gòu)Fig.9 Vortex structure around the submarine body and rudder with L/LOA=16

5.2 舵面壓力分布

圖10～圖14 給出了不同縮尺比潛艇的舵中剖面處的壓力分布。壓力采用無(wú)因次系數(shù)CP表示，CP=P/橫坐標(biāo)用無(wú)因次量x/C表示，x/C=0表示舵前緣位置，x/C=1 表示舵尾緣位置。在靠近舵前緣位置，舵面分別達(dá)到正壓和負(fù)壓峰值。在舵剖面中后區(qū)域壓力分布趨于平緩。從圖中可以看出，不同縮尺比下舵剖面壓力分布存在顯著差異，特別是舵背風(fēng)面。在L/LOA=1 時(shí)，舵背風(fēng)面壓力分布以x/C=0.2 為分界線呈現(xiàn)出不同的分布，x/C=0～0.2 區(qū)域，表現(xiàn)為負(fù)壓峰值；x/C=0.2～1.0 區(qū)域，舵面壓力趨近于0，且分布十分平緩。隨著L/LOA逐漸增加，舵背風(fēng)面壓力分布開(kāi)始發(fā)生變化，分界線逐步向尾緣移動(dòng)，且靠近尾緣部分舵面壓力開(kāi)始出現(xiàn)不規(guī)則波動(dòng)，但整體壓力差有所增大。

圖10 L/LOA=1 時(shí)舵中剖面處表面壓力分布Fig.10 Pressure distribution on the midsection of the rudder with L/LOA=1

圖11 L/LOA=2 時(shí)舵中剖面處表面壓力分布Fig.11 Pressure distribution on the midsection of the rudder with L/LOA=2

圖12 L/LOA=4 時(shí)舵中剖面處表面壓力分布Fig.12 Pressure distribution on the midsection of the rudder with L/LOA=4

圖13 L/LOA=8 時(shí)舵中剖面處表面壓力分布Fig.13 Pressure distribution on the midsection

圖14 L/LOA=16 時(shí)舵中剖面處表面壓力分布Fig.14 Pressure distribution on the midsection of the rudder with L/LOA=16

5.3 水動(dòng)力及力矩

圖15～圖17 給出了作用在不同縮尺比潛艇上的無(wú)因次縱向力、橫向力及轉(zhuǎn)首力矩的計(jì)算結(jié)果。對(duì)于縱向力，在雷諾數(shù)Re=1.4×107～31.3×107范圍內(nèi)變化最為明顯，這意味著該雷諾數(shù)區(qū)間內(nèi)尺度效應(yīng)的影響十分顯著。隨著雷諾數(shù)進(jìn)一步增加，縱向力系數(shù)的變化趨緩，尺度效應(yīng)的影響有所下降，但仍然難以忽略。對(duì)于橫向力，在雷諾數(shù)Re=1.4×107～3.9×107范圍內(nèi)變化較為明顯，隨著雷諾數(shù)進(jìn)一步增加，橫向力系數(shù)的變化顯著減緩，意味著尺度效應(yīng)的影響明顯減弱。此外，轉(zhuǎn)首力矩與橫向力的變化規(guī)律基本一致。從以上分析可以看出，與潛艇阻力性能相關(guān)的縱向力系數(shù)其尺度效應(yīng)較為明顯，且隨雷諾數(shù)的增加其影響變化較小。相較而言，與潛艇操縱性能相關(guān)的橫向力及轉(zhuǎn)首力矩系數(shù)其尺度效應(yīng)在雷諾數(shù)高于某一值后影響顯著減弱，橫向力及轉(zhuǎn)首力矩系數(shù)趨近于某一定值，這也意味著在研究潛艇操縱運(yùn)動(dòng)特性時(shí)，尺度效應(yīng)的影響不可忽略。

圖15 不同雷諾數(shù)下的縱向力Fig.15 Longitudinal force under different Reynolds numbers

圖16 不同雷諾數(shù)下的橫向力Fig.16 Lateral force under different Reynolds numbers

圖17 不同雷諾數(shù)下的轉(zhuǎn)首力矩Fig.17 Yaw moment under different Reynolds numbers

6 結(jié)語(yǔ)

本文以SUBOFF 潛艇模型為研究對(duì)象，按照1，2，4，8，16 的縮尺比對(duì)潛艇進(jìn)行等比例縮放，采用CFD 方法對(duì)不同尺度的潛艇模型的直航操舵試驗(yàn)進(jìn)行數(shù)值模擬，對(duì)獲得的繞流場(chǎng)和水動(dòng)力特性進(jìn)行分析，可得到以下結(jié)論：

1）不同雷諾數(shù)下的艇體附近流場(chǎng)存在顯著差異。舵面附近的流動(dòng)分離現(xiàn)象隨雷諾數(shù)的增大而減弱，流動(dòng)分離發(fā)生位置逐步后移，導(dǎo)致泄出渦系結(jié)構(gòu)也有所減弱，其中雷諾數(shù)較小時(shí)變化更為明顯。舵面壓力分布也隨之發(fā)生改變，舵背風(fēng)面的負(fù)壓區(qū)逐漸擴(kuò)大，導(dǎo)致舵剖面壓差隨雷諾數(shù)增大而增大。

2）潛艇縱向力、橫向力及轉(zhuǎn)首力矩系數(shù)隨雷諾數(shù)的變化而發(fā)生顯著變化。隨雷諾數(shù)增大，尺度效應(yīng)對(duì)潛艇縱向力的影響雖有所下降但依然顯著，而對(duì)橫向力和轉(zhuǎn)首力矩的影響顯著減弱。