謝長川，張利娟，劉燚，楊超

(北京航空航天大學 航空科學與工程學院，北京100191)

超長航時的無人機，由于其續(xù)航能力的要求，加之燃料動力在低速中空或高高空時效率低的原因，這類飛機一般采用高效太陽能為電動螺旋槳提供動力以保障超長甚至永久續(xù)航能力.采用螺旋槳推進系統(tǒng)的高空超長航時無人機，一方面由于具有大展弦比、柔性大的特點，這類飛機在正常飛行情況下就會產(chǎn)生顯著的彈性變形，幾何非線性氣動彈性問題突出［1］;另一方面，由于柔性機翼翼載荷較低，可以預(yù)見螺旋槳滑流對機翼的氣動彈性將產(chǎn)生強烈干擾.

針對幾何非線性氣動彈性問題，Patil等［2-4］建立了考慮失速特性的梁式機翼模型，采用動力學線化方法研究了幾何非線性效應(yīng)對單獨機翼靜、動氣動彈性行為的影響;謝長川等采用推廣的片條理論［5-6］及三維升力線理論［7］對金屬單梁式機翼進行了非線性靜氣動彈性變形計算，并采用線化方法預(yù)測了顫振臨界速度，并通過風洞試驗驗證了計算;Palacios和Cesnik采用計算結(jié)構(gòu)力學(CSD)與計算流體力學(CFD)耦合計算的方法研究了柔性機翼的非線性靜氣動彈性響應(yīng)問題［8］.針對螺旋槳滑流問題，Prandtl等建立了葉素動量理論，對螺旋槳氣動特性進行了較為合理的分析［9］;Weir將渦格法應(yīng)用于螺旋槳滑流的分析，在無黏不可壓定常流動假設(shè)下對螺旋槳滑流及其與機翼的相互干擾進行分析［10］;Agostinelli等采用葉素理論和CFD方法對螺旋槳誘導(dǎo)速度及載荷進行較為詳盡的分析，并分析了螺旋槳滑流與機翼的氣動干擾［11］.近些年國內(nèi)外學者多采用“等效盤”模型或全槳葉建模的CFD方法對該問題進行研究［12-15］.從研究現(xiàn)狀來看，針對柔性機翼幾何非線性氣動彈性研究中并未考慮螺旋槳滑流的影響，而針對螺旋槳滑流的研究又大多忽略了機翼的彈性;Christian的研究雖然考慮了機翼的柔性，但其計算依賴于大量來自CFD計算及風洞實驗的數(shù)據(jù)，對于理論分析及工程設(shè)計的指導(dǎo)意義并不大;CFD方法雖然能夠提供高精度的流場分析，但它不能提供槳葉攻角、滑流誘導(dǎo)速度等原理性的數(shù)據(jù)，而且計算效率低.因此適用于大變形的“螺旋槳/大柔性機翼”系統(tǒng)的氣動彈性分析快速方法有待進一步研究.基于這種研究背景，本文提出了一種具有工程精度、適用于初步設(shè)計階段的快速分析方法，為我國高空長航時無人機的研制做技術(shù)儲備.

1 理論基礎(chǔ)

1.1 螺旋槳滑流及“1P load”的分析方法

對于拉力螺旋槳構(gòu)型的飛機，機翼的部分處于螺旋槳滑流區(qū)，螺旋槳滑流的切向速度分量在機翼上引起當?shù)毓ソ堑脑龃蠡驕p小從而在機翼上產(chǎn)生上洗區(qū)和下洗區(qū);并且滑流軸向速度較之遠前方來流增大，處于滑流區(qū)的機翼會感受到動壓增加.這樣螺旋槳滑流對機翼繞流就產(chǎn)生了干擾.反過來，機翼渦系又會在螺旋槳槳葉處誘導(dǎo)出下洗，也就是說螺旋槳與機翼的氣動干擾是相互的.雖然螺旋槳與機翼間的干擾流場為復(fù)雜的非定常流動，但僅就靜氣動彈性分析，可忽略螺旋槳的交變載荷效應(yīng)及機翼對螺旋槳的誘導(dǎo)作用［16］.

