基于非定常曲面渦格法的氣動靈敏度和線性化分析

陳金池, 黃研昕

(中國航空工業空氣動力研究院高速高雷諾數實驗室, 遼寧沈陽 110034)

引 言

顫振分析屬于氣動彈性動穩定性問題, 是高空長航時(high altitude long endurance, HALE)無人機設計的關鍵。這類無人機普遍采用大展弦比機翼設計, 具有結構質量小、 柔性大的特點, 在氣動載荷作用下容易發生大幅彈性變形。大變形會引起結構剛度的非線性變化, 還會使氣動面發生曲面變形, 影響氣動載荷的大小和分布, 進而導致顫振特性與小變形時顯著不同[1-3]。

傳統的概念設計往往在小變形范圍內進行線性顫振分析, 采用偶極子格網法(double-lattice method, DLM)等平面氣動模型即可獲得良好的精度[4-5]。然而, 平面氣動模型不能準確描述曲面變形的影響, 在大柔性機翼的分析中存在氣動載荷失真, 顫振預測不準的問題[6]。

與平面氣動模型相比, 非定常曲面渦格法(unsteady vortex-lattice method, UVLM)通過四邊形渦格模擬氣動面及尾流區域, 并通過對渦格的更新實現對曲面氣動面、 曲面尾流區域的模擬。在此基礎上, UVLM還可以結合自由尾渦模型[7], 考慮包括翻轉在內的復雜尾流外形(見圖1)。由于建模簡單、 計算效率高、 對流場外形捕捉較準確, UVLM在大變形顫振分析中具有突出優勢[8-9], 適用于概念設計階段。近年來, 許多學者還對UVLM在陣風響應[10-11]、 顫振抑制[12]、 可變形機翼飛行器[13-14]等領域的應用進行了研究。此外, 片條理論和CFD技術也具備曲面氣動建模的能力[15-16], 但分別受限于二維假設和計算成本, 在概念設計中存在限制。

圖1 UVLM自由尾渦Fig. 1 Free wake of UVLM

標準的UVLM氣動載荷求解是時域推進形式, 在處理大型顫振計算時會涉及大量的迭代及氣動載荷更新, 這無疑會削弱UVLM在計算效率方面的優勢。為了加快求解速度, 可以在大變形平衡位置對UVLM的狀態空間方程進行線性化處理, 利用氣動載荷靈敏度實現氣動載荷的快速更新[17]。同時, 線性化的氣動載荷方程也可以與線性結構動力學方程耦合, 利用特征值進行顫振分析[18-20], 提高計算效率。盡管在平衡位置附近進行線性分析, 但UVLM模擬的尾流區域仍可表現為曲面狀態, 氣動載荷計算相比平面氣動模型更加準確。此外, UVLM氣動載荷靈敏度和線性化方程還可與其他模塊耦合, 為氣動彈性控制及優化設計提供參考[21-23]。

作為線性化處理的一部分, 早期的氣動載荷靈敏度常常通過有限差分法求解。這種方法對自變量擾動的要求較高, 計算效率較低。近年來, 基于伴隨矩陣的靈敏度求解思路備受CFD研究人員的青睞[24-26]。但是, 這種思路涉及方程的變分, 對存儲空間要求較高。考慮到UVLM方程的氣動載荷與輸入量、 狀態量存在明確的因式分解關系。通過鏈式法則即可推導UVLM氣動載荷的解析靈敏度。國內外學者針對面元法的氣動載荷靈敏度開展過一些研究[27-30]。其中, 針對UVLM的研究主要考慮了氣動載荷有關節點形狀和位置的靈敏度[21,31], 缺少對節點速度和環量強度的考慮, 對UVLM線性化方程中的Jacobi矩陣設計不夠完整。

針對這一情況, 本文將首先構造線性化的UVLM狀態空間方程, 明確Jacobi矩陣包含的元素, 并著重推導氣動載荷有關節點位置、 節點速度和環量強度的解析靈敏度。在此基礎上, 本文通過數值仿真探討了這些靈敏度矩陣的分布規律, 利用有限差分法驗證了解析靈敏度的計算效率和精度。

1 計算方法

1.1 UVLM氣動載荷模型

