詹玖榆 周興華 黃 銳
(南京航空航天大學機械結構力學及控制國家重點實驗室,南京 210016)
近年來,世界航空界正著力發展具有飛行環境(如高度、速度、氣候等) 自適應、可執行多種任務(如巡航、盤旋、機動等)的變體飛行器.這類飛行器可根據飛行任務需要自主改變結構和氣動布局,在復雜的飛行環境條件下保持良好的飛行性能[1-5].在眾多變體飛行器設計方案中,折疊翼變體飛行器因滿足多任務作戰需求和提高單任務執行效率的優勢,受到了廣泛的關注.然而,在機翼受控變體過程中,不僅非定常氣動力會隨著機翼形狀改變而發生變化,機翼的慣性、彈性和阻尼特性也同時發生變化,極易誘發不同結構模態參與的參變顫振現象,給變體飛行器的飛行安全帶來挑戰.針對不同折疊構型下變體機翼的氣動?結構耦合動力學系統的參變特性,如何建立以折疊角為參數的參數化氣動彈性模型,在全參數空間內實現氣動彈性力學行為的高效、高精度預測,是變體飛行器研制過程中亟待解決的動力學問題.
當前變體飛行器動力學建模研究可分為兩大類:一類是多剛體動力學模型(如Newton-Euler 方法[6-8]、Lagrange 方法[9-10]、Kane 方法[11-13]等);另一類是柔性體動力學模型(如浮動坐標法[14]、絕對節點坐標法[15-16]等).相較于多剛體動力學模型,柔性體動力學模型更加貼近工程實際狀況且能夠應用于折疊翼的氣動彈性分析,因此受到更廣泛的關注.例如Snyder 等[17]研究了Goland 折疊翼的顫振特性,發現折疊翼鉸鏈剛度對變體機翼的顫振頻率和顫振速度的影響較大.Selitrennik 和Karpel[18]對于快速變形的折疊翼系統提出了一種高效氣動彈性建模方法.該方法采用虛擬質量法建立系統的結構動力學方程,并通過動力學方程和CFD 技術相耦合的方式,建立了可折疊式變體機翼的氣動彈性方程.數值結果表明,瞬態變形和氣動載荷突變對折疊翼的氣動彈性穩定性有著重要的影響.對于折疊翼內外鉸鏈具有雙線性剛度的非線性氣動彈性問題,Lee 等[19]提出了一種可高效預測折疊翼亞臨界極限環振蕩的數值計算方法.
上述氣動彈性動力學建模方法雖已經成功應用到可折疊式變體機翼的氣動彈性建模中,但在氣動彈性分析過程中仍需要針對不同折疊角進行重復地結構動力學建模、非定常氣動力計算以及流固耦合建模,計算效率低下且難以分析可折疊式變體機翼在完整參數空間內的氣動彈性力學行為.參數化氣動彈性建模方法為高效、高精度分析可折疊式變體機翼的氣動彈性力學行為提供了新的解決方案.Zhao等[20]和倪迎鴿等[21]采用子結構綜合和偶極子網格法構建了折疊翼的參數化氣動彈性模型,分析了模態阻尼、折疊角、鉸鏈剛度對折疊翼動力學特性的影響.Huang 等[22]提出了一種參數化的氣動伺服彈性建模方法,高效地獲得了不同折疊角下的機翼氣動伺服彈性數學模型,實現了顫振主動抑制的控制律設計和閉環氣動伺服彈性分析.
然而,上述參數化建模方法均基于子結構綜合方法,即根據折疊過程中內外翼結構動力學模型在全局坐標系中的約束關系,采用坐標變換實現全局坐標系下折疊翼結構動力學的快速生成.該參數化建模方法并未從根本上解決不同折疊構型下模態坐標不一致的難題.本文基于流形切空間插值方法,建立以折疊角為參數的折疊翼參變結構動力學模型,并耦合基于偶極子網格法的參數化非定常氣動力模型,進而建立折疊翼的參變氣動彈性模型.為了驗證該參數化模型的準確性,本文選取一小展弦比折疊翼模型為研究對象,通過對結構自由振動分析和顫振邊界的預測,考核本文提出的參數化建模方法.
不同于常規機翼的氣動彈性分析,折疊翼的氣動彈性分析需要在每一個角度重新建立結構動力學模型和氣動模型,并將之耦合得到氣動彈性模型.為了高效、準確地分析折疊翼的氣動彈性力學行為,本節介紹基于流形切空間插值的參數化氣動彈性建模方法.模型建立過程分為參數化結構動力學建模、參數化氣動力建模和參數化氣動彈性建模3 部分.
折疊翼的幾何構型如圖1 所示,折疊翼由機身和兩個能獨立旋轉的內外翼組成.折疊翼折疊過程如圖2 所示(下標α,β 和γ 分別表示機身結構、內翼結構和外翼結構).內翼可以進行0?~120?的旋轉,同時外翼始終保持水平.機身和內外翼的材料均為鋁板,其楊氏模量為71 GPa,泊松比為0.33,密度為2.7×103kg/m3,厚度分別為2.0 mm,1.0 mm 和1.0 mm.

