王紅波， 祝小平， 周洲,*， 許曉平

1.西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

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

基于非定常面元/黏性渦粒子法的低雷諾數(shù)滑流氣動(dòng)干擾

王紅波1,2， 祝小平2， 周洲1,2,*， 許曉平1,2

1.西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

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

針對(duì)太陽(yáng)能無(wú)人機(jī)螺旋槳滑流與機(jī)翼的氣動(dòng)干擾，考慮了低雷諾數(shù)流動(dòng)下氣體黏性和壓縮性影響，并根據(jù)黎曼邊界條件和渦量等效原則建立了能夠快速計(jì)算分析螺旋槳-機(jī)翼氣動(dòng)干擾的非定常面元/黏性渦粒子的混合方法。首先使用有試驗(yàn)數(shù)據(jù)的風(fēng)洞模型以及數(shù)值模擬技術(shù)對(duì)混合方法進(jìn)行驗(yàn)證，在此基礎(chǔ)上研究了不同安裝位置與工況下螺旋槳與機(jī)翼的氣動(dòng)干擾。結(jié)果表明：螺旋槳對(duì)軸向氣流的加速以及滑流誘導(dǎo)的上洗和下洗效應(yīng)使機(jī)翼氣動(dòng)力呈現(xiàn)出增升增阻的現(xiàn)象，機(jī)翼升阻比有所下降。較大的弦向間距以及較高的垂直安裝位置在減緩機(jī)翼升阻比下降的同時(shí)也使得螺旋槳拉力有所減小。對(duì)于多個(gè)螺旋槳的氣動(dòng)干擾，不同的槳葉旋轉(zhuǎn)方向?qū)е聶C(jī)翼氣動(dòng)力不同的變化規(guī)律，當(dāng)旋轉(zhuǎn)方向與機(jī)翼翼尖渦反向時(shí)，螺旋槳滑流能夠抑制翼尖渦的強(qiáng)度，提高機(jī)翼氣動(dòng)效率。

低雷諾數(shù)滑流； 非定常面元； 黏性渦粒子； 混合方法； 氣動(dòng)干擾

隨著全球石油問題日益加劇，以清潔太陽(yáng)能為動(dòng)力能源的太陽(yáng)能飛機(jī)將是未來最有發(fā)展前景的一類飛行器。尤其是2015年陽(yáng)光動(dòng)力2號(hào)太陽(yáng)能有人飛行器開啟的首次環(huán)球飛行之旅更是給公眾留下了深刻的印象。

由于現(xiàn)有太陽(yáng)能電池片的光電轉(zhuǎn)換效率和電池能量密度都比較低，大多數(shù)太陽(yáng)能飛行器都采用大展弦比、低翼載的氣動(dòng)布局并且動(dòng)力裝置多使用分布式拉力螺旋槳，這將使得機(jī)翼大部分面積都處在螺旋槳滑流中，如美國(guó)的太陽(yáng)神無(wú)人機(jī)機(jī)翼，滑流掃掠過的機(jī)翼面積占機(jī)翼總面積的50%[1]。同時(shí)太陽(yáng)能飛機(jī)多在高空低速飛行，低雷諾數(shù)特征[2]明顯，滑流與低雷諾數(shù)流動(dòng)的共同作用將明顯改變太陽(yáng)能飛機(jī)的氣動(dòng)性能。因此在太陽(yáng)能飛機(jī)的初步設(shè)計(jì)階段，應(yīng)充分考慮分布式螺旋槳的個(gè)數(shù)、轉(zhuǎn)速、安裝位置等參數(shù)對(duì)全機(jī)氣動(dòng)力的影響從而實(shí)現(xiàn)最優(yōu)設(shè)計(jì)。然而在該階段如何快速、準(zhǔn)確計(jì)算分析螺旋槳與機(jī)翼各自的氣動(dòng)力以及它們的相互干擾成為亟待解決的問題。

關(guān)于螺旋槳的計(jì)算，現(xiàn)已發(fā)展成熟的有早期的動(dòng)量理論、葉素理論、渦流理論、升力線理論[3]。雖然這些理論方法計(jì)算效率較高，但精度比較低，并且它們只是針對(duì)單獨(dú)螺旋槳的計(jì)算，當(dāng)需要考慮螺旋槳-螺旋槳、螺旋槳-機(jī)翼之間復(fù)雜氣動(dòng)干擾時(shí)，還要借助于其他方法建立這些理論與機(jī)翼氣動(dòng)力的聯(lián)系。

隨著計(jì)算流體力學(xué)(CFD)技術(shù)的發(fā)展，對(duì)于螺旋槳、旋翼等這類運(yùn)動(dòng)部件的數(shù)值模擬也出現(xiàn)了一些較為成熟可靠的計(jì)算方法。例如忽略螺旋槳幾何外形將其簡(jiǎn)化為無(wú)限薄圓盤的激勵(lì)盤方法[4]，利用相對(duì)運(yùn)動(dòng)原理將非定常問題簡(jiǎn)化為定常問題的多重參考系[5](Multiple Reference Frames, MRF) 法，完全真實(shí)模擬螺旋槳旋轉(zhuǎn)運(yùn)動(dòng)的非定?；凭W(wǎng)格 (Sliding Mesh) 法[6-7]。雖然CFD方法具有較高的計(jì)算精度，但其計(jì)算效率卻相對(duì)較低，尤其是采用滑移網(wǎng)格法計(jì)算將更加耗時(shí)[8]，因此CFD方法也無(wú)法滿足太陽(yáng)能飛機(jī)初步設(shè)計(jì)階段準(zhǔn)確并且快速的計(jì)算需求。

