耦合分布螺旋槳的大柔性機翼靜氣彈研究

吳軒,周洲,王正平

(西北工業大學 航空學院,陜西 西安 710072)

面對未來航空器高經濟性、低能耗性、低噪聲、低污染的“綠色航空發展要求”,采用分布式螺旋槳電推進系統的高空長航時無人機正被應用于各種場景[1]。此類飛機普遍具有大展弦比、低翼載、大柔性的特點,在飛行載荷作用下,機翼會產生較大的幾何非線性變形,傳統的線性氣動彈性分析手段已經不再適用。同時,雖然單個螺旋槳相對于大展弦比機翼尺度較小,但沿展向分布的多個螺旋槳,其載荷和滑流都會對機翼的流場、載荷分布、變形分布有較大的干擾[2]。

針對大柔性機翼的幾何非線性氣彈問題,Patil[3]采用本征梁模型和有限狀態入流模型系統研究了機翼結構的幾何非線性對結構剛度、氣動力、氣動彈性的影響。Palacios等[4]采用應變梁模型與計算流體力學(CFD)相耦合的方法,研究了機翼的非線性靜氣彈響應問題。張馳[5]使用共旋轉梁模型對非線性結構進行建模,在有限元的總體和單元坐標系中重新定義了飛行力學的地軸系、體軸系和風軸系,研究了耦合飛行力學下的氣彈響應問題。針對螺旋槳對機翼氣動彈性的影響研究,Hodges等[6]將螺旋槳推力定義為具有常值的集中力,研究了推力對機翼彎扭顫振的影響,發現推力對彎扭顫振同時具有促進和抑制作用。Agostinelli等[7]采用葉素理論和CFD方法,并結合線性有限元模型,分析了螺旋槳載荷、滑流誘導速度對彈性機翼的氣動干擾。Teixeira等[8]使用黏性渦粒子法將螺旋槳滑流影響和陀螺效應引入到柔性機翼的氣彈分析中,指出螺旋槳滑流增加了機翼變形且陀螺效應的影響可以忽略。謝長川等[9]使用動量理論和蘭金渦核模擬了螺旋槳滑流對柔性機翼的氣動影響,結果表明旋轉的滑流會改變機翼當地繞流攻角,從而影響機翼氣動力分布和變形分布。

從現有研究現狀來看,對于柔性機翼在分布式螺旋槳影響下的非線性氣彈響應問題,傳統研究通常僅考慮螺旋槳的拉力,并將其作為集中力處理,并未考慮滑流的影響。文獻[10]指出對于柔性飛機,未考慮螺旋槳滑流影響是仿真結果存在誤差的重要原因之一。而對于螺旋槳滑流對機翼的影響,現有研究更注重二者之間的氣動干擾問題,忽略了柔性機翼的幾何非線性變形特征。同時,螺旋槳分布位置所帶來的氣動彈性影響也亟待研究。

對此,本文建立了一種考慮分布式螺旋槳載荷和滑流影響的柔性機翼非線性靜氣彈分析方法。首先,利用單個螺旋槳相對于大展弦比機翼尺度較小、轉速較大的特性,采用滑流管模型模擬滑流對機翼的誘導速度,避免了螺旋槳復雜的尾跡計算。螺旋槳環量分布通過渦流葉素理論獲取,機翼氣動力則使用曲面渦格法進行計算。其次,在共旋轉法的節點坐標系中,推導出展向不同位置處螺旋槳的局部柱坐標系,完成分布式螺旋槳與機翼結構在非線性大變形下的耦合建模。隨后,對螺旋槳與機翼耦合的構型進行氣動升力計算,將結果與現有方法對比,驗證了所提氣彈模型的準確性。最后,通過分布式螺旋槳與柔性機翼耦合的算例,對比分析了有無螺旋槳載荷和滑流以及螺旋槳分布位置對大柔性機翼靜氣彈響應的影響,再次表明在幾何非線性氣彈模型中引入分布式螺旋槳影響的重要性。

1 螺旋槳/機翼氣動模型

1.1 渦流葉素理論

采用結合渦流模型的片條理論對小負荷螺旋槳進行氣動建模,假定螺旋槳葉片與機翼一樣,沿槳葉展向分布著環量變化的附著渦系Γ(r),在展向每一位置處均有自由渦溢出,并在后緣附近形成一個扭曲的渦面。對于每個單元中的葉素,其半徑為r,弦長為b,受力和速度分解圖如圖1所示。

圖1 葉素受力和速度分解圖