圖1 折疊翼幾何尺寸及構型Fig.1 Configuration of folding wing

圖2 折疊翼折疊過程Fig.2 Morphing process of folding wing
通過有限元方法得到折疊翼機身、內翼和外翼的質量、剛度矩陣,然后通過子結構綜合法[23-24]得到折疊翼結構動力學方程如下

式中,θ 是折疊角,M(θ)是折疊翼的質量矩陣,K(θ)是折疊翼的剛度矩陣,x(t) ∈Rn為物理坐標,值得注意的是,本文沒有考慮折疊翼的結構阻尼.因為實際飛行器的阻尼特性難以獲得,所以建立參數化的阻尼模型更加困難.由結構動力學方程可以得到折疊翼的振型矩陣Φ(θ)=[φ1(θ) φ2(θ) ··· φn(θ)],取Φ(θ)前m列得到降階矩陣φ(θ)=[φ1(θ) φ2(θ) ··· φm(θ)],令x(t)=φ(θ)ξ(θ,t),ξ(θ,t)是廣義坐標,同時對式(1)左乘φ(θ)T,得到折疊翼的低維結構動力學方程

折疊過程中質量矩陣、剛度矩陣、振型矩陣和氣動力矩陣都發生了顯著變化,如果要對折疊翼進行氣動彈性分析,就需要對每個折疊角進行重復有限元建模和氣動力建模,建模效率低下.參數化建模即給出相應參數(如折疊角θ)能快速得到結構動力學模型和氣動力模型,這既能提高分析折疊翼氣動彈性行為的效率,也為后續折疊翼的顫振主動控制研究打下良好的基礎.
矩陣流形是指具有特殊性質(正定性、對稱性、正交性或非奇異性)的矩陣組成的流形.在流形上插值得到的矩陣仍具有原矩陣的特殊性質.正因為這一特點,有學者提出基于流形矩陣的插值方法[25-32].令表示由一系列折疊角度所組成的集合,矩陣Pi表示折疊角θi對應的系統矩陣(如質量矩陣、剛度矩陣和振型矩陣),Pi張成的空間Si可視為矩陣流形M 上的一點.參數化建模即對于一個新參數θN,獲得其對應的矩陣PN.顯然,矩陣流形插值是一個可行的方法,但是通常矩陣流形不是“線性”的,因此需要將矩陣流形上的點映射到“平直”的切空間,在切空間中進行插值,然后再從切空間映射回流形上,其推導過程如下.


圖3 給出了流形切空間插值方法的示意圖.由于在參數化結構動力學建模中涉及兩種不同流形,上述的指數映射和對數映射使用抽象的數學表達形式,下面將具體給出這兩種不同流形的插值方法.

