洪方文，張志榮，劉登成，鄭巢生，翟樹成

（1.中國(guó)船舶科學(xué)研究中心，江蘇無錫214082；2.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，江蘇無錫214082）

0 引 言

自計(jì)算機(jī)誕生以來，數(shù)值計(jì)算成為科學(xué)研究的一種主要手段。在流體力學(xué)領(lǐng)域中，經(jīng)過七十年的發(fā)展，計(jì)算流體力學(xué)（CFD）已經(jīng)成為船舶流動(dòng)分析的主要工具。船舶螺旋槳CFD 已經(jīng)被廣泛應(yīng)用于推進(jìn)器設(shè)計(jì)中，使得設(shè)計(jì)目標(biāo)實(shí)現(xiàn)的確定性大為提高。通過廣泛使用發(fā)現(xiàn)，大部分情況下螺旋槳敞水性能CFD預(yù)報(bào)的精度誤差都能達(dá)到不超過3%的水平，滿足螺旋槳設(shè)計(jì)的工程需要。不過同樣也會(huì)出現(xiàn)預(yù)報(bào)精度很差的情況，甚至誤差會(huì)達(dá)到20%的狀況，這說明螺旋槳設(shè)計(jì)中CFD的應(yīng)用還存在一定的技術(shù)風(fēng)險(xiǎn)，進(jìn)一步完善船舶CFD技術(shù)，提高其計(jì)算精度是一項(xiàng)重要的工作。

數(shù)值求解螺旋槳流動(dòng)問題始于20 世紀(jì)80 年代，最初針對(duì)非湍流控制方程、有勢(shì)流方程[1]和Euler方程[2]、低馬赫數(shù)可以壓縮流動(dòng)N-S 方程[3-5]及不可壓縮層流N-S 方程[6-7]進(jìn)行求解。90 年代初期開始出現(xiàn)了求解不可壓縮RANS方程的螺旋槳湍流流場(chǎng)CFD計(jì)算[8-10],使得計(jì)算精度有了很大的提高。在90 年代中期螺旋槳的定常流計(jì)算已相當(dāng)廣泛，并應(yīng)用到螺旋槳水動(dòng)力性能的預(yù)報(bào)中[11-13]。到90 年代末期，螺旋槳CFD 已經(jīng)比較成熟，可以較好地預(yù)報(bào)螺旋槳敞水性能。1998 年，第22 屆ITTC 推進(jìn)技術(shù)委員會(huì)在法國(guó)Grenoble 專門舉行了RANS/面元法螺旋槳性能比較計(jì)算研討會(huì)[14-16]，研討會(huì)的主要結(jié)論是RANS 和面元法都可精確地預(yù)報(bào)螺旋槳的敞水性能，能夠用于螺旋槳設(shè)計(jì)階段的性能預(yù)報(bào)。進(jìn)入本世紀(jì)后，CFD 已經(jīng)廣泛應(yīng)用到螺旋槳的精細(xì)流動(dòng)模擬和性能分析上[17-20]。現(xiàn)在求解RANS的CFD 雖然已被廣泛用于螺旋槳水動(dòng)力性能計(jì)算和流動(dòng)模擬，但總是存在一些不足的地方。水動(dòng)力計(jì)算中在重負(fù)荷的情況下誤差變大，雖然雷諾數(shù)增加計(jì)算誤差有減小的趨勢(shì)[21-23]，但仍會(huì)出現(xiàn)不可接受的水平。在梢渦和轂渦的計(jì)算中，渦心尺寸和速度耗散總是偏大[24-27]。

出現(xiàn)這些現(xiàn)象，主要是由于螺旋槳葉片的流動(dòng)狀態(tài)沒有得到正確的模擬。RANS方法一般把流動(dòng)作為全湍流進(jìn)行模擬，然而對(duì)于模型螺旋槳流動(dòng)，由于雷諾數(shù)在105量級(jí)，在槳葉上還存在較大面積的層流狀態(tài)[28]，需要使用具有轉(zhuǎn)捩能力的湍流模型進(jìn)行模擬[23，29-31]。在本研究中將利用具有轉(zhuǎn)捩能力的湍流模型計(jì)算螺旋槳模型流動(dòng)，分析其計(jì)算螺旋槳水動(dòng)力特性的精度。

1 數(shù)值方法

1.1 控制方程

控制方程使用旋轉(zhuǎn)坐標(biāo)系下不可壓縮流體雷諾平均質(zhì)量和動(dòng)量守恒方程:

1.2 數(shù)值求解