在實際飛行中,葉素運動軌跡為螺旋線,前飛速度為V1,旋轉角速度為Ω,槳盤平面內的切向速度V2=Ωr,則氣流相對于葉素的幾何合成速度Vgeo為

(1)

Vgeo與旋轉平面之間的夾角即為初始入流角φ0

(2)

葉素受軸向誘導速度Va、切向誘導速度Vt影響,其實際入流速度Vreal為

(3)

經誘導后,氣流的實際入流角φ為

(4)

令葉素的安裝角為θ(葉素弦線方向與旋轉平面之間的夾角),尾渦誘導產生的干涉角為β,則氣流相對于葉素的當地迎角α為

α=θ-β-φ0

(5)

在迎角α下,當地升阻力系數和馬赫數及雷諾數有關,即Cl=Cl(α,Re,Ma),Cd=Cd(α,Re,Ma)。根據渦流理論,經槳尖環量損失修正[11]后的當地葉素環量為

(6)

式中

式中:R為螺旋槳半徑;NB為螺旋槳槳葉數目。

根據庫塔-茹科夫斯基升力定理,當地環量分布與升力系數關系為

(7)

聯立求解(6)～(7)式,可得螺旋槳環量分布Γ(r),進一步可以得到該單元處的拉力和扭矩

(8)

將單元處拉力dT和扭矩dM沿徑向和周向積分,即可得到螺旋槳總拉力和總扭矩

(9)

式中,Rc為螺旋槳輪轂半徑。

1.2 螺旋槳滑流管尾跡模型

對于本文研究的帶有分布式螺旋槳的大展弦比機翼來說,機翼飛行速度遠遠小于螺旋槳旋轉速度,因此可以忽略螺旋槳的非定常效應,同時假設機翼附近的螺旋槳尾跡不發生變形,從而避免了復雜的尾跡更新過程,提高氣彈分析計算效率。

螺旋槳滑流對機翼的氣動影響主要分為2個部分:①滑流的軸向誘導速度分量使得機翼當地流動速度增加、動壓增大;②滑流的切向誘導速度分量引發局部上洗或者下洗流動。因此,如圖2所示,本文采用Veldhuis[12]提出的滑流管模型,將螺旋槳尾渦分成軸向渦系和切向渦系2個部分。

圖2 螺旋槳尾跡渦系模型

使用畢奧薩伐爾定律計算微元渦段對空間任一點Q(xQ,yQ,zQ)的誘導速度

(10)

式中,dVQ為微元渦段對Q點的誘導速度;Γ表示渦段環量大小;S表示點Q到渦段的向量;dl表示微元渦段向量。

圖2中的軸向渦和切向渦可以表示為

(11)

式中:Γ(r)為槳葉的環量分布,由1.1節中給出;NB為槳葉數目;Ω表示螺旋槳旋轉速度;U∞表示自由來流速度。

在螺旋槳局部柱坐標系中,對(10)式沿槳盤周向φ和縱向z積分,并結合(11)式,最終可得到切向渦對Q點的誘導速度為

(12)

軸向渦對Q點的誘導速度為

(13)

式中:b=rsinφ-xQ,c=-rcosφ-yQ。

槳盤處的誘導速度是個小量,可以忽略,因此滑流管模型對Q點的最終誘導速度為

(14)

1.3 機翼曲面渦格法氣動力計算

鑒于螺旋槳與機翼的相對尺寸較小且螺旋槳旋轉速度遠高于機翼的飛行速度,故忽略機翼對螺旋槳的誘導速度。同時,螺旋槳之間距離較遠,故忽略螺旋槳之間的氣動干擾,著重關注螺旋槳載荷對柔性機翼氣動的影響。

氣動載荷使柔性機翼的升力面發生彎曲扭轉變形從而形成空間曲面,使用曲面渦格法對其進行氣動力計算。如圖3所示,將氣動面沿展向和弦向劃分網格,共計M個單元。在每個氣動單元內布置渦環基本解,控制點位于渦格1/4弦向處,此布置形式自然滿足Kelvin條件和流體邊界條件。

圖3 考慮滑流影響的機翼曲面渦格

在當前構型下,升力面附著渦滿足物面不可穿透條件,即在每個渦格的控制點處物面法方向nj上的誘導速度分量為零,即滿足

(15)

2 螺旋槳/機翼結構模型

2.1 機翼幾何非線性建模