圖3 流形切空間插值示意圖Fig.3 Sketch map of manifold tangent space interpolation
在折疊翼低維結構動力學方程中,Mr(θ)和Kr(θ)均為SPD(symmetric positive definite)矩陣,張成的空間為SPD 流形.SPD 流形的對數映射為

式中,logm 和expm 分別指矩陣對數和矩陣指數.
在氣動彈性計算中,還有必要對振型矩陣φ(θ)進行參數化建模.但是φ(θ)是n×m維矩陣,不再適用于SPD 流形切空間插值.同時φ(θ)的列向量表示折疊翼結構的某階模態,因此參數化振型矩陣的每一列也應該表示折疊翼的某階模態.為得到參數化振型矩陣,本文引入約束矩陣?(θ),其滿足約束條件
本文基于偶極子網格法(DLM) 建立折疊翼的參數化非定常氣動力模型.偶極子網格法通過求解無量綱法洗速度與壓差系數的函數來確定空氣動力影響系數,從而獲得氣動網格上的氣動力,計算方程如下


Fa是作用在結構網格點上的等效力向量,于是我們得到了由氣動力引起的作用在個有限元節點上的參數化的空氣動力等效力,如下

模態空間下,折疊翼參數化氣動彈性方程如下



通過參數化的狀態空間方程,可求解不同風速、不同折疊角下的氣動彈性系統的響應,也可以進行氣動彈性穩定性分析.
本節以小展弦比折疊翼模型為研究對象,分別從結構動力學和氣動彈性穩定性兩個方面驗證本文所提出的參數化建模方法的準確性.結構動力學模擬方面,重點關注參數化模型預測固有頻率、固有振型等結構特性隨折疊角的變化規律.氣動彈性穩定性分析方面,預測折疊翼的顫振臨界速度和顫振頻率隨折疊角變化的規律.
參數化結構動力學模型的準確性是變體飛行器參數化氣動彈性建模的前提.本文首先建立了0?,10?,25?,50?,80?,100?和120?七個折疊角下的折疊翼有限元模型,通過流形切空間插值得到基于以上角度的參數化結構動力學模型.
如圖4 所示,參數化建模方法計算得到的前8 階模態頻率和直接法計算得到的前8 階模態頻率吻合較好.相較于其他折疊角度,基于流形切空間插值的參數化建模方法在預測115?折疊角下的第8 階固有頻率(預測頻率為101.27 Hz)與直接法(計算頻率為102.86 Hz)存在一定差異,誤差約為1.5%,仍處于工程可接受范圍.

圖4 參數化建模方法與直接法計算頻率對比Fig.4 Comparison of natural frequency between parametric model and direct method
除了結構固有頻率之外,折疊翼的結構固有振型對廣義氣動力的計算以及氣動彈性分析的精度也有重要影響.圖5 給出了5?折疊角下參數化建模方法得到的固有振型與直接法得到的振型(藍色為參數化建模計算的振型,紅色為直接法計算的振型)的對比.對比結果表明,參數化建模方法所計算的前8階固有振型與直接法的計算結果吻合較好.

圖5 5?折疊角時,各階模態對比Fig.5 Comparison of modes at 5?folding angle

圖5 5?折疊角時,各階模態對比(續)Fig.5 Comparison of modes at 5?folding angle(continued)
MAC(modal assurance criterion)值是檢驗振型函數預測誤差的重要指標.為了驗證參數化建模方法所預測的固有振型與真實振型的偏差,本文分別給出了在5?,30?,65?和105?折疊角下折疊翼前8 階固有振型的MAC 值分布.如圖6 所示,除第6 階彈性模態外,其余各階模態MAC 值均大于0.99.第6 階模態出現誤差的原因是該插值方法在兩個較小的參數值之間可能形成拐點,使得插值結果與真實值產生較大差異.

圖6 4 個典型角度下振型的MAC 值Fig.6 MAC values of mode shapes under four specific folding angles

