黃 杰， 姚衛(wèi)星,2， 姜志平， 周丹發(fā)

(1. 南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室 南京，210016)(2. 南京航空航天大學(xué)飛行器先進(jìn)設(shè)計(jì)技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室 南京，210016)(3. 中國(guó)航空工業(yè)集團(tuán)公司沈陽(yáng)飛機(jī)設(shè)計(jì)研究所 沈陽(yáng)，110035)

引 言

氣動(dòng)彈性是飛行器設(shè)計(jì)過程中必須考慮的問題，相比靜氣動(dòng)彈性，顫振對(duì)飛行器帶來(lái)的危害更嚴(yán)重，若設(shè)計(jì)不當(dāng)，飛行動(dòng)壓達(dá)到顫振臨界動(dòng)壓時(shí)，飛行器將發(fā)生機(jī)毀人亡的事故。早期的飛行器顫振分析基于經(jīng)典顫振理論，典型的分析方法有v-g法和p-k法[1-2]。經(jīng)典顫振理論為頻域分析方法，非定常氣動(dòng)力的分析通常采用偶極子網(wǎng)格法[3-4]，該氣動(dòng)力模型為線性模型，只適用于亞聲速和超聲速問題，無(wú)法分析高度非線性的跨聲速問題。當(dāng)攻角較大時(shí)，流動(dòng)分離也會(huì)導(dǎo)致氣動(dòng)力非線性，故偶極子網(wǎng)格法也無(wú)法分析大攻角問題。經(jīng)典顫振理論的結(jié)構(gòu)動(dòng)力學(xué)通常采用模態(tài)疊加法分析，同樣也為線性模型。因此，經(jīng)典顫振理論僅能分析亞聲速或超聲速、小攻角、線性結(jié)構(gòu)等條件下的顫振問題。

對(duì)于大型客機(jī)和高速無(wú)人機(jī)，巡航速度大多處于跨聲速區(qū)域，其氣動(dòng)力是高度非線性的。而大展弦比無(wú)人機(jī)機(jī)翼柔度較大，在氣動(dòng)力作用下會(huì)產(chǎn)生較大的幾何變形，機(jī)翼必然會(huì)產(chǎn)生明顯的幾何非線性效應(yīng)，而幾何非線性會(huì)造成翼面結(jié)構(gòu)在非定常氣動(dòng)力作用下發(fā)生極限環(huán)振蕩(limit cycle oscillation,簡(jiǎn)稱LCO)[5-6]，同時(shí)發(fā)生顫振問題和LCO現(xiàn)象時(shí)可稱為極限環(huán)顫振[7-9]。LCO顫振分析的關(guān)鍵是準(zhǔn)確計(jì)算LCO幅值和頻率。因此，經(jīng)典顫振理論無(wú)法分析這類高度非線性的LCO顫振問題。隨著計(jì)算流體力學(xué)(computation fluid dynamics,簡(jiǎn)稱CFD)和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)(computational structural dynamics,簡(jiǎn)稱CSD)的發(fā)展，一些學(xué)者開始采用基于CFD/CSD的耦合方法進(jìn)行LCO顫振的研究，非定常氣動(dòng)力通過有限體積法求解Euler方程或Navier-Stokes方程獲得，而結(jié)構(gòu)瞬態(tài)響應(yīng)分析通常采用非線性有限元法，因此這是一種時(shí)域方法。目前LCO顫振的研究對(duì)象主要集中在壁板，而機(jī)翼的LCO顫振研究相對(duì)較少。Schairer等[10]通過試驗(yàn)方法研究了切尖三角翼的LCO顫振問題，獲得了三角翼的模態(tài)及不同來(lái)流動(dòng)壓下的翼尖LCO幅值和LCO頻率。文獻(xiàn)[11-13]通過CFD/CSD耦合方法研究了切尖三角翼的LCO顫振，并將計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行了對(duì)比。該耦合方法具有松耦合的特點(diǎn)，耦合時(shí)間精度較低，具有時(shí)間滯后效應(yīng)，已有學(xué)者證明了松耦合方法僅具有一階耦合時(shí)間精度[14]。

