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

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

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

2.西北工業(yè)大學(xué) 無人機特種技術(shù)重點實驗室， 西安 710065

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

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

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

2.西北工業(yè)大學(xué) 無人機特種技術(shù)重點實驗室， 西安 710065

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

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

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

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

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

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

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

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

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

1 計算方法

1.1 非定常面元法

1.1.1 速度勢方程

對于已知邊界的升力體，當它在無黏、無旋、不可壓的流場中運動時，其連續(xù)方程的表達式為[10]

(1)

(2)

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

(3)

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

1.1.2 邊界條件

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

σ=n·Q∞

(4)

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

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

(5)

式中：

(6)

(7)

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

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

μw=μu-μd

(8)

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

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

(9)

1.1.3 非定常氣動力

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

(10)

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

1.1.4 壓縮性影響

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

(11)

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

1.2 黏性渦粒子法

黏性渦粒子對渦量場的描述是基于拉格朗日坐標系。對于不可壓的有黏流動，可將Navier-Stokes方程表達成速度-渦量形式[13]：

(12)

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

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

(13)

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

(14)

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

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

(15)

(16)

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

(17)

1.2.1 拉伸項求解

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

(18)

速度梯度的形式為

(19)

1.2.2 擴散項求解

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

(20)

式中：Vp、Vq分別為第p、q個渦團的體積。

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

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

(21)

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

1.2.4 渦粒子尾跡對面元的影響

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

(22)

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

2 方法驗證

2.1 低雷諾數(shù)翼型

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

表1 E387翼型計算結(jié)果與實驗值對比Table 1 Comparison of calculated results with experimental data for E387 airfoil

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

2.2 螺旋槳模型驗證

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

(23)

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

表2 拉力系數(shù)計算結(jié)果與實驗值對比

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

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

3 螺旋槳/機翼氣動干擾

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

機翼為展弦比AR=2的矩形翼，無扭轉(zhuǎn)，無安裝角；飛行高度H=20 km，雷諾數(shù)Re=4.88×105，馬赫數(shù)Ma=0.101 7。x軸正方向為自由來流方向，螺旋槳設(shè)計轉(zhuǎn)速為1 500 r/min，并沿x軸正向逆時針旋轉(zhuǎn)。

圖2 單獨螺旋槳與機翼布局
Fig.2 Configuration of single propeller and wing

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

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

3.1 單獨螺旋槳位置影響

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

圖3 螺旋槳安裝位置對機翼氣動力的影響
Fig.3 Effect of propeller positions on aerodynamic force of wing

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

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

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

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

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

3.2 雙槳對機翼的影響

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

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

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

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

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

圖7 雙槳氣動力計算結(jié)果
Fig.7 Calculated results of aerodynamic force of double propellers

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

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

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

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

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

4 螺旋槳性能

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

[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] 段中喆, 劉沛清, 屈秋林． 修正的動力盤模型與三維模擬螺旋槳滑流比較[J]．北京航空航天大學(xué)學(xué)報, 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] 徐家寬, 白俊強, 黃江濤, 等． 考慮螺旋槳滑流影響的機翼氣動優(yōu)化設(shè)計[J]． 航空學(xué)報, 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] 王紅波, 祝小平, 周洲, 等． 太陽能無人機螺旋槳滑流氣動特性分析[J]． 西北工業(yè)大學(xué)學(xué)報, 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]． 北京: 國防工業(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] 李鳳蔚．空氣與氣體動力學(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] 魏鵬． 旋翼非定常流場的黏性渦數(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ò)出版時間

時間：2016-09-091135

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

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

＊

.E-mailzhouzhou@nwpu.edu.cn

王紅波， 祝小平， 周洲， 等. 基于非定常面元/黏性渦粒子法的低雷諾數(shù)滑流氣動干擾J. 航空學(xué)報,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

(責任編輯： 李明敏)

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