圖6 4 個典型角度下振型的MAC 值(續)Fig.6 MAC values of mode shapes under four specific folding angles(continued)
為了進一步驗證該參數化氣動彈性模型的準確性,本文采用參數化建模方法計算了折疊翼的顫振邊界.如圖7 所示,參數化氣動彈性模型成功預測了折疊翼顫振邊界隨折疊角的變化規律:隨著折疊角變化,折疊翼的顫振模態發生了復雜的模態切換現象.在0?~30?折疊角區間,折疊翼的第2 階模態發生顫振.當折疊角位于35?~80?區間時,結構的第3 階模態發生顫振.在85?~105?折疊角區間,折疊翼發生顫振的是第4 階模態.當機翼折疊角進一步增大到110?~120?范圍時,折疊翼的第3 階模態再次發生顫振.
此外,本文將參數化氣動彈性模型對顫振邊界的預測結果與直接法預測的結果進行了對比研究.如圖7 所示,在大部分折疊角下,參數化折疊翼模型所預測的顫振臨界速度和顫振頻率均與直接計算結果吻合較好,頻率誤差最大值小于1%,顫振臨界速度誤差最大值低于4%.然而,圖示結果表明,在105?的折疊角處,參數化建模方法所預測的顫振臨界速度和顫振頻率與直接法計算結果存在較大差異.從顫振頻率隨折疊角變化的曲線(圖7(a)) 來看,發生較大差異的主要原因是參數化建模方法預測的顫振模態分支與直接法計算得到的顫振模態分支不同.為了深入探究引起顫振模態分支變化的原因,圖8 分別給出了105?折疊角下直接法和參數化建模方法計算得到的氣動彈性系統根軌跡分布圖.如圖所示,該折疊角下兩種計算方法所預測的根軌跡分布圖相似,主要差異在于直接法預測了兩種顫振形態.第一類顫振形態是第4 階模態發生顫振,顫振臨界速度是54.6 m/s;第二類顫振形態是第3 階模態發生小阻尼顫振,顫振臨界速度為32.3 m/s.本文提出的參數化建模方法預測的第4 階模態發生顫振的臨界速度是54.2 m/s,與直接法結果吻合較好.

圖7 參數化建模方法與直接法氣動彈性計算對比Fig.7 Comparison between parametric modeling and direct method in aeroelasticity

圖8 參數化模型與直接法計算根軌跡圖對比Fig.8 Comparison of root locus calculated by parametric modeling and direct method
和直接法相比,本文提出的基于流形切空間插值的參數化建模方法無須重復進行固有模態分析,計算效率遠高于直接法.表1 展示了使用兩種方法分別計算圖7 所示的顫振邊界所需的時間成本對比.直接法是指通過MSC.Nastran 對給定折疊角度下的折疊翼進行有限元模型的設計,并使用商業軟件對氣動和結構網格進行插值、進行非定常氣動力的計算以獲得此構型下的顫振速度;而參數化建模的方法不需要手動地重復有限元建模、結構和氣動插值,僅僅改變參數即可得到此時結構的剛度、質量矩陣,單次時間為0.7 s,可忽略不計,因此其總時間花費僅為傳統方法54.3%.

表1 兩種方法計算時間對比Table 1 Comparison of the time cost for the two methods
本文提出了一種基于流形切空間插值的折疊翼參數化氣動彈性建模方法,實現了折疊翼的氣動彈性模型的高效建立.相較于現有的參數化建模方法,本文提出的參變建模方法有效解決了變體過程中存在的模態坐標不一致性.為了驗證該參數化模型在預測折疊翼氣動彈性力學行為的準確性,本文以折疊翼為研究對象,分別從折疊翼結構固有頻率、模態振型、顫振邊界等方面進行了算例驗證.結果表明,該參數化建模方法可高效、高精度預測可折疊式變體機翼的固有頻率和顫振邊界隨折疊角的演變規律.在后續工作中,可將該方法擴展到解決快速變體過程的時變動力學建模,重點圍繞模態坐標一致性問題和時域非定常氣動力計算等方面開展研究.