高魯棒性的螺旋槳片條理論非線性修正方法

范中允，周洲,*，祝小平，王睿，王科雷

1. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072 2. 西北工業(yè)大學(xué) 無(wú)人機(jī)特種技術(shù)重點(diǎn)實(shí)驗(yàn)室，西安 710065

垂直起降無(wú)人機(jī)成為現(xiàn)代無(wú)人機(jī)領(lǐng)域的一個(gè)重要研究方向。小型垂直起降無(wú)人機(jī)以定距螺旋槳為主要推進(jìn)裝置，螺旋槳的氣動(dòng)特性優(yōu)劣將對(duì)無(wú)人機(jī)的飛行性能有較大影響。然而，垂直起降無(wú)人機(jī)螺旋槳的設(shè)計(jì)需要兼顧懸停和巡航狀態(tài)，至少包含兩個(gè)設(shè)計(jì)點(diǎn)，同時(shí)工作狀態(tài)跨度較大，對(duì)螺旋槳的分析設(shè)計(jì)工具提出了大跨度、高魯棒性的要求。

CFD時(shí)間代價(jià)過(guò)大，難以直接用于迭代設(shè)計(jì)，片條理論(也稱葉素動(dòng)量理論，BEMT)由于其計(jì)算效率高而被廣泛用于螺旋槳設(shè)計(jì)。片條理論由來(lái)已久，目前已有許多關(guān)于其修正或改進(jìn)方法的研究。2005年，Martínez 等[1]依據(jù)4組試驗(yàn)結(jié)果提出了一種無(wú)自由參數(shù)的葉素理論修正模型。2007年，Hunsaker[2]針對(duì)垂直起降無(wú)人機(jī)的設(shè)計(jì)問(wèn)題，提出將葉素理論與計(jì)算流體力學(xué)的動(dòng)量方程相結(jié)合，對(duì)螺旋槳周圍的誘導(dǎo)流場(chǎng)進(jìn)行低計(jì)算代價(jià)的數(shù)值模擬，并去掉了小迎角假設(shè)。2009年，Dossing[3]在葉素理論中加入了尾跡旋轉(zhuǎn)和擴(kuò)張修正，以改善在槳尖和槳根的預(yù)測(cè)精度。2010年，Rwigema[4]為了考慮滑流對(duì)螺旋槳的影響，在葉素動(dòng)量理論的基礎(chǔ)上將一段渦管加入到螺旋槳后部。Vaquero等[5]以葉素動(dòng)量理論結(jié)果為參考，進(jìn)行了低雷諾數(shù)槳葉特性的試驗(yàn)研究。2012年，Whitmore和Merrill[6]將葉素動(dòng)量理論推導(dǎo)中的小迎角假設(shè)去掉并進(jìn)行了理論推導(dǎo)，以得到能夠考慮大迎角的葉素動(dòng)量理論。Burdett和Van Treuren[7]對(duì)比了葉素理論和葉素動(dòng)量理論對(duì)于槳葉扭轉(zhuǎn)角設(shè)計(jì)的影響，發(fā)現(xiàn)兩種理論的設(shè)計(jì)結(jié)果在槳轂扭轉(zhuǎn)角差別約15°，功率系數(shù)相差0.08。2013年，Ding等[8]將Wilson模型與葉素理論相結(jié)合進(jìn)行了槳葉設(shè)計(jì)。2016年，El Khchine和Srili[9]提出一種將CFD方法與葉素動(dòng)量理論方法相結(jié)合的數(shù)值模擬方法，以減小CFD的計(jì)算代價(jià)。2017年，Barnes[10]為了改善螺旋槳載荷分布的精度，將葉素理論與螺旋形尾跡相結(jié)合，并將葉素理論與升力線理論相結(jié)合以考慮后掠槳葉在跨聲速的性能。Dumakude和Kamper[11]結(jié)合軸向和徑向誘導(dǎo)因子，將葉素動(dòng)量理論和CFD的計(jì)算結(jié)果進(jìn)行了對(duì)比和分析。

盡管國(guó)內(nèi)外學(xué)者已經(jīng)對(duì)提高片條理論精度和非線性修正進(jìn)行了大量研究，但這類研究多是針對(duì)螺旋槳的常規(guī)工作段，在極端條件下的計(jì)算穩(wěn)定性問(wèn)題卻鮮有研究。同時(shí)，鑒于螺旋槳與機(jī)翼/機(jī)體的耦合設(shè)計(jì)趨勢(shì)，螺旋槳分析工具要能夠提供更多的流場(chǎng)信息，如環(huán)量分布及誘導(dǎo)速度場(chǎng)，在這方面基于片條理論的研究也有所欠缺。在極端工況下，如靜止及高前進(jìn)比狀態(tài)，實(shí)踐證明片條理論的計(jì)算穩(wěn)定性和環(huán)量預(yù)測(cè)精度都將下降，這是由于理論公式中存在著解的不確定性問(wèn)題。因此，需要對(duì)片條理論計(jì)算穩(wěn)定性和環(huán)量分布的特性進(jìn)行更加深入的研究。