面元法自得到發(fā)展以來，它以其簡(jiǎn)單易用、精度較好、計(jì)算效率高的特點(diǎn)在飛行器初步設(shè)計(jì)階段得到了廣泛應(yīng)用[9]。但是傳統(tǒng)的面元法由于基于無(wú)黏不可壓的假設(shè)因此未能考慮氣體黏性影響，對(duì)于低雷諾數(shù)流動(dòng)，這將會(huì)引起較大的計(jì)算誤差。其次，在面元法計(jì)算過程中，升力面的尾跡形狀需要人為事先指定，而尾跡形狀對(duì)計(jì)算結(jié)果又會(huì)產(chǎn)生很大影響[10]，尤其是模擬尾跡的卷起時(shí)，尾跡面元相互靠近容易產(chǎn)生數(shù)值奇異問題從而導(dǎo)致計(jì)算結(jié)果誤差較大。

在拉格朗日體系下，求解渦量形式Navier-Stokes方程的黏性渦粒子法能精確模擬黏性流動(dòng)且具有無(wú)網(wǎng)格優(yōu)點(diǎn)，雖然它的計(jì)算效率整體上比基于歐拉坐標(biāo)系的CFD方法高，但是對(duì)于多個(gè)物體誘導(dǎo)速度場(chǎng)的計(jì)算其計(jì)算量仍然較大。

為避免上述問題并兼顧計(jì)算效率和計(jì)算精度，本文采用非定常面元/黏性渦粒子的混合方法[11]對(duì)螺旋槳-機(jī)翼的氣動(dòng)干擾進(jìn)行研究分析。通過非定常面元法求解機(jī)翼和螺旋槳的氣動(dòng)載荷，采用黏性渦粒子法非定常模擬升力面尾跡的發(fā)展，并根據(jù)誘導(dǎo)速度相等的原則建立尾跡面元與渦粒子的聯(lián)系從而實(shí)現(xiàn)快速、準(zhǔn)確的氣動(dòng)干擾計(jì)算?；谠摲椒?，首先選擇有實(shí)驗(yàn)數(shù)據(jù)的低雷諾數(shù)翼型和螺旋槳進(jìn)行方法驗(yàn)證，然后計(jì)算分析不同螺旋槳的安裝位置、螺旋槳個(gè)數(shù)、旋轉(zhuǎn)方向?qū)λ陨硇阅芗皺C(jī)翼氣動(dòng)力的影響規(guī)律。

1 計(jì)算方法

1.1 非定常面元法

1.1.1 速度勢(shì)方程

對(duì)于已知邊界的升力體，當(dāng)它在無(wú)黏、無(wú)旋、不可壓的流場(chǎng)中運(yùn)動(dòng)時(shí)，其連續(xù)方程的表達(dá)式為[10]

(1)

(2)

φ∞=U∞x+V∞y+W∞z

(3)

式中：U∞、V∞和W∞為來流速度分量。

1.1.2 邊界條件

在升力面的邊界上，由于物面不可穿透，因此速度沿物面的法向分量應(yīng)為零。這一條件稱為黎曼邊界條件(Neumann Boundary Condition)。根據(jù)該條件可得到源(匯)的強(qiáng)度為

σ=n·Q∞

(4)

式中：Q∞為自由來流速度。

當(dāng)在物面邊界上分布一系列源(匯)和偶極子時(shí)，連續(xù)的物面可以離散化為NB個(gè)物體面元和Nw個(gè)尾跡面元，這樣式(2)的積分形式就可以寫成多項(xiàng)式表達(dá)式從而實(shí)現(xiàn)數(shù)值求解，即

(5)

式中：

(6)

(7)

其中：Ck和Bk分別為第k個(gè)四邊形面元上(角點(diǎn)編號(hào)為1,2,3,4)單位強(qiáng)度的偶極子和源(匯)對(duì)空間任意一點(diǎn)P的影響系數(shù)；Cw的計(jì)算式與Ck相同。

為使式(5)的方程組可以求解，未知的尾跡面元的偶極子強(qiáng)度μw由未知的物面偶極子強(qiáng)度來代替，利用庫(kù)塔條件可以建立它們之間的關(guān)系：

μw=μu-μd

(8)

式中：μu和μd分別為物面尾緣上下表面的偶極子強(qiáng)度。

結(jié)合式(6)～式(8)，可將式(5)的多項(xiàng)式寫成矩陣形式從而直接求得偶極子強(qiáng)度。

(9)

1.1.3 非定常氣動(dòng)力

當(dāng)求解得到所有偶極子強(qiáng)度后便可獲得各個(gè)面元控制點(diǎn)處的誘導(dǎo)速度，再根據(jù)伯努利方程最終可計(jì)算出各面元上的壓力系數(shù)為

(10)

式中：p和pref分別為當(dāng)?shù)仂o壓和遠(yuǎn)場(chǎng)參考?jí)毫?；ρ為氣體密度；vref為機(jī)體坐標(biāo)系內(nèi)面元的運(yùn)動(dòng)速度；Q為面元運(yùn)動(dòng)速度與面元上誘導(dǎo)速度的合速度。

1.1.4 壓縮性影響

使用面元法計(jì)算旋轉(zhuǎn)部件的氣動(dòng)力時(shí)，文獻(xiàn)[10]未考慮氣體壓縮性的影響，而本文對(duì)于螺旋槳?dú)鈩?dòng)力的計(jì)算，其靠近槳尖剖面的當(dāng)?shù)伛R赫數(shù)大于0.3，氣流的壓縮性影響已經(jīng)不能忽略，為降低計(jì)算誤差，本文引入普朗特-格勞爾特法則進(jìn)行可壓縮修正[12]：

(11)

式中：Cp0為不可壓流的壓力系數(shù)；Ma∞為自由來流馬赫數(shù)。

1.2 黏性渦粒子法

