不均勻進流場下螺旋槳非定常力的數值模擬

張凱奇，葉金銘，于安斌

(海軍工程大學 艦船工程系，湖北 武漢 430033)

0 引 言

對于工作在船后或艇后的螺旋槳，由于其處在徑向和周向都不均勻的三維伴流場中，因而會產生周期性變化的非定常力，進而激發出離散譜噪聲。所以，螺旋槳非定常力的分析和研究對預報螺旋槳線譜噪聲和提高潛艇的隱身性至關重要。對螺旋槳非定常力的研究，模型試驗是比較傳統的方法，但由于試驗周期長、成本高，還會受到尺度效應、外界干擾、伴流場的模擬等因素的影響。因而模型試驗方法受到了很大的限制，這方面的文章也很少見。從20世紀70年代開始，基于勢流理論的螺旋槳非定常力數值分析得到了迅速發展。在這方面的代表學者，國外主要有Kerwin[1]，Hoshino[2]，國內有陳家棟[3]、譚延壽[4]、熊鷹[5]等。但這種方法需要處理庫塔條件和假定尾渦形狀，并且很難較好地處理粘性的影響。

隨著計算流體力學CFD技術的迅猛發展，粘性數值計算方法逐漸成為研究螺旋槳非定常力的另一有效方法。張志榮等[6]基于CFD技術，采用3種方法對斜流條件下大側斜螺旋槳的非定常水動力性能進行預報研究，分析了不同方法對計算螺旋槳非定常水動力性能的適用性。胡小菲[7]對非均勻來流下的螺旋槳所受的非定常力進行數值模擬，并對網格的收斂性進行分析，計算結果與試驗值吻合較好。黃振宇[8]采用SSTkω湍流模型和滑移網格技術對非均勻來流條件下的2個螺旋槳所受的非定常力進行數值模擬，計算結果與試驗值吻合較好。在此基礎上，對全附體潛艇模型后的大測斜螺旋槳的非定常力進行計算。劉登成等[9]基于Fluent軟件，以DTMB4119為研究對象，對速度進口到槳盤面的距離以及網格數量進行分析。并根據建立的數值方法，對HSP螺旋槳非定常力進行數值預報。馬超[10]采用滑移網格技術，結合RNG湍流模型對Suboff全附體艇體下的螺旋槳的非定常力進行數值計算，獲得螺旋槳轉動一周過程中軸承力的變化規律。

本文基于STAR-CCM+仿真軟件，使用滑移網格技術，對非均勻來流下的螺旋槳非定常力預報方法進行研究，通過與試驗值進行比較，具體分析了時間離散格式、螺旋槳每時間步旋轉的度數以及網格密度對計算結果的影響。并在此數值方法的基礎上，對艇體非均勻伴流場下的螺旋槳所受到的非定常力進行了數值分析。

1 控制方程與湍流模型的選取

1.1 控制方程

不可壓縮粘性流體的基本控制方程由連續性方程和RANS方程組成，其張量形式為：

1.2 湍流模型

顯然，上面的控制方程不封閉，因而要對流場進行求解還需要選取適當的湍流模型。通過對相關文獻的閱讀，發現對于螺旋槳非定常水動力性能的數值模擬，研究人員較多的采用SSTk-ω湍流模型。因此，本文選取SSTk-ω湍流模型對方程組中的雷諾應力項進行封閉。

該模型由Menter提出，又稱為剪切應力輸運kω模型，它綜合了邊界層外部k-ε模型獨立性和近壁kω模型穩定性的優點。在近壁面區域采用k-ω湍流模式，在邊界層外部采用k-ε湍流模式。該模型的k方程和ω方程分別為：

2 非均勻來流下螺旋槳非定常力的數值模擬

2.1 計算對象與計算域

對于非均勻來流下螺旋槳非定常力的計算，本文選取DTMB4119槳為研究對象。該槳模為3葉槳，直徑D=0.305 m，轂徑比Dh/D=0.2，設計進速系數J=0.833，無側斜和縱傾，右旋槳，幾何模型如圖1所示。