為了耦合螺旋槳模型,采用基于位移變量的共旋轉法對柔性機翼結構進行幾何非線性分析,其物理意義更加明確。如圖4所示,共旋轉法通過局部節點坐標系和單元坐標系來描述結構的剛體運動和變形位形。

圖4 共旋轉法單元與節點坐標系

(16)

式中,R(α)是廣義旋轉位移對應的旋轉矩陣。根據Rodrigues對旋轉矩陣的參數化描述方法,可以得到旋轉矩陣和旋轉偽矢量的關系式

(17)

式中:α=|α|為旋轉偽矢量的模;S(α)為斜對稱矩陣

R(γ)=UTT

(18)

(19)

(20)

所有推導的節點坐標系和單元坐標系都是結構位形的隨動坐標系,隨結構變形實時更新。

根據(20)式,三維梁在局部單元坐標系下的位移為pl,總體坐標系下的位移為p:

(21)

對(21)式進行微分處理,即可得到總體坐標系下的位移增量δp和單元坐標系下的位移增量δpl的轉換關系為

δpl=Fδp

(22)

式中,F為轉換矩陣。

令梁單元在總體坐標系下的內力為q,在單元坐標系下的內力為ql,根據虛功原理,有

(23)

在單元坐標系下,共旋轉法認為力與位移的平衡關系是線彈性的,同時根據(22)式,可以得到

q=FTql=FTKlpl

(24)

式中,Kl是線性剛度矩陣。

對(24)式進行變分處理即可得到增量平衡方程和對應的切向剛度矩陣

(25)

式中:FTKlF為材料剛度矩陣;Ktσ(ql)為幾何剛度矩陣且為單元局部內力ql的函數,KT為共旋轉法的切線剛度矩陣,具體表達式見文獻[5]。

2.2 螺旋槳機翼耦合建模

上文1.2節中推導的滑流誘導速度是在螺旋槳局部柱坐標系下進行的,然而由于機翼的大變形,螺旋槳的位置和方向隨機翼結構變形實時變化,需要考慮在結構模型中推導螺旋槳局部柱坐標系的表達形式,完成慣性坐標系、結構單元節點坐標系和螺旋槳柱坐標系的統一以及載荷變形等數據傳遞。

如圖5所示,對于一般安裝在機翼前緣的螺旋槳推進系統,在慣性系下,螺旋槳局部柱坐標系原點位于中心點Op處,平面極坐標rp(rp,φp)位于槳葉旋轉平面,zp方向垂直于旋轉平面與推力線方向OTP反向。在基于共旋轉法的機翼非線性結構模型中,螺旋槳安裝位置處機翼節點坐標系為Ui(Uix,Uiy,Uiz),節點坐標為Oi。

圖5 螺旋槳柱坐標系與結構節點坐標系

(26)

(27)

忽略推進系統旋轉產生的回轉力矩效應,螺旋槳的推力Tp和轉矩Mp在1.1節中可以通過渦流葉素理論獲得。則在慣性系中,螺旋槳安裝位置處結構單元節點所受到的力Fi和力矩Mi為

(28)

3 螺旋槳/機翼非線性靜氣彈分析

3.1 空間梁樣條插值

采用空間梁樣條完成在大柔性變形下,機翼結構節點與氣動節點之間的位移與載荷插值。在共旋轉梁模型中,結構節點平動位移和轉動位移分別為Xs和Ψs。令氣動節點到對應截面結構節點的矢量為ξa,假設機翼在變形過程中,截面保持剛性,則氣動節點位移Xa可以通過結構節點位移插值得到

Xa=Xs+Us(Ψs)ξa=GXs

(29)

式中,Us為截面轉動位移對應的轉換矩陣;G為位移插值矩陣,具體表達式見文獻[14]。

設結構界面的虛位移為δXs,氣動界面的虛位移為δXa,氣動界面的載荷為Fa,結構界面的載荷為Fs,根據虛功原理有

(30)

根據(29)式,可以得到氣動結構節點的位移增量關系式

δXa=GδXs

(31)

將(31)式代入(30)式,得

(32)

根據虛位移的任意性,可以得到氣動載荷與結構載荷之間的插值關系式為

Fs=GTFa

(33)

3.2 非線性靜氣彈分析框架

大柔性飛機從本質上來說是一個分布式的參數系統,全機不同位置處結構位移、氣動響應呈現非線性變化?？紤]柔性飛機的特點以及分布式螺旋槳的影響,本文建立如圖6所示的非線性靜氣彈分析框架并完成相應各個模塊的程序編寫。

圖6 非線性靜氣彈分析框架

