適用于非周期流固耦合問題的時間譜方法

楊體浩，白俊強，史亞云，楊一雄

西北工業大學 航空學院，西安 710072

飛行器設計是涉及氣動、結構、控制和飛行力學等其他學科在內的復雜的多學科問題。通常這些學科之間會相互耦合，尤其是氣動與結構之間會產生明顯的流固耦合現象，除了靜氣動彈性變形以外，還包括顫振等在內的動氣動彈性現象[1-2]。這些流固耦合現象會對飛行器的性能、穩定性和安全性造成顯著的影響。隨著現代飛行器展弦比的不斷增大以及機翼柔性的增加，氣動、結構之間的耦合問題也變的越來越突出[3]。這使得發展高效、準確的動氣動彈性仿真、分析工具以及考慮動氣動彈性影響(如顫振邊界)的氣動/結構多學科優化設計系統成為決定未來先進飛行器的研制能否成功的關鍵因素之一。

直接將計算流體力學(Computational Fluid Dynamics, CFD)程序與計算結構動力學(Computational Structural Dynamics, CSD)程序耦合，建立一套時域的氣動、結構仿真分析和優化設計工具是最直接的方式。但是，這種基于高精度求解器的時域方法需要花費大量的計算時間，無法滿足實際的工程應用需求。因此，針對考慮動氣動彈性影響的氣動/結構優化設計與分析問題，需要發展新的方法以提高計算效率。

目前針對非定常問題，降階模型(Reduced-Order Model, ROM)被視為一類很有應用價值的方法。在氣動彈性以及優化設計領域，本征正交分解(Proper Orthogonal Decomposition, POD)[4]是一種最為常見的降階模型。POD以一組基于最小二乘思想的正交基去近似高維問題。但是POD模型的建立需要采集樣本快照，借助高精度的CFD求解器生成樣本快照需要花費大量的計算時間，這在一定程度上降低了POD的計算效率。除此之外，包括POD在內的許多降階模型缺乏對系統參數變化的魯棒性，比如馬赫數、結構參數和機翼型面等參數的變化[5]。系統可變參數的個數越多，參數變化范圍越大，降階模型的精度和魯棒性越差。然而，許多實際的工程應用問題包含了大量的可變參數。因此，借助降階模型難以對具有大量變參的系統進行可靠的靈敏度分析。這也使得在優化設計領域，無法將降階模型與基于梯度信息的高效優化設計方法相結合，如伴隨方法[6]。

針對非定常問題，諧波平衡(Harmonic Balance, HB)[7]是另一種很有吸引力的方法。HB通過將系統控制方程的狀態變量用截斷的傅里葉系數進行替換，最終將非定常的控制方程轉化為耦合的定常控制方程。HB直接對系統控制方程進行處理，不需要樣本快照。因此，HB具有很高的計算精度和魯棒性。借助HB方法，可以準確地對整個系統進行靈敏度分析，并且很容易與伴隨方法耦合對非定常問題進行優化設計，大大地提高優化設計效率。Choi等[8]將HB與伴隨方相結合建立了非定常氣動伴隨優化設計系統，并以此對直升機旋翼進行了氣動外形優化設計。但是，HB方法只適用于周期性的非定常問題，具有一定的局限性，而很多實際的流固耦合問題都是非周期的。

除了降階模型和HB方法，采用效率更高的低精度方法是另一種選擇。這類方法控制方程形式簡單，求解速度快，非常適用于飛行器概念設計以及初始設計階段。通過耦合修正模型，針對很多非定常問題，低精度方法都可以達到期望的計算精度。Timme[9]利用耦合了全速勢方程和非定常邊界層方程的求解器，對跨聲速機翼進行了氣動彈性穩定性研究。研究結果表明，當不存在明顯的流動分離時，全速勢方程計算的結果與基于非定常雷諾平均Navier-Stokes(URANS)方程的計算結果以及實驗數據都較為接近。對于低亞聲速問題，基于勢流理論的流固耦合分析方法被廣泛的用來研究包括經典顫振問題在內的多種動氣動彈性問題。NASTRAN就是一種被廣泛應用于氣動彈性問題研究的基于勢流理論的商業軟件[10]。Hesse和Palacios[11]利用非定常渦格法和復合材料梁有限元模型對大柔性飛行器的飛行動力學問題進行了成功的研究。