圖 1 DTMB4119槳的幾何模型Fig. 1 DTMB4119 model

在螺旋槳上建立柱坐標系（x，r，θ），原點為槳盤面與槳軸的交點，x軸與槳軸重合，指向下游為正；r為徑向，指向外為正；θ為周向，從槳后向槳前看順時針為正，0°在12點鐘位置。

考慮到要使用滑移網格技術進行非定常計算，因而計算流場區域分為兩部分，外部的靜止域和內部的旋轉域。外部區域取為槳前4D、槳后7.5D、半徑為2.5D的大圓柱形，內部區域取為槳前0.23D、槳后0.4D、半徑為0.65D的小圓柱形，2個圓柱的軸線與槳軸重合，整個計算域如圖2所示。

圖 2 數值模擬計算域Fig. 2 Computation domain of numerical simulation

2.2 邊界條件與網格劃分

整個計算區域的邊界分為進流邊界、出流邊界、遠場邊界、螺旋槳以及槳軸表面。具體的邊界設置為：進流邊界取為速度入口邊界類型，出流邊界取為壓力出口邊界類型，遠場邊界取為對稱面邊界類型，螺旋槳以及槳軸表面取為無滑移壁面邊界類型。對于非均勻來流的設置，本文采用Jessup（1990）[11]在DT-RC報告中九階角頻的軸流場，在進流邊界各個半徑處（R為螺旋槳半徑）的周向速度分布如圖3所示。本文在設計工況J=0.833下，取螺旋槳轉速為25rps，則平均來流速度U0=6.35 m/s。

圖 3 DTRC 4119槳的進流條件Fig. 3 The inflow conditions of DTRC 4119

STAR-CCM+自帶的網格生成技術可以快速高效的生成質量很好的計算網格，本文采用以六面體為核心的切割體網格（Trimmer）、拉伸型網格（Extruder）和棱柱層網格（Prism Layer Mesher）對計算域進行離散。考慮到在入口邊界處要進行非均勻來流的設置，為了更好的反應各個半徑周向不同位置處速度的變化，對入口邊界的網格劃分，要從外到內逐步加密，如圖4所示。最終得到的入口邊界速度分布云圖如圖5所示，從圖中可以看出，速度分布的九階特征比較明顯，因而入口邊界的網格滿足非均勻來流速度設置的要求。另外，根據螺旋槳周圍流場的流動特點，通過體積控制（Volumetric Controls）功能模塊對靠近螺旋槳的區域分層逐步加密（見圖6），以便更好地捕捉流動細節，獲得更加精確地流場信息。最后，考慮到在螺旋槳的導邊、隨邊、葉梢以及葉根處各個物理量的變化梯度較大，對其也進行了加密（見圖7）。

圖 4 入口邊界網格劃分Fig. 4 Inlet boundary mesh

圖 5 入口邊界的速度云圖Fig. 5 velocity contour of inlet boundary

圖 6 對稱面上的網格分布Fig. 6 Mesh distribution of symmetry plane

2.3 計算結果與分析

基于上述的計算模型、網格系統，運用有限體積法對螺旋槳的非定常力進行數值模擬。通過將數值計算結果與Jessup （1990）[1]的試驗數據進行對比，分析時間離散格式與單位時間步長旋轉角度以及網格密度對計算結果的影響，得到了適用于螺旋槳非定常力計算的數值分析方法。

2.3.1 時間離散格式與單位時間步長螺旋槳旋轉角度的確定

圖 7 槳葉和槳榖表面上的網格分布Fig. 7 Mesh distribution of blades and hub

為了研究時間離散格式與單位時間步長螺旋槳旋轉角度對非定常力的影響，本文分別在1階時間離散格式和2階時間離散格式下，保持網格劃分方式以及網格密度不變，通過改變螺旋槳單位時間步長的旋轉角度（3.6°，2.7°，1.8°，0.9°，0.45°）對螺旋槳的非定常力進行數值計算。將得到的數值計算結果與試驗值進行比較，如圖8所示。圖中，0°在12點鐘位置，從槳后向槳前看，順時針方向旋轉。