流動(dòng)的數(shù)值求解使用商用軟件FLUENT，該軟件對(duì)控制方程的離散使用有限體積法。動(dòng)量守恒方程和湍流方程中的對(duì)流項(xiàng)離散使用二階迎風(fēng)格式，擴(kuò)散項(xiàng)使用二階中心差分格式,流場(chǎng)中物理量梯度計(jì)算使用基于單元的Green-Gauss 方法。離散方程求解利用SIMPLE 方法和Gauss-Seidel 迭代，同時(shí)求解過程中使用多重網(wǎng)格技術(shù)加速迭代的收斂。壓力計(jì)算松馳因子取0.3，速度計(jì)算松馳因子取0.7，湍流計(jì)算松馳因子取0.8。計(jì)算中各方程的收斂條件為10-4，但整個(gè)計(jì)算過程的收斂通過螺旋槳的推力系數(shù)變化來判斷，當(dāng)螺旋槳推力系數(shù)變化小于10-3時(shí)，可認(rèn)為計(jì)算過程收斂，結(jié)束計(jì)算。

1.3 初始邊界條件

在旋轉(zhuǎn)坐標(biāo)系下，假設(shè)流動(dòng)具有與螺旋槳葉片數(shù)相同的周期性，計(jì)算區(qū)域取一個(gè)流道的扇形區(qū)域。進(jìn)口在螺旋槳前方6D處，出口在螺旋槳后方10D處，區(qū)域半徑為6D，D是螺旋槳的直徑。為了劃分高質(zhì)量的網(wǎng)格，把計(jì)算區(qū)域劃分為6塊子區(qū)域，子區(qū)域形式顯示在輻射面內(nèi)，如圖1所示。

計(jì)算區(qū)域邊界包含進(jìn)口邊界、出口邊界、頂部圓周邊界、螺旋槳葉片、槳轂、前軸、后軸和周期邊界。進(jìn)口邊界和頂部圓周邊界施加速度進(jìn)口邊界條件，設(shè)定進(jìn)口速度和湍流相關(guān)參數(shù)；出口邊界施加壓力邊界條件，給定出口壓力；螺旋槳葉片和槳轂使用不可滑移物面邊界條件；前軸和后軸使用滑移物面邊界條件；周期邊界施加周期邊界條件；計(jì)算的速度初始條件使用絕對(duì)坐標(biāo)下的來流均勻速度。

圖1 計(jì)算區(qū)域劃分形式示意Fig.1 Division of computing domains

2 對(duì)象和網(wǎng)格

研究對(duì)象選為28 000 DWT多用途船的螺旋槳模型，其主參數(shù)如表1所示。

表1 螺旋槳模型主參數(shù)Tab.1 Main parameters of propeller model

在計(jì)算域內(nèi)，包含螺旋槳的區(qū)域使用非結(jié)構(gòu)化網(wǎng)格，其他區(qū)域使用結(jié)構(gòu)化網(wǎng)格。螺旋槳葉片導(dǎo)隨邊網(wǎng)格加密，導(dǎo)隨邊網(wǎng)格尺度以1%D 為基礎(chǔ)變化，葉片上的網(wǎng)格形式如圖2所示。葉片附近網(wǎng)格以螺旋槳葉片上的網(wǎng)格尺寸為基礎(chǔ)，向外以1.1的比例增長(zhǎng)；在螺旋槳的進(jìn)口區(qū)域，軸向布置15 個(gè)網(wǎng)格，進(jìn)口處尺寸為1.0D，與螺旋槳區(qū)域交接處網(wǎng)格尺寸為6%D；在螺旋槳的出口區(qū)域，軸向布置20個(gè)網(wǎng)格，出口處尺寸為1.0D，與螺旋槳區(qū)域交接處為6%D，徑向使用15個(gè)網(wǎng)格；在頂面處尺寸為1.0D，與螺旋槳區(qū)域交接處網(wǎng)格尺寸為12%D，周向均勻布置39 個(gè)網(wǎng)格；螺旋槳的前后區(qū)域使用桶型平推網(wǎng)格。在研究中將以這里描述的網(wǎng)格為基礎(chǔ)進(jìn)行加密和稀疏。

圖2 螺旋槳葉片上的網(wǎng)格形式Fig.2 Meshes on the propeller blade

3 結(jié)果和分析

計(jì)算中進(jìn)口速度V=1.782 45 m/s，螺旋槳轉(zhuǎn)速N=17 r/s，螺旋槳進(jìn)速系數(shù)J=0.5，計(jì)算中參考?jí)毫橐粋€(gè)大氣壓，出口相對(duì)壓力設(shè)為0.0。分別利用SST k-ω湍流模型和γ - R?eθ轉(zhuǎn)捩模型，以及4套網(wǎng)格進(jìn)行計(jì)算，計(jì)算結(jié)果如表2和表3所示。KT和10KQ的試驗(yàn)結(jié)果為0.149 1和0.192 0。從表中可以看出，對(duì)于SST k-ω 模型，KQ的計(jì)算誤差很小，并且不同的網(wǎng)格數(shù)得到的計(jì)算結(jié)果很穩(wěn)定，但KT的計(jì)算誤差較大，隨網(wǎng)格數(shù)的增加，呈增大和發(fā)散趨勢(shì)。γ - R?eθ模式計(jì)算得到的誤差，KT和KQ都在一個(gè)較好的水平，并且隨網(wǎng)格數(shù)的變化穩(wěn)定性較好。