“1P load”是指螺旋槳槳盤面內(nèi)載荷［17］.當氣流沿螺旋槳軸線流動時，螺旋槳上產(chǎn)生拉力、反扭矩，其方向均沿著螺旋槳軸向.但是，當氣流非對稱流經(jīng)螺旋槳槳盤時(如圖1所示)，螺旋槳槳葉產(chǎn)生周期性的氣動力變化，螺旋槳上的凈載荷并不沿軸向，而包含螺旋槳槳盤面內(nèi)的載荷，即“1P load”.準確地預(yù)測螺旋槳載荷對于結(jié)構(gòu)設(shè)計十分重要:“1P load”會產(chǎn)生一個平均值，對機翼升力及扭轉(zhuǎn)都有影響，同樣地，這些變化的面內(nèi)載荷會對結(jié)構(gòu)疲勞造成很大的影響.

圖1 螺旋槳葉素力和速度圖Fig.1 Load and velocity diagram of blade element-momentum theory

1.1.1 Prandtl修正的葉素動量理論

本文基于Prandtl修正的葉素動量(BEM)理論對螺旋槳滑流進行分析.考慮螺旋槳非對稱入流情況，將螺旋槳槳葉沿半徑及方位角兩個尺度離散，螺旋槳槳葉被分成若干個處于方位角ψ下的葉素，如圖1所示，在螺旋槳局部坐標系Oxyz中，原點O位于槳盤中心，Ox軸沿螺旋槳安裝軸向后，Oz軸在槳盤平面內(nèi)向上，螺旋槳沿Ox軸正方向逆時針旋轉(zhuǎn)為正轉(zhuǎn)，反之為反轉(zhuǎn).

葉素處速度的軸向分量:

旋轉(zhuǎn)速度分量:

設(shè)wa，wt分別為軸向誘導(dǎo)速度和旋轉(zhuǎn)誘導(dǎo)速度，引入軸向速度誘導(dǎo)因子a和旋轉(zhuǎn)速度誘導(dǎo)因子a′，即 wa=aV1，wt=a′V2.

根據(jù)動量定理可得葉素處拉力及扭矩為

根據(jù)葉素處氣動力關(guān)系，求得該方位角下拉力與扭矩為

其中，B為槳葉數(shù);c為葉素弦長;φ為氣流相對葉素的入流角;Cn，Ct分別為拉力系數(shù)和扭力系數(shù).

令方程(3)和方程(5)相等、方程(4)和方程(6)相等，并引入 Prandtl葉尖損失因子F［17］，那么可以得到誘導(dǎo)因子的表達式為

滑流區(qū)氣流誘導(dǎo)速度是槳盤處氣流誘導(dǎo)速度的2倍［9］.根據(jù)Prandtl修正的 BEM 理論計算得到的螺旋槳正轉(zhuǎn)時滑流誘導(dǎo)速度分布見圖2:軸向誘速在槳軸兩側(cè)對稱分布，體現(xiàn)了滑流軸向加速的效果;旋轉(zhuǎn)誘速的分布取決于螺旋槳的轉(zhuǎn)向，在槳軸兩側(cè)反對稱分布，圖示為螺旋槳正轉(zhuǎn)結(jié)果.

圖2 螺旋槳滑流誘導(dǎo)速度分布圖Fig.2 Figure of propeller slipstream’s induced velocities distribution

1.1.2 “1P load”的計算

關(guān)于“1P load”的研究還不深入，對其產(chǎn)生的原因、帶來的影響可參考的文獻極少.因此有必要對“1P load”的產(chǎn)生機理做簡要的揭示.

如圖1所示，以螺旋槳在槳盤平面內(nèi)正轉(zhuǎn)為例說明，氣流以相對螺旋槳安裝軸夾角α流過，則來流速度在 Oz軸正向有分量 V∞sinα.［0°，180°］方位角范圍內(nèi)稱為前行槳葉區(qū)(因為氣流相對槳葉的切向速度是槳葉旋轉(zhuǎn)速度和V∞sinα之和)，而在［180°，360°］方位角范圍內(nèi)情況相反.

沿Oz方向的垂直面內(nèi)載荷產(chǎn)生的原因容易揭示:相對于后行槳葉，前行槳葉有更大的氣流攻角和流速，因此產(chǎn)生更大的升力、阻力，這兩個力在旋轉(zhuǎn)平面內(nèi)的垂直分量較大，因此，整個螺旋槳面內(nèi)的凈載荷就是垂直向上的.沿Oy方向的水平面內(nèi)載荷產(chǎn)生的機理并不是很直接，因為它基于復(fù)雜的非定常氣動力，主導(dǎo)原因是誘導(dǎo)速度在槳盤平面內(nèi)的變化.