圖 8 時間離散格式及旋轉角度對非定常力的影響Fig. 8 The impact of temporal discretization and rotation angle on unsteady hydrodynamics

通過比較分析，可以得出如下結論：首先，無論是采用1階時間離散格式還是2階時間離散格式，隨著螺旋槳單位時間步長旋轉角度的減小，數值計算結果與試驗值越來越接近，精度也越來越高。另外，在螺旋槳單位時間步長旋轉角度相同的情況下，2階格式的數值計算結果與試驗值的吻合度明顯比1階格式更高；最后，對于2階時間離散格式，當螺旋槳單位時間步長的旋轉角度小于1.8°時，數值計算結果與1.8°差別很小。綜上所述，本文認為對于螺旋槳非定常力的計算，螺旋槳單位時間步長的旋轉角度應為1.8°同時選取2階時間離散格式。

2.3.2 網格密度對螺旋槳非定常力的影響

對于非均勻來流下螺旋槳非定常力的數值模擬，由于靠近槳葉表面附近的流動較為復雜，各個物理量的變化梯度較大，因而要對其周圍的網格進行加密。而在STAR-CCM+仿真軟件中，可以通過設置槳葉表面網格的邊界增長率（Boundary Growth Rate）的快慢，對其周圍網格疏密程度進行調整。為了研究槳葉表面周圍網格疏密程度對計算結果的影響，本文設計了3套計算網格（見圖9），邊界增長率分別為medium，slow，very slow，對應的旋轉域的網格數為90萬、130萬和250萬。根據2.3.1節的研究結論，在2階時間離散格式以及螺旋槳單位時間步長的旋轉角度為1.8°的情況下，對3套計算網格下的螺旋槳非定常力進行數值模擬，得到的計算結果如圖10所示。

從圖中可以看出，在給定的螺旋槳單位時間步長的旋轉角度和時間離散格式下，當槳葉表面網格的邊界增長率設定為medium時，得到的非定常力數值計算結果與試驗值吻合的較好；并且繼續對槳葉表面周圍的網格進行加密，得到的非定常力數值計算結果與medium的結果相比差別很小。因而對于螺旋槳非定常力的計算，槳葉表面網格的邊界增長率設定為medium就可以達到計算精度的要求。

圖 9 三套計算網格Fig. 9 3 sets of mesh

圖 10 網格疏密程度對非定常力的影響Fig. 10 The impact of mesh on unsteady hydrodynamics

3 全附體潛艇模型后的螺旋槳非定常力計算

基于上節的數值方法和旋轉域網格劃分形式，本節對Suboff全附體潛艇模型下3個螺旋槳（DTMB-4119，DTMB4381，DTMB4383）的非定常力進行數值模擬。Suboff全附體潛艇模型由一個水滴形軸對稱體、一個指揮臺圍殼和4個尾翼組合而成，具體的幾何參數見文獻[12]。在對自航狀態進行計算之前，本節先對無槳時的全附體潛艇尾流場進行了數值分析，得到了槳盤面r/R=0.25（R為主艇體最大半徑）處的速度分布，并與試驗值進行比較（見圖11）。從圖中可以看出，在幅值方面，各個方向上的伴流分數與試驗值吻合的較好；對于相位，從0°到180°，計算值與試驗值的相位吻合較好，而從180°到360°，數值模擬值與試驗值的相位有一定角度的偏差。由于數值模擬計算結果的對稱性較好，x方向伴流分數的峰值都與艇后舵的位置相對應，且最大值對應于指揮臺圍殼的位置，因而可以認為這二者的偏差是由試驗誤差造成的。總體而言，計算值與試驗值吻合的較好，驗證了潛艇網格劃分的合理性，為艇槳結合計算打下了基礎。

