轉捩模型在螺旋槳數值計算中的應用

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

（1.中國船舶科學研究中心，江蘇無錫214082；2.船舶振動噪聲重點實驗室，江蘇無錫214082）

0 引 言

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

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

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

1 數值方法

1.1 控制方程

控制方程使用旋轉坐標系下不可壓縮流體雷諾平均質量和動量守恒方程:

1.2 數值求解

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

1.3 初始邊界條件

在旋轉坐標系下，假設流動具有與螺旋槳葉片數相同的周期性，計算區域取一個流道的扇形區域。進口在螺旋槳前方6D處，出口在螺旋槳后方10D處，區域半徑為6D，D是螺旋槳的直徑。為了劃分高質量的網格，把計算區域劃分為6塊子區域，子區域形式顯示在輻射面內，如圖1所示。

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

圖1 計算區域劃分形式示意Fig.1 Division of computing domains

2 對象和網格

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

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

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

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

3 結果和分析

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

表2 KT和KQ計算結果Tab.2 Calculating results of KT and KQ

表3 KT和KQ計算誤差Tab.3 Calculating error of KT and KQ

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

圖3 槳葉表面流動狀態Fig.3 Flow traces on blade surface

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

表4 不同雷諾數計算結果Tab.4 Calculating results of KT and KQ at different Re numbers

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

4 結 語

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