另一方面，在片條理論中，葉素的氣動(dòng)特性將直接影響所得螺旋槳的氣動(dòng)特性，螺旋槳的非線性特性主要由葉素的非線性體現(xiàn)，因此在計(jì)算時(shí)需要考慮葉素的大跨度、多重非線性的因素。例如，當(dāng)在懸停狀態(tài)時(shí)，螺旋槳工作在很低前進(jìn)比下，部分葉素的當(dāng)?shù)赜呛艽笊踔烈呀?jīng)發(fā)生分離；大多數(shù)小型無(wú)人機(jī)的螺旋槳都工作在低雷諾數(shù)狀態(tài)下，而在平飛過(guò)程中螺旋槳的槳尖馬赫數(shù)也到了需要考慮壓縮性的范圍。

對(duì)于葉素氣動(dòng)模型，Carroll和Marcum[12]針對(duì)葉素分離狀態(tài)進(jìn)行了葉素動(dòng)量理論與CFD方法的結(jié)果對(duì)比。Morgado等[13]提出了一種根據(jù)翼型前緣半徑拓展大迎角翼型氣動(dòng)特性的方法，可將翼型大迎角模型范圍提高到90°。Dorfling和Rokhsaz[14]為了改進(jìn)葉素理論在非線性區(qū)域的計(jì)算精度，提出了一種用于葉素理論的翼型非線性氣動(dòng)模型，預(yù)測(cè)迎角可以達(dá)到90°，同時(shí)進(jìn)行了壓縮性修正。Kown等[15]使用二維CFD方法在特定雷諾數(shù)下建立了截面翼型升阻系數(shù)隨迎角變化的氣動(dòng)數(shù)據(jù)庫(kù),并用于葉素動(dòng)量理論。Macneill和Verstraete[16]針對(duì)小尺寸無(wú)人機(jī)螺旋槳的低雷諾數(shù)和槳葉失速問(wèn)題，對(duì)翼型建立了大范圍的迎角和雷諾數(shù)數(shù)據(jù)庫(kù)，并進(jìn)行了旋轉(zhuǎn)效應(yīng)的處理。

可以發(fā)現(xiàn)，葉素動(dòng)量理論的葉素氣動(dòng)數(shù)據(jù)，主要由兩種方法獲得：一是用簡(jiǎn)單的翼型氣動(dòng)力模型，這種方法拓展性較好，而適應(yīng)性有限，雖然大迎角特性可以較好體現(xiàn)，但是壓縮性、低雷諾數(shù)特性反映有限；二是建立翼型數(shù)據(jù)庫(kù)，可以使用CFD或者風(fēng)洞數(shù)據(jù)，這種方法在數(shù)據(jù)點(diǎn)充足的情況下準(zhǔn)確性較好，在維度較低時(shí)使用廣泛，但是倘若參數(shù)較多、非線性較強(qiáng)時(shí)，對(duì)未知函數(shù)形狀的高維插值或擬合存在困難，同時(shí)逐漸增大的數(shù)據(jù)量，使得數(shù)據(jù)存儲(chǔ)和使用也越來(lái)越復(fù)雜。因此，在考慮翼型多維度的多重非線性時(shí)，氣動(dòng)模型需要進(jìn)行改進(jìn)。

針對(duì)以上兩方面的問(wèn)題，作為翼型數(shù)據(jù)庫(kù)方法的延伸，本文首先建立人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network, ANN)對(duì)翼型不同狀態(tài)下的非線性效應(yīng)進(jìn)行特征提取，解決了翼型多維度的多重非線性氣動(dòng)特性快速預(yù)測(cè)的問(wèn)題；然后基于上述模型，分析了片條理論在極端狀態(tài)下的局限性，提出了一種環(huán)量迭代修正方法，解決了極端環(huán)境下片條理論計(jì)算穩(wěn)定性問(wèn)題和解的奇異性問(wèn)題。最后，通過(guò)試驗(yàn)數(shù)據(jù)和CFD方法，對(duì)本文提出的片條理論修正方法進(jìn)行準(zhǔn)確性及魯棒性驗(yàn)證。

1 基于神經(jīng)網(wǎng)絡(luò)的翼型氣動(dòng)分析模型

葉素氣動(dòng)特性的準(zhǔn)確性將直接影響片條理論的準(zhǔn)確性。本文利用神經(jīng)網(wǎng)絡(luò)優(yōu)異的無(wú)先驗(yàn)函數(shù)非線性擬合能力，來(lái)兼顧多種翼型非線性氣動(dòng)特性和較高的預(yù)測(cè)響應(yīng)速度。擬建立BP(Back Propagation)人工神經(jīng)網(wǎng)絡(luò)，作為翼型氣動(dòng)特性預(yù)測(cè)模型，并對(duì)其性能進(jìn)行驗(yàn)證分析。

1.1 神經(jīng)網(wǎng)絡(luò)模型

神經(jīng)網(wǎng)絡(luò)的基本原理已經(jīng)有諸多論述[17]，本文不再贅述。典型的BP神經(jīng)網(wǎng)絡(luò)構(gòu)型如圖1所示，包括一個(gè)輸入層x1,x2,…,xn、一個(gè)輸出層o1,o2，…，on和若干隱含層。適當(dāng)增加隱含層或增加隱含層結(jié)點(diǎn)可以提高神經(jīng)網(wǎng)絡(luò)對(duì)數(shù)據(jù)非線性的捕捉能力。考慮到螺旋槳葉素氣動(dòng)模型有多重非線性問(wèn)題，本文選用雙隱含層的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。