如圖1螺旋槳槳葉半徑r處的任一葉素，其氣流及氣動力分布情況:dFt和dT分別為葉素切向力及垂直拉力的微分量.由幾何關(guān)系，可以推導(dǎo):

根據(jù)上述分析，求得在V∞=14 m/s，螺旋槳轉(zhuǎn)速ns=135 r/s，αflight=5°時螺旋槳面內(nèi)載荷在一個旋轉(zhuǎn)周期隨方位角的變化情況，如圖3所示.

圖3 “1P load”一周期內(nèi)隨方位角變化圖Fig.3 1P force distribution vs azimuth angle

可見，沿y軸的載荷的平均值基本為0，而沿z軸方向的面內(nèi)載荷有較大的平均值，是“1P load”的主要部分，這部分載荷不僅對機翼升力和扭轉(zhuǎn)變形有貢獻，并且這種交變載荷很可能會影響機翼的穩(wěn)定性，在后續(xù)穩(wěn)定性研究中應(yīng)給予重視.圖4給出沿z向的面內(nèi)載荷隨入流攻角的變化情況，隨著入流攻角的增加，面內(nèi)垂直載荷不斷增加.

圖4 “1P load”中z方向的載荷隨飛行攻角變化圖Fig.4 Vertical force and moment of“1P load”vs flight angle of attack

1.2 滑流對大柔性機翼誘導(dǎo)的模擬

雖然螺旋槳直徑相對機翼展長只是一個小數(shù)，但螺旋槳滑流對整個機翼都有誘導(dǎo)效果.本文計算中簡化認為滑流區(qū)為一個從槳盤向后拖出的以槳盤為底面的圓柱體.旋轉(zhuǎn)的滑流會帶動滑流區(qū)邊界外的氣流旋轉(zhuǎn)，這樣螺旋槳滑流就在整個機翼上誘導(dǎo)出下洗或上洗.處于滑流區(qū)內(nèi)的機翼，其上的誘導(dǎo)速度已通過BEM理論計算得到;對于滑流區(qū)外的機翼，滑流的誘導(dǎo)作用，本文采用蘭金(Rankine)渦核模型［18］模擬，如圖5所示，螺旋槳滑流簡化成一條由螺旋槳旋轉(zhuǎn)中心向后延伸的一根有限渦段.根據(jù)工程估算經(jīng)驗，渦段長度可取為6倍機翼根弦長.根據(jù)畢奧-薩法爾定理，蘭金渦核在滑流區(qū)外任一點的誘導(dǎo)速度為

其中，R為螺旋槳槳盤半徑;r為誘導(dǎo)點到渦核軸線的距離;Γ為蘭金渦核的渦強，由滑流旋轉(zhuǎn)速度確定.

圖5 Rankine渦對空間一點誘導(dǎo)示意圖Fig.5 Geometry induced schematic diagram of the Rankine vortex

1.3 螺旋槳/大柔性機翼靜氣動彈性分析流程

針對大柔性機翼非線性氣動彈性問題，前期研究工作可參考文獻［19］.這里沿用對柔性機翼的非線性氣彈分析方法，在此基礎(chǔ)上加入螺旋槳推進裝置，建立了“螺旋槳/大柔性機翼”系統(tǒng)的靜氣動彈性分析流程，如圖6所示.通過BEM理論計算得到的螺旋槳滑流通過蘭金渦作用到機翼的空間馬蹄渦系上，接著根據(jù)三維升力線方法得到機翼載荷，將螺旋槳載荷及機翼載荷一并加載到結(jié)構(gòu)有限元模型上，利用NASTRAN非線性靜力分析，得到結(jié)構(gòu)變形.根據(jù)結(jié)構(gòu)變形更新模型的氣動構(gòu)型，并開始新一輪計算.迭代計算以機翼主梁翼尖的垂向變形為收斂條件，直到滿足收斂條件則停止計算.

圖6 螺旋槳/大柔性機翼靜氣動彈性分析流程圖Fig.6 Flow chart for nonlinear static aeroelastic analysis of the propeller/wing system

