楊體浩，白俊強(qiáng)，史亞云，楊一雄

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

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

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

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

針對(duì)非定常問題，諧波平衡(Harmonic Balance, HB)[7]是另一種很有吸引力的方法。HB通過將系統(tǒng)控制方程的狀態(tài)變量用截?cái)嗟母道锶~系數(shù)進(jìn)行替換，最終將非定常的控制方程轉(zhuǎn)化為耦合的定?？刂品匠?。HB直接對(duì)系統(tǒng)控制方程進(jìn)行處理，不需要樣本快照。因此，HB具有很高的計(jì)算精度和魯棒性。借助HB方法，可以準(zhǔn)確地對(duì)整個(gè)系統(tǒng)進(jìn)行靈敏度分析，并且很容易與伴隨方法耦合對(duì)非定常問題進(jìn)行優(yōu)化設(shè)計(jì)，大大地提高優(yōu)化設(shè)計(jì)效率。Choi等[8]將HB與伴隨方相結(jié)合建立了非定常氣動(dòng)伴隨優(yōu)化設(shè)計(jì)系統(tǒng)，并以此對(duì)直升機(jī)旋翼進(jìn)行了氣動(dòng)外形優(yōu)化設(shè)計(jì)。但是，HB方法只適用于周期性的非定常問題，具有一定的局限性，而很多實(shí)際的流固耦合問題都是非周期的。

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

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

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

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

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

1 Chebyshev譜方法

1.1 標(biāo)準(zhǔn)的Chebyshev譜方法

Chebyshev譜方法是一種以Chebyshev多項(xiàng)式插值為理論基礎(chǔ)的時(shí)間譜方法[18]。理論上，任意一個(gè)定義在時(shí)間區(qū)間[-1,1]內(nèi)的光滑函數(shù)u(t)，在tm時(shí)刻的函數(shù)值可以表示為

(1)

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

(2)

(3)

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

進(jìn)一步，函數(shù)u(t)在Chebyshev控制點(diǎn)處的時(shí)間導(dǎo)數(shù)可以推導(dǎo)為

(4)

(5)

(6)

(7)

1.2 Chebyshev算子的投影變換

Chebyshev控制點(diǎn)在時(shí)間區(qū)間內(nèi)的分布情況會(huì)對(duì)其計(jì)算結(jié)果的精度產(chǎn)生一定的影響。標(biāo)準(zhǔn)的Chebyshev譜方法的控制點(diǎn)在時(shí)間區(qū)間的兩端分布較為密集，這種分布形態(tài)使得標(biāo)準(zhǔn)的Chebyshev譜方法需要更多的控制點(diǎn)個(gè)數(shù)去達(dá)到期望的計(jì)算精度[19]。因此，采用Kosloff和Tal-Ezer[20]提出的投影變換函數(shù)改變Chebyshev譜方法的控制點(diǎn)分布，使其分布更為均勻。投影變換的表達(dá)式為

(8)

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

(9)

(10)

1.3 驗(yàn)證算例

經(jīng)過投影變換的Chebyshev算子被稱為投影的Chebyshev算子(Mapped Chebyshev Operator，MCO)。利用非周期的三角函數(shù)驗(yàn)證投影的Chebyshev算子的計(jì)算精度，揭示投影的Chebyshev譜方法的數(shù)學(xué)特性。三角函數(shù)的表達(dá)式為

(11)

(12)

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

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

圖2 投影的Chebyshev控制點(diǎn)分布Fig.2 Distribution of mapped Chebyshev control points

計(jì)算結(jié)果表明，即便對(duì)于非周期函數(shù)，Chebyshev算子也可以用較少的控制點(diǎn)的個(gè)數(shù)得到滿足精度要求的計(jì)算結(jié)果。

2 基于投影的Chebyshev譜方法的非定常氣動(dòng)力模型

2.1 非定常面元法

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

(13)

(14)

(15)

(16)

圖3 尾跡模型Fig.3 Wake model

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

(17)