黏性渦粒子對(duì)渦量場(chǎng)的描述是基于拉格朗日坐標(biāo)系。對(duì)于不可壓的有黏流動(dòng)，可將Navier-Stokes方程表達(dá)成速度-渦量形式[13]：

(12)

式中:u為自由來流速度;ω為渦量;ν為運(yùn)動(dòng)黏性系數(shù)。

使用渦粒子方法進(jìn)行計(jì)算時(shí)首先需要對(duì)連續(xù)的渦量場(chǎng)進(jìn)行離散化，本文使用可攜帶渦量、體積的渦團(tuán)[14]對(duì)渦量場(chǎng)離散化，因此空間某一點(diǎn)在時(shí)刻t的渦量可表達(dá)為

(13)

式中:N為渦元的數(shù)量；xq和aq為第q個(gè)渦團(tuán)的位置和渦矢量；ζε為高斯分布函數(shù)：

(14)

式中：Rv=|x-xq|/ε為任意一點(diǎn)到渦元的無(wú)量綱距離，ε為光滑參數(shù)，其作用是在使用畢奧-薩伐爾定律計(jì)算誘導(dǎo)速度時(shí)避免產(chǎn)生數(shù)值奇異性。

當(dāng)連續(xù)的渦量場(chǎng)被離散化后，控制方程式(12)可分化為

(15)

(16)

式中：u∞為來流速度；誘導(dǎo)速度uind的計(jì)算可通過畢奧-薩伐爾定律求得

(17)

1.2.1 拉伸項(xiàng)求解

本文采用直接法求解拉伸項(xiàng)，即直接構(gòu)建速度梯度矩陣，并用該梯度矩陣與渦元的渦量相乘，其表達(dá)式為

(18)

速度梯度的形式為

(19)

1.2.2 擴(kuò)散項(xiàng)求解

求解擴(kuò)散項(xiàng)的方法主要有粒子強(qiáng)度交換 (Particle Strength Exchange, PSE) 法[16]、隨機(jī)走步法[14]等。最常用的為粒子強(qiáng)度交換法，它的求解思想是將拉普拉斯算子替換為一個(gè)積分函數(shù)：

(20)

式中：Vp、Vq分別為第p、q個(gè)渦團(tuán)的體積。

1.2.3 尾跡面元與渦粒子等效原則

當(dāng)非定常面元的偶極子強(qiáng)度已知時(shí)，需使用黏性渦粒子方法非定常模擬升力面尾跡的發(fā)展，因此需要建立兩者的聯(lián)系。在此過程中，本文在升力面的尾緣建立了兩行尾跡面元，第一行尾跡面元用來滿足庫(kù)塔條件，第二行尾跡面元為前一個(gè)時(shí)間步的第一組面元，其偶極子強(qiáng)度已知。對(duì)于該行面元，使用Hess等效原則[17]建立兩種計(jì)算方法之間的聯(lián)系，即令渦粒子和尾跡面元對(duì)空間任意一點(diǎn)的誘導(dǎo)速度相同：

(21)

式中：各個(gè)變量的定義見文獻(xiàn)[17]。將尾跡面元等效成渦粒子微團(tuán)后升力面尾跡的發(fā)展則利用式(15)和式(16)通過求解計(jì)及黏性的速度-渦量形式的Navier-Stokes方程來確定。

1.2.4 渦粒子尾跡對(duì)面元的影響

當(dāng)考慮渦粒子尾跡對(duì)升力面氣動(dòng)力的影響時(shí)，最直觀的解決方法為計(jì)算所有渦粒子對(duì)升力面面元的誘導(dǎo)速度[18]。此時(shí)，式(4)用來計(jì)算源(匯)強(qiáng)度的表達(dá)式應(yīng)變化為

(22)

綜合式(21)和式(22)可以看出，尾跡面元的偶極子強(qiáng)度決定了渦粒子的渦量大小，而渦粒子強(qiáng)度對(duì)升力面面元的誘導(dǎo)速度又直接影響了升力面面元上偶極子的強(qiáng)度，因此基于非定常面元/黏性渦粒子的混合計(jì)算方法是一種強(qiáng)耦合算法。

2 方法驗(yàn)證

2.1 低雷諾數(shù)翼型

由于太陽(yáng)能飛機(jī)的飛行狀態(tài)多為高空低雷諾數(shù)流動(dòng)，氣體的黏性影響已經(jīng)不能忽略。為校驗(yàn)本文采用的非定常面元/黏性渦粒子混合方法對(duì)低雷諾數(shù)流動(dòng)的計(jì)算精度，選取有實(shí)驗(yàn)數(shù)據(jù)的低雷諾數(shù)翼型E387[19]進(jìn)行對(duì)比驗(yàn)證。實(shí)驗(yàn)雷諾數(shù)Re=4.6×105，馬赫數(shù)Ma=0.13。表1給出了實(shí)驗(yàn)、CFD方法以及本文的混合方法三者計(jì)算結(jié)果的對(duì)比，表中：CL為升力系數(shù)，CD為阻力系數(shù)，Cm為俯仰力矩系數(shù)，α為迎角。其中CFD的結(jié)果是采用成熟的商業(yè)軟件Fluent計(jì)算所得，湍流模型為k-kL-ω三方程轉(zhuǎn)捩模型。需要說明的是，文中所有的CFD方法驗(yàn)證都采用該湍流模型。

表1 E387翼型計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比Table 1 Comparison of calculated results with experimental data for E387 airfoil

通過對(duì)比發(fā)現(xiàn)，計(jì)算得到的升力和力矩系數(shù)與實(shí)驗(yàn)結(jié)果、CFD數(shù)值模擬結(jié)果都十分接近，其誤差均在5%以內(nèi)。阻力系數(shù)計(jì)算結(jié)果的誤差相對(duì)較大，然而對(duì)于太陽(yáng)能飛機(jī)的初步設(shè)計(jì)分析，其誤差也在可以接受的范圍內(nèi)。所以本文采用的非定常面元/黏性渦粒子的混合計(jì)算方法整體上對(duì)低雷諾數(shù)流動(dòng)的計(jì)算具有較高的精度。

