陳金池, 黃研昕

(中國(guó)航空工業(yè)空氣動(dòng)力研究院高速高雷諾數(shù)實(shí)驗(yàn)室, 遼寧沈陽 110034)

引 言

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

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

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

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

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

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

針對(duì)這一情況, 本文將首先構(gòu)造線性化的UVLM狀態(tài)空間方程, 明確Jacobi矩陣包含的元素, 并著重推導(dǎo)氣動(dòng)載荷有關(guān)節(jié)點(diǎn)位置、 節(jié)點(diǎn)速度和環(huán)量強(qiáng)度的解析靈敏度。在此基礎(chǔ)上, 本文通過數(shù)值仿真探討了這些靈敏度矩陣的分布規(guī)律, 利用有限差分法驗(yàn)證了解析靈敏度的計(jì)算效率和精度。

1 計(jì)算方法

1.1 UVLM氣動(dòng)載荷模型

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

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

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

(1)

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

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

(2)

其中,Abv和Awk是尾渦渦強(qiáng)演化矩陣;Bbv和Bwk是尾渦節(jié)點(diǎn)演化矩陣, 均由1和0組成;ζbv和ζwk分別表示附著渦和尾渦節(jié)點(diǎn)的位置,vwk表示尾渦節(jié)點(diǎn)的局部速度, 由流場(chǎng)的誘導(dǎo)速度及外部擾動(dòng)組成。為避免尾渦數(shù)量過多帶來的額外計(jì)算成本, 可在尾流區(qū)域達(dá)到一定長(zhǎng)度時(shí)設(shè)置截?cái)? 忽略遠(yuǎn)處尾渦的影響。

在每個(gè)控制點(diǎn)處, UVLM的氣動(dòng)載荷升力分量L垂直于非環(huán)量速度向上, 誘導(dǎo)阻力D沿非環(huán)量速度方向。L和D的矩陣表達(dá)式如下

(3)

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

(4)

1.2 線性化狀態(tài)空間方程

將方程(1)～(4)整合, 可得到UVLM氣動(dòng)載荷在離散時(shí)刻的狀態(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)動(dòng)力學(xué)方程耦合, 形成氣動(dòng)彈性方程。

圖3 UVLM氣動(dòng)力線性化過程的參數(shù)關(guān)系圖Fig. 3 Parameter relationship of aerodynamic linearization by the UVLM

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

(7)

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

(8)

其中, Δt為時(shí)間步長(zhǎng)。

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

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

(9)

1.3 氣動(dòng)載荷解析靈敏度

(10)

(11)

其中

(12)

(13)

(14)

其中, 下標(biāo)i表示附著渦渦格節(jié)點(diǎn)索引, 每個(gè)元素表示節(jié)點(diǎn)i發(fā)生擾動(dòng)對(duì)附著渦渦格k的影響。

按類似的方法可推導(dǎo)方程(12)～(13)中其余偏導(dǎo)項(xiàng), 并獲得方程(11)中氣動(dòng)載荷解析靈敏度的表達(dá)式, 在此不一一贅述。相比其他計(jì)算方法, 采用鏈?zhǔn)椒▌t推導(dǎo)的解析靈敏度可直接使用矩陣運(yùn)算符, 避免了數(shù)值計(jì)算中的循環(huán)過程。

2 計(jì)算結(jié)果

將長(zhǎng)直機(jī)翼置于均勻分布的流場(chǎng)中, 忽略彈性變形和重力的影響, 通過氣動(dòng)載荷系數(shù)驗(yàn)證UVLM建模的準(zhǔn)確性, 再對(duì)平衡狀態(tài)下的機(jī)翼施加小擾動(dòng), 驗(yàn)證解析靈敏度的計(jì)算精度及效率。

2.1 氣動(dòng)載荷驗(yàn)證

構(gòu)建5個(gè)攻角為5°, 弦長(zhǎng)為1 m的長(zhǎng)直機(jī)翼, 這些機(jī)翼的展弦比分別為4, 8, 12, 20, 100。將這些機(jī)翼依次置于速度50 m/s的均勻流場(chǎng)中, 通過UVLM模型獲得升力系數(shù)與誘導(dǎo)阻力系數(shù), 并與Katz等[32]提供的參考值進(jìn)行對(duì)比。計(jì)算中, 對(duì)UVLM尾流區(qū)域在9 m處進(jìn)行截?cái)? 忽略遠(yuǎn)處尾渦的影響。