式中：Cpk為第k個(gè)物面面元控制點(diǎn)處的壓力系數(shù)；pk和pref分別為作用在第k個(gè)物面面元控制點(diǎn)處的靜壓和無窮遠(yuǎn)處自由來流的靜壓；ρ為大氣密度；Qk為第k個(gè)物面面元控制點(diǎn)處的當(dāng)?shù)亓黧w速度；Φk為第k個(gè)物面面元控制點(diǎn)處的速度勢(shì)。

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

根據(jù)Chebyshev譜方法的理論可知，將Chebyshev譜方法與非定常面元法結(jié)合以后，只有在Chebyshev控制點(diǎn)處的狀態(tài)變量為待求未知變量。顯然，在給定的時(shí)間區(qū)內(nèi)的任意時(shí)刻，物面的運(yùn)動(dòng)速度和位置形狀可表示為

(18)

(19)

根據(jù)尾跡形狀的遞推關(guān)系式(式(16))可知，在ti=t0+ndt時(shí)刻，尾跡的形狀為

μwake,μsurface,σsurface,Vgust)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

顯然，當(dāng)初始條件給定后，通過迭代求解上述控制方程可以獲得在選定的Chebyshev控制點(diǎn)處的面元的強(qiáng)度。進(jìn)而借助Chebyshev算子可以計(jì)算得到定義的時(shí)間區(qū)間內(nèi)任一時(shí)刻下的面元強(qiáng)度。聯(lián)立式(10)與式(17)可計(jì)算得到對(duì)應(yīng)狀態(tài)下所有物面面元的壓力系數(shù)：

(27)

2.3 驗(yàn)證算例

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

(28)

式中:迎角α0=0°;運(yùn)動(dòng)周期T=0.4 s。等直機(jī)翼以過50%弦點(diǎn)的平行于前緣的直線為軸做俯仰運(yùn)動(dòng)。

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

圖4 機(jī)翼表面面元Fig.4 Surface panels of wing

圖5 Chebyshev譜方法與時(shí)間推進(jìn)方法計(jì)算得到的非定常氣動(dòng)力對(duì)比Fig.5 Comparison of unsteady aerodynamic forces calculated by Chebyshev spectral method and time-marching method

3 基于Chebyshev譜方法的結(jié)構(gòu)動(dòng)力學(xué)模型

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

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

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

圖6 坐標(biāo)系的定義Fig.6 Definition of coordinate systems

t時(shí)刻梁模型上任意一點(diǎn)，相對(duì)于未變形構(gòu)型(即t=0時(shí)刻的構(gòu)型)的應(yīng)變可表示為

(29)

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

根據(jù)Hamilton原理有：

(30)

(31)

(32)

(33)

式中：M、C和K分別為切向質(zhì)量矩陣、切向阻尼矩陣和切向剛度矩陣?？紤]幾何非線性時(shí)，M、C和K不是常值矩陣，而是與狀態(tài)變量相關(guān)的函數(shù)。基于增量形式的結(jié)構(gòu)動(dòng)力學(xué)方程，利用相應(yīng)的數(shù)值求解方法(如Newton-Raphson方法)進(jìn)行迭代求解，即可獲得梁模型的結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)。

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

幾何非線性梁有限元模型可以寫成結(jié)構(gòu)動(dòng)力學(xué)方程的一般形式：

(34)

通過等價(jià)變換將二階微分方程轉(zhuǎn)化為一階微分方程：

(35)

(36)

(37)

(38)

(39)

顯然，借助Chebyshev譜方法，非定常形式的結(jié)構(gòu)動(dòng)力學(xué)模型被轉(zhuǎn)化為在選定的Chebyshev控制點(diǎn)處的定常方程。當(dāng)給定初始條件后，通過迭代求解控制方程可以獲得結(jié)構(gòu)狀態(tài)變量在所有Chebyshev控制點(diǎn)處的值。之后，利用Chebyshev算子可以計(jì)算得到在給定時(shí)間區(qū)間內(nèi)的結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)。

3.3 驗(yàn)證算例

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