2.2 螺旋槳模型驗(yàn)證

對(duì)于螺旋槳的計(jì)算，本文采用的計(jì)算方法計(jì)及了壓縮性的影響。本算例選用有詳細(xì)試驗(yàn)數(shù)據(jù)的風(fēng)洞模型[20]進(jìn)行對(duì)比驗(yàn)證。該試驗(yàn)?zāi)Ｐ陀蓛蓚€(gè)矩形槳葉組成，槳葉剖面翼型為NACA0012，直徑為2.286 m，展弦比為6，槳葉剖面沿徑向方向無(wú)幾何扭轉(zhuǎn)和氣動(dòng)扭轉(zhuǎn)，槳葉總距為8°，實(shí)驗(yàn)轉(zhuǎn)速為1 250 r/min，大氣計(jì)算條件為標(biāo)準(zhǔn)海平面大氣，槳葉葉尖馬赫數(shù)為Ma=0.439，該實(shí)驗(yàn)中，槳葉拉力系數(shù)的定義[21]為

(23)

式中：Ω為轉(zhuǎn)速；R為槳葉半徑；T為拉力。表2同樣給出了計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比。

表2 拉力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比

由表2可知，CFD方法的結(jié)果與實(shí)驗(yàn)幾乎相同，而本文計(jì)及壓縮性影響的混合方法的計(jì)算結(jié)果與實(shí)驗(yàn)值也十分接近，誤差只有3.1%，同時(shí)，圖1給出了不同徑向位置r/R的槳葉剖面的壓力系數(shù)Cp的分布對(duì)比。圖中：“Reference”曲線為文獻(xiàn)[10]使用面元法而未考慮壓縮性影響的計(jì)算結(jié)果。顯然，本文計(jì)及壓縮性影響后其計(jì)算精度有了比較明顯的提升。

圖1 槳葉不同徑向位置的弦向壓力分布對(duì)比
Fig.1 Chordwise pressure distribution at different radial sections on blades

3 螺旋槳/機(jī)翼氣動(dòng)干擾

本節(jié)研究螺旋槳對(duì)機(jī)翼氣動(dòng)力的干擾，首先從單獨(dú)螺旋槳安裝位置的影響開始著手。單獨(dú)螺旋槳與機(jī)翼的氣動(dòng)布局構(gòu)型如圖2所示。

機(jī)翼為展弦比AR=2的矩形翼，無(wú)扭轉(zhuǎn)，無(wú)安裝角；飛行高度H=20 km，雷諾數(shù)Re=4.88×105，馬赫數(shù)Ma=0.101 7。x軸正方向?yàn)樽杂蓙砹鞣较?，螺旋槳設(shè)計(jì)轉(zhuǎn)速為1 500 r/min，并沿x軸正向逆時(shí)針旋轉(zhuǎn)。

圖2 單獨(dú)螺旋槳與機(jī)翼布局
Fig.2 Configuration of single propeller and wing

為進(jìn)一步證明非定常面元/黏性渦粒子混合方法的計(jì)算效率和計(jì)算精度，表3給出了該方法與CFD方法對(duì)機(jī)翼和螺旋槳?dú)鈩?dòng)力計(jì)算結(jié)果的對(duì)比。由此可見，所有的計(jì)算值與CFD結(jié)果都十分接近從而再次校驗(yàn)了混合方法的精度。然而在計(jì)算效率上，本文的方法則具有比較明顯的優(yōu)勢(shì)。對(duì)于太陽(yáng)能飛機(jī)的初步設(shè)計(jì)，其重心在于快速獲得螺旋槳對(duì)機(jī)翼宏觀氣動(dòng)力的影響而對(duì)全機(jī)流場(chǎng)的細(xì)節(jié)關(guān)注較少，因此使用本文的方法則能夠在保證計(jì)算精度的同時(shí)顯著提高計(jì)算效率。

表3 混合法與CFD計(jì)算結(jié)果及計(jì)算效率對(duì)比

3.1 單獨(dú)螺旋槳位置影響

考慮單獨(dú)螺旋槳位置對(duì)機(jī)翼的影響主要分為弦向、垂直安裝位置兩個(gè)影響因素。圖3為2° 迎角下不同螺旋槳位置對(duì)機(jī)翼升阻力的影響結(jié)果，其計(jì)算結(jié)果與文獻(xiàn)[22]中的變化規(guī)律相一致。同時(shí)，在圖3中，針對(duì)一組螺旋槳位置，本文也給出了CFD數(shù)值模擬結(jié)果以便對(duì)比。圖中：Δx/R為螺旋槳旋轉(zhuǎn)中心與機(jī)翼前緣的無(wú)量綱水平距離，Δz/R為螺旋槳中心相對(duì)于機(jī)翼弦平面的無(wú)量綱垂直距離，其值為正時(shí)表明螺旋槳高于弦平面。

圖3 螺旋槳安裝位置對(duì)機(jī)翼氣動(dòng)力的影響
Fig.3 Effect of propeller positions on aerodynamic force of wing

對(duì)于拉力螺旋槳布局，它對(duì)機(jī)翼氣動(dòng)力影響的第一個(gè)作用是螺旋槳通過旋轉(zhuǎn)為槳盤后方的來流注入能量，提高了軸向氣流的總壓和動(dòng)壓。總壓的增加導(dǎo)致機(jī)翼的阻力隨之增大，而動(dòng)壓的升高直接體現(xiàn)為軸向氣流的加速進(jìn)而引起升力的增加。隨著弦向間距的增大，滑流流管逐漸收縮，當(dāng)滑流遇到后方的機(jī)翼時(shí)，其流管收縮與加速的趨勢(shì)進(jìn)一步加強(qiáng)，滑流“增升”的效果也愈加明顯。因此螺旋槳滑流的存在導(dǎo)致機(jī)翼氣動(dòng)力出現(xiàn)了“增升增阻”的現(xiàn)象，例如當(dāng)螺旋槳位置為Δx/R=0.5,Δz/R=0.5時(shí)， 升力增量為2.0%，阻力增量為13.3%，升阻比下降10%。

