帶前置定子導(dǎo)管槳非定常力特性計(jì)算分析

王 磊，唐登海，劉登成

（中國(guó)船舶科學(xué)研究中心，江蘇無錫214082）

0 引 言

本文研究的帶前置定子導(dǎo)管槳由環(huán)狀導(dǎo)管、定子和轉(zhuǎn)子組成，具有推進(jìn)效率高、激振力低、臨界航速高等優(yōu)勢(shì)[1]，該類型推進(jìn)器常用于潛器、AUV等水下航行器、魚雷和潛艇等。

與敞開式螺旋槳相比，帶前置定子導(dǎo)管槳由于轉(zhuǎn)定子時(shí)刻處于相對(duì)運(yùn)動(dòng)之中，因而不可避免會(huì)產(chǎn)生非定常力，加之推進(jìn)器運(yùn)行于船艉的非均勻尾流之中，這使推進(jìn)器受到的非定常力更加復(fù)雜。推進(jìn)器非定常力一方面通過軸系激勵(lì)船體艉部產(chǎn)生振動(dòng)，另一方面還直接輻射噪聲，干擾航行體有效執(zhí)行任務(wù)。因此，有必要對(duì)帶前置定子導(dǎo)管槳的非定常力進(jìn)行研究，以便為設(shè)計(jì)低激振的帶前置定子導(dǎo)管槳提供依據(jù)。

針對(duì)帶前置定子導(dǎo)管槳的研究，目前大多集中于其定常水動(dòng)力性能，非定常力研究尚有一定的局限性。在勢(shì)流理論研究方面，王國(guó)強(qiáng)[2]使用升力面-面元耦合方法計(jì)算了帶前置定子導(dǎo)管槳的定常水動(dòng)力性能；劉小龍[3]使用面元法計(jì)算了帶前置定子導(dǎo)管槳的非定常性能，但未對(duì)非均勻來流對(duì)推進(jìn)器非定常力特性影響規(guī)律進(jìn)行計(jì)算分析。隨著計(jì)算機(jī)性能的提高和CFD 技術(shù)的進(jìn)步，采用求解RANS方程等粘流方法對(duì)帶前置定子導(dǎo)管槳進(jìn)行性能預(yù)報(bào)和流場(chǎng)計(jì)算成為可能。潘光[4]使用RANS 方法對(duì)后置定子的導(dǎo)管螺旋槳進(jìn)行了定常數(shù)值模擬，得到了其敞水特性曲線；洪方文[5]對(duì)帶前置定子導(dǎo)管槳進(jìn)行了單流道定常計(jì)算，得到了轉(zhuǎn)定子表面的定常壓力分布特性；劉登成[6]對(duì)帶前置定子導(dǎo)管槳進(jìn)行了定常數(shù)值預(yù)報(bào)，研究了導(dǎo)管間隙對(duì)導(dǎo)管槳時(shí)均水動(dòng)力性能的影響；饒志強(qiáng)[7]對(duì)均勻來流下泵噴推進(jìn)器非定常力頻率進(jìn)行了理論推導(dǎo)，分析了轉(zhuǎn)定子葉數(shù)與非定常力頻率的關(guān)系，并采用RANS 方法進(jìn)行了非定常計(jì)算，結(jié)果與理論推導(dǎo)吻合良好；張志榮等[8]對(duì)帶前置定子導(dǎo)管槳的梢隙流動(dòng)進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)梢隙流動(dòng)中存在4 種類型的渦結(jié)構(gòu)。在實(shí)驗(yàn)研究方面，舒禮偉等[9]在循環(huán)水槽中對(duì)帶前置定子導(dǎo)管槳導(dǎo)管脈動(dòng)壓力進(jìn)行了實(shí)驗(yàn)測(cè)量，得到了脈動(dòng)壓力在導(dǎo)管周向和軸向的分布規(guī)律；Suranarayana[10-11]對(duì)軸對(duì)稱水下航行體后泵噴推進(jìn)器開展了水動(dòng)力和空泡性能實(shí)驗(yàn)研究。

