基于CFD降階的非線性氣動彈性穩(wěn)定性分析

周　強， 陳　剛， 李躍明

(西安交通大學 機械結構強度與振動國家重點實驗室，西安　710049)

?

基于CFD降階的非線性氣動彈性穩(wěn)定性分析

周強， 陳剛， 李躍明

(西安交通大學 機械結構強度與振動國家重點實驗室，西安710049)

為了快速尋找基于CFD/CSD的非線性氣彈系統(tǒng)的顫振邊界，根據Lyapunov穩(wěn)定性理論對非線性流固耦合系統(tǒng)進行穩(wěn)定性分析。首先通過微擾理論對非線性流固耦合系統(tǒng)處理建立近似線性化狀態(tài)空間方程，再通過POD(Proper Orthogonal Decomposition)方法將高維狀態(tài)空間方程降階為便于分析的降階系統(tǒng)，根據系統(tǒng)所有特征值即可判定原始非線性系統(tǒng)穩(wěn)定性。Lyapunov穩(wěn)定性理論主要針對非線性系統(tǒng)，在實現過程中采用了POD降階的方法，與大多數對降階系統(tǒng)穩(wěn)定性判定不同，其數學理論上是反映原始非線性流固耦合系統(tǒng)穩(wěn)定性。POD降階方法從系統(tǒng)內部流場出發(fā)，可以較好反映系統(tǒng)內部特性。二維三維算例仿真結果驗證了該方法預測顫振邊界的正確性，分析發(fā)現，在亞音速階段，穩(wěn)定性主要由于結構模態(tài)主導；而在跨音速和超音速流動階段，氣彈穩(wěn)定性主要由受流體特性影響。

流固耦合；李雅普洛夫；正則正交分解；狀態(tài)空間；顫振

飛行器氣動彈性是研究飛行器結構與氣動力相互耦合關系，是典型的流固耦合振動問題。傳統(tǒng)的線性化方法很難處理跨音速階段的激波、渦脫落等流體非線性影響[1]。而基于CFD/CSD耦合計算飛行器氣動彈性問題，可以很好的處理絕大多數流體非線性問題并得到廣泛應用。但是CFD/CSD耦合計算氣動彈性問題，由于其巨大的計算量，即使氣動彈性中線性小擾動問題——顫振，分析一種狀態(tài)下的氣彈響需數小時至幾天。近年來，基于CFD氣動彈性降階模型(ROM)被認為解決這些關鍵問題的主要途徑之一[2]。即利用ROM來代替原始非線性CFD/CSD系統(tǒng)，再通過對ROM穩(wěn)定性分析來判斷非線性氣彈系統(tǒng)穩(wěn)定性。

目前應用于氣動彈性分析的ROM主要有兩大類：一類是系統(tǒng)辨識方法，利用系統(tǒng)的輸入和輸出之間關系來建立一個低階狀態(tài)空間模型代替原始氣動彈性系統(tǒng)，有Volterra 級數法[3-5]，ARMA模型[6]以及神經網絡[7]等系統(tǒng)識別方法，其缺點相當于把原始系統(tǒng)當成一個黑箱子，所以此類方法可以得到與原系統(tǒng)一樣的輸入輸出關系，但是并不能反映系統(tǒng)的本質特性。另一類是基于流場特征模態(tài)的方法[8]，從流固耦合系統(tǒng)方程出發(fā)，利用當地CFD流場數據，通過微擾理論建立耦合系統(tǒng)線性化空間方程，通過尋找的特征模態(tài)來描述流場變化特性，將全階方程投影到特征模態(tài)空間中得到降階系統(tǒng)。POD降階是目前應用比較廣泛的尋找最優(yōu)特征模態(tài)方法之一，與大多數系統(tǒng)識別方法不一樣，POD方法需建立在流體系統(tǒng)內部流場上，所以能保留CFD流場中激波等非線性特性，更能反映流體的內部結構特征。自從POD方法首次成功被用于氣動彈性問題中后[9]，POD降階不斷發(fā)展推廣應用到簡單翼型[10]、機翼、渦輪葉片[11]以及全機模型[12]中，所解決的問題也從簡單響應分析擴展到參數變化的自適應降階，控制[13]等問題中。目前廣泛采用根軌跡法對ROM進行顫振邊界的尋找[14]，但是該方法反映的實際上是ROM系統(tǒng)的穩(wěn)定性，并不能在數學上證明原始非線性系統(tǒng)的穩(wěn)定性，只能在數值響應上驗證所尋找顫振點的正確性。