圖1 BP神經(jīng)網(wǎng)絡(luò)框架Fig.1 Framework of BP neural network

構(gòu)建形如圖1所示的4層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，網(wǎng)絡(luò)輸入量為雷諾數(shù)Re、馬赫數(shù)Ma、迎角α，輸出量為升力系數(shù)CL、阻力系數(shù)CD、力矩系數(shù)Cm。選取第一隱含層結(jié)點(diǎn)數(shù)12個(gè)，第二隱含層結(jié)點(diǎn)數(shù)16個(gè)。隱含層的傳輸函數(shù)選用雙正切Sigmoid函數(shù)[18]，輸出層為線性傳輸函數(shù)。盡管在本文的片條理論計(jì)算中沒(méi)有用到葉素力矩系數(shù)，但考慮到將來(lái)研究螺旋槳變距特性時(shí)的需要，本文仍然將力矩系數(shù)包含在神經(jīng)網(wǎng)絡(luò)結(jié)果內(nèi)進(jìn)行分析。

作為初步研究，以螺旋槳常用的經(jīng)典翼型Clark Y作為樣本翼型，均布計(jì)算樣本共616個(gè)，在迎角維度上插值加密，最終樣本點(diǎn)數(shù)目共計(jì)1 463個(gè)。選取樣本參數(shù)范圍為：雷諾數(shù)Re=0.6×105～40×105，馬赫數(shù)Ma=0.05～0.55，迎角α=-3°～15°。采用CFD方法進(jìn)行樣本特性計(jì)算，使用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格量為2.7萬(wàn)，網(wǎng)格質(zhì)量大于0.7，正交性大于0.8，選擇剪切應(yīng)力輸運(yùn)(SST)k-ω湍流模型。

1.2 神經(jīng)網(wǎng)絡(luò)性能驗(yàn)證

本節(jié)就氣動(dòng)特性預(yù)測(cè)、非線性特征提取和包線拓展3個(gè)方面的性能，對(duì)所構(gòu)建的Clark Y翼型氣動(dòng)特性分析模型進(jìn)行驗(yàn)證。

1.2.1 氣動(dòng)特性預(yù)測(cè)

通過(guò)交叉驗(yàn)證檢查發(fā)現(xiàn)，相對(duì)誤差較大的點(diǎn)一般在數(shù)值接近0或迎角大于10°(非線性較強(qiáng))的位置。為了直觀反映神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的準(zhǔn)確性，選擇雷諾數(shù)Re=3.0×105、馬赫數(shù)Ma=0.4 的狀態(tài)，預(yù)測(cè)結(jié)果隨迎角的變化曲線如圖2所示。

圖2 氣動(dòng)系數(shù)隨迎角的變化曲線 (Re=3.0×105, Ma=0.4)Fig.2 Change curves of aerodynamic coefficients with angle of attack (Re=3.0×105, Ma=0.4)

由圖2可知，升力系數(shù)和阻力系數(shù)的預(yù)測(cè)值與CFD吻合度較高，盡管神經(jīng)網(wǎng)絡(luò)對(duì)力矩系數(shù)的預(yù)測(cè)在數(shù)值上有一定誤差，但是其變化規(guī)律是一致的。

1.2.2 非線性特征提取

圖2已經(jīng)展示了神經(jīng)網(wǎng)絡(luò)對(duì)迎角非線性特征的良好捕捉能力，圖3為不同雷諾數(shù)和不同馬赫數(shù)下的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)曲線。

圖3 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)升阻極曲線Fig.3 Prediction of CL vs CD polar curves of neural network

由圖3(a)可以看出，雷諾數(shù)降低，零升阻力增大，最大升力系數(shù)減小，升阻比減小，升力線斜率略有降低。神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果與Silverstein[19]對(duì)Clark Y翼型的研究結(jié)果一致。由圖3(b)可以看出，馬赫數(shù)增加，零升阻力略有增加，升阻比降低，最大升力系數(shù)減小，隨著迎角增大馬赫數(shù)的影響變大。神經(jīng)網(wǎng)絡(luò)可以反映翼型氣動(dòng)性能在不同馬赫數(shù)下的變化特征。

1.2.3 包線拓展

在小型螺旋槳計(jì)算分析時(shí)，不可避免地會(huì)觸及計(jì)算邊界。許多小尺寸螺旋槳槳葉雷諾數(shù)常在6.0×104以下，這一邊界是最容易超越的。本節(jié)以雷諾數(shù)下限為例，說(shuō)明神經(jīng)網(wǎng)絡(luò)對(duì)外插的適應(yīng)性問(wèn)題。圖4為Re=1×104～6×104下的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)曲線。

圖4 超過(guò)雷諾數(shù)下限時(shí)的氣動(dòng)系數(shù)曲線(Ma=0.1)Fig.4 Aerodynamic coefficient curves for Re over the lower limit (Ma=0.1)

在迎角小于約8°時(shí)，預(yù)測(cè)結(jié)果趨勢(shì)仍然合理。迎角大于8°以后，大迎角非線性特征過(guò)強(qiáng)，外插預(yù)測(cè)結(jié)果出現(xiàn)震蕩，因此在本文中大迎角的超包線使用將被嚴(yán)格限制。

2 螺旋槳片條理論修正方法