圖 11 槳盤面處的速度分布Fig. 11 Velocity distribution at propeller disk

4381和4383均為5葉槳，轂徑比均為0.2，但前者無側斜，后者的側斜角為72°。3部螺旋槳的直徑統一取為，DH為Suboff主艇體最大直徑，螺旋槳轉速統一設置為720 r/min。將計算得到的單個槳葉非定常力變化規律與劉志華[13]的計算結果進行比較如圖12所示。從圖中可以看出，在相位方面，本文計算值與其計算結果吻合較好，非定常力的變化規律與槳盤面處軸向伴流分數的分布相對應，伴流分數的低點對應非定常力的高點，伴流分數的高點對應非定常力的低點。在幅值方面，本文計算的非定常力變化幅度與劉的計算結果相比稍微有點偏大，在最低點和最高點處這種差別更明顯。因而，對于艇槳結合下螺旋槳的非定常力數值計算，迫切需要相應的試驗結果進行驗證。另外，從圖中還可以看出，在計算條件和網格劃分都相同的條件下，大側斜螺旋槳可以明顯降低單個槳葉非定常力的變化幅度。

圖13給出了3部螺旋槳無因次化的全槳非定常力隨旋轉角度的變化情況，從圖中可以直觀看出，推力的變化幅度隨著槳葉葉數的增加以及側斜的增大而降低。為了對非定常力的變化進行定量研究，本文將計算得到的全槳非定常力進行傅里葉分析，得到了1階葉頻的軸向軸承力，并與擬定常方法的結果進行比較如圖14所示。對于CFD方法，通過傅里葉分析得到的（4119，4381，4383）3部螺旋槳的軸向軸承力的1階葉頻幅值分別為2.267 N，1.56 N，0.467 N。 其中，5葉槳4381的幅值相比3葉槳4119降低了30% ，大側斜螺旋槳4383的幅值相比無側斜槳4381降低了70%。而采用擬定常經驗公式得到的不同螺旋槳的軸向軸承力變化不大，基本相當，因而擬定常方法不適應多槳葉螺旋槳尤其是大側斜螺旋槳的軸承力計算，無法滿足工程需要。

圖 12 單個槳葉的非定常力Fig. 12 Unsteady hydrodynamics of single blade

圖 13 全槳的非定常力Fig. 13 Unsteady hydrodynamics of blades

圖 14 1階葉頻軸向軸承力Fig. 14 Axial bearing force of 1-st order blade frequency

4 結 語

本文基于STAR-CCM+對非均勻來流下的螺旋槳非定常力進行了數值模擬，研究了時間離散格式、單位時間步長螺旋槳旋轉角度與網格密度對計算結果的影響。并在此基礎上，對Suboff全附體潛艇模型下3個螺旋槳的非定常力進行數值模擬。得到結論主要有：

1）時間離散格式和單位時間步長螺旋槳旋轉角度對非均勻來流下螺旋槳非定常力的影響較大，通過與試驗結果的比較分析，認為螺旋槳單位時間步長的旋轉角度應為1.8°同時選取2階時間離散格式。

2）在網格密度方面，槳葉表面網格的邊界增長率設定為medium就可以達到計算精度的要求。

3）對于全附體潛艇模型下螺旋槳的非定常力計算，槳葉側斜以及槳葉個數對非定常力的變化幅度影響較大，隨著側斜角和槳葉個數的增加，非定常力的變化幅度會大幅降低。

[ 1 ]KERWIN J E, LEE C S. Prediction of steady and unsteady marine propeller performance by numerical lifting-surface theory[J]. Societyof Naval Architects and Marine Engineers Transactions, 1978, 86: 218–253.

[ 2 ]HOSHINO T. Hydrodynamic analysis of propeller in unsteady flow using a surface panel method[J]. Journal of the Society of Naval Architects of Japan, 1993, 74: 71–87.

