載人潛器阻力的數(shù)值計算方法分析

李 佳，黃德波，鄧 銳

（哈爾濱工程大學(xué)多體船技術(shù)國防重點學(xué)科實驗室，哈爾濱 150001）

1 引 言

潛器的種類繁多，用途非常廣泛，即使作為同一用途的潛器其三維幾何形狀也各不相同，而且常常帶有附體及支架等復(fù)雜結(jié)構(gòu)，因此對潛器水動力性能的研究有一定的困難，針對不同問題使用的研究方法也很多。

計算流體力學(xué)（CFD）方法周期短、費用低，是研究潛器水動力性能的有效手段，將CFD方法的計算結(jié)果與試驗數(shù)據(jù)進(jìn)行比較，可以驗證CFD方法的有效并實現(xiàn)物理試驗所不容易完成的模擬。Le等[1]用數(shù)值方法求解了有制導(dǎo)推進(jìn)器的AUV的定常流動，用有限差分法求解了流函數(shù)的軸對稱微分方程，并將結(jié)果用于計算速度場和壓力場。Fuglestad等[2]用CFD方法計算了HUGIN 3000原型的阻力，并與模型試驗數(shù)據(jù)進(jìn)行了比較，驗證了CFD計算的有效性和高雷諾數(shù)阻力預(yù)報的可能性。張懷新等[3]用有限體積法求解了碟狀水下航行器三維粘性繞流場的不可壓縮RANS方程，得出了航行器阻力性能的數(shù)值計算結(jié)果并與試驗數(shù)據(jù)進(jìn)行了比較。Nishi等[4]用商業(yè)CFD軟件STAR－CD計算了水下航行器的阻力，并與拖曳水池試驗結(jié)果進(jìn)行了比較，考慮阻力和推進(jìn)的平衡，得出航行器的推進(jìn)性能。Phillips等[5]用CFD方法計算了三個不同形狀和尺寸的AUV的阻力，并將結(jié)果與試驗數(shù)據(jù)進(jìn)行了比較，驗證了使用RANS流動求解器完成概念設(shè)計的可能性。劉玉秋等[6]采用混合均質(zhì)流理論、Realizable k-ε湍流模型對某水下非流線型航行體模型進(jìn)行了自然狀態(tài)和通氣產(chǎn)生超空泡狀態(tài)下阻力的數(shù)值模擬，得到了自然和通氣狀態(tài)下非流線型水下航行體阻力的變化規(guī)律。Listak等[7]用CFD方法得出了水下仿生機(jī)器人的水動力特性，并將數(shù)值模擬得出的阻力值與實際物理測量的阻力值進(jìn)行了比較。Seo等[8]為了測量水下滑翔物的性能，用CFD處理器得到水動力系數(shù)從而對俯仰進(jìn)行了控制仿真。

本文所研究的潛器是載人潛器，載人潛器在設(shè)計上形式多樣，關(guān)于水動力性能方面的資料在國內(nèi)外并不多見。本文用數(shù)值模擬方法對一艘載人潛器的縮尺模型的水動力性能進(jìn)行了研究，得到了可行的計算方法和有參考價值的結(jié)論。作者首先對潛器計算模型的繞流場進(jìn)行區(qū)域離散，得到適于計算的多塊結(jié)構(gòu)化網(wǎng)格。然后采用不同的湍流模型，對潛器周圍的繞流場進(jìn)行數(shù)值模擬，并將結(jié)果與試驗值進(jìn)行對比。改變?nèi)肟诘耐牧鲝?qiáng)度，探討入口湍流強(qiáng)度對數(shù)值計算結(jié)果的影響。

2 區(qū)域離散

本文所研究的載人潛器模型的主要尺度參數(shù)如表1所示，表中給出了潛器計算模型主體和底部主要附體的主尺度參數(shù)，縮尺比為1:6。