針對(duì)松耦合方法的時(shí)間滯后效應(yīng)及耦合時(shí)間精度較低的缺點(diǎn)，筆者發(fā)展了一種跨聲速翼面極限環(huán)顫振的CFD/CSD緊耦合分析方法，其特點(diǎn)是在傳統(tǒng)松耦合方法的基礎(chǔ)上增加了內(nèi)部偽迭代分析過程。采用Schairer等的試驗(yàn)?zāi)Ｐ瓦M(jìn)行了緊耦合和松耦合分析，并將計(jì)算結(jié)果和試驗(yàn)結(jié)果進(jìn)行了對(duì)比，以驗(yàn)證緊耦合方法的耦合時(shí)間精度。

1 氣動(dòng)力和結(jié)構(gòu)動(dòng)力學(xué)控制方程

不考慮體積力和內(nèi)熱源情況下，直角坐標(biāo)系下的流體動(dòng)力學(xué)Navier-Stokes控制方程的積分形式為

(1)

其中：W為守恒向量；Fc為對(duì)流通量；Fv為黏性通量；?V為控制體V的邊界面；n為?V的外法線單位向量。

(2)

(3)

(4)

(5)

其中：ρ為密度；p為壓強(qiáng)；u，v和w分別為直角坐標(biāo)系下的速度分量；qx，qy和qz分別為直角坐標(biāo)系下的熱流分量；e為單位質(zhì)量氣體的總能量，其表達(dá)式為p/((γ-1)ρ)+(u2+v2+w2)/2；γ為氣體比熱比；τij為黏性應(yīng)力分量。

將式(1)按有限體積法進(jìn)行空間離散可得

(6)

其中：Wi和Vi分別為控制體i的守恒向量和體積；NF為控制體邊界面的數(shù)目；ΔSN為第N個(gè)邊界面的面積。

由于對(duì)流通量Fc具有高度非線性特點(diǎn)，并且集中體現(xiàn)了流場(chǎng)的對(duì)流特征，筆者采用Roe格式[15]對(duì)其進(jìn)行空間離散，為獲得單調(diào)解，采用完全迎風(fēng)的二階MUSCL格式[16]離散分裂后的無(wú)黏通量,并采用minmod限制器使空間離散格式達(dá)到空間二階精度，且非定常問題的求解采用雙時(shí)間步長(zhǎng)法[17]。

以上是氣動(dòng)力分析的控制方程和數(shù)值計(jì)算方法，而通過有限元法離散后結(jié)構(gòu)動(dòng)力學(xué)運(yùn)動(dòng)方程可表示為

(7)

本研究的機(jī)翼極限環(huán)非線性顫振僅考慮結(jié)構(gòu)的幾何非線性，不考慮其材料非線性，且采用Newmark法進(jìn)行結(jié)構(gòu)非線性瞬態(tài)動(dòng)力學(xué)的求解，其屬于隱式求解方法。將式(7)寫為結(jié)構(gòu)非線性動(dòng)力學(xué)的平衡方程為

(8)

Fi體現(xiàn)了結(jié)構(gòu)的幾何非線性，需要將其線性化以方便求解

Fi,n+1=Fi,n+KT,nΔu

(9)

其中：Δu=un+1-un。

將式(9)帶入式(8)，即可得到線性化的結(jié)構(gòu)動(dòng)力學(xué)的平衡方程

(10)

由于式(9)在一定程度上屬于一階線性化，故一階線性化的式(10)求解精度較差，可引入子迭代Newmark方法來(lái)求解式(10)，其平衡方程為

(11)

基于子迭代Newmark方法的結(jié)構(gòu)瞬態(tài)非線性動(dòng)力學(xué)求解的詳細(xì)過程可參考文獻(xiàn)[14]。

2 全隱式緊耦合分析策略

本研究非定常氣動(dòng)力采用CFD方法求解，而結(jié)構(gòu)非線性動(dòng)力學(xué)采用基于有限元的CSD方法求解，且CFD和CSD均采用隱式時(shí)間推進(jìn)，相比顯示格式，隱式格式能采用更大的時(shí)間步長(zhǎng)且具有更好的穩(wěn)定性。