除了對(duì)自由來流的軸向加速作用外，螺旋槳滑流還將誘導(dǎo)均勻的自由來流產(chǎn)生周向速度。對(duì)于螺旋槳的上行槳葉和下行槳葉，其周向速度在螺旋槳轉(zhuǎn)軸兩側(cè)等值反向。因此對(duì)于處在滑流區(qū)內(nèi)的翼段，周向速度與自由來流速度的矢量疊加將會(huì)增大部分機(jī)翼前緣的有效迎角同時(shí)減小另一部分機(jī)翼的有效迎角，此即為螺旋槳滑流引起的上洗和下洗效應(yīng)。從圖4給出的機(jī)翼展向升力分布可以清楚看到因上洗和下洗效應(yīng)改變了滑流區(qū)內(nèi)機(jī)翼的有效迎角，機(jī)翼的展向升力分布在螺旋槳轉(zhuǎn)軸位置(圖中豎直線所示)兩側(cè)同時(shí)出現(xiàn)了升力“波峰”和“波谷”的現(xiàn)象。螺旋槳的旋轉(zhuǎn)方向決定了“波峰”和“波谷”出現(xiàn)的位置，而轉(zhuǎn)速以及安裝位置則最終決定了機(jī)翼展向升力分布的形狀。

圖5給出了在Δx/R=1.0,Δz/R=-0.5安裝位置下螺旋槳對(duì)機(jī)翼上下表面壓力分布的影響，圖中的豎直線代表了螺旋槳轉(zhuǎn)軸。壓力等值線在機(jī)翼前緣的分布較為集中表明滑流對(duì)機(jī)翼氣動(dòng)力的影響主要集中在這一區(qū)域。由于螺旋槳的垂直安裝位置低于機(jī)翼弦平面，并且自由來流迎角為正，因此機(jī)翼下表面壓力受到的滑流影響比機(jī)翼上表面較為明顯。在下翼面前緣以及螺旋槳轉(zhuǎn)軸右側(cè)，上洗效應(yīng)導(dǎo)致的壓力分布稍大于左側(cè)下洗區(qū)內(nèi)的壓力分布。

圖4 機(jī)翼展向升力分布曲線
Fig.4 Spanwise lift distribution curves of wing

圖5 機(jī)翼上下表面壓力云圖
Fig.5 Pressure counter of upper and lower wing surface

3.2 雙槳對(duì)機(jī)翼的影響

為保證影響因素的單一性，本節(jié)將在螺旋槳同一弦向位置、垂直位置以及固定轉(zhuǎn)速這一前提下研究分析雙槳展向間距、旋轉(zhuǎn)方向?qū)C(jī)翼氣動(dòng)力的影響。兩個(gè)螺旋槳的弦向安裝位置均為Δx/R=1.0，垂直位置Δz/R=-0.5，螺旋槳轉(zhuǎn)速仍為1 500 r/min，機(jī)翼來流迎角仍為2°。

圖6為雙槳布局的后視圖，圖中的箭頭指向規(guī)定了兩個(gè)螺旋槳的旋轉(zhuǎn)正方向。當(dāng)兩個(gè)螺旋槳都以正方向旋轉(zhuǎn)時(shí)則標(biāo)記為(P1+, P2+)。由于在無(wú)側(cè)滑條件下，螺旋槳(P1+, P2+)和(P1-, P2-)兩個(gè)轉(zhuǎn)向在理論上對(duì)機(jī)翼氣動(dòng)力的干擾結(jié)果是相同的，因此兩個(gè)螺旋槳4種旋轉(zhuǎn)方向的組合此時(shí)可以縮減為3種情況，即(P1+, P2+)，(P1-, P2+)，(P1+, P2-)。圖中：Δy表示兩個(gè)螺旋槳旋轉(zhuǎn)中心的間距，并以Δy/R對(duì)中心間距進(jìn)行無(wú)量綱化處理。計(jì)算中Δy/R最大值為6，此時(shí)兩個(gè)螺旋槳的轉(zhuǎn)軸正好位于左右翼尖位置。

圖7為兩個(gè)螺旋槳的旋轉(zhuǎn)方向以及間距變化對(duì)機(jī)翼氣動(dòng)力的影響結(jié)果。顯然，雙槳的不同旋轉(zhuǎn)方向使機(jī)翼氣動(dòng)力表現(xiàn)出不同的變化特征。

當(dāng)雙槳轉(zhuǎn)向?yàn)?P1-,P2+)時(shí)，滑流的增升效果最明顯，并隨著間距的增大呈現(xiàn)進(jìn)一步增加的趨勢(shì)。這是由于在該轉(zhuǎn)向下，滑流誘導(dǎo)的上洗區(qū)域位于兩個(gè)螺旋槳之間，當(dāng)槳間距增大時(shí)，上洗區(qū)域隨之增大，這就意味著機(jī)翼有效迎角增加的范圍越大。當(dāng)兩個(gè)螺旋槳的展向位置逐漸靠近翼尖尤其是當(dāng)Δy/R=6，即兩個(gè)螺旋槳轉(zhuǎn)軸位于左右翼尖時(shí)，(P1-,P2+)的旋轉(zhuǎn)方向正好與兩個(gè)翼尖渦方向相反從而削弱了翼尖渦的強(qiáng)度，增大了機(jī)翼的有效展弦比。所以在此位置下，機(jī)翼升力明顯提升，阻力進(jìn)一步下降，并且機(jī)翼的阻力小于干凈機(jī)翼的情況。