Veldhuis[20]及劉沛清[21]均已經(jīng)對(duì)螺旋槳的片條理論方法進(jìn)行了詳細(xì)的推導(dǎo)。本節(jié)針對(duì)理論在應(yīng)用時(shí)的奇異性和不穩(wěn)定性，結(jié)合渦流理論進(jìn)行修正方法研究，并對(duì)修正方法有效性進(jìn)行驗(yàn)證。

2.1 基本片條理論及其局限性

片條理論是一種由二維到三維的計(jì)算方法。圖5為螺旋槳葉素的受力及速度分解圖。圖中:V0為來(lái)流速度，β為干涉角，Vt為環(huán)向誘導(dǎo)速度，Va為軸向誘導(dǎo)速度，α為實(shí)際迎角,V′為誘導(dǎo)螺距增量，φ和φ分別為幾何入流角和實(shí)際入流角，Ω為角速度，r為葉素半徑，Vgeo、Vreal分別為幾何入流速度、實(shí)際入流速度，γ為阻升角，dR、dL、dD、dT、dF分別為葉素剖面合力、升力、阻力、拉力、扭力。如圖5所示，葉素的實(shí)際迎角受誘導(dǎo)速度影響，減小了β。在給定計(jì)算條件下，圖中的β、Vt、Va是未知的。

為了求解β，結(jié)合葉素理論和動(dòng)量理論可以推導(dǎo)出隱式方程[21]：

(1)

(2)

(3)

(4)

式中：σ為實(shí)度；NB為槳葉數(shù)；b為葉素弦長(zhǎng)；ns為轉(zhuǎn)速。

由隱式方程式(1)求得干涉角后，可由文獻(xiàn)[21]中的公式求得螺旋槳的拉力、扭矩和功率。

在非線性不強(qiáng)、前進(jìn)比適中的條件下，傳統(tǒng)片條理論方法能夠較好地處理絕大多數(shù)的螺旋槳分析計(jì)算問(wèn)題。然而，極端狀態(tài)下片條理論的計(jì)算穩(wěn)定性不足。

一方面，在片條理論中，求解拉力系數(shù)、扭矩系數(shù)都是基于一個(gè)軸向速度干擾系數(shù)a：

(5)

但是當(dāng)升力系數(shù)趨于0時(shí)(CL→0)，有

tanφtan(φ+γ)→tanφtan(φ+90°)=-1

圖5 螺旋槳葉素受力分析Fig.5 Blade element forces analysis of propeller

那么軸向速度干擾系數(shù)將會(huì)趨于無(wú)窮大(a→∞)， 這將導(dǎo)致在高前進(jìn)比、拉力很小時(shí)，誘導(dǎo)速度計(jì)算結(jié)果容易出現(xiàn)奇異性。因此，若直接由傳統(tǒng)片條理論獲得誘導(dǎo)速度，其只能用于低前進(jìn)比狀態(tài)。

另一方面，片條理論在極端狀態(tài)下存在非物理解。例如，圖6為很低和很高前進(jìn)比下片條理論計(jì)算得到的環(huán)量Γ分布。其中，螺旋槳模型均為JXF16×8，高度H=400 m。

圖6(a)為來(lái)流速度V0=0.1m/s、轉(zhuǎn)速ns=948.65 r/min狀態(tài)下的環(huán)量分布，圖6(b)為來(lái)流速度V0=11.11 m/s、轉(zhuǎn)速ns=2 200 r/min狀態(tài)下的環(huán)量分布。圖中r/R表示槳葉剖面半徑與螺旋槳半徑的比值，表征剖面的徑向無(wú)量綱站位。圖6(a)所示狀態(tài)幾乎是靜止?fàn)顟B(tài)，葉素的幾何迎角很大，已經(jīng)進(jìn)入非線性區(qū)；圖6(b)所示狀態(tài)則是高來(lái)流速度、低轉(zhuǎn)速狀態(tài)，葉素的幾何迎角較小，產(chǎn)生小量的負(fù)拉力。兩圖環(huán)量分布均在某一徑向位置突變?yōu)榱?，而其他徑向位置?jì)算正常。然而該處葉素的幾何迎角并非零升迎角，會(huì)產(chǎn)生有限值的升力或環(huán)量，因此計(jì)算所得環(huán)量或升力系數(shù)與葉素所處狀態(tài)不符，更與實(shí)際的槳葉環(huán)量分布不符。

圖6 螺旋槳槳葉環(huán)量分布Fig.6 Circulation distribution of propeller blade

分析原因可以發(fā)現(xiàn)，干涉角方程式(1)的解是不唯一的。當(dāng)升力系數(shù)趨于0(CL→0)時(shí)，式(1)右邊的分母會(huì)趨于無(wú)窮((1-tanγtanβ)→∞)。這使得實(shí)際迎角為零升迎角(CL==0)的干涉角恒為方程的解之一。因此，當(dāng)方程的解使得實(shí)際迎角為零升迎角(CL==0)時(shí)，當(dāng)?shù)厝~素的環(huán)量也計(jì)算得零，而這與葉素所處的實(shí)際流動(dòng)狀態(tài)不符。這表明，該隱式方程的解是不唯一的，存在沒(méi)有物理意義的解。當(dāng)計(jì)算得到的解是非物理解時(shí)，所得到的槳葉環(huán)量分布會(huì)在某些剖面突然降為零，即得到與實(shí)際不符的結(jié)果。

