不均勻進(jìn)流場下螺旋槳非定常力的數(shù)值模擬

張凱奇，葉金銘，于安斌

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

0 引 言

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

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

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

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

1.1 控制方程

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

1.2 湍流模型

顯然，上面的控制方程不封閉，因而要對流場進(jìn)行求解還需要選取適當(dāng)?shù)耐牧髂Ｐ汀Ｍㄟ^對相關(guān)文獻(xiàn)的閱讀，發(fā)現(xiàn)對于螺旋槳非定常水動力性能的數(shù)值模擬，研究人員較多的采用SSTk-ω湍流模型。因此，本文選取SSTk-ω湍流模型對方程組中的雷諾應(yīng)力項進(jìn)行封閉。

該模型由Menter提出，又稱為剪切應(yīng)力輸運(yùn)kω模型，它綜合了邊界層外部k-ε模型獨(dú)立性和近壁kω模型穩(wěn)定性的優(yōu)點(diǎn)。在近壁面區(qū)域采用k-ω湍流模式，在邊界層外部采用k-ε湍流模式。該模型的k方程和ω方程分別為：

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

2.1 計算對象與計算域

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

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

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

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

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

2.2 邊界條件與網(wǎng)格劃分

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

圖 3 DTRC 4119槳的進(jìn)流條件Fig. 3 The inflow conditions of DTRC 4119

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

圖 4 入口邊界網(wǎng)格劃分Fig. 4 Inlet boundary mesh

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

圖 6 對稱面上的網(wǎng)格分布Fig. 6 Mesh distribution of symmetry plane

2.3 計算結(jié)果與分析

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

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

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

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

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

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

2.3.2 網(wǎng)格密度對螺旋槳非定常力的影響

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

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

圖 9 三套計算網(wǎng)格Fig. 9 3 sets of mesh

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

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

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

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

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

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

圖 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 結(jié) 語

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

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

2）在網(wǎng)格密度方面，槳葉表面網(wǎng)格的邊界增長率設(shè)定為medium就可以達(dá)到計算精度的要求。

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

[ 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 ]陳家棟, 董世湯. 非定常螺旋槳表面壓力面元法預(yù)報計算[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 ]熊鷹. 非均勻流中螺旋槳空泡和脈動壓力的數(shù)值和試驗研究[D]. 武漢: 武漢理工大學(xué), 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 ]張志榮, 洪方文. 斜流中螺旋槳性能數(shù)值分析方法研究[C]//船舶力學(xué)學(xué)術(shù)會議暨《船舶力學(xué)》創(chuàng)刊十周年紀(jì)念學(xué)術(shù)會議. 銀川, 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 ]胡小菲, 黃振宇, 洪方文. 螺旋槳非定常力的粘性數(shù)值分析[J]. 水動力學(xué)研究與進(jìn)展, 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]//第九屆全國水動力學(xué)學(xué)術(shù)會議暨第二十二屆全國水動力學(xué)研討會文集. 成都, 2009.HUANG Zhen-yu. Numerical simulation of unsteady force of marine propeller[C]// Proceedings of 9th Hydrodynamics Meeting, 2009.

[ 9 ]劉登成, 洪方文, 張志榮. 等. 伴流中螺旋槳非定常力粘性數(shù)值方法研究[C]//第二十三屆全國水動力學(xué)研討會暨第十屆全國水動力學(xué)學(xué)術(shù)會議文集. 西安, 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]. 船舶力學(xué), 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.