本文采用CFD 商用軟件Fluent 對(duì)帶前置定子導(dǎo)管槳非定常力特性開展研究。首先計(jì)算9 階來流下P4119槳的非定常力，驗(yàn)證所用計(jì)算方法的準(zhǔn)確性，在此基礎(chǔ)上以非均勻來流周期數(shù)和轉(zhuǎn)定子葉數(shù)為自變量，系統(tǒng)地制定了計(jì)算方案，使用RANS 方法和滑移網(wǎng)格技術(shù)對(duì)帶前置定子導(dǎo)管槳進(jìn)行了非定常性能模擬，并使用快速傅里葉變換的手段提取轉(zhuǎn)子非定常力的頻域信息，分析轉(zhuǎn)子的非定常力與非均勻來流及轉(zhuǎn)定子數(shù)目的關(guān)系。

1 數(shù)值計(jì)算方法及驗(yàn)證

本文采用RANS 方法計(jì)算推進(jìn)器非定常性能，為了對(duì)所采用數(shù)值計(jì)算方法的準(zhǔn)確性進(jìn)行驗(yàn)證，選取P4119 槳為計(jì)算模型，如圖1 所示。P4119 槳是ITTC 確定的一個(gè)螺旋槳標(biāo)準(zhǔn)模型，被廣泛用于螺旋槳性能研究和CFD計(jì)算結(jié)果驗(yàn)證。美國(guó)泰勒水池的Jessup[12]曾通過實(shí)驗(yàn)測(cè)量了P4119 槳在3、6、9、12周期的非均勻來流下所受到的非定常力。本文對(duì)9 周期非均勻來流下P4119 槳的非定常力進(jìn)行了計(jì)算，所用的P4119 槳的尺寸與實(shí)驗(yàn)相同，其直徑D為304.8 mm。

本文采用商業(yè)軟件Fluent 進(jìn)行數(shù)值計(jì)算，動(dòng)量方程對(duì)流項(xiàng)采用二階迎風(fēng)格式離散，擴(kuò)散項(xiàng)采用中心差分格式離散，壓力-速度耦合采用SIMPLE算法，湍流模式選擇SSTk-ω模型。

計(jì)算域?yàn)橐粋€(gè)與螺旋槳共軸的圓柱體，進(jìn)口位于槳盤面上游3D處，出口位于槳盤面下游10D處，圓柱體直徑為10D。計(jì)算域劃分為外部的靜止域和螺旋槳周圍的旋轉(zhuǎn)域。網(wǎng)格劃分如圖2所示，靜止域采用結(jié)構(gòu)網(wǎng)格劃分，網(wǎng)格數(shù)為277 萬；旋轉(zhuǎn)域采用4 面體非結(jié)構(gòu)網(wǎng)格劃分，網(wǎng)格數(shù)為131 萬，第一層網(wǎng)格的y+在20-200范圍內(nèi)。

邊界條件設(shè)定為：進(jìn)口和外圍圓柱面均定義為速度入口，出口定義為壓力出口，物面為無滑移壁面，靜止域與旋轉(zhuǎn)域的交界面設(shè)置為interface，采用滑移網(wǎng)格技術(shù)處理槳的轉(zhuǎn)動(dòng)，螺旋槳轉(zhuǎn)速為600 r/min。

圖1 P4119槳Fig.1 Propeller P4119

圖2 靜止域（左）和旋轉(zhuǎn)域（右）網(wǎng)格Fig.2 Grids of stationary domain(left)and rotational domain(right)

非均勻來流在進(jìn)口處通過PROFILE定義，其平均速度為2.54 m/s。非均勻來流只有軸向分量，徑向分量和切向分量為零。軸向分量具有明顯的九周期特性，圖3 給出了9 周期非均勻來流的無量綱速度Vx/Vin（以平均來流速度無量綱化）沿周向的分布。

計(jì)算時(shí)先采用多參考系（MRF）方法進(jìn)行定常計(jì)算，收斂后以定常計(jì)算得到的結(jié)果為初值，進(jìn)行非定常計(jì)算，時(shí)間步長(zhǎng)為0.000 25 s，即每個(gè)周期包含400 個(gè)時(shí)間步，每個(gè)時(shí)間步螺旋槳旋轉(zhuǎn)0.9°，共計(jì)算1 200 步，即螺旋槳旋轉(zhuǎn)3 周。圖4 給出了計(jì)算結(jié)果與Jessup實(shí)驗(yàn)結(jié)果的對(duì)比，可見計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。