除此之外，基于勢流理論的流固耦合分析方法還被廣泛的與其他優化設計技術相結合，進行考慮動氣動彈性影響的優化設計研究。Haghighat等[12]利用渦格法、梁有限元模型和智能優化算法針對大展弦比柔性機翼進行了考慮陣風減緩的多學科優化設計研究。Stanford和Beran[13]利用Theodorsen氣動力模型、非線性梁有限元模型和ONERA動失速修正模型，借助設計手段，通過優化結構質量和剛度分布去改變大展弦比柔性機翼的顫振邊界和極限環振蕩(Limit Cycle Oscillation, LCO)特性。Mallik等[14]利用非線性結構有限元模型和考慮了壓縮性修正的Theodorsen氣動力模型進行顫振邊界的計算，并借助遺傳算法將其與傳統針對靜力學問題的氣動/結構多學科優化設計系統相結合，針對支撐翼(Truss-Braced Wing, TBW)布局飛機探討了顫振約束對設計結果的影響。其研究結果表明，針對TBW，是否考慮顫振約束對整個優化設計結果的性能有較為明顯的影響。

目前國際上，在考慮了動氣動彈性影響(尤其是顫振邊界)的優化設計領域方面，主流的做法還是利用頻域法或者直接時域法進行動氣動彈性特性的計算分析，然后借助智能算法將其引入到優化設計過程中。這種借助智能算法的優化設計方法計算量很大，尤其是當整個系統的設計變量個數較多，或者優化設計系統的靜力學問題部分采用高精度求解器的時候。除此之外，也有一部份學者直接把伴隨方法與傳統的時域流固耦合分析方法進行結合，將動氣動彈性的影響引入到氣動結構多學科優化設計中。如，Zhang等[15]基于Euler和復雜有限元模型求解器，建立了針對非定常流固耦合問題的伴隨優化設計系統，并對M6機翼進行了考慮顫振約束的氣動結構多學科優化設計研究。雖然伴隨方法的采用克服了基于智能算法的優化設計系統的缺點，但是直接將伴隨方法與傳統的時域流固耦合求解方法相結合的處理方式，需要在每一個時間推進步上求解伴隨方程。通常非定常流固耦合問題的時域推進求解過程需要成百上千甚至更多的時間推進步，這意味著采取直接耦合的方式需要進行成百上千甚至更多次的伴隨方程的求解，整個求解過程計算量很大，大大削弱了基于伴隨方法的氣動/結構多學科優化設計系統的計算效率。因此，需要發展新的方法，以顯著提針對非定常流固耦合多學科設計問題的伴隨優化設計系統的計算效率。

在飛行器設計領域，針對非定常流固耦合問題建立高效、可靠的設計、分析工具十分必要。然而，現有的方法都無法很好地滿足這一需求。針對這一問題，本文的主要研究目的是針對一般運動形式的低亞聲速問題，將Chebyshev譜方法與基于非定常面元法、幾何非線性梁有限元模型的流固耦合分析方法進行耦合，建立面向考慮動氣動彈性影響(主要考慮經典顫振問題)的氣動/結構一體化設計的，可與伴隨方法相結合的高效、可靠的流固耦合分析方法。Chebyshev譜方法直接對流固耦合控制方程進行處理，將非定常問題轉化為在選取的Chebyshev控制點處耦合的定常問題。因此，與HB方法類似，Chebyshev譜方法具有較高的計算精度和可靠性，并且很容易與伴隨方法結合建立針對非定常流固耦合問題的伴隨優化設計系統。由于Chebyshev譜方法將非定常問題轉化為在選取的Chebyshev控制點處耦合的定常問題進行求解，因此在與伴隨方法進行結合時，只需在選取的Chebyshev控制點處進行伴隨方程的求解即可，以此大大地節省計算量提高優化設計效率。但是與HB方法不同的是，Chebyshev譜方法不僅適用于周期性問題，還適用于非周期性問題。Choi等[16-17]將Chebyshev譜方法與URANS求解器相結合進行了非定常氣動力的計算分析，得到了令人滿意的求解精度。雖然，他成功地利用Chebyshev譜方法進行了非定常氣動力的計算分析，但是目前其研究成果具有比較明顯的局限性，距離實現Chebyshev譜方法與伴隨方法相結合進行非定常流固耦合問題優化設計的目標還較遠。一方面，Choi等引入的是周期性邊界條件，這意味著雖然采用的是Chebyshev譜方法，但是處理的依舊是達到了穩定狀態的周期性非定常問題，而自然界中很多現象是非周期性的；另一方面，Choi等僅僅將Chebyshev譜方法與非定常氣動模型進行了耦合，還沒有將其與結構動力學模型進行耦合建立相應的流固耦合分析方法。相比之下，本文將Chebyshev譜方法同時與非定常氣動力模型以及結構動力學模型進行耦合，建立了相應的流固耦合分析方法。同時，向系統中引入初始條件而非周期性邊界條件，使得所建立的流固耦合分析方法可以計算非定常問題的瞬態響應，即適用于周期性非定常問題也適用于非周期性的非定常問題。