三維復(fù)雜體的控制域進(jìn)行區(qū)域離散時，生成非結(jié)構(gòu)化網(wǎng)格所用的時間少，而且操作也比較容易，但是往往生成的非結(jié)構(gòu)化網(wǎng)格數(shù)量很大，進(jìn)行數(shù)值求解時是一般的PC機(jī)或工作站所無法實現(xiàn)的。而且對于生成的非結(jié)構(gòu)化網(wǎng)格很難進(jìn)行人為控制，所生成的網(wǎng)格會產(chǎn)生局部區(qū)域過于密集或系數(shù)的現(xiàn)象，網(wǎng)格數(shù)量無法預(yù)測。

本文所研究的三維載人潛器外型復(fù)雜，并且在主體上帶有許多附體，因此將控制域離散為非結(jié)構(gòu)網(wǎng)格的數(shù)量將會很大。結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)按照順序存儲，根據(jù)數(shù)組的下標(biāo)可以方便地索引和查找。結(jié)構(gòu)網(wǎng)格可以減少相應(yīng)的存儲開銷，而且由于網(wǎng)格具有貼體性，流場的計算精度可以大大提高，因此作者花費了大量的時間將控制域離散為結(jié)構(gòu)網(wǎng)格。

潛器周圍流場的控制域離散如圖1所示，入口邊界到潛器最前端距離為L（L為潛器計算模型主體的總長），出口邊界到潛器最后端距離為2L，控制域半徑為L。由于潛器周圍繞流場的數(shù)值模擬為不可壓縮外流問題，入口邊界設(shè)置為速度入口，出口邊界設(shè)置為自由出流，來流方向與控制域軸向平行，周向邊界設(shè)置為對稱邊界。

表1 潛器主要尺度參數(shù)Tab.1 Main parameters of a manned submersible

控制域采用多塊搭接結(jié)構(gòu)化網(wǎng)格進(jìn)行離散，將控制域根據(jù)需要分為96個子塊，在不同子塊內(nèi)分別生成結(jié)構(gòu)網(wǎng)格。多塊搭接結(jié)構(gòu)網(wǎng)格中相鄰兩塊公共邊界上的網(wǎng)格點為兩塊公共的網(wǎng)格點，不同子塊的網(wǎng)格共同充滿整個流場，相互之間沒有重疊區(qū)，相鄰塊之間以公共邊界相連接。在流場計算中，不同塊之間的物理信息通過邊界點進(jìn)行交換。作者盡可能描述潛器型體的具體構(gòu)形，避免過多簡化，只在潛器實際幾何體的基礎(chǔ)上略去較簡單的支桿以及較小的推進(jìn)器。控制域的結(jié)構(gòu)化網(wǎng)格采用O型和H型拓?fù)浣Y(jié)構(gòu)，潛器主體以及主要附體的表面網(wǎng)格分別見圖2～3。

3 基本控制方程

在研究載人潛器的水動力性能時，流體可壓縮性的影響是被忽略的。假設(shè)流體不可壓，則潛器粘性繞流場中流體運動的連續(xù)性方程為[9]：

對于不可壓縮流體，雷諾平均RANS方程的張量形式為[9]：

其中：ui為速度分量時均值，ui′為速度分量脈動值，為雷諾應(yīng)力項。

壓力和速度用SIMPLEC算法耦合，壓力方程離散采用二階格式。

4 湍流模型

為了得到適用于所研究潛器的數(shù)值計算的湍流模型，作者采用了RNG k-ε湍流模型、Realizable k-ε湍流模型、標(biāo)準(zhǔn)k-ω湍流模型以及SST k-ω湍流模型分別對潛器周圍繞流場進(jìn)行數(shù)值計算，并根據(jù)計算結(jié)果改變?nèi)肟谶吔绲耐牧鲝?qiáng)度從而得到適用于潛器數(shù)值計算的方法。

RNG k-ε湍流模型輸運方程[10]：

Yakhot等[11]提出重整化群（Renormalization Group）的理論，RNG k-ε模型有效改善了精度，考慮到了湍流漩渦，為湍流Prandtl數(shù)提供了一個解析公式。

Realizable k-ε湍流模型輸運方程[10]：

Realizable k-ε模型對于平板和圓柱射流的發(fā)散比率預(yù)測更精確，對于旋轉(zhuǎn)流動、強(qiáng)逆壓梯度的邊界層流動、流動分離和二次流有很好的表現(xiàn)。