如圖2所示, UVLM將機翼的平均氣動面離散為一系列結構網格(黑色實線), 在結構網格下游1/4弦長處布置附著渦渦格(紅色虛線), 通過對附著渦的更新模擬氣動面曲面效應。隨著計算過程的推進, 附著渦會從后緣脫落形成尾渦(黑色虛線), 集合附著渦與尾渦的位置即可較準確地反映氣動面及流場的幾何特征。UVLM在每個渦格內設置渦環, 通過環量強度完成氣動載荷計算。

圖2 UVLM渦格劃分示意圖Fig. 2 Discrete vortex panels of UVLM

附著渦渦格的幾何中心可視為氣動載荷的控制點。在計算過程中, 需要在附著渦渦格的控制點處施加不可穿透邊界條件, 用以獲得附著渦的環量強度。與平面氣動模型不同, UVLM附著渦沿曲面分布, 邊界條件中的誘導速度并非豎直向下, 而是沿著渦格的法向量向下。以t+1時刻為例, 不可穿透邊界條件可表示為

(1)

其中,Γbv和Γwk為兩組列向量, 分別表示附著渦和尾渦的環量強度; 矩陣Wbb和Wbw分別表示附著渦和尾渦關于氣動載荷作用點的影響系數, 用以計算誘導速度,Wbb和Wbw中的元素由渦格節點位置決定, 可通過Biot-Savart定律獲得;Unc表示附著渦控制點處的非環量速度, 包括剛體運動速度、 彈性振動速度、 空氣擾動速度等。上標n表示渦格的法向量方向。

尾流區域的演化包括兩部分, 即附著渦的脫落和已有尾渦的擴散。演化方程如下

(2)

其中,Abv和Awk是尾渦渦強演化矩陣;Bbv和Bwk是尾渦節點演化矩陣, 均由1和0組成;ζbv和ζwk分別表示附著渦和尾渦節點的位置,vwk表示尾渦節點的局部速度, 由流場的誘導速度及外部擾動組成。為避免尾渦數量過多帶來的額外計算成本, 可在尾流區域達到一定長度時設置截斷, 忽略遠處尾渦的影響。

在每個控制點處, UVLM的氣動載荷升力分量L垂直于非環量速度向上, 誘導阻力D沿非環量速度方向。L和D的矩陣表達式如下

(3)

其中,E1(2)是由0, 1, -1組成的稀疏矩陣, 用以計算有效渦強。G1(2)是對角矩陣, 用以計算投影面積, 由附著渦的幾何信息和局部攻角決定。U1(2)(3)是誘導速度對角矩陣。根據本文需要,G1(2)及U1(2)(3)的表達形式如下

(4)

1.2 線性化狀態空間方程

將方程(1)～(4)整合, 可得到UVLM氣動載荷在離散時刻的狀態空間形式

fA(xt+1,ut+1)=gA(xt,ut)
yt+1=hA(xt+1,ut+1)

(5)

其中,u,y,x分別表示方程的輸入量、 輸出量、 狀態量, 分別為

(6)

方程(5)可以與結構動力學方程耦合, 形成氣動彈性方程。

圖3 UVLM氣動力線性化過程的參數關系圖Fig. 3 Parameter relationship of aerodynamic linearization by the UVLM

將方程(5)中的變量表示為在平衡位置x0處的基準量與擾動量之和

(7)

再將方程(7)代入方程(1)～(3), 基于上述參數關系可得到線性化的狀態空間方程

(8)

其中, Δt為時間步長。

將方程(8)中的Jacobi矩陣整合, 線性化的狀態空間方程可表示為

EΔxt+1=AΔxt+BΔut+1, Δyt+1=CΔxt+1+DΔut+1

(9)

1.3 氣動載荷解析靈敏度

(10)

(11)

其中

(12)

(13)

(14)

其中, 下標i表示附著渦渦格節點索引, 每個元素表示節點i發生擾動對附著渦渦格k的影響。

按類似的方法可推導方程(12)～(13)中其余偏導項, 并獲得方程(11)中氣動載荷解析靈敏度的表達式, 在此不一一贅述。相比其他計算方法, 采用鏈式法則推導的解析靈敏度可直接使用矩陣運算符, 避免了數值計算中的循環過程。

2 計算結果