圖6 雙槳布局后視圖
Fig.6 Configuration of double propellers (back view)

圖7 雙槳?dú)鈩?dòng)力計(jì)算結(jié)果
Fig.7 Calculated results of aerodynamic force of double propellers

當(dāng)雙槳以(P1+,P2-)的方向旋轉(zhuǎn)時(shí)，其情況則與(P1-,P2+)的結(jié)果正好相反。即在兩個(gè)螺旋槳之間的翼段內(nèi)，滑流減小了機(jī)翼的有效迎角，隨著中心間距的增大，機(jī)翼迎角減小的范圍進(jìn)一步擴(kuò)大導(dǎo)致升力逐漸降低；與此同時(shí)，在Δy/R=6位置上，兩個(gè)螺旋槳的旋轉(zhuǎn)方向與左右翼尖渦同向，這一作用降低了機(jī)翼有效展弦比從而使得機(jī)翼升力明顯減小，誘導(dǎo)阻力迅速增大。所以在該轉(zhuǎn)向下，機(jī)翼氣動(dòng)效率的損失也最為嚴(yán)重，升阻比相對(duì)于干凈機(jī)翼下降了10.37%。

圖8給出了雙槳間距Δy/R=4，螺旋槳轉(zhuǎn)向分別為(P1+, P2-)、(P1+, P2+)時(shí)機(jī)翼展向升力分布的對(duì)比，其中，y/b=-1代表機(jī)翼左翼尖位置，y/b=1代表機(jī)翼右翼尖位置。從圖中可以看到雙槳的存在(圖中豎直虛線所示)將機(jī)翼表面的流動(dòng)分成了3個(gè)區(qū)域(Ⅰ、Ⅱ、Ⅲ)。在(P1+, P2-)轉(zhuǎn)向下，兩個(gè)槳葉在區(qū)域Ⅱ內(nèi)共同誘導(dǎo)出了下洗效應(yīng)導(dǎo)致該區(qū)域內(nèi)的升力分布整體上小于干凈機(jī)翼，而Ⅰ、Ⅲ區(qū)域則為上洗區(qū)，其升力略高于干凈構(gòu)型。由于2個(gè)螺旋槳的旋轉(zhuǎn)方向關(guān)于機(jī)翼中心面對(duì)稱，所以機(jī)翼的展向升力也呈現(xiàn)對(duì)稱分布。

圖8 雙槳展向升力分布對(duì)比
Fig.8 Spanwise lift distribution affected by double propellers

對(duì)于(P1+, P2+)的旋轉(zhuǎn)方向，其結(jié)果將有所變化。P2旋轉(zhuǎn)方向的改變使上洗區(qū)和下洗區(qū)的位置有所互換。由于P1的轉(zhuǎn)向未改變，所以在P1轉(zhuǎn)軸兩側(cè)，兩種轉(zhuǎn)向下的升力分布幾乎相同。

由此可見，相比于單獨(dú)螺旋槳對(duì)機(jī)翼的氣動(dòng)干擾，分布式螺旋槳滑流對(duì)機(jī)翼的展向升力影響更為復(fù)雜也更為顯著。

4 螺旋槳性能

在低速流動(dòng)條件下，由于螺旋槳滑流與機(jī)翼氣動(dòng)力的干擾是相互的，因此滑流對(duì)機(jī)翼產(chǎn)生影響的同時(shí)，其自身的性能也會(huì)有一定的變化。

圖9為單槳-機(jī)翼氣動(dòng)布局下螺旋槳弦向、垂直位置對(duì)自身拉力的影響，此外本文也給出了一組CFD數(shù)值模擬結(jié)果，兩者的最大計(jì)算誤差只有3.5%，說明本文的方法在考慮機(jī)翼氣動(dòng)干擾下對(duì)螺旋槳拉力的計(jì)算也有良好的精度。

圖9 螺旋槳拉力系數(shù)隨弦向位置變化關(guān)系
Fig.9 Effect of single propeller streamwise locations on propeller thrust coefficient

計(jì)算結(jié)果表明隨著弦向間距的增大，其拉力逐漸減小并在Δx/R>2之后逐漸趨于平緩，這是由于機(jī)翼的存在對(duì)槳盤后方的滑流產(chǎn)生了一定的阻滯作用，提高了槳盤后方的壓力。當(dāng)螺旋槳距離機(jī)翼前緣較遠(yuǎn)時(shí)，這種阻滯作用有所減弱，螺旋槳的拉力逐漸接近于無(wú)機(jī)翼干擾下的單獨(dú)螺旋槳的情況。另外，在同一弦向位置下，螺旋槳較低的垂直安裝位置相對(duì)而言有利于提高自身的拉力，這與螺旋槳對(duì)機(jī)翼的干擾情況正好相反。

圖10為(P1-, P2+)轉(zhuǎn)向下雙槳拉力、效率η與其展向間距的變化關(guān)系。同時(shí)文中也給出了CFD計(jì)算結(jié)果。顯然，兩者的變化趨勢(shì)相一致，拉力系數(shù)的計(jì)算誤差約為4%，螺旋槳效率的誤差約為6%。從圖中可知，雙槳之間的氣動(dòng)干擾對(duì)各自的性能影響不大，兩個(gè)螺旋槳的拉力和效率幾乎重合。而展向間距的變化則使兩個(gè)螺旋槳的性能逐漸下降，在Δy/R=6.0位置處螺旋槳的推進(jìn)效率下降了2%。由此可見，機(jī)翼對(duì)螺旋槳性能的影響也是不能忽視的。