這一現(xiàn)象降低了片條理論對(duì)極端狀態(tài)的計(jì)算魯棒性和可靠性。編程實(shí)踐發(fā)現(xiàn)，螺旋槳在很低或很高前進(jìn)比下更容易出現(xiàn)上述問(wèn)題。分析原因可能如下：

1) 很低前進(jìn)比下剖面翼型處于大迎角狀態(tài)，非線性較強(qiáng)，升力系數(shù)在真實(shí)干涉角解附近不單調(diào)，導(dǎo)致數(shù)值求解方法很難找到真實(shí)解，相比之下非物理解更容易找到。

2) 很高前進(jìn)比下，剖面翼型的實(shí)際升力系數(shù)很小，與零升迎角對(duì)應(yīng)的非物理解極為接近，導(dǎo)致數(shù)值求解時(shí)無(wú)法分辨，從而解出非物理解。

非物理解的干涉角正好使得當(dāng)?shù)厝~素升力為零，而升力為零對(duì)應(yīng)的誘導(dǎo)速度也為零，誘導(dǎo)速度引起的干涉角也為零，那么葉素所處的實(shí)際迎角便不為零，會(huì)產(chǎn)生有限值的升力。由此可見，非物理解在其本身是自相矛盾的，其對(duì)應(yīng)的干涉角、升力系數(shù)、誘導(dǎo)速度、實(shí)際迎角等參數(shù)之間互不匹配。

2.2 片條理論修正方法

由于非物理解導(dǎo)致各參數(shù)不匹配，考慮采用各參數(shù)相互迭代計(jì)算的方法以使之逐漸協(xié)調(diào)。因此，本節(jié)結(jié)合渦流理論相關(guān)方法，構(gòu)建了如圖7所示的環(huán)量迭代框架。首先將片條理論得到的環(huán)量作為環(huán)量初值，由渦流理論求得誘導(dǎo)速度，再由幾何關(guān)系求得干涉角，進(jìn)而得到實(shí)際迎角及升力系數(shù)，由升力系數(shù)又可以得到槳葉環(huán)量。如此反復(fù)迭代直到槳葉環(huán)量的殘差小于給定值。

圖7 環(huán)量迭代流程Fig.7 Iteration flow falt of circulation

在環(huán)量迭代框架中，渦流理論與片條理論相互結(jié)合、互為補(bǔ)充，框架仍然采用片條理論以干涉角為核心的計(jì)算思路，但計(jì)算公式借用渦流理論方法進(jìn)行推導(dǎo)。為了使框架更加完整、清晰，下文給出在迭代過(guò)程中所需的環(huán)量、誘導(dǎo)速度、誘導(dǎo)速度修正、拉力和扭矩的求解方法。

1) 環(huán)量

在不考慮葉素阻力的情況下，環(huán)量計(jì)算式為

(6)

式中：T(r)為半徑r處葉素單位展長(zhǎng)的拉力；ρ為空氣密度。由于拉力中包含阻力分量，因此當(dāng)拉力較小即剖面升阻比較低時(shí)，會(huì)引起較大的誤差。本文直接使用當(dāng)?shù)仄拭嬉硇蜕ο禂?shù)來(lái)得到環(huán)量，即

(7)

式中：Vlocal為當(dāng)?shù)貋?lái)流速度(包含誘導(dǎo)速度)。

2) 誘導(dǎo)速度

國(guó)內(nèi)外對(duì)渦流理論已有諸多論述[22]。按照渦流理論方法的推導(dǎo)，槳盤處的誘導(dǎo)速度與當(dāng)?shù)丨h(huán)量Γ的關(guān)系為[22]

軸向誘導(dǎo)速度：

(8)

環(huán)向誘導(dǎo)速度：

(9)

在公式推導(dǎo)中，假設(shè)V′沿徑向不變。這一假設(shè)僅在貝茲的螺旋槳最高效率條件[21]下才成立。為了得到變化的V′，本文根據(jù)圖5進(jìn)行了幾何分析和公式推導(dǎo)。由葉素速度分量的幾何關(guān)系可得

(10)

V′-Va=Vttanφ

(11)

代入誘導(dǎo)速度式(8)和式(9)可得

(12)

化簡(jiǎn)后得到

(13)

當(dāng)沿槳葉的環(huán)量分布已知時(shí)，即可迭代求解該隱式方程得到V′。值得一提的是，當(dāng)給定V′時(shí)，即可得到環(huán)量Γ沿半徑r的分布，而這個(gè)環(huán)量分布即為無(wú)限槳葉數(shù)、不計(jì)葉素阻力時(shí)的螺旋槳貝茲最佳環(huán)量分布。

若假設(shè)合成誘導(dǎo)速度(Va+Vt)與真實(shí)入流速度Vreal相互垂直[23]，則可以得到

(14)

由于按式(9)求解環(huán)向誘導(dǎo)速度較為簡(jiǎn)單，因此可根據(jù)式(14)求出軸向誘導(dǎo)速度Va，但是，真實(shí)入流角φ事先是未知的，式(14)無(wú)法直接使用，因此本文將用隱式方程式(13)的解作為初值求出干涉角，得到φ的初值后再使用式(14)求解。

3) Prandtl槳尖/槳根修正