將長直機翼置于均勻分布的流場中, 忽略彈性變形和重力的影響, 通過氣動載荷系數驗證UVLM建模的準確性, 再對平衡狀態下的機翼施加小擾動, 驗證解析靈敏度的計算精度及效率。

2.1 氣動載荷驗證

構建5個攻角為5°, 弦長為1 m的長直機翼, 這些機翼的展弦比分別為4, 8, 12, 20, 100。將這些機翼依次置于速度50 m/s的均勻流場中, 通過UVLM模型獲得升力系數與誘導阻力系數, 并與Katz等[32]提供的參考值進行對比。計算中, 對UVLM尾流區域在9 m處進行截斷, 忽略遠處尾渦的影響。

在第1步迭代計算中, 流場從靜止開始起風, 機翼會受到一個較大的脈動載荷, 導致升力系數出現較大的初始值。隨著尾渦的生成, 這一脈動載荷帶來的影響逐漸減小。為便于比較, 本文將初始升力系數設置為0.5。對比結果如圖4所示, UVLM升力在迭代過程中逐漸達到與參考值接近的穩定狀態, 并且迭代步數隨著展弦比的增大而增加, 穩定時的升力系數也隨著展弦比的增大而接近5°攻角下的二維翼型計算值0.55。

圖4 UVLM升力系數與展弦比的關系Fig. 4 Variation of lift coefficient with aspect ratio

在誘導阻力系數的對比中, 選取展弦比為8的機翼, 分析了誘導阻力系數及分量的變化情況。圖5顯示計算結果與Katz和Poltkin的結果吻合, 進一步驗證了UVLM建模的準確性。

圖5 UVLM誘導阻力系數及分量瞬態變化(展弦比8)Fig. 5 Variation of induced drag coefficient and components at an aspect ratio of 8

2.2 氣動載荷靈敏度驗證

將機翼的平均氣動面劃分為20個展向單元和2個弦向單元。在UVLM時間推進計算中, 選取0.025 s作為UVLM時間推進步長。計算結果表明, 算例機翼的氣動載荷在起風1.5 s后達到穩定。因此, 選取2 s時刻穩定狀態作為靈敏度計算的平衡狀態。在機翼向下游方向30 m處設置截斷, 忽略30 m外的尾渦影響。尾流區域如圖6所示。計算中, 附著渦及尾渦的編號順序如圖7所示。

圖6 UVLM尾渦在2 s時刻的渦格情況Fig. 6 UVLM wake at 2 s

圖7 UVLM渦格及節點編號Fig. 7 Numbering scheme of UVLM vortex panels and nodes

(a) Analytical sensitivity

表2 升力關于渦格節點位置靈敏度結果及對比(y)

表3 升力關于渦格節點位置靈敏度結果及對比(z)

表4 升力關于渦格節點速度靈敏度結果及對比(y)

表5 升力關于渦格節點速度靈敏度結果及對比(z)

解析靈敏度與有限差分法的計算時間對比見表6, 采用解析靈敏度的計算效率顯著高于有限差分法。本文最大僅涉及40×63的靈敏度矩陣計算, 在這一情況下, 解析解靈敏度的計算效率大約是有限差分法的5～6倍。這種優勢在進行更大型的計算時會更加明顯, 這是兩種方法的理論基礎決定的。解析靈敏度在計算時直接使用矩陣運算符, 不需進行任何循環或迭代。而有限差分法作為一種數值計算方法, 必須對一個變量進行擾動, 再通過多次循環計算才能獲得靈敏度結果, 過程中的計算消耗遠大于解析解。

表6 解析解與有限差分法計算時間對比

3 結論

本文詳細描述了UVLM的狀態空間建模過程, 并推導了UVLM升力關于附著渦強度以及渦格節點運動信息的解析靈敏度。經驗證, 解析靈敏度的數值大小、 分布情況均與有限差分法一致, 各靈敏度項的計算誤差不超過0.4%。相比有限差分法, 解析靈敏度直接進行矩陣運算, 避免了擾動量選取不當帶來的收斂性問題, 也避免了大量數值計算帶來時間消耗, 體現出高精度、 高效、 高魯棒性的優勢。

氣動力解析靈敏度可進一步與結構響應方程耦合, 構建氣動彈性分析系統。除了可大幅提高大型系統的穩定性分析效率, 解析靈敏度在氣動彈性優化設計及流動控制領域都具有良好的應用價值。