基于CFD的地面顫振模擬試驗非定常氣動力重構方法研究

宋巧治， 王彬文， 李曉東

(中國飛機強度研究所 飛機結構振動技術研究室，西安 710065)

飛行器顫振是由結構與空氣相互耦合而發生的一種氣動彈性失穩現象，顫振問題對飛行器的飛行安全構成很大威脅。雖然顫振主動抑制算法研究較多，但是距離工程實用尚有一定的距離，因此目前防止顫振發生的有效途徑是通過分析和試驗獲得準確的顫振邊界并嚴格控制飛行速度在顫振包線以內。目前的顫振的驗證手段主要有風洞試驗和飛行試驗[1]，兩者都具有周期長、成本高、風險高的特點，因此國外學者提出了一種地面顫振模擬試驗驗證方法，通過在地面重構結構的非定常氣動力實現顫振邊界的測試。

地面顫振模擬的概念最早由Kearns[2]提出，1974年法國的Rajagopal研究了二自由度舵面模型的地面顫振模擬試驗，但是由于缺乏先進的數字計算機技術，因此只能成為一個設想。隨著硬件條件的發展，地面顫振模擬試驗得以實現，Karkle等[3]研究了地面顫振模擬試驗方案，并進行了舵面顫振和導彈氣動伺服彈性特性測試； Zeng等[4]系統地研究了氣動力重構和激勵力控制方法，基于平板結構對地面顫振模擬試驗進行了驗證。Liseykin等[5]在非線性顫振驗證進行了初步應用嘗試。Dhital等[6-7]研究了氣動力等效降階的方法。

國內潘樹祥等[8]在20世紀80年代就開始了地面熱顫振模擬試驗技術研究，通過石英燈加熱的方式模擬氣動加熱環境，對熱環境下結構的顫振特性進行了測試，為熱顫振驗證提供了一種可行手段。Wu等[9-11]研究了氣動力重構以及氣動力模擬加載控制等問題，并利用導彈模型等結構的對方法進行了驗證。胡巍等[12]實現了帶有舵面結構的氣動力重構方法。邵崇暉[13]研究了地面顫振模擬試驗在壁板結構的顫振測試中的應用。

1 非定常氣動力重構

影響地面顫振模擬試驗的因素主要為重構氣動力的精度，當前國內外大多數研究都是基于頻域氣動力模型，氣動力計算的精度在跨音速范圍內較差，因而限制了地面顫振模擬試驗應用馬赫數的范圍。基于計算流體動力學(computational fluid dynamics,CFD)的時域方法具有較高的計算精度，但是計算效率較低。由于實際飛行過程中的非定常氣動力作用是實時的，因此要求重構的氣動力模型具有實時性，為了將CFD方法用于地面顫振模擬試驗中，需要研究模型降階方法，提高氣動力的計算效率。

跨音速飛行中由于激波的出現，氣動力呈現較強的非線性特性。Dowell[14]經過研究后發現在跨音速范圍內，雖然氣動力整體是非線性的，但是非線性主要體現在靜壓背景，而非定常氣動力在結構作小幅值振蕩時，仍可以認為是線性的，仿真分析也證明了Dowell假設的準確性，該研究為后續降階模型(reduced order model，ROM)奠定了基礎。

近年來，模型降階方法開始應用于非定常氣動力計算中。鑒于氣動力降階能夠大幅提升計算效率，因此國內外學者對降階方法開展了大量研究，并形成了眾多的成果，當前方法按照原理主要分為三類：基于系統辨識的方法、基于流場模態的方法以及正交投影分解(proper orthogonal decomposition，POD)方法[15]。

基于系統辨識的方法由于原理簡單，是其中較為常用的一種方法，該方法在降階過程中需要對分析模型進行激勵，以便獲得較寬頻帶范圍的系統輸入和輸出信號，然后利用輸入和輸出建立氣動力降階模型。Gupta等[16-17]提出了以顫振飛行試驗中常用的3211信號作為激勵形式的非定常氣動力降階方法，后期經過發展逐漸成熟。張偉偉等[18]基于3211信號對氣動力降階模型進行了研究，認為該信號激勵頻帶較寬，但是需要反復調節信號的周期以充分激勵期望的頻帶。

本文在氣動力模型降階的基礎上，提出坐標轉換方法，建立物理坐標系下的氣動力模型，同時對插值點進行縮聚，滿足地面顫振模擬試驗插值點的數量要求，最終建立了地面顫振模擬試驗可用的非定常氣動力模型。