標(biāo)準(zhǔn)k-ω湍流模型輸運方程[10]：

Kolmogorov首先提出了k-ω湍流模型的概念，Launder和Spalding對k-ω模型進(jìn)行了改進(jìn)，Wilcox進(jìn)一步對k-ω湍流模型進(jìn)行應(yīng)用研究。標(biāo)準(zhǔn)k-ω模型基于Wilcox k-ω模型，是為考慮低雷諾數(shù)、可壓縮性和剪切流傳播而修改的。標(biāo)準(zhǔn)k-ω湍流模型能夠很好地處理近壁面流動問題，但對自由來流有很大的依賴性。

SST k-ω湍流模型輸運方程[10]：

SST k-ω模型是由Menter發(fā)展的，考慮了湍流剪應(yīng)力的輸運特性，能夠準(zhǔn)確地預(yù)報由于逆壓梯度導(dǎo)致的流動分離點和分離區(qū)域。

5 自由來流湍流度

如果兩種流體運動具有幾何相似的邊界，并且它們的速度場是幾何相似的，也就是說，如果它們有幾何相似的流線，則這兩種流體運動是動力學(xué)相似的。根據(jù)動力學(xué)相似原理，如果兩個物體幾何相似，其它相關(guān)條件如相對粗糙度、自由來流湍流度、空泡數(shù)等相同，則物體的阻力系數(shù)就可表示為：

上式中的四個變量中可以略去次要因素而考慮主要因素來實現(xiàn)流動相似，對于載人潛器阻力性能的研究，我們考慮的主要因素是雷諾數(shù)Re。

當(dāng)動力學(xué)相似所考慮的其它條件相同時，物體在靜水中（湍流度為0）航行和自由來流流經(jīng)靜止物體時物體的阻力系數(shù)是不同的[12]，這是因為自由來流的湍流度不同。也就是說載人潛器的計算模型的數(shù)值模擬無法做到與實際情況相同，自由來流的湍流度不能設(shè)置為0，因此入口邊界的湍流度的取值會對模擬結(jié)果產(chǎn)生影響。

充分發(fā)展的管內(nèi)流動可以根據(jù)半經(jīng)驗公式給出入口邊界處的湍流強(qiáng)度和湍流長度尺度，而不可壓縮外部繞流只能根據(jù)經(jīng)驗或試驗給出入口的湍流強(qiáng)度I和湍動粘度比μt/μ，湍流強(qiáng)度應(yīng)為很小的量，而湍動粘度比應(yīng)在1～10之間。為了研究適于潛器流場數(shù)值計算的入口邊界條件，入口處湍流強(qiáng)度和湍動粘度比的設(shè)置首先根據(jù)文獻(xiàn)[13]取I=0.5%，μt/μ=5，然后改變?nèi)肟谔幍耐牧鲝?qiáng)度I，用不同的湍流模型進(jìn)行潛器周圍流場的數(shù)值計算，將計算結(jié)果與試驗值相比較，得出適于計算的入口湍流強(qiáng)度。

6 計算結(jié)果及分析

為了探討自由來流湍流強(qiáng)度對潛器阻力數(shù)值模擬結(jié)果的影響，以及不同湍流模型對自由來流湍流度的依賴程度，首先設(shè)置入口湍流強(qiáng)度I=0.5%以及I=0.1%，分別采用RNG k-ε湍流模型、Realizable k-ε湍流模型、標(biāo)準(zhǔn)k-ω湍流模型以及SST k-ω湍流模型對潛器周圍繞流場進(jìn)行數(shù)值模擬。阻力計算結(jié)果如圖4～5所示，其中Vm為來流速度，Rt為潛器總阻力。對比可見，改變?nèi)肟谕牧鲝?qiáng)度時，采用標(biāo)準(zhǔn)k-ω湍流模型得到的潛器阻力計算結(jié)果差別很大，而采用其它三種湍流模型得到的計算結(jié)果基本沒有差別。下面對標(biāo)準(zhǔn)k-ω湍流模型對自由來流湍流強(qiáng)度的敏感程度作簡要分析。