整個分析框架分為2個部分:螺旋槳機翼氣動計算模塊和非線性結構計算模塊。在氣動計算模塊,首先通過2.2節中非線性結構與螺旋槳的耦合建模,可以確定不同位置處螺旋槳的實時隨動位形,并對空間氣動曲面進行渦格劃分。其次,通過1.1節中介紹的渦流葉素理論可以獲取各個螺旋槳的載荷和環量分布,并結合1.2節和1.3節中的滑流管模型和曲面渦格法,完成考慮滑流影響的大柔性機翼氣動特性計算。在結構計算模塊,首先通過空間梁樣條將螺旋槳載荷和氣動載荷轉換加載到結構有限元節點上。同時,利用迭代計算的方法求解2.1節中的增量平衡公式(25),直至結構達到變形收斂狀態,最終完成考慮分布式螺旋槳影響的柔性機翼非線性靜氣彈分析并輸出結果。

4 方法驗證及算例分析

4.1 氣彈模型驗證

采用密歇根大學的高空長航時驗證無人機[15](X-HALE)為研究對象,該機翼分為內翼段和外翼段,展長分別為2 m,1 m,弦長為0.2 m。機翼在翼根處有5°的安裝角,外翼段上反10°,其他特性參數見表1。如圖7所示,5個螺旋槳分布在中心點以及機翼半展長的±1/3,±2/3位置處,螺旋槳距離機翼前緣前伸0.2 m、下移0.028 m,相對翼根安裝角為0°,螺旋槳型號為APC11×5.5E,參數見文獻[16]。

圖7 X-HALE機翼示意圖

設定自由來流速度v∞=14 m/s,來流攻角為2°,機翼翼根處固支。設定螺旋槳轉速為6 000 r/min,其拉力為1.70 N。根據收斂性研究,機翼弦向劃分12個渦格,展向劃分96個渦格。根據3.2節中的靜氣彈分析框架,計算機翼在剛性不變性和柔性大變形下的收斂位移和升力分布,計算結果如圖8所示。

圖8 X-HALE機翼位移與升力分布

圖8a)為機翼受來流氣動載荷和螺旋槳載荷共同作用下的垂向位移分布圖,翼尖最大位移約占半展長的35%。本文結果和文獻[15]結果基本一致,靜變形誤差基本可以忽略,驗證了本文非線性靜氣彈分析模型的準確性。

圖8b)為機翼的展向升力系數分布圖,可以看出,滑流軸向誘導速度可以引起翼面加速流動從而使得局部升力變大,切向誘導速度在機翼附近引發上洗和下洗流動也會改變機翼的展向升力分布。機翼的整體升力系數如表2所示,變形后機翼的升力系數相對未變形構型降低17%左右,本文和文獻[15]的結果相對誤差在2%以內。

表2 機翼升力系數對比

升力分布以及升力系數的計算結果表明,本文采用的氣動模型以及編寫的程序是合理準確的,可以有效獲取螺旋槳對機翼的氣動影響。同時,計算時間相對于CFD等高階方法大大降低,僅需30 s左右,本文氣動模型可應用于此類飛機初始設計階段的快速分析。

4.2 柔性機翼耦合分布式螺旋槳算例分析

選用一廣泛被用作基準模型的大展弦比柔性機翼[17]作為研究對象,使用3.2節中的求解框架探究螺旋槳的載荷、滑流、分布位置對大柔性機翼的靜氣彈影響。圖9所示為機翼半模示意圖,10個螺旋槳沿展向對稱布置在機翼半展長的3/8,4/8,5/8,6/8,7/8位置處。所有螺旋槳弦向距離前緣0.5m,垂向安裝距離為0,槳軸與機翼共面,轉向為正(機翼外側為上洗區域)。

圖9 耦合分布式螺旋槳的柔性機翼(半模)示意圖

螺旋槳具體幾何參數見文獻[18],機翼的具體參數見表3。自由來流速度為25 m/s,空氣密度為0.088 9 kg/m3,所有螺旋槳轉速皆為4 400 r/min,機翼翼根處固支,安裝角為4°。

表3 大柔性機翼參數

基于收斂性研究,機翼結構有限元模型共分為32個梁單元,機翼升力面弦向劃分12個單元,展向劃分128個單元,結構和氣動的計算收斂誤差皆在1%以內。

4.2.1 螺旋槳對結構變形的影響