在第1步迭代計(jì)算中, 流場(chǎng)從靜止開始起風(fēng), 機(jī)翼會(huì)受到一個(gè)較大的脈動(dòng)載荷, 導(dǎo)致升力系數(shù)出現(xiàn)較大的初始值。隨著尾渦的生成, 這一脈動(dòng)載荷帶來的影響逐漸減小。為便于比較, 本文將初始升力系數(shù)設(shè)置為0.5。對(duì)比結(jié)果如圖4所示, UVLM升力在迭代過程中逐漸達(dá)到與參考值接近的穩(wěn)定狀態(tài), 并且迭代步數(shù)隨著展弦比的增大而增加, 穩(wěn)定時(shí)的升力系數(shù)也隨著展弦比的增大而接近5°攻角下的二維翼型計(jì)算值0.55。

圖4 UVLM升力系數(shù)與展弦比的關(guān)系Fig. 4 Variation of lift coefficient with aspect ratio

在誘導(dǎo)阻力系數(shù)的對(duì)比中, 選取展弦比為8的機(jī)翼, 分析了誘導(dǎo)阻力系數(shù)及分量的變化情況。圖5顯示計(jì)算結(jié)果與Katz和Poltkin的結(jié)果吻合, 進(jìn)一步驗(yà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

2.2 氣動(dòng)載荷靈敏度驗(yàn)證

將機(jī)翼的平均氣動(dòng)面劃分為20個(gè)展向單元和2個(gè)弦向單元。在UVLM時(shí)間推進(jìn)計(jì)算中, 選取0.025 s作為UVLM時(shí)間推進(jìn)步長(zhǎng)。計(jì)算結(jié)果表明, 算例機(jī)翼的氣動(dòng)載荷在起風(fēng)1.5 s后達(dá)到穩(wěn)定。因此, 選取2 s時(shí)刻穩(wěn)定狀態(tài)作為靈敏度計(jì)算的平衡狀態(tài)。在機(jī)翼向下游方向30 m處設(shè)置截?cái)? 忽略30 m外的尾渦影響。尾流區(qū)域如圖6所示。計(jì)算中, 附著渦及尾渦的編號(hào)順序如圖7所示。

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

圖7 UVLM渦格及節(jié)點(diǎn)編號(hào)Fig. 7 Numbering scheme of UVLM vortex panels and nodes

(a) Analytical sensitivity

表2 升力關(guān)于渦格節(jié)點(diǎn)位置靈敏度結(jié)果及對(duì)比(y)

表3 升力關(guān)于渦格節(jié)點(diǎn)位置靈敏度結(jié)果及對(duì)比(z)

表4 升力關(guān)于渦格節(jié)點(diǎn)速度靈敏度結(jié)果及對(duì)比(y)

表5 升力關(guān)于渦格節(jié)點(diǎn)速度靈敏度結(jié)果及對(duì)比(z)

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

表6 解析解與有限差分法計(jì)算時(shí)間對(duì)比

3 結(jié)論

本文詳細(xì)描述了UVLM的狀態(tài)空間建模過程, 并推導(dǎo)了UVLM升力關(guān)于附著渦強(qiáng)度以及渦格節(jié)點(diǎn)運(yùn)動(dòng)信息的解析靈敏度。經(jīng)驗(yàn)證, 解析靈敏度的數(shù)值大小、 分布情況均與有限差分法一致, 各靈敏度項(xiàng)的計(jì)算誤差不超過0.4%。相比有限差分法, 解析靈敏度直接進(jìn)行矩陣運(yùn)算, 避免了擾動(dòng)量選取不當(dāng)帶來的收斂性問題, 也避免了大量數(shù)值計(jì)算帶來時(shí)間消耗, 體現(xiàn)出高精度、 高效、 高魯棒性的優(yōu)勢(shì)。

氣動(dòng)力解析靈敏度可進(jìn)一步與結(jié)構(gòu)響應(yīng)方程耦合, 構(gòu)建氣動(dòng)彈性分析系統(tǒng)。除了可大幅提高大型系統(tǒng)的穩(wěn)定性分析效率, 解析靈敏度在氣動(dòng)彈性優(yōu)化設(shè)計(jì)及流動(dòng)控制領(lǐng)域都具有良好的應(yīng)用價(jià)值。