根據(jù)上面的計算結(jié)果，入口湍流強(qiáng)度的變化范圍取為I=0.1%-0.5%，不同來流速度下采用適當(dāng)?shù)淖杂蓙砹魍牧鲝?qiáng)度，用標(biāo)準(zhǔn)k-ω湍流模型得出的潛器阻力計算結(jié)果如圖6所示。相應(yīng)速度下設(shè)置的入口湍流強(qiáng)度與雷諾數(shù)的關(guān)系曲線見圖7，其中Re=VmL/ν為雷諾數(shù)，運動粘性系數(shù)ν=1.12434×10-6m2/s，圖中還給出了該關(guān)系曲線的擬合趨勢曲線。由圖6可知，設(shè)置適當(dāng)?shù)娜肟谕牧鲝?qiáng)度采用標(biāo)準(zhǔn)k-ω湍流模型對潛器模型的阻力計算的結(jié)果與試驗值是非常吻合的。對于不可壓縮外流而言，入口湍流強(qiáng)度的取值并不像內(nèi)流一樣有經(jīng)驗公式可循，我們只能根據(jù)試驗值找出一定的規(guī)律。圖7中自由來流湍流強(qiáng)度取值的趨勢曲線的函數(shù)為I=2×107×Re-1.2298（%），此函數(shù)關(guān)系是對比載人潛器模型的阻力試驗值，采用標(biāo)準(zhǔn)k-ω湍流模型進(jìn)行載人潛器模型的阻力計算時對于不同來流速度下入口湍流強(qiáng)度的設(shè)置依據(jù)。

采用RNG k-ε湍流模型、Realizable k-ε湍流模型以及SST k-ω湍流模型對潛器的阻力進(jìn)行數(shù)值計算時，自由來流的湍流強(qiáng)度分別取I=0.5%和I=0.1%，阻力的計算結(jié)果變化很小。因此自由來流湍流強(qiáng)度的變化范圍在I=0.1%-0.5%時，采用這三種湍流模型對潛器進(jìn)行阻力計算的結(jié)果將不會有過大差別。

圖8～9給出了潛器模型速度為Vm＝2.6m/s，采用Realizable k-ε湍流模型，入口湍流強(qiáng)度設(shè)置為I＝0.5%時潛器主體和主要附體的壓力云圖，由圖中可以看到附體前后壓差很大，并且附體之間的流場產(chǎn)生了相互干擾。

7 結(jié) 論

為了得到適于潛器周圍繞流場進(jìn)行數(shù)值模擬的計算方法，本文對某載人潛器周圍控制域進(jìn)行多塊搭接結(jié)構(gòu)化網(wǎng)格離散，分別采用RNG k-ε湍流模型、Realizable k-ε湍流模型、標(biāo)準(zhǔn)k-ω湍流模型以及SST k-ω湍流模型對潛器周圍繞流場進(jìn)行數(shù)值模擬，并改變自由來流的湍流強(qiáng)度的設(shè)置，比較了阻力計算結(jié)果并將其與試驗值進(jìn)行了對比。計算結(jié)果分析表明：

（1）采用標(biāo)準(zhǔn)k-ω湍流模型對潛器周圍繞流場進(jìn)行數(shù)值計算時，自由來流湍流強(qiáng)度在I=0.1%-0.5%范圍內(nèi)變化，得到的潛器阻力計算結(jié)果有所不同。根據(jù)阻力試驗值，可以給出自由來流湍流強(qiáng)度I=2×107×Re-1.2298（%），采用標(biāo)準(zhǔn)k-ω湍流模型進(jìn)行潛器阻力的數(shù)值計算時不同速度下入口湍流強(qiáng)度的設(shè)置要有變化，可以參照此函數(shù)關(guān)系進(jìn)行設(shè)置。