對于非線性系統(tǒng)，本文所采用的Lyapunov穩(wěn)定性判據主要針對非線性氣彈系統(tǒng)穩(wěn)定性分析，基本思路為將非線性系統(tǒng)運動方程在足夠小的鄰域內進行泰勒展開，導出一次近似線性化系統(tǒng)，再根據線性化系統(tǒng)特征值在復平面的分布判斷非線性鄰域內的穩(wěn)定性。若線性化系統(tǒng)特征值均有負實部，則原非線性系統(tǒng)在鄰域內穩(wěn)定；若線性化系統(tǒng)包含正實部特征值，則原非線性系統(tǒng)在領域內不穩(wěn)定。該方法原理上對非線性系統(tǒng)在鄰域內泰勒展開，剛好與目前廣泛使用的POD降階理論在數學本質上一致的。而目前所見文獻中POD降階大多關于效率與精度的分析[15-16]，采用響應分析判斷系統(tǒng)穩(wěn)定性，快速尋找顫振邊界的方法還不多見。本文根據微擾動方法和POD方法建立了降階的氣動彈性系統(tǒng)狀態(tài)空間方程，分析系統(tǒng)隨著動壓變化的穩(wěn)定性特性變化。由于本文所采用的方法是針對非線性系統(tǒng)，建立的降階系統(tǒng)由非線性流固耦合系統(tǒng)方程導出的，所以更能反映系統(tǒng)內部的非線性特性。

1　理論基礎

對于氣動彈性穩(wěn)定性分析，一種常見的方法是對所建立的ROM模型采用根軌跡法，即變化動壓或者VF的數值，觀察結構自由度對應的特征值變化，繪制根軌跡變化曲線。若實軸值穿過零點，則為顫振臨界點。而本文采用Lyapunov第一穩(wěn)定性理論對非線性系統(tǒng)穩(wěn)定性分析，所以本文首先將非線性流固耦合方程在足夠小鄰域的泰勒展開得到的，再通過POD降階方法得到低維度的氣動彈性系統(tǒng)，觀察降階系統(tǒng)的所有特征值正負即可判斷穩(wěn)定性。具體流程首先將基于CFD的Euler求解器線性化(1.2節(jié))，并通過POD降階方法將高維線性化狀態(tài)空間方程降階(1.3節(jié))，最后統(tǒng)一流體和結構方程的時間量綱，并耦合構成顫振分析系統(tǒng)的低維狀態(tài)空間模型(1.4節(jié))，通過不同方法分析其穩(wěn)定性。

1.1非線性氣動彈性控制方程

基于CFD/CSD耦合求解氣動彈性方法是分別求解流體、結構動力學方程，再通過交界面的數據插值，如此反復迭代，觀察響應曲線趨勢判斷發(fā)散或收斂來分析其穩(wěn)定性?；贑FD的氣動彈性系統(tǒng)控制方程如下：

(1)

1.2基于CFD的線性化方程

根據Lyapunov第一穩(wěn)定性理論，需對非線性系統(tǒng)在鄰域內泰勒展開，所以本文對非線性流體動力學方程在當地CFD流場下對流場變量和結構位移泰勒展開[17]。對于氣動彈性系統(tǒng)，選取結構未變形的定常狀態(tài)流場作為初始流場，此時

(2)

(3)

整理成狀態(tài)空間形式為

(4)

式(4)為流體系統(tǒng)狀態(tài)空間模型，為多輸入多輸出系統(tǒng)，系統(tǒng)輸入為模態(tài)位移和速度xs，而系統(tǒng)輸出為廣義模態(tài)力fext(t)，流體守恒變量w為系統(tǒng)的狀態(tài)變量。由于該線性化模型是基于當地CFD流場微擾動得到的，可以保留激波位置等流場非線性特性，所以該模型具有流體非線性特性。若CFD網格數為n，對于三維問題，w為5n1維矩陣，而A為5n5n維矩陣，對于一般的三維CFD問題n往往為105～107量級，如此大的矩陣處理大大耗費了計算內存，數學上的絕大多數算法對如此巨大維數矩陣都將失效，尤其求其特征值是非常困難的。所以本文通過降階方法將高階系統(tǒng)式(4)降階成低階系統(tǒng)，以便于響應快速系統(tǒng)分析以及設計。

1.3流體動力學降階系統(tǒng)

