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

李 佳，黃德波，鄧 銳

（哈爾濱工程大學(xué)多體船技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，哈爾濱 150001）

1 引 言

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

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

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

2 區(qū)域離散

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

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

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

潛器周?chē)鲌?chǎng)的控制域離散如圖1所示，入口邊界到潛器最前端距離為L(zhǎng)（L為潛器計(jì)算模型主體的總長(zhǎng)），出口邊界到潛器最后端距離為2L，控制域半徑為L(zhǎng)。由于潛器周?chē)@流場(chǎng)的數(shù)值模擬為不可壓縮外流問(wèn)題，入口邊界設(shè)置為速度入口，出口邊界設(shè)置為自由出流，來(lái)流方向與控制域軸向平行，周向邊界設(shè)置為對(duì)稱(chēng)邊界。

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

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

3 基本控制方程

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

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

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

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

4 湍流模型

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

RNG k-ε湍流模型輸運(yùn)方程[10]：

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

Realizable k-ε湍流模型輸運(yùn)方程[10]：

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

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

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

SST k-ω湍流模型輸運(yùn)方程[10]：

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

5 自由來(lái)流湍流度

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

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

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

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

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

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

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

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

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

7 結(jié) 論

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

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

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

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

[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é)報(bào),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]劉玉秋,于開(kāi)平,張嘉鐘.水下非流線型航行體減阻的數(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]王福軍.計(jì)算流體動(dòng)力學(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ù)與實(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.