2 算例分析

2.1 計算模型

以文獻［19］中所用某單梁式柔性機翼的右半翼展模型為算例，運用本文所建立的幾何非線性分析方法，對其在螺旋槳滑流作用下的靜氣動彈性進行分析.機翼外形如圖7所示，具體參數(shù)見文獻［19］.電機選為XM2815A，最大功率95 W，重20 g.螺旋槳選擇5X5E，重4 g，兩葉定距，槳葉翼型為NACA0016.螺旋槳安裝在機翼主梁靠根部143.5 mm處向前伸出48 mm，螺旋槳及其安裝軸材料的剛度很大，在結(jié)構(gòu)分析中忽略螺旋槳部分的彈性效應(yīng).

圖7 大柔性機翼模型結(jié)構(gòu)示意圖Fig.7 Structure of the very flexible rectangular wing

根據(jù)實際結(jié)構(gòu)建立的結(jié)構(gòu)非線性有限元模型如圖8所示.機翼主梁剛度特性用梁單元模擬，質(zhì)量特性用分散于翼面的多個集中質(zhì)量模擬，螺旋槳系統(tǒng)簡化為翼根剛性連接于主梁的短梁，電機質(zhì)量以集中質(zhì)量平均分配到短梁的3個節(jié)點上，螺旋槳質(zhì)量加載于短梁前端點上.初始氣動面模型如圖9所示，沿展向按內(nèi)翼段、滑流區(qū)、外翼段3部分共劃分30個片條.

圖8 模型初始有限元結(jié)構(gòu)圖Fig.8 Initial finite element model of the structure

圖9 模型初始氣動構(gòu)型圖Fig.9 Initial aerodynamic model of the structure

2.2 幾何非線性靜氣彈分析結(jié)果

下面具體討論幾何非線性計算情況下“螺旋槳/大柔性機翼”系統(tǒng)的靜氣動彈性特性，主要包括機翼氣動特性和變形分布等情況.

2.2.1 螺旋槳滑流對機翼氣動力的影響

圖10給出了風速V=14 m/s，飛行攻角5°，螺旋槳正、反轉(zhuǎn)(轉(zhuǎn)速分別為 0，130，135，140 r/s)情況下，機翼根部固支約束，氣動升力、誘導(dǎo)阻力、側(cè)力沿機翼展向的分布圖.當螺旋槳正轉(zhuǎn)時，滑流在內(nèi)側(cè)機翼區(qū)誘導(dǎo)下洗，在外側(cè)機翼區(qū)誘導(dǎo)上洗，使得內(nèi)側(cè)機翼有效攻角減小，外側(cè)機翼有效攻角增加，內(nèi)側(cè)機翼氣動力減小，外側(cè)機翼氣動力增加，當螺旋槳反轉(zhuǎn)時，情況與之相反;側(cè)力是氣動升力沿片條側(cè)向的分量，大小隨升力而變化;另外，滑流區(qū)機翼動壓的增加使得滑流區(qū)內(nèi)升力顯著增加.誘導(dǎo)阻力的分布形式主要取決于機翼氣動片條的下洗角，是機翼渦系和滑流蘭金渦線共同誘導(dǎo)的效果，在滿足物面邊界條件情況下，機翼的馬蹄渦系分布發(fā)生變化，使下洗角分布呈“M”型，因此誘導(dǎo)阻力分布如圖10(b)所示.為進一步說明，文中還給出了轉(zhuǎn)速為135 r/s時氣動模型各片條的蘭金渦誘導(dǎo)攻角、馬蹄渦誘導(dǎo)攻角、合成有效攻角的分布圖，如圖11所示.可見，繞流下洗角在滑流區(qū)受馬蹄渦、滑流及動壓增加的共同影響，在滑流區(qū)外主要受馬蹄渦的影響.

圖10 機翼氣動力沿翼展分布圖Fig.10 Distribution of aerodynamic force in the three orientations along wingspan

圖11 氣動片條攻角沿翼展變化圖Fig.11 Distribution of strip’s attack angle distribution along the wingspan

2.2.2 螺旋槳滑流對機翼變形的影響