由于前面得到的式(4)的高維度，很難求得該系統(tǒng)的特征值，本文采用數據驅動的POD降階方法。POD降階方法是通過選擇一組數據樣本集來構造降階過程中所需要的變換矩陣來得到降階系統(tǒng)。給予式(4)系統(tǒng)結構模態(tài)位移脈沖輸入作為激勵，產生一系列樣本數據{wk}, wk∈Rn為n維向量，則m個n維樣本數據構造形成一個大型矩陣，記為：

式中：數據{wk}是給系統(tǒng)一個激勵信號而得到的響應數據，POD方法就是尋找一個最優(yōu)化的基底向量Ψr=(Ψ1,Ψ2,…,Ψr),(r?n)，使得投影到該基底子空間的降階系統(tǒng)，能夠很好的保持原系統(tǒng)的動力學特性，因此向量w可以用以下向量表示,

(5)

(6)

式(6)為流體降階系統(tǒng)，和式(4)所代表的全階系統(tǒng)具有同樣的性質。該系統(tǒng)的階數為r，一般為10～102量級，遠小于原始系統(tǒng)的階數。如果對于某些問題(控制設計等)若階數仍然較高，可采用平衡截斷法[15]將該系統(tǒng)進一步降階以達到要求。由于平衡截斷法可以保持系統(tǒng)的可觀性和可控性，同時保持原系統(tǒng)的穩(wěn)定性，所以對其進一步降階能保持系統(tǒng)特性。

1.4顫振系統(tǒng)狀態(tài)空間模型

根據前面方法得到低維度的流體降階狀態(tài)空間方程式(6)后，耦合結構動力學方程，構成一個低階的顫振系統(tǒng)狀態(tài)空間方程。和流體動力學線性化方法類似，結構動力學方程需在靜平衡狀態(tài)泰勒展開得到，結構動力學方程寫為

(7)

式中：M為質量陣，K為剛度陣，F為節(jié)點氣動力向量，X為節(jié)點位移向量。根據模態(tài)疊加原理，可得

X=Φu

(8)

式中：Ф為機翼模態(tài)質量歸一化后的前N階模態(tài)矩陣，代入式(7)后兩邊同時左乘ФT，根據固有模態(tài)正交性可得模態(tài)坐標下的微擾動結構動力學方程為：

(9)

(10)

式中：

一般流體與結構時間量綱不一致，流體求解是采用無量綱化可知t=τ·Vf/Lref=τ/ω，τ為有量綱時間即結構動力學方程時間量綱形式，Lref為特征長度，式(6)統(tǒng)一表達成結構時間量綱形式為

(11)

將式(11)和結構動力學方程式(10)耦合在一起，得到氣動彈性系統(tǒng)狀態(tài)空間方程

(12)

式中：

xflutter=[wr,xs]T，

式(12)為顫振分析系統(tǒng)的低維狀態(tài)空間模型，對于其系統(tǒng)穩(wěn)定性，采用根軌跡法和Lyapunov第一穩(wěn)定性原理來分別尋找氣彈系統(tǒng)顫振邊界，驗證所建立顫振分析系統(tǒng)的準確性，并分析比較不同穩(wěn)定性分析方法對原始非線性流固耦合系統(tǒng)穩(wěn)定性判斷的合理性。

2　數值分析

2.1二維兩自由度模型

對于二維模型，本文選擇Isogai Wing，該模型是一個典型的二維兩自由度跨音速氣動彈性標準算例，大致表征一個大展弦比跨音速后掠機翼外側翼面的彈性運動方式。翼型為NACA 64A010，其動力學方程可參考文獻[18]，結構參數見表1。

表1　Isogai Wing結構參數

首先采用根軌跡方法對Ma=0.75狀態(tài)下的穩(wěn)定性判斷，從圖1的沉浮模態(tài)在復平面的根軌跡可以得到在無因次顫振速度VF=1.36時軌跡線穿過虛軸，可知VF=1.36為該狀態(tài)的臨界顫振點，對降階模型在VF=1.36下，給予沉浮速度初始值作為初始擾動，觀察其在該狀態(tài)下的自由響應。從圖2可知，響應曲線呈等幅振動狀態(tài)，時域響應上判斷該狀態(tài)確實是臨界顫振點，驗證了根軌跡法對此狀態(tài)顫振點尋找的正確性。

圖1　沉浮模態(tài)根軌跡圖(Ma=0.75)Fig.1 Root locus of plunge model(Ma=0.75)

圖2　VF=1.36下降階系統(tǒng)的位移響應(Ma=0.75)Fig.2 System displacement response at VF=1.36