圖3 九周期非均勻來流的速度分布Fig.3 Velocity distribution of the ninth-harmonic incoming flow

圖4 P4119槳非定常推力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的比較Fig.4 Comparison between calculated and measured unsteady forces of P4119 propeller

2 帶前置定子導(dǎo)管槳非定常力特性分析

2.1 研究對(duì)象

圖5 帶前置定子導(dǎo)管槳Fig.5 Ducted propeller with a pre-swirl stator

本文的研究對(duì)象是一個(gè)帶前置定子導(dǎo)管槳，如圖5所示。導(dǎo)管長(zhǎng)為119 mm，定子弦長(zhǎng)為39.4 mm，轉(zhuǎn)子直徑D為200 mm，0.7 半徑處螺距比（P/D）0.7R為1.274 9，轉(zhuǎn)子與導(dǎo)管之間的間隙為1 mm。

帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力主要是由于船艉非均勻尾流和定子尾流的影響而產(chǎn)生的，對(duì)船艉非均勻流場(chǎng)進(jìn)行傅里葉分析，發(fā)現(xiàn)尾流場(chǎng)中包含不同周期數(shù)的脈動(dòng)量，加之定子尾流的影響，導(dǎo)致轉(zhuǎn)子進(jìn)流場(chǎng)周向不均勻。這樣，在轉(zhuǎn)子葉片旋轉(zhuǎn)過程中，進(jìn)流大小和攻角均是周期性變化的，結(jié)果導(dǎo)致轉(zhuǎn)子受力也是周期性變化的。本文主要研究帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力與非均勻來流周期數(shù)和轉(zhuǎn)定子葉數(shù)之間的關(guān)系，為此設(shè)計(jì)了包括多個(gè)工況的計(jì)算方案，如表1 所示。對(duì)于不同工況，通過改變來流平均速度來保證帶前置定子導(dǎo)管槳的總推力系數(shù)相同。

表1 帶前置定子導(dǎo)管槳計(jì)算方案Tab.1 Calculation proposals with different parameters

2.2 計(jì)算方法

計(jì)算域?yàn)榕c導(dǎo)管槳共軸的圓柱體，進(jìn)口位于導(dǎo)管進(jìn)口上游3D處，出口位于導(dǎo)管進(jìn)口下游10D處，圓柱體直徑為10D。計(jì)算域劃分為3 部分：外部流場(chǎng)域、定子域和轉(zhuǎn)子域，前兩者為靜止域，轉(zhuǎn)子域?yàn)樾D(zhuǎn)域。

網(wǎng)格采用混合網(wǎng)格。對(duì)于外部流場(chǎng)域和定子域劃分H 型結(jié)構(gòu)網(wǎng)格，如圖6所示，對(duì)于轉(zhuǎn)子域采用四面體非結(jié)構(gòu)網(wǎng)格劃分。定子域和轉(zhuǎn)子域均劃分單流道網(wǎng)格。定子、轉(zhuǎn)子和導(dǎo)管表面均劃分邊界層網(wǎng)格，第一層距壁面約0.2 mm，y+在10-110范圍內(nèi)。推進(jìn)器各部分網(wǎng)格如圖7所示。

圖6 外部流場(chǎng)域網(wǎng)格分布 Fig.6 Grid distribution of outer domain

數(shù)值離散方法、湍流模式及邊界條件設(shè)置與P4119 槳非定常力計(jì)算一致，其中轉(zhuǎn)子轉(zhuǎn)速設(shè)定為1 800 r/min，各階軸向非均勻進(jìn)流分量以周向正弦分布形式定義，具體表達(dá)式為

圖7 帶前置定子導(dǎo)管槳表面網(wǎng)格Fig.7 Surface grid of ducted propeller with stator