傳統(tǒng)的機(jī)翼非線性顫振分析采用基于CFD/CSD的串行松耦合方法，其在相同的時(shí)間步進(jìn)行氣動(dòng)力和節(jié)點(diǎn)位移的數(shù)據(jù)交換，并進(jìn)行迭代推進(jìn)求解。松耦合方法簡(jiǎn)單明了，實(shí)現(xiàn)起來(lái)較容易，但在任意n到n+1時(shí)間步內(nèi)CFD(或CSD)求解過程中節(jié)點(diǎn)位移(或氣動(dòng)力)不變，即凍結(jié)邊界條件，這樣會(huì)造成時(shí)間滯后的效應(yīng)，且隨著分析的進(jìn)行，累積誤差會(huì)越來(lái)越大。文獻(xiàn)[14]已經(jīng)證明了即使CFD和CSD分析均達(dá)到了二階時(shí)間精度，由于耦合邊界上氣動(dòng)力和節(jié)點(diǎn)位移在交換時(shí)間上的滯后，傳統(tǒng)松耦合方法只能達(dá)到一階時(shí)間精度。這明顯不滿足機(jī)翼極限環(huán)顫振對(duì)耦合方法時(shí)間精度的要求。

基于傳統(tǒng)松耦合方法時(shí)間精度較低的缺點(diǎn)，筆者采用全隱式緊耦合方法進(jìn)行機(jī)翼極限環(huán)顫振的研究。其在任意一個(gè)時(shí)間步內(nèi)CFD與CSD進(jìn)行反復(fù)的內(nèi)迭代，通常將其稱之為偽迭代。當(dāng)流體與結(jié)構(gòu)滿足精度要求，跳出偽迭代進(jìn)行下一個(gè)時(shí)間步的求解，即進(jìn)入物理迭代的分析。相比松耦合方法，緊耦合方法增加了偽迭代，并要求偽迭代收斂后才進(jìn)行物理迭代的計(jì)算，這在一定程度上能消除松耦合方法在時(shí)間推進(jìn)過程中累積的誤差，即降低松耦合方法存在的時(shí)間滯后效應(yīng)，具有更高的耦合時(shí)間分析精度。由于偽迭代完全收斂后再進(jìn)入物理迭代，造成了計(jì)算量大大增加。因此在實(shí)際的緊耦合分析中通常需要設(shè)置一個(gè)最大偽迭代次數(shù)，即當(dāng)偽迭代達(dá)到最大次數(shù)時(shí)即使其還未完全收斂，也跳出偽迭代進(jìn)入物理迭代。

如圖1所示，緊耦合方法的主要分析步驟如下：

圖1 緊耦合方法分析流程Fig.1 Analysis process of tightly coupled method

1) 進(jìn)行CFD和CSD建模，并進(jìn)行CFD定常流場(chǎng)的計(jì)算，計(jì)算得到的定常氣動(dòng)力結(jié)果作為耦合分析的初始條件；

2) 進(jìn)行ti時(shí)刻的非定常氣動(dòng)力的計(jì)算，并將壁面壓力熱Fi傳遞給CSD模型，進(jìn)行CFD與CSD的偽迭代分析，直到偽迭代達(dá)到最大迭代次數(shù)后結(jié)束該時(shí)間步內(nèi)的偽迭代分析；

3) 進(jìn)行下一個(gè)時(shí)間步的分析，即進(jìn)入物理迭代的計(jì)算；

4) 若物理迭代分析時(shí)間ti到達(dá)分析總時(shí)間ttotal，結(jié)束分析，否者返回步驟2，直至ti=ttotal為止。