圖3　結構自由度特征值實部和系統(tǒng)最大特征值實部隨VF的變化(Ma=0.8)Fig.3 Real part of eigenvalue of strutural model and system’s maximum versus VF’s change(Ma=0.8)

如圖3和圖4結果所示，在Ma=0.8狀態(tài)下，分別采用根軌跡法和Lyapunov第一穩(wěn)定性判據來尋找該狀態(tài)的顫振點，根軌跡圖采用沉浮和俯仰模態(tài)對應特征值的實部隨著VF的變化來表示，而Lyapunov第一穩(wěn)定性判據則采用特征值實部最大值隨著VF變化表示。Lyapunov穩(wěn)定性判據是包含根軌跡方法的，比根軌跡方法描述系統(tǒng)穩(wěn)定性更加完整全面。

圖4　VF=0.91下降階系統(tǒng)的位移響應(Ma=0.8)Fig.4 System displacement response at VF=0.91(Ma=0.8)

當Ma=0.875狀態(tài)下，采用上述的方法得到結果見圖5?？芍^察結構自由度所對應的特征值變化并不能得到該馬赫數下的顫振點，而采用觀察系統(tǒng)最大特征值實部正負可以判斷系統(tǒng)穩(wěn)定性，很快得到該馬赫數下的顫振點。由圖中曲線可知，該馬赫數下有三個顫振點，分別是VF=0.61，VF=1.70，VF=2.36，由其值的正負可知，在VF<0.61狀態(tài)下，沒有正實部，系統(tǒng)是穩(wěn)定的；在0.612.36狀態(tài)下，實部最大值大于零，系統(tǒng)又不穩(wěn)定。

圖5　結構自由度特征值實部和系統(tǒng)最大特征值實部隨VF的變化(Ma=0.875)Fig.5 Real part of eigenvalue of strutural model and system’s maximum versus VF’s change(Ma=0.875)

為了驗證是否在該狀態(tài)有三個顫振點，采用降階系統(tǒng)和CFD對不同VF條件下的響應進行計算，圖6結果表明在所得到的顫振點，響應都是呈等幅振動，圖7則對CFD和降階系統(tǒng)響應的頻域比較，二者結果非常吻合，說明采用Lyapunov判據所得到的顫振點是非常準確完備的。

由于該狀態(tài)有多個顫振點，這就導致了S型顫振邊界的出現，結果與參考文獻[18-19]趨勢大小接近，見圖8。而根軌跡法對此狀態(tài)的預測則失效，說明在該狀態(tài)下，對POD法所建立的ROM來觀察結構模態(tài)對應的根軌跡并不能找到顫振點，其穩(wěn)定性主要由跨音速流體非線性特性決定。

圖6　不同VF下降階系統(tǒng)的位移響應(Ma=0.875)Fig.6 System displacement response at different(Ma=0.875)

圖7　不同VF下CFD和降階系統(tǒng)頻域比較(Ma=0.875)Fig.7 Frequency domain comparsion of CFD and ROM’s system(Ma=0.875)

圖8　Isogai Wing隨馬赫數變化的顫振邊界Fig.8 The flutter boundary of Isogai Wing

2.2三維機翼模型

本文的三維POD/ROM是基于三維CFD/CSD求解器基礎上發(fā)展而來的，計算模型采用國際上用于檢驗顫振計算方法的標準模型——AGARD 445.6[20]。由于CFD/CSD求解器對于該模型的驗證已在文獻[21]中詳細討論，所以本文直接采用ROM方法尋找隨Ma變化的顫振邊界并和原始CFD/CSD系統(tǒng)計算比較。流體網格分為22塊，總數為202 568，見圖9，其全階自由度為202 568×5=1 012 840。對于三維模型采用POD降階后的階數仍然較高，本文采用平衡截斷法將該系統(tǒng)進一步降階以達到要求。首先通過POD方法將全階線性化方程降階為500階，然后通過平衡截斷法降階成50階，結構動力學采用前4階模態(tài)代替，耦合結構動力學方程，式(12)的維數為58。

圖9　AGARD 445.6機翼物面和對稱面網格Fig.9 Fluid grids of AGARD 445.6 wing’s surface and symmetry plane

采用觀察系統(tǒng)狀態(tài)空間方程式(12)是否包含正實部特征值方法，來判斷非線性流固耦合系統(tǒng)是否穩(wěn)定。圖10表示不同馬赫數下系統(tǒng)特征值實部最大值隨遠場速度的變化。采用該方法可以很有效方便的找到對應馬赫數下的顫振邊界。