式中，V0為軸向速度分量均值，A（A= 0.1）表示軸向速度脈動(dòng)幅值，B表示指定的進(jìn)流周期數(shù)，θ為周向角度，范圍為0°～360°。非均勻進(jìn)流的徑向和切向分量為零。

計(jì)算時(shí)首先使用MRF 方法進(jìn)行定常計(jì)算，收斂后以定常計(jì)算的結(jié)果為初值進(jìn)行非定常計(jì)算，非定常計(jì)算的時(shí)間步為8.333 33E-5 s，即每個(gè)周期包含400 個(gè)時(shí)間步，轉(zhuǎn)子每步旋轉(zhuǎn)0.9°，計(jì)算1 200 步即轉(zhuǎn)子旋轉(zhuǎn)3周。

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

2.3.1 不同周期數(shù)來流下帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力特性分析

為了研究非均勻來流周期數(shù)對(duì)帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力的影響，選擇轉(zhuǎn)子數(shù)為6，定子數(shù)為8 的導(dǎo)管槳進(jìn)行研究，計(jì)算了4、5、6、7、8、11、12、13 周期非均勻來流下轉(zhuǎn)子的非定常力。當(dāng)來流周期數(shù)為6 時(shí)，轉(zhuǎn)子推力系數(shù)在一個(gè)周期內(nèi)的變化如圖8 所示，呈現(xiàn)明顯的6 周期特性。對(duì)其進(jìn)行傅里葉分析，結(jié)果如圖9所示，縱軸為相對(duì)幅值，即轉(zhuǎn)子推力脈動(dòng)量幅值與推進(jìn)器總推力平均值的千分比。結(jié)果表明，轉(zhuǎn)子推力在葉頻處出現(xiàn)明顯的峰值，其它頻率處的幅值可以忽略不計(jì)。

圖8 六周期來流下轉(zhuǎn)子非定常推力系數(shù)計(jì)算結(jié)果Fig.8 Calculated unsteady thrust coefficient of rotor in sixth-harmonic incoming flow

圖9 六周期來流下轉(zhuǎn)子非定常推力FFT分析結(jié)果Fig.9 Harmonic components of rotor thrust in sixth-harmonic incoming flow

對(duì)不同周期數(shù)非均勻來流下的計(jì)算結(jié)果，提取轉(zhuǎn)子的三向非定常力，并進(jìn)行快速傅里葉分析，表2為轉(zhuǎn)子三向非定常力相對(duì)幅值。

表2 不同周期數(shù)來流下轉(zhuǎn)子非定常力相對(duì)幅值（單位：‰）Tab.2 Unsteady forces of rotor with different inflows(unit:‰)

可以發(fā)現(xiàn)，對(duì)于推力，6周期來流產(chǎn)生1倍葉頻的非定常推力，12周期來流產(chǎn)生2倍葉頻的非定常推力，其它工況下非定常推力可以忽略不計(jì)；對(duì)于側(cè)向力，5 周期和7 周期來流產(chǎn)生1 倍葉頻的非定常側(cè)向力，11 周期和13 周期來流產(chǎn)生2 倍葉頻的非定常側(cè)向力，其他工況下非定常側(cè)向力可忽略不計(jì)，這與螺旋槳的非定常力特性是一致的。由表2 還發(fā)現(xiàn)，7 周期伴流產(chǎn)生的轉(zhuǎn)子非定常側(cè)向力要小于5周期伴流產(chǎn)生的轉(zhuǎn)子非定常側(cè)向力，13周期伴流產(chǎn)生的轉(zhuǎn)子非定常側(cè)向力要小于11周期伴流產(chǎn)生的轉(zhuǎn)子非定常側(cè)向力，這是因?yàn)榉蔷鶆騺砹髦芷跀?shù)越多，相鄰伴流峰周向間距就越小，轉(zhuǎn)子葉片相對(duì)于伴流峰的相對(duì)側(cè)斜角就越大，導(dǎo)致非定常力幅值偏小。由于同樣的原因，12周期來流產(chǎn)生的轉(zhuǎn)子2倍葉頻非定常推力小于6周期來流產(chǎn)生的轉(zhuǎn)子1倍葉頻非定常推力，11周期和13周期來流產(chǎn)生的轉(zhuǎn)子2倍葉頻非定常側(cè)向力也小于5 周期和7 周期來流產(chǎn)生的轉(zhuǎn)子1 倍葉頻非定常側(cè)向力。此外，還注意到，5、7、11、13周期來流各自產(chǎn)生的轉(zhuǎn)子Y、Z兩個(gè)方向非定常力幅值相同。