本文面向大展弦比飛行器的非定常流固耦合多學科優化設計問題，采用將Chebyshev譜方法與基于非定常面元法和幾何非線性梁有限元模型的流固耦合分析方法進行結合，建立了針對低亞音速流動的高效、可靠的，可與伴隨方法相結合的流固耦合分析方法。所建立的分析方法主要針對不考慮氣動非線性的動氣動彈性問題，如經典的顫振問題。通過Goland機翼顫振速度的計算對建立的基于Chebyshev譜方法的流固耦合分析系統的精度、可靠性和收斂性等特點進行了探討研究。

1 Chebyshev譜方法

1.1 標準的Chebyshev譜方法

Chebyshev譜方法是一種以Chebyshev多項式插值為理論基礎的時間譜方法[18]。理論上，任意一個定義在時間區間[-1,1]內的光滑函數u(t)，在tm時刻的函數值可以表示為

(1)

Ti(tm)=cos(icos-1(tm))i=0,1,…,N

(2)

(3)

式中：Ti為Chebyshev多項式;ai為Chebyshev系數;N為在時間區間[-1,1]內選取的Chebyshev控制點的個數。顯然，一旦Chebyshev控制點處的函數值已知，那么在定義的時間區間內，任意時刻的函數值可以通過式(1)求得。

進一步，函數u(t)在Chebyshev控制點處的時間導數可以推導為

(4)

(5)

(6)

(7)

1.2 Chebyshev算子的投影變換

Chebyshev控制點在時間區間內的分布情況會對其計算結果的精度產生一定的影響。標準的Chebyshev譜方法的控制點在時間區間的兩端分布較為密集，這種分布形態使得標準的Chebyshev譜方法需要更多的控制點個數去達到期望的計算精度[19]。因此，采用Kosloff和Tal-Ezer[20]提出的投影變換函數改變Chebyshev譜方法的控制點分布，使其分布更為均勻。投影變換的表達式為

(8)

通常對于實際的非定常問題，時間區間是[t0,tt]而不是[-1,1]，因此需要通過等價變換將時間區間從[t0,tt]變化到[-1,1]。采用線性變換函數：

(9)

(10)

1.3 驗證算例

經過投影變換的Chebyshev算子被稱為投影的Chebyshev算子(Mapped Chebyshev Operator，MCO)。利用非周期的三角函數驗證投影的Chebyshev算子的計算精度，揭示投影的Chebyshev譜方法的數學特性。三角函數的表達式為

(11)

(12)

圖1為具有3種不同Chebyshev控制個數的投影的Chebyshev算子的計算結果與解析解的對比，其中圖1(a)為計算得到的函數值u的對比，圖1(b)為計算得到的函數值u相對于解析解uAnalytic的絕對誤差對比。圖2為Chebyshev控制點在定義的時間區間內的分布示意圖。顯然，投影的Chebyshev算子的計算精度隨著Chebyshev控制點個數的增加而不斷增大。當控制點的個數增加到14，投影的Chebyshev算子的計算結果幾乎與解析解重合。

圖1 投影的Chebyshev算子的計算結果與解析解的對比Fig.1 Comparisons between results obtained from mapped Chebyshev operator and analytic solution

圖2 投影的Chebyshev控制點分布Fig.2 Distribution of mapped Chebyshev control points

計算結果表明，即便對于非周期函數，Chebyshev算子也可以用較少的控制點的個數得到滿足精度要求的計算結果。