圖10　不同馬赫數下的顫振點的尋找Fig.10 The flutter point searching at different Mach number

為了驗證所得到的顫振邊界是否準確，圖11比較了降階模型、CFD/CSD耦合計算以及實驗的結果，可以看出結果非常吻合，說明該方法快速尋找顫振邊界和CFD/CSD耦合計算具有同樣的精度。計算效率上，采用降階模型對于指定馬赫數下尋找顫振邊界，計算時間只需幾秒鐘，而CFD/CSD耦合計算需要反復試算尋找顫振邊界，指定馬赫數下顫振邊界的尋找往往需要數天。更為重要的是，該低維狀態(tài)空間方程可以很方便有效的與控制系統(tǒng)耦合，為顫振抑制提供一個高精度數學模型。

圖11　AGARD 446.5 顫振邊界比較Fig.11 The flutter boundary of AGARD 445.6 wing by different method

3　結　論

本文根據Lyapunov穩(wěn)定性理論對非線性流固耦合系統(tǒng)進行穩(wěn)定性分析，采用微擾動理論和POD降階建立了基于CFD的氣動彈性系統(tǒng)低維狀態(tài)空間模型。采用根軌跡法和Lyapunov第一穩(wěn)定性原理對二維兩自由度顫振模型在不同馬赫數下，進行顫振邊界的尋找。結果分析表明：

(1) 若對POD法建立的ROM采用根軌跡法尋找顫振邊界，在亞音速階段可以與Lyapunov得到的結果一致，而在跨音速區(qū)域，觀察結構模態(tài)對應的根軌跡并不能找到顫振點，根軌跡法失效。說明跨音速階段其穩(wěn)定性并不是由某一結構模態(tài)決定的，而是主要由跨音速流體非線性特性決定。

(2) Lyapunov穩(wěn)定性理論針對原始非線性系統(tǒng)，對系統(tǒng)最大特征值實部正負分析，其理論上包含根軌跡法所得結果，因而更全面反映非線性系統(tǒng)的穩(wěn)定性。

(3) 最后采用Lyapunov第一穩(wěn)定性原理對三維氣彈模型進行顫振邊界的尋找，所得結果很好的與CFD/CSD耦合計算所得結果一致，大大提高了尋找顫振邊界的效率。

[1] 楊國偉. 計算氣動彈性若干研究進展[J]. 力學進展, 2009, 39(4): 406-420.

YANG Guowei. Recent progress on computational aeroelasticity[J]. Advances in Mechanics, 2009, 39(4): 406-420.

[2] 陳剛, 李躍明. 非定常流場降階模型及應用研究進展與展望[J]. 力學進展, 2011, 41(6): 686-701.

CHEN Gang, LI Yueming. Advances and prospects of the reduced order model for unsteady flow and its application[J]. Advances in Mechanics, 2011, 41(6): 686-701.

[3] SILVA W A, BARTELS R E. Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code[J]. Journal of Fluids and Structures, 2004, 19(6): 729-745.

[4] 劉曉燕，楊超，吳志剛. 適用于氣動彈性的小波非定常氣動力降階方法[J]. 航空學報,2010,6(31): 1149-1155.

LIU Xiaoyan, YANG Chao, WU Zhigang. Wavelet-based reduced-order method for unsteady aerodynamics applicale to aeroelasticity[J]. Acta Astronautica Sinica, 2010,6(31): 1149-1155.

[5] 聶雪媛，劉中玉，楊國偉. 基于CFD氣動力辨識模型的氣動彈性數值計算[J]. 振動與沖擊,2014,33(20):20-25.

NIE Xueyuan, LIU Zhongyu, YANG Guowei. Identification of unsteady model CFD-based for aeroelastic numerical computation [J]. Journal of Vibration and Shock,2014,33(20): 20-25.

[6] RAVEH D E, ZAIDE A. Numerical simulation and reduced-order modeling of airfoil gust response[J]. AIAA Journal, 2006,44(8):1826-1834.

[7] HUANG R, HU H, ZHAO Y. Nonlinear reduced-order modeling for multiple-input/multiple-output aerodynamic systems[J]. AIAA Journal,2014,52(6):1219-1231.

[8] BADCOCK K J, SEBASTIAN T, SIMAO M, et al. Transonic aeroelastic simulation for instability searches and uncertainty analysis[J]. Progress in Aerospace Sciences, 2011, 47(5): 392-423.