2 氣動力模型重構原理

常用的氣動力降階模型大都基于廣義坐標建立，而地面顫振模擬試驗中需要的是物理坐標下的氣動力模型，為了完成氣動力模型由廣義坐標向物理坐標的轉換，推導了兩個坐標系的相互轉換關系。同時基于廣義力等效的方式，完成了物理坐標系下的氣動力縮聚模型的建立，具體流程如圖1所示。

圖1 非定常氣動力重構方案Fig.1 Scheme of unsteady aerodynamic reproduction

3 氣動力模型降階

以跨音速標準模型AGARD445.6機翼作為驗證對象，該機翼采用NACA65A004對稱翼型，具有明顯的跨音速氣動特征，半展長762 mm，1/4弦處后掠角45，展弦比(展長與平均弦長)1.652 5，跟梢比0.66，模型尺寸如圖2所示。

圖2 AGARD445.6機翼模型尺寸Fig.2 Model size of AGARD445.6 wing

為了對模型進行分析，首先建立了結構的有限元模型，通過調節材料特性參數，使分析模態頻率與試驗結果吻合，表1為最終分析結果與試驗結果對比，可以看出模型與試驗結果前四階模態頻率誤差都在1%以內，證明分析模型能夠反映真實結構的動態特性。

表1 分析結構模態頻率與試驗結果對比Tab.1 Comparison of model frequency between analysis and test

為了獲得結構的顫振邊界，基于CFD/CSD(computational fluid dynamic/computational solid dynamic)耦合的方法進行分析，求解采用無黏Euler方程、理想氣體模型，空間離散采用二階迎風格式，物理時間步長設置為0.002 s，通過給定大氣來流條件，分析不同動壓條件下的結構響應趨勢，得到了結構的顫振邊界，計算來流條件如表2所示，結構和氣動網格如圖3和圖4所示。

表2 顫振計算大氣狀態條件Tab.2 Airflow condition of flutter analysis

圖3 結構網格劃分Fig 3 Mesh generation on structure

圖4 結構表面及對稱面氣動網格Fig.4 Fluid mesh on wing surface and symmetry

設計了用于模態激勵的3211信號，以結構前四階模態振型作為強迫位移分布，各階模態的幅值均設置為0.001 5，四階模態的廣義位移幅值見圖5。通過流固耦合分析，獲得結構節點的非定常氣動力向量，并利用坐標變換計算獲得結構的廣義氣動力。

圖5 前四階模態強迫位移幅值Fig.5 Forced displacement magnitude of first four modes

以廣義位移為輸入，廣義氣動力為輸出，采用自回歸滑動平均ARMA模型對氣動力進行系統建模，該模型離散時間形式為[19]

(1)

式中：y為CFD計算輸出廣義氣動力向量；u為CFD系統輸入廣義位移信號；na，nb分別為輸出、輸入延遲階數，本例中均取6；k為離散時間變量；A，B為系統模型待辨識參數。通過系統辨識建立了氣動力降階模型，在Ma=0.96的狀態下，前四階模態的廣義力與降階模型輸出結果對比情況，如圖6所示。由圖6可知，降階模型與原始模型的吻合度較高，同時也可以看出給定的3211激勵信號激勵下，模型的各階廣義力響應都比較明顯，因此可以證明設定的激勵信號滿足模態激勵的要求。

圖6 結構前四階模態廣義氣動力Fig.6 The generalized force of first four modes

經過以上降階過程，建立了廣義坐標下的氣動力影響系數矩陣Qhh，但試驗過程中使用的物理量必須是在物理坐標下描述的，為此需要完成廣義坐標下物理量與物理坐標下物理量的相互轉換。

4 坐標變換

由結構模態分析理論可知，結構響應和激勵力的廣義坐標和物理坐標定義為

x=Φq

(2)

F=ΦTf

(3)

式中：x為物理坐標下結構位移響應；q為廣義坐標，Φ為模態矩陣；F為廣義力；f為物理坐標下力向量。

由于模態矩陣本身不具備正交性，模態變換方程不能直接通過前乘振型矩陣完成，同時考慮到結構有限元節點往往數量眾多，求解獲得所有階次模態矩陣較為困難，因此無法對模態矩陣求逆以完成上述轉換，為了解決上述轉換問題，利用了模態矩陣關于質量陣正交(假設模態矩陣按照質量陣進行歸一化，見式(4))這一特性完成坐標變換。

ΦTMΦ=I