由于片條理論及渦流理論均使用了槳盤平均假設(shè)，其計(jì)算所得誘導(dǎo)速度是槳盤平均值。若將Prandtl修正的機(jī)理看作修正由槳盤平均引起的誤差，那么Prandtl的槳尖/槳根修正系數(shù)可以作用在槳葉誘導(dǎo)速度上[23-24]。修正公式為[23]

(15)

(16)

(17)

(18)

4) 拉力和扭矩

根據(jù)葉素受力分析可得剖面拉力dT和扭矩dM為

dT=(ρΓVt-dDsinφ)dr

(19)

dM=(ρΓVa+dDcosφ)rdr

(20)

為了將式(19)和式(20)中的阻力分量項(xiàng)換成與環(huán)量項(xiàng)相同的形式，作如下假設(shè)：參照環(huán)量的定義，定義假想矢量“阻量”ΓD。環(huán)量矢量Γ的方向垂直于來(lái)流速度矢量V，而阻量ΓD的方向則沿著矢量V：

(21)

(22)

(23)

阻量同樣可以進(jìn)行矢量運(yùn)算，給定環(huán)量下的剖面拉力和扭矩分別為

(24)

(25)

獲得剖面的拉力和扭矩后沿槳葉積分即可得到螺旋槳拉力和扭矩。

2.3 環(huán)量迭代方法驗(yàn)證分析

為說(shuō)明環(huán)量迭代修正方法對(duì)片條理論非物理解、不穩(wěn)定性的修正效果，本節(jié)根據(jù)計(jì)算實(shí)踐中遇到的錯(cuò)誤點(diǎn)，以JXF16×8螺旋槳為對(duì)象，將槳葉由0.2R～R間隔0.05等分17個(gè)截面，選取極端狀態(tài)進(jìn)行傳統(tǒng)方法和改進(jìn)方法的對(duì)比。圖8中，BEMT是傳統(tǒng)片條理論計(jì)算結(jié)果，pBEMT是僅僅加入Prandtl修正的片條理論計(jì)算結(jié)果，rBEMT是本文提出的環(huán)量迭代修正后的結(jié)果。

圖8 槳葉環(huán)量分布對(duì)比Fig.8 Comparisons for blade circulation distribution

2.3.1 很低前進(jìn)比計(jì)算驗(yàn)證

選取很低前進(jìn)比的計(jì)算狀態(tài)：V0=0.1 m/s、轉(zhuǎn)速ns=948.65 r/min，前進(jìn)比J=0.016，高度H=400 m。圖8(a)為該狀態(tài)下修正前后的環(huán)量分布圖。在此狀態(tài)下，葉素的非線性效應(yīng)較強(qiáng)，傳統(tǒng)片條理論的葉素干涉角求解方法更容易找不到合理解。圖示狀態(tài)在徑向位置r=0.5R處得到了非物理的解，使得此處的升力系數(shù)為零，從而環(huán)量在此處突變?yōu)榱?，顯然不符合事實(shí)。通過(guò)前述環(huán)量迭代方法，干涉角逐漸與其他參數(shù)相協(xié)調(diào)，徹底消除了環(huán)量分布的突變。

2.3.2 很高前進(jìn)比計(jì)算驗(yàn)證

選取很高前進(jìn)比的計(jì)算狀態(tài)：V0=40 km/h、轉(zhuǎn)速ns=2 200 r/min，前進(jìn)比J=0.745 6，高度H=400 m。圖8(b)為該狀態(tài)下修正前后的環(huán)量分布圖。在此狀態(tài)下，剖面翼型的當(dāng)?shù)赜呛苄?，真?shí)升力系數(shù)很小，隱式方程式(1)的真實(shí)解與非物理解比較接近，數(shù)值計(jì)算更加難以分辨真實(shí)解與非物理解。圖示環(huán)量在r=0.35R處突變?yōu)榱?，而環(huán)量迭代使得干涉角與安裝角、升力系數(shù)相容，從而消除了環(huán)量的突變。

2.3.3 環(huán)量迭代對(duì)比分析

選擇計(jì)算狀態(tài)：V0=0.1 m/s、轉(zhuǎn)速ns=4 000 r/min，前進(jìn)比J=0.003 69，高度H=400 m。圖8(c)為該狀態(tài)下修正前后的環(huán)量分布圖。從修正前后的對(duì)比可以看出，傳統(tǒng)方法在槳根處由于非線性較強(qiáng)產(chǎn)生了輕微波動(dòng)，而修正后環(huán)量變得更加平滑；迭代修正后環(huán)量比修正前環(huán)量整體偏高，尤其在槳尖比單純使用Prandtl修正時(shí)偏高明顯。本文和文獻(xiàn)[15]的計(jì)算實(shí)踐表明，片條理論計(jì)算所得拉力和功率一般低于試驗(yàn)值(小拉力狀態(tài)除外)，由于環(huán)量迭代修正使得環(huán)量增大，因此修正后的計(jì)算結(jié)果更加接近實(shí)際值。

3 修正片條理論方法驗(yàn)證