[9] DOWELL E H, HALL K C, ROMANOWSKI M C. Eigenmode analysis in unsteady aerodynamics: reduced order models[J]. Applied Mechanics Reviews,1997,50(6):371-386.

[10] WILLCOX K, PERAIRE J. Balanced model reduction via the proper orthogonal decomposition[J]. AIAA Journal,2002,40(11):2323-2330.

[11] EPUREANU B. A parametric analysis of reduced order models of viscous flows in turbomachinery[J]. Journal of Fluids and Structures, 2003, 17(7): 971-982.

[12] AMSALLEM D, FARHAT C. Interpolation method for adapting reduced-order models and application to aeroelasticity[J]. AIAA Journal,2008,46(7):1803-1813.

[13] CHEN G, SUN J, LI Y M. Adaptive reduced-order model-based control-law design for active flutter suppression[J]. Journal of Aircraft, 2012, 49(4): 973-980.

[14] 張偉偉，葉正寅. 基于非定常氣動力辨識技術的氣動彈性數值模擬[J]. 航空學報, 2006, 27(4): 579-583.

ZHANG Weiwei, YE Zhengyin. Numerical simulation of aeroelasticity basing on indentification technology of unsteady aerodynamic loads[J].Acta Astronautica Sinica,2006,27(4):579-583.

[15] 陳剛，李躍明，閆桂榮,等. 基于POD降階模型的氣動彈性快速預測方法研究[J]. 宇航學報,2009,30(5):1765-1770.

CHEN Gang, LI Yueming, YAN Guirong, et al. A fast aeroelastic response prediction method based on proper orthogonal decomposition reduced order model[J]. Journal of Astronautics,2009,30(5):1765-1770.

[16] LIEU T, FARHAT C, LESOINNE M. Pod-based aeroelastic analysis of a complete f-16 conguration: rom adaptation and demonstration[C]//46th AIAA/ASME/ACSE/AHS/ASC Structure, Structural Dynamics & Materials conference. Austin,TX, 2005.

[17] LESOINNE M, SARKIS M, HETMANIUK U, et al. A linearized method for the frequency analysis of three-dimensional fluid/structure interaction problems in all flow regimes[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(24): 3121-3146.

[18] FENG L, CAI J, ZHU Y, et al. Calculation of wing flutter by a coupled fluid-structure method[J]. Journal of Aircraft, 2001, 38(2): 334-342.

[19] ALONSO J J, JAMESON A. Fully-implicit time-marching aeroelastic solutions[C]//32th AIAA Aerospace Sciences Meeting and Exhibit. Reno,NV, 1994.

[20] CARSON Y E. Agard standard aeroelastic configurations for dynamic response I-wing 445.6[S].AGARD-R-756,1988.

[21] 周強, 陳剛, 李躍明. 復雜外形飛行器跨音速氣動彈性CFD/CSD耦合數值模擬研究[C]//第十三屆全國空氣彈性學術交流會會議.哈爾濱, 2013.

Nonlinear aeroelastic stability analysis based on CFD reduced order models

ZHOU Qiang，CHEN Gang，LI Yueming

(State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China)

In order to quickly find the flutter boundary of a nonlinear aeroelstic system based on CFD/CSD, this paper used the Lyapunov stability theory to analyze the stability of a nonlinear fluid-structure coupling system. Firstly, a linearized state space model was obtained based on a nonlinear fluid-structure coupling system through the small perturbation theory; then a reduced order model was achieved by reducing the high dimensional linearized model through a POD(proper orthogonal decomposition) method. According to the eigenvalues of the ROM system, the stability of the original nonlinear system could be determined. Lyapunov stability theory was mainly aimed at the nonlinear system, and the POD method was adopted to realize the process. Different from other stability methods, the mathematical theory reflected the stability of the original nonlinear fluid-structure coupling system. As the POD/ROM used in this paper was based on the CFD fluid field, so it could better reflect the internal characteristics of the original nonlinear coupling system. Numerical cases including two-dimensional aerofoil and three-dimensional wing’s models were used to validate the accuracy of the stability method. The results show that in the subsonic domain, the stability is determined by the structural model; but in transonic and supersonic domain, the aeroelastic stability is mainly affected by fluid characteristics.

fluid-structure interaction; Lyapunov; proper orthogonal decomposition (POD); state-space; flutter

國家自然科學基金(11272005;11472206;11371288)

2015-07-14修改稿收到日期：2015-08-23

周強 男，博士生，1989年生

陳剛 男，教授，博士生導師，1979年生

V211；V215.3

A

10.13465/j.cnki.jvs.2016.16.004