2.3.2 不同定子葉片數(shù)目時(shí)轉(zhuǎn)子非定常力特性分析

為了研究轉(zhuǎn)定子葉片數(shù)目對(duì)帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力的影響，選擇轉(zhuǎn)子葉片數(shù)為7，定子葉片數(shù)分別為6、7、8、9的導(dǎo)管槳進(jìn)行研究。

由于前置定子會(huì)影響轉(zhuǎn)子的實(shí)際進(jìn)流場(chǎng),進(jìn)而影響轉(zhuǎn)子非定常力，因此首先對(duì)轉(zhuǎn)子進(jìn)流場(chǎng)進(jìn)行分析。以7 葉定子為例，計(jì)算均勻來流下導(dǎo)管+7 葉定子模型流場(chǎng)，取定子近后方流場(chǎng)作為轉(zhuǎn)子進(jìn)流場(chǎng)，對(duì)0.7R半徑處的軸向速度進(jìn)行傅里葉分析（R為轉(zhuǎn)子半徑），結(jié)果如圖10 所示，縱軸為軸向相對(duì)速度的脈動(dòng)幅值。由圖可知，由于7葉定子的影響，轉(zhuǎn)子進(jìn)流場(chǎng)存在7、14、21等7的整數(shù)倍周期的脈動(dòng)，其中7周期脈動(dòng)幅度最大。除7的整數(shù)倍周期外，其它周期數(shù)的脈動(dòng)幅度都很小。

對(duì)轉(zhuǎn)子葉片數(shù)為7，定子葉片數(shù)分別為6、7、8、9 的導(dǎo)管槳，計(jì)算其在均勻來流下的非定常力，表3為轉(zhuǎn)子非定常力相對(duì)幅值。

圖10 轉(zhuǎn)子進(jìn)流場(chǎng)0.7R半徑處軸向速度FFT分析結(jié)果Fig.10 Harmonic components of axial velocity of rotor inflow at 0.7R

表3 不同定子葉片數(shù)目時(shí)轉(zhuǎn)子非定常力葉頻分量相對(duì)幅值（單位：‰）Tab.3 Unsteady forces of rotor blade frequency components with different stator blade numbers(unit:‰)

對(duì)于推力，當(dāng)定子數(shù)為7 葉時(shí)，轉(zhuǎn)子在1 倍、2 倍和3 倍葉頻處均產(chǎn)生非定常推力，這是由7 葉定子產(chǎn)生的7、14、21周期尾流導(dǎo)致的，其它工況下轉(zhuǎn)子非定常推力可忽略不計(jì);對(duì)于側(cè)向力，只有當(dāng)定子葉數(shù)為6 葉和8 葉時(shí)，轉(zhuǎn)子產(chǎn)生非定常側(cè)向力，這是由6 葉和8 葉定子產(chǎn)生的6 周期和8 周期尾流所導(dǎo)致的。同時(shí)，由表3 中還可以發(fā)現(xiàn)，定子數(shù)為8 葉時(shí)轉(zhuǎn)子的非定常側(cè)向力要小于定子數(shù)為6 葉時(shí)轉(zhuǎn)子的非定常側(cè)向力，而這兩種工況下各自的轉(zhuǎn)子Y、Z兩個(gè)方向非定常力的幅值是相同的。

2.3.3 非均勻來流和轉(zhuǎn)定子葉片數(shù)目共同作用下轉(zhuǎn)子非定常力特性分析

