陳金池, 黃研昕
(中國航空工業(yè)空氣動力研究院高速高雷諾數(shù)實驗室, 遼寧沈陽 110034)
顫振分析屬于氣動彈性動穩(wěn)定性問題, 是高空長航時(high altitude long endurance, HALE)無人機設(shè)計的關(guān)鍵。這類無人機普遍采用大展弦比機翼設(shè)計, 具有結(jié)構(gòu)質(zhì)量小、 柔性大的特點, 在氣動載荷作用下容易發(fā)生大幅彈性變形。大變形會引起結(jié)構(gòu)剛度的非線性變化, 還會使氣動面發(fā)生曲面變形, 影響氣動載荷的大小和分布, 進(jìn)而導(dǎo)致顫振特性與小變形時顯著不同[1-3]。
傳統(tǒng)的概念設(shè)計往往在小變形范圍內(nèi)進(jìn)行線性顫振分析, 采用偶極子格網(wǎng)法(double-lattice method, DLM)等平面氣動模型即可獲得良好的精度[4-5]。然而, 平面氣動模型不能準(zhǔn)確描述曲面變形的影響, 在大柔性機翼的分析中存在氣動載荷失真, 顫振預(yù)測不準(zhǔn)的問題[6]。
與平面氣動模型相比, 非定常曲面渦格法(unsteady vortex-lattice method, UVLM)通過四邊形渦格模擬氣動面及尾流區(qū)域, 并通過對渦格的更新實現(xiàn)對曲面氣動面、 曲面尾流區(qū)域的模擬。在此基礎(chǔ)上, UVLM還可以結(jié)合自由尾渦模型[7], 考慮包括翻轉(zhuǎn)在內(nèi)的復(fù)雜尾流外形(見圖1)。由于建模簡單、 計算效率高、 對流場外形捕捉較準(zhǔn)確, UVLM在大變形顫振分析中具有突出優(yōu)勢[8-9], 適用于概念設(shè)計階段。近年來, 許多學(xué)者還對UVLM在陣風(fēng)響應(yīng)[10-11]、 顫振抑制[12]、 可變形機翼飛行器[13-14]等領(lǐng)域的應(yīng)用進(jìn)行了研究。此外, 片條理論和CFD技術(shù)也具備曲面氣動建模的能力[15-16], 但分別受限于二維假設(shè)和計算成本, 在概念設(shè)計中存在限制。
圖1 UVLM自由尾渦Fig. 1 Free wake of UVLM
標(biāo)準(zhǔn)的UVLM氣動載荷求解是時域推進(jìn)形式, 在處理大型顫振計算時會涉及大量的迭代及氣動載荷更新, 這無疑會削弱UVLM在計算效率方面的優(yōu)勢。為了加快求解速度, 可以在大變形平衡位置對UVLM的狀態(tài)空間方程進(jìn)行線性化處理, 利用氣動載荷靈敏度實現(xiàn)氣動載荷的快速更新[17]。同時, 線性化的氣動載荷方程也可以與線性結(jié)構(gòu)動力學(xué)方程耦合, 利用特征值進(jìn)行顫振分析[18-20], 提高計算效率。盡管在平衡位置附近進(jìn)行線性分析, 但UVLM模擬的尾流區(qū)域仍可表現(xiàn)為曲面狀態(tài), 氣動載荷計算相比平面氣動模型更加準(zhǔn)確。此外, UVLM氣動載荷靈敏度和線性化方程還可與其他模塊耦合, 為氣動彈性控制及優(yōu)化設(shè)計提供參考[21-23]。
作為線性化處理的一部分, 早期的氣動載荷靈敏度常常通過有限差分法求解。這種方法對自變量擾動的要求較高, 計算效率較低。近年來, 基于伴隨矩陣的靈敏度求解思路備受CFD研究人員的青睞[24-26]。但是, 這種思路涉及方程的變分, 對存儲空間要求較高??紤]到UVLM方程的氣動載荷與輸入量、 狀態(tài)量存在明確的因式分解關(guān)系。通過鏈?zhǔn)椒▌t即可推導(dǎo)UVLM氣動載荷的解析靈敏度。國內(nèi)外學(xué)者針對面元法的氣動載荷靈敏度開展過一些研究[27-30]。其中, 針對UVLM的研究主要考慮了氣動載荷有關(guān)節(jié)點形狀和位置的靈敏度[21,31], 缺少對節(jié)點速度和環(huán)量強度的考慮, 對UVLM線性化方程中的Jacobi矩陣設(shè)計不夠完整。
針對這一情況, 本文將首先構(gòu)造線性化的UVLM狀態(tài)空間方程, 明確Jacobi矩陣包含的元素, 并著重推導(dǎo)氣動載荷有關(guān)節(jié)點位置、 節(jié)點速度和環(huán)量強度的解析靈敏度。在此基礎(chǔ)上, 本文通過數(shù)值仿真探討了這些靈敏度矩陣的分布規(guī)律, 利用有限差分法驗證了解析靈敏度的計算效率和精度。
如圖2所示, UVLM將機翼的平均氣動面離散為一系列結(jié)構(gòu)網(wǎng)格(黑色實線), 在結(jié)構(gòu)網(wǎng)格下游1/4弦長處布置附著渦渦格(紅色虛線), 通過對附著渦的更新模擬氣動面曲面效應(yīng)。隨著計算過程的推進(jìn), 附著渦會從后緣脫落形成尾渦(黑色虛線), 集合附著渦與尾渦的位置即可較準(zhǔn)確地反映氣動面及流場的幾何特征。UVLM在每個渦格內(nèi)設(shè)置渦環(huán), 通過環(huán)量強度完成氣動載荷計算。
圖2 UVLM渦格劃分示意圖Fig. 2 Discrete vortex panels of UVLM
附著渦渦格的幾何中心可視為氣動載荷的控制點。在計算過程中, 需要在附著渦渦格的控制點處施加不可穿透邊界條件, 用以獲得附著渦的環(huán)量強度。與平面氣動模型不同, UVLM附著渦沿曲面分布, 邊界條件中的誘導(dǎo)速度并非豎直向下, 而是沿著渦格的法向量向下。以t+1時刻為例, 不可穿透邊界條件可表示為
(1)
其中,Γbv和Γwk為兩組列向量, 分別表示附著渦和尾渦的環(huán)量強度; 矩陣Wbb和Wbw分別表示附著渦和尾渦關(guān)于氣動載荷作用點的影響系數(shù), 用以計算誘導(dǎo)速度,Wbb和Wbw中的元素由渦格節(jié)點位置決定, 可通過Biot-Savart定律獲得;Unc表示附著渦控制點處的非環(huán)量速度, 包括剛體運動速度、 彈性振動速度、 空氣擾動速度等。上標(biāo)n表示渦格的法向量方向。
尾流區(qū)域的演化包括兩部分, 即附著渦的脫落和已有尾渦的擴散。演化方程如下
(2)
其中,Abv和Awk是尾渦渦強演化矩陣;Bbv和Bwk是尾渦節(jié)點演化矩陣, 均由1和0組成;ζbv和ζwk分別表示附著渦和尾渦節(jié)點的位置,vwk表示尾渦節(jié)點的局部速度, 由流場的誘導(dǎo)速度及外部擾動組成。為避免尾渦數(shù)量過多帶來的額外計算成本, 可在尾流區(qū)域達(dá)到一定長度時設(shè)置截斷, 忽略遠(yuǎn)處尾渦的影響。
在每個控制點處, UVLM的氣動載荷升力分量L垂直于非環(huán)量速度向上, 誘導(dǎo)阻力D沿非環(huán)量速度方向。L和D的矩陣表達(dá)式如下
(3)
其中,E1(2)是由0, 1, -1組成的稀疏矩陣, 用以計算有效渦強。G1(2)是對角矩陣, 用以計算投影面積, 由附著渦的幾何信息和局部攻角決定。U1(2)(3)是誘導(dǎo)速度對角矩陣。根據(jù)本文需要,G1(2)及U1(2)(3)的表達(dá)形式如下
(4)
將方程(1)~(4)整合, 可得到UVLM氣動載荷在離散時刻的狀態(tài)空間形式
fA(xt+1,ut+1)=gA(xt,ut)
yt+1=hA(xt+1,ut+1)
(5)
其中,u,y,x分別表示方程的輸入量、 輸出量、 狀態(tài)量, 分別為
(6)
方程(5)可以與結(jié)構(gòu)動力學(xué)方程耦合, 形成氣動彈性方程。
圖3 UVLM氣動力線性化過程的參數(shù)關(guān)系圖Fig. 3 Parameter relationship of aerodynamic linearization by the UVLM
將方程(5)中的變量表示為在平衡位置x0處的基準(zhǔn)量與擾動量之和
(7)
再將方程(7)代入方程(1)~(3), 基于上述參數(shù)關(guān)系可得到線性化的狀態(tài)空間方程
(8)
其中, Δt為時間步長。
將方程(8)中的Jacobi矩陣整合, 線性化的狀態(tài)空間方程可表示為
EΔxt+1=AΔxt+BΔut+1, Δyt+1=CΔxt+1+DΔut+1
(9)
(10)
(11)
其中
(12)
(13)
(14)
其中, 下標(biāo)i表示附著渦渦格節(jié)點索引, 每個元素表示節(jié)點i發(fā)生擾動對附著渦渦格k的影響。
按類似的方法可推導(dǎo)方程(12)~(13)中其余偏導(dǎo)項, 并獲得方程(11)中氣動載荷解析靈敏度的表達(dá)式, 在此不一一贅述。相比其他計算方法, 采用鏈?zhǔn)椒▌t推導(dǎo)的解析靈敏度可直接使用矩陣運算符, 避免了數(shù)值計算中的循環(huán)過程。
將長直機翼置于均勻分布的流場中, 忽略彈性變形和重力的影響, 通過氣動載荷系數(shù)驗證UVLM建模的準(zhǔn)確性, 再對平衡狀態(tài)下的機翼施加小擾動, 驗證解析靈敏度的計算精度及效率。
構(gòu)建5個攻角為5°, 弦長為1 m的長直機翼, 這些機翼的展弦比分別為4, 8, 12, 20, 100。將這些機翼依次置于速度50 m/s的均勻流場中, 通過UVLM模型獲得升力系數(shù)與誘導(dǎo)阻力系數(shù), 并與Katz等[32]提供的參考值進(jìn)行對比。計算中, 對UVLM尾流區(qū)域在9 m處進(jìn)行截斷, 忽略遠(yuǎn)處尾渦的影響。
在第1步迭代計算中, 流場從靜止開始起風(fēng), 機翼會受到一個較大的脈動載荷, 導(dǎo)致升力系數(shù)出現(xiàn)較大的初始值。隨著尾渦的生成, 這一脈動載荷帶來的影響逐漸減小。為便于比較, 本文將初始升力系數(shù)設(shè)置為0.5。對比結(jié)果如圖4所示, UVLM升力在迭代過程中逐漸達(dá)到與參考值接近的穩(wěn)定狀態(tài), 并且迭代步數(shù)隨著展弦比的增大而增加, 穩(wěn)定時的升力系數(shù)也隨著展弦比的增大而接近5°攻角下的二維翼型計算值0.55。
圖4 UVLM升力系數(shù)與展弦比的關(guān)系Fig. 4 Variation of lift coefficient with aspect ratio
在誘導(dǎo)阻力系數(shù)的對比中, 選取展弦比為8的機翼, 分析了誘導(dǎo)阻力系數(shù)及分量的變化情況。圖5顯示計算結(jié)果與Katz和Poltkin的結(jié)果吻合, 進(jìn)一步驗證了UVLM建模的準(zhǔn)確性。
圖5 UVLM誘導(dǎo)阻力系數(shù)及分量瞬態(tài)變化(展弦比8)Fig. 5 Variation of induced drag coefficient and components at an aspect ratio of 8
將機翼的平均氣動面劃分為20個展向單元和2個弦向單元。在UVLM時間推進(jìn)計算中, 選取0.025 s作為UVLM時間推進(jìn)步長。計算結(jié)果表明, 算例機翼的氣動載荷在起風(fēng)1.5 s后達(dá)到穩(wěn)定。因此, 選取2 s時刻穩(wěn)定狀態(tài)作為靈敏度計算的平衡狀態(tài)。在機翼向下游方向30 m處設(shè)置截斷, 忽略30 m外的尾渦影響。尾流區(qū)域如圖6所示。計算中, 附著渦及尾渦的編號順序如圖7所示。
圖6 UVLM尾渦在2 s時刻的渦格情況Fig. 6 UVLM wake at 2 s
圖7 UVLM渦格及節(jié)點編號Fig. 7 Numbering scheme of UVLM vortex panels and nodes
(a) Analytical sensitivity
表2 升力關(guān)于渦格節(jié)點位置靈敏度結(jié)果及對比(y)
表3 升力關(guān)于渦格節(jié)點位置靈敏度結(jié)果及對比(z)
表4 升力關(guān)于渦格節(jié)點速度靈敏度結(jié)果及對比(y)
表5 升力關(guān)于渦格節(jié)點速度靈敏度結(jié)果及對比(z)
解析靈敏度與有限差分法的計算時間對比見表6, 采用解析靈敏度的計算效率顯著高于有限差分法。本文最大僅涉及40×63的靈敏度矩陣計算, 在這一情況下, 解析解靈敏度的計算效率大約是有限差分法的5~6倍。這種優(yōu)勢在進(jìn)行更大型的計算時會更加明顯, 這是兩種方法的理論基礎(chǔ)決定的。解析靈敏度在計算時直接使用矩陣運算符, 不需進(jìn)行任何循環(huán)或迭代。而有限差分法作為一種數(shù)值計算方法, 必須對一個變量進(jìn)行擾動, 再通過多次循環(huán)計算才能獲得靈敏度結(jié)果, 過程中的計算消耗遠(yuǎn)大于解析解。
表6 解析解與有限差分法計算時間對比
本文詳細(xì)描述了UVLM的狀態(tài)空間建模過程, 并推導(dǎo)了UVLM升力關(guān)于附著渦強度以及渦格節(jié)點運動信息的解析靈敏度。經(jīng)驗證, 解析靈敏度的數(shù)值大小、 分布情況均與有限差分法一致, 各靈敏度項的計算誤差不超過0.4%。相比有限差分法, 解析靈敏度直接進(jìn)行矩陣運算, 避免了擾動量選取不當(dāng)帶來的收斂性問題, 也避免了大量數(shù)值計算帶來時間消耗, 體現(xiàn)出高精度、 高效、 高魯棒性的優(yōu)勢。
氣動力解析靈敏度可進(jìn)一步與結(jié)構(gòu)響應(yīng)方程耦合, 構(gòu)建氣動彈性分析系統(tǒng)。除了可大幅提高大型系統(tǒng)的穩(wěn)定性分析效率, 解析靈敏度在氣動彈性優(yōu)化設(shè)計及流動控制領(lǐng)域都具有良好的應(yīng)用價值。