為說(shuō)明本文方法在正常和極端狀態(tài)下都具有較好的準(zhǔn)確性和計(jì)算魯棒性，本節(jié)結(jié)合項(xiàng)目試驗(yàn)和CFD方法，以JXF16×8櫸木螺旋槳為主，對(duì)比說(shuō)明本文修正后的片條理論在零來(lái)流、高空高速等極端狀態(tài)下都具有較好的魯棒性。同時(shí)，為了說(shuō)明本文修正方法與幾何模型無(wú)關(guān)，對(duì)其他螺旋槳和流動(dòng)狀態(tài)差異較大的旋翼均具有普適性，選取Falcon16×10碳纖維螺旋槳和某1.5 m直徑旋翼，輔助說(shuō)明本文計(jì)算方法對(duì)不同對(duì)象均具有同樣的魯棒性和準(zhǔn)確性。

3.1 JXF16×8螺旋槳算例

1) 螺旋槳車載試驗(yàn)

JXF16×8螺旋槳是一種成熟的櫸木螺旋槳商品，結(jié)合項(xiàng)目需要進(jìn)行了螺旋槳車載試驗(yàn)。圖9為螺旋槳車載試驗(yàn)系統(tǒng)的示意圖。試驗(yàn)采用六分量天平測(cè)量螺旋槳的拉力和扭矩。數(shù)據(jù)采集器與六分量天平和計(jì)算機(jī)直接相連。電機(jī)固定在六分量天平上，天平則固定在由車廂伸出的撐桿上。試驗(yàn)來(lái)流速度范圍約為0～50 km/h，功率范圍約為0～300 W，海拔450 m。圖10和圖11為天平、數(shù)據(jù)采集器及車載試驗(yàn)的實(shí)物圖。

將試驗(yàn)獲得的離散點(diǎn)擬合為一條前進(jìn)比曲線，使用時(shí)換算得到對(duì)應(yīng)速度下的曲線。由于靜止?fàn)顟B(tài)的轉(zhuǎn)速曲線無(wú)法由前進(jìn)比曲線換算得到，因此在進(jìn)行靜止?fàn)顟B(tài)對(duì)比時(shí)直接使用試驗(yàn)數(shù)據(jù)。圖12為螺旋槳效率η和拉力系數(shù)CT試驗(yàn)離散點(diǎn)擬合示意圖。

圖9 螺旋槳車載試驗(yàn)系統(tǒng)示意圖Fig.9 Schematic of vehicle testing system for propeller

圖10 六分量天平(左)和數(shù)據(jù)采集器(右)Fig.10 Six-component balance (left) and data collector (right)

圖11 車載試驗(yàn)實(shí)物圖Fig.11 Ground vehicle test photos

2) 螺旋槳計(jì)算及結(jié)果分析

為了與商用CFD軟件進(jìn)行對(duì)比，本文同時(shí)使用Fluent對(duì)部分狀態(tài)進(jìn)行了計(jì)算。計(jì)算采用SSTk-ω湍流模型，基于結(jié)構(gòu)-非結(jié)構(gòu)混合網(wǎng)格，采用多重參考系方法進(jìn)行計(jì)算。計(jì)算網(wǎng)格中，旋轉(zhuǎn)域非結(jié)構(gòu)網(wǎng)格量為95.598 1萬(wàn)，外部靜止域結(jié)構(gòu)化網(wǎng)格量為281.147 8萬(wàn)。計(jì)算狀態(tài)選?。核俣萔0=0.406 3 m/s(靜止?fàn)顟B(tài)計(jì)算收斂預(yù)設(shè)流速)，高度H=400 m，螺旋槳轉(zhuǎn)速ns=1 000～5 000 r/min(接近試驗(yàn)狀態(tài))。圖13為螺旋槳幾何模型示意圖。

在片條理論計(jì)算中，將0.2R內(nèi)的部分作為槳轂(圖13中深色部分)，不參與片條理論計(jì)算，槳葉等分為17段。但在CFD計(jì)算中則包含整個(gè)模型。

圖12 JXF16×8試驗(yàn)數(shù)據(jù)示例Fig.12 Examples of JXF16×8 test data

圖13 JXF16×8 螺旋槳幾何模型Fig.13 Geometry model of JXF16×8 propeller

圖14為靜止?fàn)顟B(tài)下螺旋槳拉力、功率、力效的片條理論、CFD及試驗(yàn)曲線。圖15為40 km/h前進(jìn)速度下螺旋槳拉力和功率的片條理論及試驗(yàn)曲線。

在全狀態(tài)下，CFD結(jié)果及修正的片條理論與試驗(yàn)結(jié)果較為吻合，片條理論的拉力及功率相對(duì)誤差約在5%以內(nèi)。出現(xiàn)誤差的原因可能有以下幾點(diǎn)：

圖14 靜止?fàn)顟B(tài)JXF16×8螺旋槳特性曲線Fig.14 Performance curves of JXF16×8 propeller at 0 km/h

圖15 40 km/h時(shí)的JXF16×8螺旋槳特性曲線Fig.15 Performance curves of JXF16×8 propeller at 40 km/h

1) 幾何模型誤差：螺旋槳槳尖形狀和槳轂的建模與實(shí)物有一定差別，引起微小誤差。

2) 片條理論中，僅有槳葉參與計(jì)算，槳轂產(chǎn)生的影響被忽略。在轉(zhuǎn)速不大時(shí)，槳轂產(chǎn)生的阻力相對(duì)于槳葉拉力變大，因此片條理論在小轉(zhuǎn)速時(shí)的預(yù)測(cè)結(jié)果偏樂(lè)觀。