(4)

式中，M為結構的質量矩陣。

ΦTMx=ΦTMΦq

(5)

整理可得

q=ΦTMx

(6)

式(6)完成了響應由物理坐標向模態廣義坐標的轉換。

力的坐標變換不能直接通過前乘矩陣的方式完成變換，為此采用如下假設，即如果兩種物理坐標下力向量分布轉換為廣義力等效，則結構的模態響應也等效，同時考慮到對顫振影響的主要模態為低階模態，因此低階模態廣義力等效即可保證顫振邊界的等效，假設節點力向量服從以下分布

f=MΦf0

(7)

式中，f0為假設的任意向量，則相應的廣義力為

F=ΦTf=ΦTMΦf0=f0

(8)

代入式(7)可以得到

f=MΦf0=MΦF

(9)

由于MΦ為滿秩矩陣，MΦf0可以生成所有向量，上述假設不失一般性。因此由式(9)可以實現模態廣義坐標下的力到物理坐標下的節點力向量轉換，物理坐標下的氣動力影響系數矩陣可以描述為

Qkk=MΦQhhΦTM

(10)

通過式(10)操作建立了物理坐標下降階氣動力模型，該模型以結構節點響應為輸入，以結構節點氣動載荷為輸出。由于模型的輸入是所有節點的響應，輸出是全部節點的氣動載荷，而在地面顫振模擬試驗中，受設備數量和空間位置的限制，能夠記錄的結構節點響應和模擬的氣動載荷加載點數量有限。為了能夠滿足試驗要求，需要對結構插值節點進行縮聚處理。

羅扎諾夫的寫作方法簡直快趕上當下的某些“身體寫作”了。 我們似乎不應該誤解其寫作的“身體性”，而是應該把“身體性”理解為極度的“任性”。 他說：

5 插值節點縮聚

顫振控制方程可以通過式(11)描述

(11)

式中：M，C，K分別為結構的質量陣、阻尼陣和剛度陣；x為節點響應列向量；f為氣動力向量，可以通過式(12)計算獲得

f=qdMΦQhhΦTMx

(12)

式中：qd為來流動壓；x為所有節點的響應；f為所有節點的氣動力。為了使結構的顫振特性不發生改變，節點縮聚需要滿足條件：①根據縮聚后節點能夠獲得與全部節點等效的變形；②縮聚后節點力能夠等效所有節點的節點力。

為了滿足上述兩個條件，引入了插值方法，由于本文采用的結構模型為三維模型，因此插值采用薄板插值(thin panel spline, TPS)方法，薄板插值方法變形控制方程為

(13)

無限遠處邊界條件定義為

(14)

利用插值點構造滿足邊界條件的變形場，由此可以得到插值矩陣

h=GTPSx

(15)

式中，GTPS為插值矩陣，若要第一個條件滿足，則應有

‖x-Grxr‖<ε

(16)

式中：x為所有節點響應向量；xr為縮聚后節點響應向量；Gr為r集到s集節點的插值矩陣；ε為任意小量。引入模態變換，則有

Δ=‖Φs-GrΦr‖<ε

(17)

式中：Φs為模型所有節點模態振型矩陣；GrΦr為縮聚后節點振型矩陣在所有節點上的插值。

僅考慮低階模態，如果縮聚后節點振型矩陣在所有節點上插值與所有節點模態振型矩陣等效，則條件1滿足；

在條件1滿足情況下，根據虛功原理，則有式(18)成立，條件2滿足

(18)

為了使條件1滿足，采用了優化算法對縮聚節點的位置進行了優化，采用4激振點和4個拾振點的配置方案，優化方法采用遺傳算法，采用縮聚后節點前兩階模態振型在所有節點上的插值振型與原始振型的差異最小作為優化的目標函數，由于結構的振動響應主要沿Z方向，因此本文中計算振型只考慮Z方向的分量，其他方向不做考慮，優化變量取縮聚后節點位置，由于實際試驗中邊界處的約束條件自動滿足，因此在縮聚過程中，增加了邊界處的虛擬插值點以對應實際的邊界條件。

由于結構顫振發生時，顫振的主參與模態集中在低階模態，因此采用模態疊加原理，若要使式(17)成立，則應該有式(19)成立

Δ1=‖Φs1-GrΦr1‖<ε

(19)

式中：Φs1為關心的結構模態的振型；GrΦr1為選擇的測點處的振型分量在整個結構節點上的插值，通過優化算法進行求解，使兩者的偏差Δ1最小，即可利用少量測點代替所有結構節點的目標。