[ 3 ]陳家棟, 董世湯. 非定常螺旋槳表面壓力面元法預報計算[J].中國造船, 1998(1): 21–26.CHEN Jia-dong, DONG Shi-tang. Prediction of pressure distributions on unsteady propeller blades using potential-based panel method[J]. Ship Buiding of China, 1998(1): 21–26.

[ 4 ]譚延壽, 賀偉. 螺旋槳非定常軸承力計算[J]. 船海工程, 2006,171(2): 42–45.TAN Yan-shou, HE Wei. Calculation of unsteady shaft forces of propeller[J]. Ship & Ocean Engineering, 2006, 171(2):42–45.

[ 5 ]熊鷹. 非均勻流中螺旋槳空泡和脈動壓力的數值和試驗研究[D]. 武漢: 武漢理工大學, 2002.XIONG Ying. Numerical and experimental research on propeller-induced pressure fluctuations and cavitation in nonuniform flow[D]. Wuhan: Wuhan University of Technology,2002.

[ 6 ]張志榮, 洪方文. 斜流中螺旋槳性能數值分析方法研究[C]//船舶力學學術會議暨《船舶力學》創刊十周年紀念學術會議. 銀川, 2007.ZHANG Zhi-rong, HONG Fang-wen. Research on performance analysis method for propeller in oblique flow[C]//Proceedings of Ship Dynamics Meeting, 2007, 9.

[ 7 ]胡小菲, 黃振宇, 洪方文. 螺旋槳非定常力的粘性數值分析[J]. 水動力學研究與進展, 2009, 24(6): 734–739.HU Xiao-fei, HUANG Zhen-yu, HONG Fang-wen. Unsteady hydrodynamics forces of propeller predicted with viscous CFD[J]. Chinese Journal of hydrodynamics, 2009, 24(6):734–739.

[ 8 ]黃振宇. 螺旋槳的非定常力計算[C]//第九屆全國水動力學學術會議暨第二十二屆全國水動力學研討會文集. 成都, 2009.HUANG Zhen-yu. Numerical simulation of unsteady force of marine propeller[C]// Proceedings of 9th Hydrodynamics Meeting, 2009.

[ 9 ]劉登成, 洪方文, 張志榮. 等. 伴流中螺旋槳非定常力粘性數值方法研究[C]//第二十三屆全國水動力學研討會暨第十屆全國水動力學學術會議文集. 西安, 2011.LIU Deng-cheng, HONG Fang-wen, ZHANG Zhi-rong, et al.Research on viscous numerical method of propeller unsteady force in wake[C]// Proceedings of 10th Hydrodynamics Meeting, 2011.

[10]馬超, 肖鋒, 張志誼. 等. 艇體非均勻伴流場下螺旋槳壓力場[J]. 噪聲與振動控制, 2013, 35(2): 49–53.MA Chao, XIAO Feng, ZHANG Zhi-yi, et al. Numerical study on pressure pulsation near the propeller with non-uniform wake field[J]. Journal of Noise and Vibration Control, 2013,35(2): 49–53.

[11]JESSUP S. Measurement of multiple blade rate unsteady propeller force [R]. David Taylor Research Center, Ship Hydromechanics Department Research and Development Report, DTRC-90/015 May 1990.

[12]HUANG T, LIU H L, GROVES N, et al. Measurements of flows over an axisymmetric body with various appendages in a wind tunnel: the DARPA SUBOFF experimental program[C]//Proceeding of 19th Symposium on Naval Hydrodynamics. Seoul, Korea, 1992.

[13]劉志華, 熊鷹, 韓寶玉. 消渦整流片對潛艇馬蹄渦的控制及其與輔翼效果的比較[J]. 船舶力學, 2011, 15(10):1102–1105.LIU Zhi-hua, XIONG Ying, HAN Bao-yu. Comparison on the submarine horseshoe vortex control effects by vortex control bafflers and fillets[J]. Journal of Ship Mechanics, 2011, 15(10):1102–1105.