表2 KT和KQ計(jì)算結(jié)果Tab.2 Calculating results of KT and KQ

表3 KT和KQ計(jì)算誤差Tab.3 Calculating error of KT and KQ

計(jì)算顯示，使用帶有轉(zhuǎn)捩的湍流模型計(jì)算結(jié)果與試驗(yàn)更為吻合，這說明槳葉表面有很大的區(qū)域保持著層流狀態(tài)。在RANS 方法中，默認(rèn)物體表面流動(dòng)全處于湍流狀態(tài)，而實(shí)際中物體表面起始流動(dòng)總是處于層流狀態(tài)，湍流是隨著層流的發(fā)展轉(zhuǎn)捩后形成的。在流動(dòng)雷諾數(shù)很高時(shí)，物面開始的層流區(qū)域較小，絕大部分區(qū)域都處于湍流狀態(tài)，對(duì)于RANS 方法中默認(rèn)引起的誤差是很小的。在流動(dòng)雷諾數(shù)較小的情況，層流區(qū)域占有物面相當(dāng)?shù)谋壤?，如果全部?dāng)作湍流處理，可能會(huì)引起較大的誤差。本文計(jì)算對(duì)象在17 r/s 的轉(zhuǎn)速下0.7R 處的流動(dòng)雷諾數(shù)為3.4×105。一般情況下，對(duì)于物面的臨界雷諾數(shù)為2.0×105-6.0×105。這樣看來，槳葉表面可能還有50%的區(qū)域處于層流狀態(tài)。圖3給出了兩種湍流模式計(jì)算的槳葉表面流態(tài)，同時(shí)給出了另一螺旋槳的試驗(yàn)結(jié)果。從流態(tài)看非轉(zhuǎn)捩模型計(jì)算的結(jié)果，壓力面分離區(qū)較小，且吸力面和壓力面沿半徑方向的流動(dòng)更明顯，呈現(xiàn)湍流的流動(dòng)特性。轉(zhuǎn)捩模型計(jì)算的結(jié)果與試驗(yàn)結(jié)果更相似，在吸力面有較大的分離區(qū)，吸力面和壓力面徑向流動(dòng)更明顯，呈現(xiàn)層流的流動(dòng)特性，這與試驗(yàn)結(jié)果一致。在數(shù)值模擬時(shí)，如果使用全湍流模式將對(duì)計(jì)算結(jié)果有較大的影響。

圖3 槳葉表面流動(dòng)狀態(tài)Fig.3 Flow traces on blade surface

轉(zhuǎn)捩湍流模型應(yīng)有更寬的適用范圍，它既可以處理層流占比大的問題，也可以處理層流占比小的流動(dòng)問題，也就預(yù)示著在高雷諾數(shù)時(shí)，轉(zhuǎn)捩湍流模型計(jì)算結(jié)果將與非轉(zhuǎn)捩湍流模型一致，以及在低雷諾數(shù)時(shí)，對(duì)于不同的對(duì)象，其計(jì)算結(jié)果的精度都將得到改善。表4 列出了雷諾數(shù)提高10倍的計(jì)算結(jié)果，可以看出在高雷諾數(shù)的情況下兩種湍流模式計(jì)算的螺旋槳水動(dòng)力系數(shù)相差很小，在0.1%以內(nèi)。表5 給出了另外4 個(gè)算例的計(jì)算誤差，從表中可以看出，使用轉(zhuǎn)捩湍流模型計(jì)算結(jié)果的精度都得到了不同程度的改善，說明轉(zhuǎn)捩湍流模型具有更寬的適用范圍。

表4 不同雷諾數(shù)計(jì)算結(jié)果Tab.4 Calculating results of KT and KQ at different Re numbers

表5 系列算例的計(jì)算結(jié)果誤差Tab.5 Calculating errors of KT and KQ for a few propellers

4 結(jié) 語

在模型試驗(yàn)中，螺旋槳葉片流動(dòng)雷諾數(shù)Re=5.0×105左右，雖然已經(jīng)大于相關(guān)試驗(yàn)規(guī)程中要求的臨界雷諾數(shù)，但螺旋槳葉片表面仍然有相當(dāng)區(qū)域的流動(dòng)為層流狀態(tài)。一般情況下，螺旋槳模型敞水的流動(dòng)計(jì)算都使用基于雷諾平均方程的湍流模式，即把螺旋槳葉片表面的流動(dòng)全部作為湍流處理。但作為全湍流處理，并不符合實(shí)際的流動(dòng)情況，所以計(jì)算結(jié)果有時(shí)會(huì)出現(xiàn)較大的偏差，必須使用帶轉(zhuǎn)捩的湍流模式進(jìn)行螺旋槳模型的敞水流動(dòng)計(jì)算，才能保證計(jì)算精度。