以上耦合流程中涉及到氣動(dòng)力和節(jié)點(diǎn)位移在耦合面上的數(shù)據(jù)交換，而CSD的網(wǎng)格尺寸通常遠(yuǎn)大于CFD網(wǎng)格尺寸，故在耦合面上不能直接通過節(jié)點(diǎn)一一對(duì)應(yīng)的方式進(jìn)行以上耦合變量的傳遞，需要采用插值算法完成數(shù)據(jù)的傳遞。筆者采用三維薄板樣條插值方法進(jìn)行氣動(dòng)力和節(jié)點(diǎn)位移的插值計(jì)算。此外，結(jié)構(gòu)變形會(huì)導(dǎo)致流場(chǎng)網(wǎng)格的重新調(diào)整，筆者采用基于無(wú)限插值方法的動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)流體網(wǎng)格的變形，其詳細(xì)計(jì)算方法可參考文獻(xiàn)[18]。

3 跨聲速切尖三角翼極限環(huán)顫振分析

3.1 數(shù)值模型

筆者采用跨聲速切尖三角翼試驗(yàn)?zāi)Ｐ蚚11]進(jìn)行研究，通過試驗(yàn)可知，該模型的非線性因素來(lái)自于氣動(dòng)載荷作用下的結(jié)構(gòu)大變形和跨聲速非線性氣動(dòng)力。此剪切三角翼為等厚度薄板，其厚度為0.889mm，材料為鋼，彈性模量為200GPa，泊松比為0.25，密度為7 850kg/m3，其平面形狀和幾何尺寸如圖2(a)所示，且翼根為固支邊界條件。來(lái)流馬赫數(shù)在0.86和0.879之間，攻角為0°。為了與試驗(yàn)來(lái)流參數(shù)相匹配，來(lái)流動(dòng)壓的單位采用psi，其與Pa的變換關(guān)系為1psi=6 895Pa。圖2為CFD和CSD數(shù)值模型，由于采用松耦合方法分析該模型極限環(huán)顫振的文獻(xiàn)中未考慮流場(chǎng)的黏性效應(yīng)，為了對(duì)比分析，本研究CFD模型僅求解Euler方程。此外CSD模型采用殼單元模擬，且考慮結(jié)構(gòu)大變形帶來(lái)的幾何非線性效應(yīng)。

圖2 流體和結(jié)構(gòu)模型網(wǎng)格Fig.2 Fluid and structural meshes

筆者進(jìn)行了結(jié)構(gòu)模態(tài)分析，結(jié)構(gòu)振型的試驗(yàn)和分析結(jié)果如圖3所示，機(jī)翼前3階振型分別為一階彎曲、一階扭轉(zhuǎn)和二階彎曲，且分析結(jié)果與試驗(yàn)結(jié)果吻合得很好。表1列出了機(jī)翼前3階固有頻率的試驗(yàn)和分析結(jié)果[11]，固有頻率的計(jì)算值與試驗(yàn)值最大相對(duì)誤差為1.82%，這說(shuō)明本研究的CSD模型能夠反映切尖三角翼的基本動(dòng)力學(xué)性能。

圖3 機(jī)翼振型對(duì)比Fig.3 Comparison of modal shape of wing

表1 機(jī)翼固有頻率的對(duì)比

3.2 機(jī)翼幾何非線性的影響

為了研究結(jié)構(gòu)幾何非線性對(duì)顫振的影響，筆者采用緊耦合方法同時(shí)分析了來(lái)流動(dòng)壓q=2.58psi下考慮和不考慮幾何非線性時(shí)的翼尖瞬態(tài)位移響應(yīng)，其中耦合時(shí)間步長(zhǎng)Δt取0.000 4s，每個(gè)時(shí)間步內(nèi)的最大為迭代為10。圖4為線性結(jié)構(gòu)(不考慮結(jié)構(gòu)幾何非線性)和非線性結(jié)構(gòu)下的結(jié)果對(duì)比圖，可知考慮結(jié)構(gòu)幾何非線性時(shí)機(jī)翼在初始擾動(dòng)下振動(dòng)幅度逐漸增加，并且在最終呈等幅振蕩狀態(tài)。這是因?yàn)橐砻鏋閯偠葷u硬的系統(tǒng)，即隨著翼面變形的增加，翼面幾何剛度也逐漸增加，機(jī)翼抵抗變形的能力得到了加強(qiáng)。線性剛度的機(jī)翼在分析過程中剛度保持不變，當(dāng)受初始擾動(dòng)時(shí)，機(jī)翼在初期呈現(xiàn)較小幅度的振動(dòng)，隨后進(jìn)入劇烈的動(dòng)態(tài)發(fā)散階段，此階段機(jī)翼的變形遠(yuǎn)大于考慮幾何非線性情況下的結(jié)果。據(jù)此可知機(jī)翼幾何非線性對(duì)時(shí)域分析結(jié)果產(chǎn)生了重要影響，必須考慮其影響才能獲得正確的顫振分析結(jié)果。