機翼非線性彎曲、扭轉變形計算結果對比如圖10～11所示,在無螺旋槳的情況下,本文計算得到的垂向彎曲變形和Murua等[17]使用SHARP程序計算的結果基本一致,翼尖最大垂向位移誤差為1.2%,表明本文氣彈模型和代碼具有一定的準確性。

圖10 柔性機翼垂向彎曲變形圖

圖11 柔性機翼扭轉變形圖

在只考慮螺旋槳拉力扭矩載荷不考慮滑流影響的情況下,翼尖最大彎曲變形相對于干凈機翼減小9.6%,翼尖最大扭轉變形減小34.3%。這主要是由于較大的彎曲變形使螺旋槳拉力產生機翼結構負扭轉,從而造成機翼局部迎角減小、升力減小,進一步使彎曲變形減小。

當引入滑流的影響時,翼尖最大彎曲變形相較于只考慮螺旋槳載荷的情況增加3.2%,翼尖最大扭轉變形減小7.9%?；鲙頇C翼局部入流速度增加,從而引起升力增加、機翼彎曲變形增加。同時,滑流引起氣動俯仰力矩減小,產生低頭力矩,從而減小扭轉變形。

4.2.2 螺旋槳對氣動特性的影響

圖12給出了在不同迎角下,螺旋槳載荷和滑流對機翼升力系數曲線的影響。由于拉力產生結構負扭轉,機翼局部迎角變小,升力系數相較于干凈機翼降低10%左右?；魇箼C翼局部來流速度增加,導致機翼升力系數相對于只考慮螺旋槳載荷的情形增加約2.5%左右。

圖12 升力系數隨迎角變化曲線

圖13為機翼縱向靜穩定裕度Cmα/Clα隨迎角變化曲線,可以看出螺旋槳的載荷使靜穩定裕度絕對值減小20%～35%,滑流使靜穩定裕度絕對值增加2%～8%,差值隨迎角增大而增大。

圖13 靜穩定裕度隨迎角變化曲線

4.2.3 螺旋槳不同分布位置的影響

分別只保留機翼半展長3/8,4/8,5/8,6/8,7/8位置處對稱的2個螺旋槳工作,關閉其他8個螺旋槳,計算條件不變,圖14～15給出了螺旋槳位于不同位置處機翼升力系數和俯仰力矩系數隨轉速變化曲線。

圖14 螺旋槳位于不同位置處升力系數隨轉速曲線

如圖14所示,當螺旋槳的位置靠近翼根時,滑流的影響更大,導致升力增加;靠近翼尖時,拉力的影響更大,導致升力減小;當螺旋槳位于半展長的3/8到4/8之間時,載荷和滑流對升力系數的影響相互抵消。從曲線斜率可知,越靠近翼尖,轉速對升力線斜率Clα影響越大。如圖15所示,螺旋槳在各個位置均增大了俯仰力矩系數,造成機翼縱向靜穩定裕度減小。同樣,越靠近翼尖,轉速對力矩系數影響越大。

圖15 螺旋槳位于不同位置處俯仰力矩系數隨轉速曲線

5 結 論

1) 本文在基于渦流理論的滑流管模型以及共旋轉法非線性結構模型的基礎上,建立了考慮分布式螺旋槳載荷和滑流影響的大柔性機翼靜氣彈快速分析求解方法。將機翼升力分布和靜氣彈變形等結果與文獻[15]對比,驗證了本文模型具有一定的準確性,能夠有效體現大柔性機翼的非線性變形特征以及分布式螺旋槳與機翼的耦合氣動特性。

2) 分布式螺旋槳的載荷和滑流會對柔性機翼的彎曲扭轉變形、氣動力分布、靜穩定裕度等產生較大的影響,使機翼的非線性效應更加顯著。螺旋槳拉力會產生機翼結構負扭轉,使當地有效迎角減小,造成氣動升力損失和靜穩定裕度降低;滑流則通過增加當地流速和改變繞流攻角,增加了升力和靜穩定裕度。在此類飛機的初始設計階段,必須考慮這些影響,以免造成設計偏差。

3) 螺旋槳在展向不同位置對大柔性機翼的性能影響不同,靠近翼根時增升,靠近翼尖時減升,且越靠近翼尖螺旋槳轉速變化帶來的氣彈影響越顯著,文獻[19]也指出了分布式螺旋槳可以控制機翼升力分布以及抑制機翼彈性變形的潛力。在后續研究中,可以根據本文模型方法和分析結果,通過調節螺旋槳的轉速、安裝位置、分布形式等參數完成大柔性機翼與分布式螺旋槳的耦合優化設計。