圖10 雙槳性能隨螺旋槳間距變化關(guān)系
Fig.10 Effect of double propellers spanwise locations on their performance

5 結(jié) 論

1) 計(jì)及氣體黏性、壓縮性影響的非定常面元/黏性渦粒子的混合方法能夠比較準(zhǔn)確地計(jì)算分析低雷諾數(shù)流動(dòng)下螺旋槳與機(jī)翼之間氣動(dòng)力的相互干擾，同時(shí)該方法的計(jì)算效率相比于CFD數(shù)值模擬方法具有比較明顯的優(yōu)勢(shì)，可以用于太陽(yáng)能飛機(jī)初步設(shè)計(jì)階段快速的氣動(dòng)力分析。

2) 螺旋槳滑流既對(duì)自由來流進(jìn)行軸向加速又誘導(dǎo)其產(chǎn)生上洗和下洗效應(yīng)從而同時(shí)改變了機(jī)翼展向、弦向的升力分布，這些作用最終使機(jī)翼氣動(dòng)力出現(xiàn)了明顯的增升增阻現(xiàn)象，機(jī)翼升阻比有所下降。

3) 螺旋槳弦向和垂直安裝位置對(duì)機(jī)翼氣動(dòng)力的影響都比較明顯。較大的弦向間距和較高的垂直位置對(duì)機(jī)翼的升力、阻力相對(duì)有利。對(duì)于多槳干擾的情況，不同的螺旋槳旋轉(zhuǎn)方向使機(jī)翼氣動(dòng)力呈現(xiàn)出不同的變化特征。在合適的轉(zhuǎn)向下，螺旋槳滑流能夠削弱翼尖渦的強(qiáng)度從而提高機(jī)翼的有效展弦比，增大升力，降低誘導(dǎo)阻力。

4) 螺旋槳之間的相互干擾對(duì)自身性能的影響較弱；而展向間距的增大則使螺旋槳本身拉力和推進(jìn)效率都有所降低。

綜上，螺旋槳與機(jī)翼之間氣動(dòng)力的相互干擾對(duì)各自的性能都產(chǎn)生了不可忽略的影響，所以在太陽(yáng)能飛機(jī)初步設(shè)計(jì)階段，應(yīng)充分考慮分布式螺旋槳相對(duì)于全機(jī)合理的安裝位置，同時(shí)應(yīng)充分利用分布式螺旋槳滑流帶來的有利影響，通過調(diào)整螺旋槳間距、轉(zhuǎn)速、旋轉(zhuǎn)方向等措施進(jìn)行滑流影響下的增升減阻的氣動(dòng)優(yōu)化設(shè)計(jì)。

[1] NOLL T E, ISHMAEL S D, HENWOOD B, et al. Technical findings, lessons learned, and recommendations resulting from the helios prototype vehicle mishap: NASA-20070022260[R]. Washington， D.C.: National Aeronautics and Space Admin Langley Research Center, 2007.

[2] WINDTE J, SCHOLZ U, RADESPIEL R. Validation of the RANS-simulations of laminar separation bubbles on airfoils[J]. Aerospace Science and Technology, 2006, 10(6): 484-494.

[3] 劉沛清． 空氣螺旋槳理論及其應(yīng)用[M]． 北京: 北京航空航天大學(xué)出版社, 2006: 55-82．

LIU P Q. Air propeller theory and application[M]. Beijing: Beihang University Press， 2006: 55-82 (in Chinese).

[4] 段中喆, 劉沛清, 屈秋林． 修正的動(dòng)力盤模型與三維模擬螺旋槳滑流比較[J]．北京航空航天大學(xué)學(xué)報(bào), 2013, 39(5): 585-589．

DUAN Z Z, LIU P Q, QU Q L. Slipstream characters comparison of improved actuator disk model and 3D propeller numerical simulation[J]. Journal of Beijing University of Aeronautics and Astronautics， 2013， 39(5): 585-589 (in Chinese).

[5] 徐家寬, 白俊強(qiáng), 黃江濤, 等． 考慮螺旋槳滑流影響的機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)[J]． 航空學(xué)報(bào), 2014, 35(11): 2910-2920．