圖4 翼尖線性和非線性位移響應(yīng)(動(dòng)壓q=2.58psi)Fig.4 Linear and nonlinear displacement responses of wing tip(q=2.58psi)

3.3 非線性臨界顫振動(dòng)壓

筆者采用緊耦合方法進(jìn)行了切尖三角翼的臨界顫振動(dòng)壓分析，同時(shí)考慮了非線性跨聲速氣動(dòng)力和結(jié)構(gòu)幾何非線性。圖5為不同動(dòng)壓下機(jī)翼翼尖的法向位移時(shí)域分析結(jié)果，根據(jù)結(jié)果可知：動(dòng)壓q為2.42psi時(shí)，在外界擾動(dòng)下翼尖的振動(dòng)呈收斂的趨勢(shì)；動(dòng)壓q為2.44psi時(shí)，翼尖的振動(dòng)逐漸發(fā)散；動(dòng)壓q為2.43psi時(shí)，翼尖呈近似等幅振動(dòng)的狀態(tài)，故本研究的切尖三角翼模型非線性臨界顫振動(dòng)壓為2.43psi。

圖5 不同動(dòng)壓下翼尖位移響應(yīng)Fig.5 Displacement responses of wing tip under different dynamic pressure

表2為p-k法[13]與本研究緊耦合方法計(jì)算得到的臨界顫振動(dòng)壓與試驗(yàn)結(jié)果的對(duì)比情況，可知緊耦合方法獲得的臨界顫振動(dòng)壓值2.43psi與試驗(yàn)值2.40psi吻合良好，而p-k法計(jì)算得到的臨界顫振動(dòng)壓為2.75psi，結(jié)果較差。這是因?yàn)閜-k法為工程算法，其氣動(dòng)力和結(jié)構(gòu)分析均為線性，無(wú)法考慮切尖三角翼的非線性跨聲速氣動(dòng)力和結(jié)構(gòu)幾何非線性。表2中的數(shù)據(jù)同樣也驗(yàn)證了緊耦合方法在求解非線性臨界顫振動(dòng)壓中的高精度。

表2 非線性顫振動(dòng)壓比較

3.4 LCO分析結(jié)果和討論

在切尖三角翼的顫振試驗(yàn)中觀察到了LCO現(xiàn)象，其中的非線性來(lái)自于非線性跨聲速氣動(dòng)力和結(jié)構(gòu)幾何非線性，且結(jié)構(gòu)幾何非線性是主要影響因素。筆者分析了不同來(lái)流動(dòng)壓下的翼尖法向位移的時(shí)域結(jié)果和相應(yīng)的相位圖，如圖6所示。從圖中結(jié)構(gòu)可知，在初始擾動(dòng)下翼尖振動(dòng)在初始階段呈現(xiàn)輕微的發(fā)散趨勢(shì)，隨后翼尖振動(dòng)進(jìn)入劇烈發(fā)散階段，翼尖位

圖6 不同動(dòng)壓下翼尖位移響應(yīng)及相位圖Fig.6 Displacement responses and phase diagrams of wing tip under different dynamic pressure

移迅速增加，并且最終作等幅振蕩的狀態(tài)，即表現(xiàn)出LCO現(xiàn)象。此外隨著來(lái)流動(dòng)壓的增加，翼尖會(huì)更劇烈且更快地進(jìn)入LCO狀態(tài)。