2 基于投影的Chebyshev譜方法的非定常氣動力模型

2.1 非定常面元法

非定常面元法由Laplace方程推導演化而來。本文采用具有定常面元強度的非定常面元法[21]，并且忽略黏性等氣動非線性的影響。其在(n+1)dt時刻的控制方程為

(13)

(14)

(15)

(16)

圖3 尾跡模型Fig.3 Wake model

通過控制方程可以求出物面面元的強度，進而作用在物面面元上的壓力系數可以通過式(17)獲得：

(17)

式中：Cpk為第k個物面面元控制點處的壓力系數；pk和pref分別為作用在第k個物面面元控制點處的靜壓和無窮遠處自由來流的靜壓；ρ為大氣密度；Qk為第k個物面面元控制點處的當地流體速度；Φk為第k個物面面元控制點處的速度勢。

2.2 非定常面元法與Chebyshev譜方法的耦合

根據Chebyshev譜方法的理論可知，將Chebyshev譜方法與非定常面元法結合以后，只有在Chebyshev控制點處的狀態變量為待求未知變量。顯然，在給定的時間區內的任意時刻，物面的運動速度和位置形狀可表示為

(18)

(19)

根據尾跡形狀的遞推關系式(式(16))可知，在ti=t0+ndt時刻，尾跡的形狀為

μwake,μsurface,σsurface,Vgust)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

顯然，當初始條件給定后，通過迭代求解上述控制方程可以獲得在選定的Chebyshev控制點處的面元的強度。進而借助Chebyshev算子可以計算得到定義的時間區間內任一時刻下的面元強度。聯立式(10)與式(17)可計算得到對應狀態下所有物面面元的壓力系數：

(27)

2.3 驗證算例

選取做非周期強迫俯仰運動的大展弦比等直機翼作為驗證耦合了Chebyshev譜方法的非定常面元法計算精度的算例。機翼由NACA0012翼型成型，展弦比為36，物面面元的劃分如圖4所示。強迫的俯仰運動迎角為

(28)

式中:迎角α0=0°;運動周期T=0.4 s。等直機翼以過50%弦點的平行于前緣的直線為軸做俯仰運動。

圖5為Chebyshev譜方法(Chebyshev Spec-tral Method, CSM)與傳統的時間推進方法計算得到的非定常氣動力系數對比圖。Chebyshev譜方法的控制點的個數為18。計算結果顯示，無論是升力系數CL還是力矩系數Cm，Chebyshev譜方法的計算結果與傳統的時間推進方法所得到的結果都有很好的吻合度。顯然，針對非周期問題，耦合了Chebyshev譜方法的非定常面元法能夠以較少的控制點個數獲得具有足夠精度的計算結果。

圖4 機翼表面面元Fig.4 Surface panels of wing

圖5 Chebyshev譜方法與時間推進方法計算得到的非定常氣動力對比Fig.5 Comparison of unsteady aerodynamic forces calculated by Chebyshev spectral method and time-marching method

3 基于Chebyshev譜方法的結構動力學模型

3.1 幾何非線性梁有限元模型

采用基于旋轉矢量的幾何非線性梁有限元模型進行結構動力學建模[22]。相比于復雜有限元模型，梁模型自由度很少，形式簡單，非常適用于飛行器概念以及初始設計階段的多學科問題的研究。對于傳統的大展弦比機翼(如B737機翼)，線性梁有限元模型能夠較為準確地捕捉機翼結構的靜力學和動力學行為。但是對于展弦比更大的柔性機翼(如美國NASA的“Helios”太陽能無人機)，在典型的使用工況下機翼結構便表現出明顯的幾何非線性特性。這種情況下，線性梁有限元模型已經不再適用，需要借助幾何非線性梁有限元模型。這類大展弦比、大柔性機翼雖然發生大變形，但是主要是彎曲大變形，而扭轉變形相對較小，對于諸如經典顫振問題這類動氣動彈性問題而言，線性氣動力模型依然適用。因此，為使所建立的方法更具一般性，本文采用可考慮幾何非線性的梁有限模型進行結構動力學建模，并將其與線性非定常氣動力模型耦合，針對大展弦比機翼進行不考慮氣動非線性的流固耦合問題研究，如經典的顫振問題。