XU J K，BAI J Q，HUANG J T, et al. Aerodynamic optimization design of wing under the interaction of propeller slipstream[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(11): 2910-2920 (in Chinese).

[6] NICOLAS T, CHRISTIAN B, NIKOLAUS A A. Numerical and experimental analysis of a generic fan-in-wing configuration[J]. Journal of Aircraft, 2009, 46(2): 656-666.

[7] TSUNG L L, PAN K C. Application of the sliding mesh technique for helicopter rotor flow simulation[J]. Journal of Aeronautics, Astronautics and Aviation, 2012, 44(3): 201-210.

[8] 王紅波, 祝小平, 周洲, 等． 太陽(yáng)能無(wú)人機(jī)螺旋槳滑流氣動(dòng)特性分析[J]． 西北工業(yè)大學(xué)學(xué)報(bào), 2015, 33(6): 913-920．

WANG H B， ZHU X P， ZHOU Z, et al. Aerodyanmic investigation on propeller slipstream flows for solar powered airplanes[J]. Journal of Northwestern Polytechnical University, 2015, 33(6): 913-920 (in Chinese).

[9] 徐華舫． 亞、超音速常位流的面元法[M]． 北京: 國(guó)防工業(yè)出版社, 1981: 159-190．

XU H F. Panel methods of subsonic and supersonic flows[M]. Beijing: National Defense Industry Press, 1981: 159-190 (in Chinese).

[10] KATZ J, PLOTKIN A. Low speed aerodynamics[M]. Cambridge: Cambridge University Press, 2000: 207-217.

[11] WILLIS D J, JAIME P, JACOB K W. A combined pFFT-multipole tree code, unsteady panel method with vortex particle wakes[J]. International Journal for Numerical Methods in Fluids, 2007, 53(8): 1399-1422.

[12] 李鳳蔚．空氣與氣體動(dòng)力學(xué)引論[M]． 西安: 西北工業(yè)大學(xué)出版社, 2007: 253-256．

LI F W. Air and gas dynamic theory[M]. Xi’an: Northwestern Polytechnical University Press, 2007:253-256 (in Chinese).

[13] 魏鵬． 旋翼非定常流場(chǎng)的黏性渦數(shù)值模擬方法及其混合方法的研究[D]． 南京: 南京航空航天大學(xué), 2012: 11-28．

WEI P. Research on viscous vortex numerical algorithm and its hybrid method for rotor unsteady flowfield[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2012: 11-28 (in Chinese).

[14] CHORIN A J. Numerical study of slightly viscous flow[J]. Journal of Fluid Mechanics, 1973, 57(4): 785-796.

[15] HE C J, ZHAO J G. Modeling rotor wake dynamics with viscous vortex particle method[J]. AIAA Journal, 2009, 47(4): 902-915.

[16] ELDREDGE J D, LEONARD A, COLONIUS T. A general deterministic treatment of derivatives in particle methods[J]. Journal of Computational Physics, 2002, 180(2): 686-709.

[17] JOIN L H. Calculation of potential flow about arbitrary three dimensional lifting bodies: AD755840[R]. California: McDonnell Douglas Corporation, 1972.

[18] CALABRETTA J. A three dimensional vortex particle-panel code for modeling propeller-airframe interaction[D]. California: California Polytechnic State University, 2010: 134-139.

[19] MCGHEE R J, WALKER B S, MILLARD B F. Experimental results for the eppler 387 airfoil at low Reynolds numbers in the langley low-turbulence pressure tunnel: NASA TM 4062[R]. Washington, D.C.: NASA Langley Research Center, 1988.

[20] CARADONNA F X, TUNG C. Experimental and analytical studies of a model helicopter rotor in hover: NASA TM-81232[R]. Washington, D.C.: NASA, 1981.

[21] BRAND A G, MCMAHON H M, KOMERATH N M. Surface pressure measurements on a body subject to vortex wake interaction[J]. AIAA Journal, 1989, 27(5): 569-574.

[22] VELDHUIS L M. Propeller wing aerodynamic interference[D]. Delft: Delft University of Technology, 2005: 187-191.

AerodynamicinteractionsatlowReynoldsnumberslipstreamwithunsteadypanel/viscousvortexparticlemethod

WANGHongbo1,2,ZHUXiaoping2,ZHOUZhou1,2，*,XUXiaoping1,2

1.SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China2.ScienceandTechnologyonUAVLaboratory,NorthwesternPolytechnicalUniversity,Xi’an710065,China

Anunsteadypanel/viscousvortexparticlehybridmethod,withtheconsiderationofairviscousandcompressibilityeffectsatlowReynoldnumber,isdevelopedbaseonequivalentvorticityprincipleandNeumannboundaryconditiontorapidlycalculatetheaerodynamicinteractionbetweenthewingandthepropellerofthesolar-poweredairplane.Experimentaldataarecomparedwithcomputationmethodtovalidatethehybridmethodproposed.Theaerodynamicinteractionsbetweenthepropellerandthewingareinvestigatedatdifferentinstallationpositionsandworkingconditions.Calculatedresultsindicatethatthedistributionofthespanwiseandthechordwisepressureareapparentlychangedbytheincreasedaxialvelocityandupwashanddownwasheffectsinducedbythepropellerslipstreamtoleadtoadecreaseoflift-tot-dragratio.Alargerchordwisedistanceandahigherverticalinstallationpositioncanreducepropellerthrusts,andcanalsodeceleratelift-to-dragratioofthewing.Forthecaseofmultipropellerinteractions,differentrotationdirectionscausedifferentaerodynamiccharacteristicsofthewing.Whenthepropellerrotationdirectionisoppositetothewingtipvortexdirection,thepropellerslipstreamscancounteractvortexstrengthsatthewingtiptoinduceanaugmentoflift-to-dragratioofthewing.

lowReynoldsnumberslipstream;unsteadypanel;viscousvortexparticle;hybridmethod;aerodynamicinteractions

2016-05-09；Revised2016-07-25；Accepted2016-09-06；Publishedonline2016-09-091135

URL：www.cnki.net/kcms/detail/11.1929.V.20160909.1135.002.html

s：NationalHigh-techResearchandDevelopmentProgramofChina(2014AA7052002);CivilAircraftSpecificProject(MIZ-2015-F-009);ShaanxiProvinceScienceandTechnologyCo-ordinationProject(2015KTCQ01-78)

2016-05-09；退修日期2016-07-25；錄用日期2016-09-06； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-09-091135

www.cnki.net/kcms/detail/11.1929.V.20160909.1135.002.html

國(guó)家“863”計(jì)劃 (2014AA7052002); 民機(jī)專項(xiàng) (MIZ-2015-F-009); 陜西省科技統(tǒng)籌項(xiàng)目 (2015KTCQ01-78)

＊

.E-mailzhouzhou@nwpu.edu.cn

王紅波， 祝小平， 周洲， 等. 基于非定常面元/黏性渦粒子法的低雷諾數(shù)滑流氣動(dòng)干擾J. 航空學(xué)報(bào),2017,38(4):120412.WANGHB,ZHUXP,ZHOUZ,etal.AerodynamicinteractionsatlowReynoldsnumberslipstreamwithunsteadypanel/viscousvortexparticlemethodJ.ActaAeronauticaetAstronauticaSinica,2017,38(4):120412.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0251

V211

A

1000-6893(2017)04-120412-11

(責(zé)任編輯： 李明敏)

＊Correspondingauthor.E-mailzhouzhou@nwpu.edu.cn