優化選取的插值點位置，如圖7所示。為了對插值點縮聚效果進行說明，計算了縮聚后的插值振型并與原始振型進行了對比，如圖8所示。采用縮聚后插值振型與原始振型的MAC值作為插值的等效的衡量指標，基于優化插值點位置獲得的前兩階模態振型與原始振型的MAC值分別為0.98和0.94，證明縮聚模型具有較高的精度。

圖7 優化獲得的激振點(拾振點)位置Fig.7 Shaker location after optimization

經過縮聚可以實現利用少量激振點和拾振點即可滿足氣動力等效模擬的要求，最終的氣動力計算模型如式(20)所示，xr為縮聚后節點響應向量，freduced為縮聚點上的等效氣動力，該氣動力模型即可在地面顫振模擬試驗中使用，氣動力模型的建立整個過程如圖9所示。

(20)

圖8 基于縮聚點插值后振型與原始振型對比Fig.8 Comparison between modal shape interpolated from condensed nodes and original shape

圖9 物理坐標下氣動力建模流程Fig.9 Aerodynamic modeling process in physical coordinates

表3 坐標變換及插值點縮聚精度分析Tab.3 Analysis of accuracy for coordinate transformation and condensation of interpolation nodes

6 結構建模與仿真

為了對本文提出的氣動力模型重構方法精度進行驗證，利用直接耦合法計算結果作為標準，采用仿真的方式進行對比驗證，以表1中Ma=0.96作為驗證狀態點。采用有限元計算獲得的質量陣、剛度陣以及模態振型矩陣，建立了結構狀態空間模型，并與氣動力模型構成閉環系統模擬結構的顫振特性，通過在不同來流條件下進行測試，獲得了不同測試速度條件下的系統響應情況，通過調節來流條件，使系統的響應幅值達到穩定，從而獲得系統的臨界穩定點，對應的測試條件即為結構的顫振邊界。圖10給出了亞臨界(來流速度326 m/s)、臨界(來流速度330 m/s)和超臨界(來流速度340 m/s)，系統仿真響應輸出情況。

圖10 不同風速條件下系統響應(亞臨界、臨界、超臨界)Fig.10 System response under different test condition(subcritical, critical and supercritical)

將仿真獲得的顫振邊界與直接耦合方法及風洞試驗結果進行了對比分析(無量綱化后的結果)，三種方式獲得的顫振速度邊界對比情況，如圖11所示。

圖11 顫振邊界仿真結果與直接耦合及試驗結果對比Fig.11 Comparison of flutter boundary simulation result with CFD/CSD coupling method and experiment

由圖11可以看出，基于重構非定常氣動力模型獲得的顫振邊界與直接耦合法獲得的顫振邊界吻合較好，證明了本文提出方法的有效性，同時為了驗證建立的模型是否滿足實時性要求，設計了半物理試驗，試驗采用NI半物理仿真設備，將建立的氣動力模型在設備上運行，設置仿真模型采樣頻率為200 Hz，試驗中未出現運行丟幀的情況(如果出現丟幀證明在一個采樣周期內模型未能完成運算)，證明了模型的運行時長小于0.005 s，從而可以證明效率滿足試驗的要求。至此已經實現氣動力的重構，并且模型在實時性和插值點數量方面均滿足地面顫振模擬試驗的要求。

7 結 論

本文依據地面顫振模擬試驗對非定常氣動力模型的要求，提出了一種基于CFD降階模型的非定常氣動力重構方法，建立了模態坐標下的非定常氣動力模型。通過提出的坐標變換方法將氣動力模型轉換至物理坐標下，采用標準模型AGARD445.6作為驗證對象，基于建立的結構模型進行了聯合仿真證實了提出方法的精度和有效性，具體結論入如下：

(1) 引入了氣動力降階模型方法，基于CFD技術建立了廣義坐標系下的氣動力影響系數矩陣，提高了氣動力計算效率。

(2) 提出了模態坐標與物理坐標量之間的轉換關系，將氣動力模型轉換至物理坐標下，避免了不完備模態振型矩陣直接求逆帶來的精度損失，提高了坐標變換的精度。

(3) 以AGARD445.6標準模型對氣動力重構方法進行了驗證，證明了提出的降階重構方法在保證氣動力精度損失較小的情況下，提高氣動力計算的效率，滿足地面顫振模擬試驗對氣動力模型計算實時性的要求。