圖12給出了相同計算情況下機翼主梁沿垂向、側(cè)向及縱向的變形.在螺旋槳拉力、扭矩及機翼氣動載荷、結(jié)構(gòu)重力的共同作用下結(jié)構(gòu)發(fā)生變形.當螺旋槳正轉(zhuǎn)時，隨著轉(zhuǎn)速的增加，由于機翼氣動升力、側(cè)力的增加，機翼主梁沿垂向、側(cè)向的變形顯著增大.總的來說，在算例構(gòu)型下，螺旋槳正轉(zhuǎn)對機翼為氣動加載，并且轉(zhuǎn)速越大，升力增加的效果越明顯;螺旋槳反轉(zhuǎn)與此相反.而機翼沿縱向變形受氣動阻力、螺旋槳拉力的共同影響，隨著轉(zhuǎn)速增加，螺旋槳拉力增大，在一定轉(zhuǎn)速后，拉力增大程度使得螺旋槳安裝處的機翼縱向向前變形.

圖12 機翼主梁3個方向變形示意圖Fig.12 Wing spar’s deformations in three orientations along wingspan

圖13還給出螺旋槳轉(zhuǎn)速135 r/s時，機翼主梁扭轉(zhuǎn)變形情況.可見，螺旋槳正轉(zhuǎn)還使得機翼扭轉(zhuǎn)變形增大.螺旋槳轉(zhuǎn)速越大，即前進比越小，螺旋槳拉力、扭矩及滑流均增強，對機翼氣動力、變形的干擾作用越強.

圖13 機翼主梁扭轉(zhuǎn)Fig.13 Wing spar’s torsional deformation

由以上計算結(jié)果可見，螺旋槳滑流對機翼氣動載荷有很大的干擾作用，使氣動力在機翼上重新分布，進而影響結(jié)構(gòu)變形特點;另外，螺旋槳轉(zhuǎn)動引起的氣動加載和減載效果，與螺旋槳的安裝位置、旋轉(zhuǎn)方向等因素有關(guān)，是一個多因素共同決定的結(jié)果.在工程應(yīng)用中若不考慮螺旋槳滑流效應(yīng)，勢必會造成結(jié)構(gòu)設(shè)計偏差，進而影響飛機的總體性能.因此，這一問題必須加以重視.

3 結(jié)論

1)螺旋槳在非對稱入流條件下會產(chǎn)生槳盤面內(nèi)載荷，其中沿槳盤垂直向上的載荷為主要載荷，它對機翼升力、扭轉(zhuǎn)變形及結(jié)構(gòu)穩(wěn)定性都有影響，應(yīng)該在設(shè)計初期就給予足夠重視.

2)由算例結(jié)果知，拉力螺旋槳在機翼流場中增加了一個螺旋渦并增加了機翼當?shù)貋砹鲃訅?，對柔性機翼靜氣彈特性產(chǎn)生影響，對機翼縱向變形的影響主要取決于螺旋槳載荷，而對機翼氣動力、垂向及展向變形的影響主要取決于螺旋槳滑流.

3)螺旋槳對柔性機翼的彈性扭轉(zhuǎn)、氣動力分布、翼尖位移和扭轉(zhuǎn)角等均產(chǎn)生較大的影響，使機翼的幾何非線性效應(yīng)更加顯著，對機翼的性能造成了一定影響，需要在設(shè)計初期就給予重視.

References)

［1］ 謝長川，楊超.大展弦比飛機的幾何非線性氣彈問題［C］//第七屆全國空氣彈性學術(shù)交流會論文集.北京:中國空氣動力學會，2001:97-102.Xie C C，Yang C.Geometric nonlinear static aeroelasticity problems of high-aspect-ratio planes［C］//7th National Conference on Aeroelastic Academic Communication.Beijing:Aerodynamic Society of China，2001:97-102(in Chinese).

［2］ Patil M J，Hodges D H.On the importance of aerodynamic and structural nonlinearities in aeroelastic behavior of high-aspect-ratio wings，AIAA-2000-1448［R］.Altanta:AIAA，2000.

［3］ Patil M J，Hodges D H，Cesnik C E S，et al.Limit cycle oscillations in high-aspect-ratio wings［C］//40th Structures Structural Dynamics，and Materials Conferernce and Exhibit.St.Louis:AIAA，1999，3:2184-2194.