分析獲得了不同來(lái)流動(dòng)壓下的翼尖LCO幅值和頻率，如圖7所示，圖7(a)中的縱坐標(biāo)為L(zhǎng)CD與機(jī)翼半展長(zhǎng)b(其值為203.2mm)之比。為了對(duì)比分析，圖中還給出了一些采用松耦合方法獲得的結(jié)果。為了保證對(duì)比分析的合理性，本研究緊耦合方法和文獻(xiàn)中松耦合方法的氣動(dòng)力分析控制方程均為Euler方程，且耦合時(shí)間步長(zhǎng)均為0.000 4s。由圖可知，緊耦合和松耦合方法獲得的LCO幅值在較低來(lái)流動(dòng)壓下吻合較好，隨著來(lái)流動(dòng)壓的增加，計(jì)算值與試驗(yàn)值誤差越來(lái)越大，尤其在來(lái)流動(dòng)壓q=3.45psi時(shí)LCO幅值的試驗(yàn)值遠(yuǎn)大于所有計(jì)算值，但可以明顯觀察到采用緊耦合方法獲得的LCO幅值更靠近試驗(yàn)值。隨著來(lái)流動(dòng)壓的增加，LCO頻率的試驗(yàn)值逐漸增加，緊耦合方法計(jì)算得到的LCO頻率能清晰地反映這一規(guī)律，而松耦合方法卻不能。此外，緊耦合方法獲得的LCO頻率也與試驗(yàn)值吻合得更好。

圖7 不同動(dòng)壓下翼尖LCO幅值和頻率Fig.7 LCO amplitude and frequency of wing tip under different dynamic pressure

以上分析說(shuō)明了緊耦合方法比傳統(tǒng)松耦合方法能更好地分析LCO幅值和頻率，具有更高的分析精度。但是在較高來(lái)流動(dòng)壓下，緊耦合方法獲得的LCO幅值及頻率與試驗(yàn)值還存在較大的誤差。主要原因有：①三角翼的扭轉(zhuǎn)變形會(huì)造成局部大攻角，三角翼在大攻角下會(huì)形成前緣分離渦和前緣低壓區(qū)，如圖8所示，這難以通過Euler方程進(jìn)行數(shù)值模擬，即存在氣動(dòng)力不精確的問題；②試驗(yàn)過程中機(jī)翼的振動(dòng)容易引起翼根的松動(dòng)[11]，降低約束剛度，造成振動(dòng)幅值增加；③翼根部為高應(yīng)力區(qū)域，有部分結(jié)構(gòu)已進(jìn)入塑性階段，而本研究的分析不考慮材料非線性，這也將引起振幅增加。

圖8 翼面壓力系數(shù)分布Fig.8 Pressure coefficient distribution of wing

4 結(jié) 論

1) 傳統(tǒng)松耦合方法分析過程中存在嚴(yán)重的時(shí)間滯后問題，會(huì)影響耦合時(shí)間精度。筆者發(fā)展了一種機(jī)翼極限環(huán)顫振的CFD/CSD全隱式緊耦合研究方法，其特點(diǎn)是在傳統(tǒng)松耦合方法的基礎(chǔ)上增加了內(nèi)部偽迭代分析過程。

2) 進(jìn)行了切尖三角翼跨聲速極限環(huán)顫振的分析，結(jié)果表明結(jié)構(gòu)幾何非線性對(duì)LCO顫振具有重要影響，且緊耦合方法獲得的機(jī)翼臨界顫振動(dòng)壓與試驗(yàn)結(jié)果吻合得很好，從而驗(yàn)證了緊耦合方法的分析精度。

3) 緊耦合方法獲得的翼尖LCO幅值和LCO頻率均優(yōu)于傳統(tǒng)松耦合方法，更靠近試驗(yàn)結(jié)果，故緊耦合方法在一定程度上能消除時(shí)間推進(jìn)累積的誤差，具有更高的耦合時(shí)間精度。

4) 本研究的緊耦合方法不僅能應(yīng)用于機(jī)翼極限環(huán)顫振的分析，還可應(yīng)用于對(duì)時(shí)間精度要求較高的機(jī)翼突風(fēng)響應(yīng)以及其他類型的瞬態(tài)耦合分析中，如氣動(dòng)熱與結(jié)構(gòu)傳熱的瞬態(tài)耦合分析，因此緊耦合方法具有廣泛地工程應(yīng)用價(jià)值。