基于旋轉矢量的幾何非線性梁有限元模型采用矢端向量RG(s,t)和笛卡爾旋轉矢量(Cartesian Rotation Vector，CRV)ψ(s,t)[23]作為描述梁模型上任意一點的位移以及扭轉變形自由度,其中s為在梁參考線上定義的弧長坐標，t表示時間。坐標系的定義如圖6所示。

圖6 坐標系的定義Fig.6 Definition of coordinate systems

t時刻梁模型上任意一點，相對于未變形構型(即t=0時刻的構型)的應變可表示為

(29)

式中：γ和κ分別為線應變和角應變；CBG表示從慣性坐標系G到材料坐標B的坐標變換矩陣；(·)′表示相對于弧長坐標s的導數；KB(s,t)=T(ψ)ψ′，T(ψ)為切向算子[22]。

根據Hamilton原理有：

(30)

(31)

(32)

(33)

式中：M、C和K分別為切向質量矩陣、切向阻尼矩陣和切向剛度矩陣。考慮幾何非線性時，M、C和K不是常值矩陣，而是與狀態變量相關的函數。基于增量形式的結構動力學方程，利用相應的數值求解方法(如Newton-Raphson方法)進行迭代求解，即可獲得梁模型的結構動力學響應。

3.2 幾何非線性梁有限元模型與Chebyshev譜方法的耦合

幾何非線性梁有限元模型可以寫成結構動力學方程的一般形式：

(34)

通過等價變換將二階微分方程轉化為一階微分方程：

(35)

(36)

(37)

(38)

(39)

顯然，借助Chebyshev譜方法，非定常形式的結構動力學模型被轉化為在選定的Chebyshev控制點處的定常方程。當給定初始條件后，通過迭代求解控制方程可以獲得結構狀態變量在所有Chebyshev控制點處的值。之后，利用Chebyshev算子可以計算得到在給定時間區間內的結構動力學響應。

3.3 驗證算例

選取端部承受強迫載荷的懸臂梁(見圖7)作為驗證耦合了Chebyshev譜方法的幾何非線性梁有限元模型的計算精度的算例。懸臂梁長l=16 m，結構屬性見表1。強迫載荷F的方向豎直向上，在整個時間歷程中方向保持不變，力的大小隨著時間呈現如圖8所示的變化(圖中Famp為強迫載荷F允許的最大幅值)。

圖7 施加強迫載荷的懸臂梁Fig.7 Cantilever beam with forced load

表1 梁結構屬性Table 1 Properties of beam structure

圖8 強迫載荷的時間歷程Fig.8 Time history of forced load

圖9為Chebyshev譜方法與傳統的時間推進方法計算得到的結構動力學響應(梁端部的撓度變形Z以及縱向速度Vz響應)對比圖。其中Chebyshev譜方法的控制點的個數為18，時間推進方法的時間步長為0.001 s。計算結果表明，無論是位移響應還是速度響應，Chebyshev譜方法的計算結果都幾乎與時間推進方法的計算結果重合。顯然，耦合了Chebyshev譜方法的幾何非線性梁有限元模型，借助有限的控制點個數，便可以以足夠高的精度準確地捕捉整個結構系統的動力學特性。

圖9 Chebyshev譜方法與時間推進方法計算得到的結構動力學響應對比Fig.9 Comparison of structural dynamic response calculated by Chebyshev spectral method and time-marching method

4 基于Chebyshev譜方法的流固耦合分析系統

利用松耦合迭代策略求解基于Chebyshev譜方法的流固耦合控制方程。由于Chebyshev譜方法分別與非定常氣動力模型以及結構動力學模型相結合，將非定常的流固耦合問題轉化為在選定的Chebyshev控制點處耦合的定常問題。因此，基于Chebyshev譜方法的動氣動彈性問題的求解過程，類似于傳統采用松耦合策略的靜氣動彈性問題的求解，只不過所有Chebyshev控制點處的控制方程是同時求解的。這種不同Chebyshev控制點處，控制方程的耦合實際上是整個時間響應歷程內不同Chebyshev控制點處的狀態變量之間時間依賴性的直接體現。

圖10 基于Chebyshev譜方法的流固耦合問題求解流程Fig.10 Calculation flow chart of fluid-structure coupling problems based on Chebyshev spectral method