圖7 施加強(qiáng)迫載荷的懸臂梁Fig.7 Cantilever beam with forced load

表1 梁結(jié)構(gòu)屬性Table 1 Properties of beam structure

圖8 強(qiáng)迫載荷的時(shí)間歷程Fig.8 Time history of forced load

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

圖9 Chebyshev譜方法與時(shí)間推進(jìn)方法計(jì)算得到的結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)對(duì)比Fig.9 Comparison of structural dynamic response calculated by Chebyshev spectral method and time-marching method

4 基于Chebyshev譜方法的流固耦合分析系統(tǒng)

利用松耦合迭代策略求解基于Chebyshev譜方法的流固耦合控制方程。由于Chebyshev譜方法分別與非定常氣動(dòng)力模型以及結(jié)構(gòu)動(dòng)力學(xué)模型相結(jié)合，將非定常的流固耦合問題轉(zhuǎn)化為在選定的Chebyshev控制點(diǎn)處耦合的定常問題。因此，基于Chebyshev譜方法的動(dòng)氣動(dòng)彈性問題的求解過程，類似于傳統(tǒng)采用松耦合策略的靜氣動(dòng)彈性問題的求解，只不過所有Chebyshev控制點(diǎn)處的控制方程是同時(shí)求解的。這種不同Chebyshev控制點(diǎn)處，控制方程的耦合實(shí)際上是整個(gè)時(shí)間響應(yīng)歷程內(nèi)不同Chebyshev控制點(diǎn)處的狀態(tài)變量之間時(shí)間依賴性的直接體現(xiàn)。

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

5 數(shù)值仿真與討論

對(duì)標(biāo)準(zhǔn)模型Goland機(jī)翼[24]進(jìn)行顫振速度的計(jì)算，以驗(yàn)證本文搭建的基于Chebyshev譜方法的流固耦合計(jì)算、分析方法的收斂特性及計(jì)算精度。Goland機(jī)翼半展長(zhǎng)6.096 m，弦長(zhǎng)1.828 8 m。采用國(guó)際上針對(duì)Goland機(jī)翼計(jì)算分析使用的梁有限元相關(guān)參數(shù)[24]，具體的結(jié)構(gòu)屬性如表2所示。圖11為Goland機(jī)翼表面面元的劃分。

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

表2 Goland機(jī)翼屬性Table 2 Properties of Goland wing

圖11 Goland機(jī)翼表面面元Fig.11 Surface panels of Goland wing

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

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

圖12 不同自由來流速度下的Goland機(jī)翼翼尖時(shí)間歷程響應(yīng)對(duì)比Fig.12 Comparison of time history of Goland wing tip at various free-stream velocities

表3 Goland機(jī)翼顫振速度對(duì)比Table 3 Comparison of flutter speed of Goland wing

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

6 結(jié) 論

1) Chebyshev譜方法利用Chebyshev算子替換系統(tǒng)的狀態(tài)變量及其導(dǎo)數(shù)，將非定常問題轉(zhuǎn)化為選取的Chebyshev控制點(diǎn)處耦合的定常問題進(jìn)行處理。由于Chebyshev譜方法直接對(duì)流固耦合的控制方程進(jìn)行處理，因此所建立的基于Chebyshev譜方法的流固耦合計(jì)算、分析方法具有較強(qiáng)的魯棒性。

2) Chebyshev譜方法的計(jì)算精度隨著Chebyshev控制點(diǎn)個(gè)數(shù)的增加而不斷提高，測(cè)試算例表明，只需少許有限的Chebyshev控制點(diǎn)，Chebyshev譜方法便可以達(dá)到期望的計(jì)算精度。

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

4) 建立的基于Chebyshev譜方法的流固耦合分析方法能夠準(zhǔn)確地捕捉流體與結(jié)構(gòu)之間的耦合作用，可用于忽略了氣動(dòng)非線性的流固耦合問題的計(jì)算、分析與設(shè)計(jì)研究，如經(jīng)典的顫振問題。

參 考 文 獻(xiàn)

[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.