3) 試驗(yàn)狀態(tài)的槳轂后面連接有支撐設(shè)備和儀器，而CFD槳轂部分后面是懸空狀態(tài)，氣流分離產(chǎn)生了較大的壓阻，使得其計(jì)算結(jié)果較為保守。

為了進(jìn)一步驗(yàn)證本文計(jì)算方法對(duì)不同高度及高前進(jìn)比狀態(tài)的計(jì)算仍然穩(wěn)定且可靠，結(jié)合項(xiàng)目計(jì)算需求選取狀態(tài)：高度500 m(12 m/s)、3 000 m(15 m/s)、10 000 m(20、33 m/s)及轉(zhuǎn)速3 000、5 000、7 000 r/min，其他CFD計(jì)算條件與本節(jié)靜止?fàn)顟B(tài)計(jì)算條件類似。CFD與rBEMT的計(jì)算對(duì)比結(jié)果如圖16所示。

圖16 JXF16×8螺旋槳不同高度下的拉力和功率Fig.16 Thrust and power of JXF16×8 propeller at different height

經(jīng)過(guò)修正的片條理論結(jié)果與CFD的計(jì)算結(jié)果吻合較好，相對(duì)誤差普遍低于5%。僅在10 000 m高度、33m/s狀態(tài)下的計(jì)算偏離較大，這是由于在該狀態(tài)下螺旋槳已經(jīng)處于深度失速的風(fēng)車狀態(tài)，部分葉素處于-7°以下的負(fù)迎角狀態(tài)，其氣動(dòng)數(shù)據(jù)是神經(jīng)網(wǎng)絡(luò)外插得到的，因此造成了一定誤差。

3.2 其他螺旋槳模型算例

為了說(shuō)明本文修正方法的幾何無(wú)關(guān)性，本節(jié)選用特性不同的Falcon16×10碳纖維螺旋槳和流動(dòng)狀態(tài)差別較大的平直旋翼來(lái)驗(yàn)證本文的計(jì)算模型對(duì)不同螺旋槳模型具有同樣的計(jì)算穩(wěn)定性和準(zhǔn)確性。

1) Falcon16×10螺旋槳算例

Falcon16×10碳纖維螺旋槳是一種著名的高效碳纖維螺旋槳商品，其車載試驗(yàn)方法與3.1節(jié)所述相同。計(jì)算狀態(tài)為：前進(jìn)速度V0=40 km/h，高度H=400 m，槳葉等分為17個(gè)截面，0.2R以內(nèi)作為槳轂不參與計(jì)算。圖17為片條理論與車載試驗(yàn)的拉力和功率曲線。

在n<6 500 r/min(槳盤載荷約為230 N/m2)時(shí)，計(jì)算結(jié)果仍然較為吻合，拉力相對(duì)誤差保持在5%以內(nèi)，功率相對(duì)誤差在8%以內(nèi)。

2) 靜止?fàn)顟B(tài)旋翼算例

旋翼幾何外形、工作狀態(tài)、流動(dòng)狀態(tài)與螺旋槳有較大差異，構(gòu)建一個(gè)直徑D=1.5 m、弦長(zhǎng)c=0.14 m、安裝角i=7°、翼型為DAE-21的簡(jiǎn)單旋翼，0.2R以內(nèi)不包含槳轂，如圖18所示。計(jì)算狀態(tài)為：來(lái)流速度V0=0.4 m/s，高度H=400 m。CFD湍流模型選為層流動(dòng)能轉(zhuǎn)捩(kkl)模型，其他計(jì)算條件與前文靜止?fàn)顟B(tài)JXF16×8相似。

圖19為片條理論與CFD計(jì)算拉力和功率曲線對(duì)比。由圖可知，片條理論計(jì)算結(jié)果與CFD吻合良好，二者之間的拉力和功率相對(duì)誤差均低于5%。

圖17 40 km/h下的Falcon16×10螺旋槳的 拉力和功率Fig.17 Thrust and power of Falcon16×10 propeller at 40 km/h

圖18 DAE-21旋翼幾何模型Fig.18 Geometry model of DAE-21 rotor

圖19 DAE-21旋翼的拉力和功率曲線Fig.19 Thrust and power curves of of DAE-21 rotor

4 結(jié) 論

1) 應(yīng)用BP神經(jīng)網(wǎng)絡(luò)建立翼型的氣動(dòng)分析模型，可以為片條理論提供高效的非線性翼型氣動(dòng)特性預(yù)測(cè)。

2) 傳統(tǒng)片條理論方程存在非物理解，導(dǎo)致在極端狀態(tài)下的計(jì)算魯棒性下降。本文提出的環(huán)量迭代修正方法彌補(bǔ)了這一不足，消除了環(huán)量分布突變錯(cuò)誤，極大提高計(jì)算穩(wěn)定性和魯棒性。

3) 修正后片條理論適用于多種螺旋槳模型和工況，具有普遍意義。針對(duì)本文計(jì)算模型和狀態(tài)，計(jì)算拉力和功率的相對(duì)誤差在常規(guī)工作段低于5%，在很高/很低前進(jìn)比等極端狀態(tài)下本文方法仍能給出相對(duì)誤差為10%以內(nèi)的可靠結(jié)果。