5 數值仿真與討論

對標準模型Goland機翼[24]進行顫振速度的計算，以驗證本文搭建的基于Chebyshev譜方法的流固耦合計算、分析方法的收斂特性及計算精度。Goland機翼半展長6.096 m，弦長1.828 8 m。采用國際上針對Goland機翼計算分析使用的梁有限元相關參數[24]，具體的結構屬性如表2所示。圖11為Goland機翼表面面元的劃分。

圖12為Chebyshev譜方法與時間推進方法計算得到的不同自由來流速度V∞下的翼尖撓度變形Z的動力學響應對比圖。Chebyshev譜方法的控制點個數為22。Goland機翼的初始迎角為0°，在t=0 s時刻，給Goland機翼施加一個大小為0.055°的小迎角擾動。計算結果顯示，在172 m/s的自由來流速度下，Goland機翼的振動成發散的趨勢。當自由來流速度減小到167.5 m/s時，Goland機翼維持等幅地周期性振動。在兩種不同的自由來流速度下，Chebyshev譜方法的計算結果都幾乎與時間推進方法(圖中“Time-marching Method”表示采用傳統的時間推進方法，計算得到的流固耦合系統在給定初始條件下的瞬態響應)的所得結果完全重合。因此，本文建立的基于Chebyshev譜方法的流固耦合計算、分析方法能夠比較準確地捕捉流體與固體之間的瞬態響應。

表2 Goland機翼屬性Table 2 Properties of Goland wing

圖11 Goland機翼表面面元Fig.11 Surface panels of Goland wing

表3為不同模型計算得到的Goland機翼顫振速度Vf的對比。顯然，相比于Patil等[25]利用片條理論計算得到的Goland機翼的顫振速度，本文的計算結果與采用幾何非線性梁有限元模型和非定常渦格法的SHARP程序[22]的計算結果更為接近，顫振速度相差不超過1%。基于片條理論得到的Goland機翼顫振速度偏低的主要原因在于經典的片條理論無法捕捉三維效應。對比結果顯示，本文基于Chebyshev譜方法、非定常面元法和幾何非線性梁有限元模型所建立的流固耦合計算方法具有足夠的可信度。

圖13為采用了松耦合求解策略的Chebyshev譜方法的收斂歷程。其中Iterarion=0表示給定的系統初值。計算結果顯示，經過大約8次迭代以后，基于Chebyshev譜方法的流固耦合計算方法便達到了收斂。收斂歷程顯示，整個迭代過程呈現出從初始時刻向終了時刻逐步收斂的趨勢。其原因在于本文所處理的流固耦合問題為初值問題。初始條件所包含的信息隨著迭代計算逐步向時間下游傳遞。

圖12 不同自由來流速度下的Goland機翼翼尖時間歷程響應對比Fig.12 Comparison of time history of Goland wing tip at various free-stream velocities

表3 Goland機翼顫振速度對比Table 3 Comparison of flutter speed of Goland wing

圖13 Chebyshev譜方法的收斂歷程Fig.13 Convergence process of Chebyshev spectral method

6 結 論

1) Chebyshev譜方法利用Chebyshev算子替換系統的狀態變量及其導數，將非定常問題轉化為選取的Chebyshev控制點處耦合的定常問題進行處理。由于Chebyshev譜方法直接對流固耦合的控制方程進行處理，因此所建立的基于Chebyshev譜方法的流固耦合計算、分析方法具有較強的魯棒性。

2) Chebyshev譜方法的計算精度隨著Chebyshev控制點個數的增加而不斷提高，測試算例表明，只需少許有限的Chebyshev控制點，Chebyshev譜方法便可以達到期望的計算精度。

3) Chebyshev譜方法具有較強的適用性，不僅適用于周期性非定常問題還適用于非周期性非定常問題。

4) 建立的基于Chebyshev譜方法的流固耦合分析方法能夠準確地捕捉流體與結構之間的耦合作用，可用于忽略了氣動非線性的流固耦合問題的計算、分析與設計研究，如經典的顫振問題。

參 考 文 獻

[1] LIVNE E, WEISSHAARW T A. Aeroelasticity of nonconventional airplane configurations-past and future[J]. Journal of Aircraft, 2003, 40(6): 1047-1065.