為了研究非均勻來流和轉(zhuǎn)定子葉片數(shù)目共同作用下帶前置定子導(dǎo)管槳轉(zhuǎn)子非定常力特性，選擇轉(zhuǎn)子葉片數(shù)為7、定子葉片數(shù)為8的導(dǎo)管槳，計(jì)算其在5、6、7、8、9周期非均勻來流和均勻來流下的非定常力。表4為轉(zhuǎn)子非定常力在葉頻處的相對(duì)幅值。

表4 不同周期數(shù)來流下轉(zhuǎn)子非定常力葉頻分量相對(duì)幅值（單位：‰）Tab.4 Unsteady forces of rotor blade frequency components with different inflows(unit:‰)

對(duì)于推力，只有7 周期來流下轉(zhuǎn)子產(chǎn)生非定常推力；而對(duì)于側(cè)向力，表4 中各個(gè)工況下，轉(zhuǎn)子均產(chǎn)生非定常側(cè)向力，這是由于8葉定子產(chǎn)生的8周期尾流的影響。此外，5、7、9周期來流下轉(zhuǎn)子非定常側(cè)向力幅值與均勻來流下相同，說明此時(shí)轉(zhuǎn)子側(cè)向力是由于定子的作用產(chǎn)生的。對(duì)于6周期和8周期來流，此時(shí)轉(zhuǎn)子同時(shí)受非均勻來流和8 葉定子產(chǎn)生的8 周期尾流共同作用，導(dǎo)致轉(zhuǎn)子非定常側(cè)向力與均勻來流下轉(zhuǎn)子非定常側(cè)向力不同。

3 結(jié) 論

本文以非均勻來流下的帶前置定子導(dǎo)管槳為研究對(duì)象，計(jì)算了多種轉(zhuǎn)定子葉數(shù)、不同周期數(shù)來流下的轉(zhuǎn)子非定常力，并重點(diǎn)分析了轉(zhuǎn)子非定常力的影響因素，得到了以下結(jié)論：

（1）CFD 計(jì)算結(jié)果表明，帶前置定子導(dǎo)管槳轉(zhuǎn)子的非定常力由轉(zhuǎn)子進(jìn)流決定，這與敞開式螺旋槳一致，但轉(zhuǎn)子進(jìn)流同時(shí)受到非均勻來流和定子尾流的影響；

（2）來流經(jīng)過定子后，轉(zhuǎn)子進(jìn)流中會(huì)增加定子葉數(shù)整數(shù)倍周期的脈動(dòng)；

（3）對(duì)于葉數(shù)為ZR的轉(zhuǎn)子，設(shè)k為正整數(shù)，進(jìn)流中的kZR階脈動(dòng)分量引起kZR階非定常推力；

（4）對(duì)于葉數(shù)為ZR的轉(zhuǎn)子，進(jìn)流中的kZR±1 階分量使轉(zhuǎn)子產(chǎn)生kZR階非定常側(cè)向力，且kZR+ 1階分量引起的側(cè)向力幅度小于kZR- 1階分量引起的側(cè)向力幅度；

（5）對(duì)于葉數(shù)為ZR的轉(zhuǎn)子，當(dāng)進(jìn)流中只存在kZR±1 階分量中的一個(gè)時(shí)，轉(zhuǎn)子Y、Z兩個(gè)方向的非定常力幅度相同；

（6）為減小非定常力，設(shè)計(jì)帶前置定子導(dǎo)管槳時(shí)，除了需要考慮非均勻進(jìn)流的影響外，還要考慮轉(zhuǎn)定子葉數(shù)相互適配，避免mZS=kZR和mZS±1 =kZR的轉(zhuǎn)定子葉數(shù)組合出現(xiàn)（ZS、ZR為定子和轉(zhuǎn)子葉數(shù)，k、m為正整數(shù)）。

本文計(jì)算了一系列工況下帶前置定子導(dǎo)管槳轉(zhuǎn)子的非定常力，但僅分析了轉(zhuǎn)子的非定常力特性。對(duì)于定子和導(dǎo)管，二者實(shí)際運(yùn)行于轉(zhuǎn)子轉(zhuǎn)動(dòng)產(chǎn)生的周期性變化的流場(chǎng)之中，也會(huì)受到非定常力作用，其非定常力特性將在后續(xù)文章中進(jìn)行分析。