［4］ Patil M J，Hodges D H，Cesnik C E S.Characterizing the effects of geometrical nonlinearities on aeroelastic behavior of high-aspect-ratio wings［C］//International Forum on Aeroelasticity and Structural Dynamics.Williamsburg，VA:［s.n.］，1999:501-510.

［5］ 謝長川.飛行器氣動彈性穩(wěn)定性靜/動耦合理論與試驗研究［D］.北京:北京航空航天大學，2009.Xie C C.Static/dynamic coupling theory and test study of aircraft aeroelastic stability［D］.Beijing:Beihang University，2009(in Chinese).

［6］ Xie C C，Leng J Z，Yang C.Geometrical nonlinear aeroelastic stability analysis of a composite high-aspect-ratio wing［J］.Shock and Vibration，2008，15(3-4):325-333.

［7］ 劉燚，謝長川.大展弦比機翼的幾何非線性靜氣動彈性分析方法［C］//第十二屆全國空氣彈性學術(shù)交流會會議論文集.北京:中國空氣動力學會，2011:227-233.Liu Y，Xie C C.Anlysis method for geometric nonlinear static aeroelasticity of high-aspect-ratio wings［C］//12th National Conference on Aeroelastic Academic Communication.Beijing:Aerodynamic Society of China，2011:227-233(in Chinese).

［8］ Palacios R，Cesnik C E S.Nonlinear aeroelastic modeling and experiments of flexible wings，AIAA-2006-2186［R］.Newport:AIAA，2006.

［9］ 劉沛清.空氣螺旋槳理論及其應(yīng)用［M］.北京:北京航空航天大學出版社，2006:56-63.Liu P Q.Air propeller theory and application［M］.Beijing:Beihang University Press，2006:56-63(in Chinese).

［10］ Weir D S.Wing loads induced by a propeller wake，AIAA-86-1967［R］.Seattle:AIAA，1986.

［11］ Agostinelli C，Liu C H，Allen C B，et al.Propeller-flexible wing interaction using rapid computational methods，AIAA-2013-2418［R］.San Diego:AIAA，2013.

［12］ 段中喆，劉沛清，屈秋林.某輕載螺旋槳滑流區(qū)三維流場特性數(shù)值研究［J］.控制工程，2012，19(5):836-841.Duan Z Z，Liu P Q，Qu Q L.Numerical research on 3-D flow field characteristics within the slipstream of a low loaded propeller［J］.Control Engineering of China，2012，19(5):836-841(in Chinese).

［13］ 鄂秦，楊國偉，李鳳蔚，等.螺旋槳滑流對飛機氣動特性影響的數(shù)值分析［J］.西北工業(yè)大學學報，1997，15(4):511-516.E Q，Yang G W.Li F W，et al.On coupling effect of two vortex systems of Chinese aircraft with turbo-propellers［J］.Journal of Northwestern Polytechnical University，1997，15(4):511-516(in Chinese).

［14］ 左歲寒，楊永.螺旋槳滑流對帶后緣襟翼機翼氣動特性影響的數(shù)值分析［J］.航空計算技術(shù)，2007，17(1):54-57.Zuo S H，Yang Y.Numercal simulation of propeller/high-lift system interaction［J］.Aeronautical Computing Technique，2007，17(1):54-57(in Chinese).

［15］ 段義乾，史愛明.一種新型的螺旋槳滑流激勵盤模型的研究方法［J］.西北工業(yè)大學學報，2012，30(6):841-846.Duan Y Q，Shi A M.A new and effective actuator disk model approach for the simulation of propeller slipstream［J］.Journal of Northwestern Polytechnical University，2012，30(6):841-846(in Chinese).

［16］ Kroo I.Propeller-wing integration for minimum induced loss［J］.Journal of Aircraft，1986，23(7):561-565.

［17］ Ortun B，Boisard R.In-plane airloads of a propeller with inflow angle:prediction vs experiment，AIAA-2012-2778［R］.New Orleans:AIAA，2012.

［18］ 黃國偉.非定常漩渦空氣東流額理論及應(yīng)用［M］.上海:上海交通大學出版社，1994:56-61.Huang G W.Unsteady vortex aerodynamics theory and application［M］.Shanghai:Shanghai Jiao Tong Unversity Press，1994:56-61(in Chinese).

［19］ Xie C C.Theoretic analysis and experiment on aeroelasticity of very flexible wing［J］.Science China，2012，55(9):2489-2500.