[2] XIANG J, YAN Y, LI D. Recent advance in nonlinear aeroelastic analysis and control of the aircraft[J]. Chinese Journal of Aeronautics, 2014, 27(1): 12-22.

[3] SU W, CESNIK C E S. Strain-based geometrically nonlinear beam formulation for modeling very flexible aircraft[J]. International Journal of Solids and Structures, 2011, 48(16): 2349-2360.

[4] LIEU T, FARHAT C. POD-based aeroelastic analysis of a complete F-16 configuration: ROM adaptation and demonstration: AIAA-2005-2295[R]. Reston, VA: AIAA, 2005.

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

[6] LYU Z, KENWAY G K, PAIGE C, et al. Automatic differentiation adjoint of the reynolds-averaged navier-stokes equations with a turbulence model: AIAA-2013-2581[R]. Reston, VA: AIAA, 2013.

[7] HALL K C, THOMAS J P, CLARK W S. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique[J]. AIAA Journal, 2002, 40(5): 879-886.

[8] CHOI S, LEE K, POTSDAM M M, et al. Helicopter rotor design using a time-spectral and adjoint-based method[J]. Journal of Aircraft, 2014, 51(2): 412-423.

[9] TIMME S. Transonic aeroelastic instability searches using a hierarchy of aerodynamic models[D]. Liverpool: University of Liverpool, 2010.

[10] PERERA M, GUO S. Structural and dynamic analysis of a seamless aeroelastic wing: AIAA-2010-2878[R]. Reston, VA: AIAA, 2010.

[11] HESSE H, PALACIOS R. Reduced-order aeroelastic models for dynamics of maneuvering flexible aircraft[J]. AIAA Journal, 2014, 52(8): 1717-1732.

[12] HAGHIGHAT S, RA MATRINS J R, LIU H H. Aeroservoelastic design optimization of a flexible wing[J]. Journal of Aircraft, 2012, 49(2): 432-443.

[13] STANFORD B, BERAN P. Direct flutter and limit cycle computations of highly flexible wings for efficient analysis and optimization[J]. Journal of Fluids and Structures, 2013, 36: 111-123.

[14] MALLIK W, KAPANIA R K, SCHETZ J A. Effect of flutter on the multidisciplinary design optimization of truss-braced-wing aircraft[J]. Journal of Aircraft, 2015, 52(6): 1858-1872.

[15] ZHANG Z, CHEN P C, WANG Q, et al. Adjoint based structure and shape optimization with flutter constraints: AIAA-2016-1176[R]. Reston, VA: AIAA, 2016.

[16] IM D K, CHOI S, MCCLURE J E, et al. Mapped Chebyshev pseudospectral method for unsteady flow analysis[J]. AIAA Journal, 2015, 53(12): 3805-3820.

[17] CHOI J Y, CHOI S, PARK J, et al. Prediction of dynamic Stability using mapped Chebyshev pseudospectral method: AIAA-2016-1347[R]. Reston, VA: AIAA, 2016.

[18] MOIN P. Fundamentals of engineering numerical analysis[M]. Cambridge: Cambridge University Press, 2010.

[19] BAYLISS A, TURKEL E. Mappings and accuracy for Chebyshev pseudo-spectral approximations[J]. Journal of Computational Physics, 1992, 101(2): 349-359.

[20] KOSLOFF D, TAL-EZER H. A modified Chebyshev pseudospectral method with anO(N-1) time step restriction[J]. Journal of Computational Physics, 1993, 104(2): 457-469.

[21] KATZ J, PLOTKINP A. Low-speed aerodynamics[M]. Cambridge: Cambridge University Press, 2001.

[22] MURUA J. Flexible aircraft dynamics with a geometrically-nonlinear description of the unsteady aerodynamics[D]. London: Imperial College London, 2012.

[23] GERADIN M, CARDONA A. Flexible multibody dynamics: A finite element approach[M]. 2001.

[24] GOLAND M. The flutter of a uniform cantilever wing[J]. Journal of Applied Mechanics, 1945, 12(4):197-208.

[25] PATIL M J, HODEGES D H, CESNIK C E S. Nonlinear aeroelastic analysis of complete aircraft in subsonic flow[J]. Journal of Aircraft, 37(5): 753-760.