（2）采用RNG k-ε湍流模型、Realizable k-ε湍流模型以及SST k-ω湍流模型對潛器的阻力進(jìn)行數(shù)值計算時，自由來流湍流強(qiáng)度在I=0.1%-0.5%范圍內(nèi)變化對潛器阻力的計算結(jié)果影響很小。采用這三種湍流模型進(jìn)行潛器阻力的數(shù)值計算時，可以將入口湍流強(qiáng)度設(shè)置為I=0.5%。

本文對采用標(biāo)準(zhǔn)k-ω湍流模型進(jìn)行潛器阻力數(shù)值計算情況時，是在有試驗值可依的情況下得出的自由來流湍流強(qiáng)度的設(shè)置隨雷諾數(shù)的變化規(guī)律，此規(guī)律還要經(jīng)過大量的計算以及試驗驗證。

[1]Le Page Y G,Holappa K W.Hydrodynamics of an autonomous underwater vehicle equipped with a vectored thruster[C]//Oceans Conference Record(IEEE).Oceans 2000,Providence RI USA,2000:2135-2140.

[2]Fuglestad A L,Grahl-Madsen M.Computational fluid dynamics applied on an autonomous underwater vehicle[C]//Proceedings of the 23rd International Conference on Offshore Mechanics and Arctic Engineering.Vancouver BC Canada,2004,USA:American Society of Mechanical Engineers,2004:447-451.

[3]張懷新,潘雨村.圓碟形潛水器阻力性能研究[J].上海交通大學(xué)學(xué)報,2006,40(6):978-982,987.Zhang Huaixin,Pan Yucun.The resistance performance of a dish-shaped underwater vehicle[J].Journal of Shanghai Jiaotong University,2006,40(6):978-982,987.

[4]Nishi Y,Kashiwagi M,Koterayama W,et al.Resistance and propulsion performance of an underwater vehicle estimated by a CFD method and experiment[C]//Proceedings of the 17th 2007 International Offshore and Polar Engineering Conference.Lisbon Portugal,2007,Cupertino,CA 95015-0189,United States:International Society of Offshore and Polar Engineers,2007:2045-2052.

[5]Phillips A,Furlong M,Turnock S R.The use of computational fluid dynamics to assess the hull resistance of concept autonomous underwater vehicles[C]//OCEANS 2007-Europe.Aberdeen Scotland United kingdom,2007.

[6]劉玉秋,于開平,張嘉鐘.水下非流線型航行體減阻的數(shù)值模擬與比較[J].工程力學(xué),2007,24(2):178-182.Liu Yuqiu,Yu Kaiping,Zhang Jiazhong.Numerical simulation and comparison on drag reduction of non-streamline submerged body[J].Engineering Mechanics,2007,24(2):178-182.

[7]Listak M,Pugal D,Kruusmaa M.CFD simulations and real world measurements of drag of biologically inspired underwater robot[C]//US/EU-Baltic International Symposium:Ocean Observations,Ecosystem-Based Management and Forecasting.Tallinn Estonia,2008.

[8]Seo D C,Jo G,Choi H S.Pitching control simulations of an underwater glider using CFD analysis[C]//Oceans’08 MTS/IEEE Kobe Technology Ocean’08 Voyage toward the Future.Kobe Japan,2008.

[9]王福軍.計算流體動力學(xué)分析—CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社,2004.Wang Fujun.Computational fluid dynamics analysis－principle and application of CFD software[M].Beijing:Tsinghua U-niversity Press,2004.

[10]Fluent Inc.,FLUENT User’s Guide[K].Fluent Inc,2003.

[11]Yakhot V,Orzag S A.Renormalization group analysis of turbulence.Ⅰ.Basic theory[J].Journal of Science Computing,1986,1(1):3-51.

[12]Douglas J F,Gasiorek J M,Swaffield J A.Fluid Mechanics[M].Great Britain:The Pitman Press,1979.

[13]韓占忠,王 敬,蘭小平.FLUENT—流體工程仿真技術(shù)與實例應(yīng)用[M].北京:北京理工大學(xué)出版社,2004.Han Zhanzhong,Wang Jing,Lan Xiaoping.FLUENT－Emulational technique and instance application of fluid engnieering[M].Beijing:Beijing Insitute of Technology Press,2004.