基于氣動(dòng)力分段表達(dá)的顫振求解方法與模態(tài)跟蹤技術(shù)的應(yīng)用

任智毅，金海波，丁運(yùn)亮

（1.上海航天技術(shù)研究院 第八設(shè)計(jì)部，上海 201109；2.南京航空航天大學(xué) 航空宇航學(xué)院，南京 210016）

0 引 言

顫振是一種典型的氣動(dòng)彈性動(dòng)不穩(wěn)定現(xiàn)象，氣動(dòng)彈性的運(yùn)動(dòng)方程通過振型線性化假設(shè)、模態(tài)降階等理論將其轉(zhuǎn)換成質(zhì)量歸一化的廣義結(jié)構(gòu)運(yùn)動(dòng)方程，通過求解該方程的非線性二次特征值問題來確定結(jié)構(gòu)顫振的臨界穩(wěn)定條件。目前大多數(shù)顫振分析方法求解其特征根都需要在指定的折合頻率區(qū)間上進(jìn)行廣義氣動(dòng)力的插值和不斷迭代兩個(gè)過程［1－2］。由于采用迭代的方法求解顫振方程不僅增加計(jì)算時(shí)間而且可能存在收斂性問題，因此Goodman［3］提出了一種不需要迭代的顫振求解方法即采用n個(gè)折合頻率將折合頻率區(qū)間分成n－1段，并通過分段插值函數(shù)（PQI）來近似表達(dá)廣義氣動(dòng)力與折合頻率之間的函數(shù)關(guān)系，從而將顫振方程的求解轉(zhuǎn)化為分段線性二次特征值問題。由于Goodman在分段邊界上只進(jìn)行了前后分段插值函數(shù)的廣義氣動(dòng)力匹配，沒有保證段與段之間插值函數(shù)一階導(dǎo)數(shù)的連續(xù)性，導(dǎo)致特征根在分段邊界上出現(xiàn)間斷和重疊現(xiàn)象，從而發(fā)生模態(tài)竄支。Eller［4］在此基礎(chǔ)上采用n－1個(gè)折合頻率將折合頻率區(qū)間分成n－2段，在分段邊界上增加了前后分段插值函數(shù)的一階導(dǎo)數(shù)連續(xù)性條件，從而保證了顫振方程特征根的連續(xù)性，因此Eller認(rèn)為該方法在顫振計(jì)算過程中不會(huì)出現(xiàn)模態(tài)竄支現(xiàn)象，故不需要模態(tài)跟蹤。但是由于分段插值函數(shù)的構(gòu)建依賴于折合頻率的選擇，將構(gòu)造的插值函數(shù)嵌入顫振方程計(jì)算過程中仍然會(huì)產(chǎn)生模態(tài)竄支現(xiàn)象，其原因主要是在顫振方程分段求解中，在同一插值段內(nèi)可能獲得其中兩支或者更多模態(tài)的特征根，由于在該段顫振方程特征值的計(jì)算中這些特征根是沒有順序的，從而出現(xiàn)各階模態(tài)對(duì)應(yīng)的特征值隨參數(shù)變化的前后對(duì)應(yīng)關(guān)系混亂即“模態(tài)交叉”現(xiàn)象，如果不采用有效的模態(tài)跟蹤技術(shù)來跟蹤這些特征根，這將給顫振模態(tài)識(shí)別帶來困難。

對(duì)于模態(tài)跟蹤技術(shù)，目前主要采用相似排序法［5］、正交檢驗(yàn)法［6］、預(yù)測(cè)跟蹤法［7］等。由于相似排序法和正交檢驗(yàn)法需要較小的步長(zhǎng)，其參數(shù)變化前后特征值的相似性和特征向量的正交性條件才近似成立，且對(duì)于復(fù)雜的系統(tǒng)可能出現(xiàn)模態(tài)跟蹤失敗現(xiàn)象，因此本文針對(duì)顫振求解PQI法，發(fā)展了一種基于特征值攝動(dòng)理論的模態(tài)跟蹤技術(shù)－預(yù)測(cè)跟蹤法，它利用特征值及左、右特征向量信息求解下一參數(shù)點(diǎn)的特征值估計(jì)量，以估計(jì)量為參考來跟蹤模態(tài)的變化情況，解決了其存在的模態(tài)竄支問題。

本文介紹了帶有模態(tài)跟蹤的PQI顫振分析的基本方法，并用兩個(gè)標(biāo)準(zhǔn)算例驗(yàn)證了該方法的計(jì)算精度和有效性。

1 顫振方程的直接求解方法

1.1 顫振的運(yùn)動(dòng)方程

在小變形和線彈性材料假設(shè)前提下，通過拉普拉斯變換將顫振分析的運(yùn)動(dòng)方程寫成式（1）［8］：

式中是廣義質(zhì)量矩陣，是廣義阻尼矩陣，V是來流速度，L是參考弦長(zhǎng)的一半，ρ是大氣密度；Qs（P）是拉普拉斯域廣義氣動(dòng)力矩陣，P是拉普拉斯參數(shù)，P＝g＋ik，g是阻尼，k為折合頻率，k＝ωL／V，ω為結(jié)構(gòu)的固有頻率。

在常規(guī)的顫振方程（1）的求解中，由于氣動(dòng)力Qs（P）矩陣是關(guān)于P的非線性非對(duì)稱矩陣，所以需要采用一種迭代的方法求解非線性矩陣特征值問題的特征根，一般的做法是給定，計(jì)算其對(duì)應(yīng)的氣動(dòng)力Q（），并結(jié)合尋優(yōu)的過程來尋找線性特征值問題的解中與相等的特征根P即為非線性特征值問題式（1）的解。

本文擬采用分段表達(dá)的方式把氣動(dòng)力Qs（P）矩陣表示為關(guān)于p的二次表達(dá)形式，如式（2）：

式中A、B、C是多項(xiàng)式的非對(duì)稱系數(shù)矩陣，j是分段數(shù)，b為折合頻率區(qū)間分段斷點(diǎn)。

采用這樣的表達(dá)方式，其優(yōu)點(diǎn)在于：① 采用這樣的表達(dá)將顫振方程（1）的非線性特征值問題轉(zhuǎn)化為分段線性特征值問題；② 由于式（2）滿足柯西－黎曼方程（3）［9］，因此整個(gè)復(fù)平面的分段二次氣動(dòng)力插值函數(shù)可以僅僅沿著虛軸來構(gòu)建，同時(shí)可以看出采用二次插值函數(shù)對(duì)阻尼的變化進(jìn)行了二階近似，由于氣動(dòng)力是在復(fù)平面上構(gòu)造的，因此其氣動(dòng)阻尼的計(jì)算更為精確。

雖然采用二次形式可以解決復(fù)平面氣動(dòng)力矩陣的近似表達(dá)問題，但由于顫振計(jì)算中非定常氣動(dòng)力的非線性特性，因此本文采用分段形式來提高近似表達(dá)精度，在分段表達(dá)中本文采用的分段方法為：首先在折合頻率求解空間構(gòu)造一系列離散點(diǎn)ki（i∈［1，n］），并在該離散點(diǎn)上計(jì)算所對(duì)應(yīng)的廣義氣動(dòng)力Q（iki），然后分段插值函數(shù)的斷點(diǎn)位置bj（j∈［1，n－1］）按如下原則構(gòu)造：

對(duì)于每一折合頻率分段j（k∈［bj，bj＋1］），滿足三個(gè)條件：如式（7）的斷點(diǎn)連續(xù)性；式（8）的斷點(diǎn)符合性；式（9）的斷點(diǎn)處分段函數(shù)一階導(dǎo)數(shù)連續(xù)性：

其中第一段與最后一段條件需要做適當(dāng)?shù)男薷?，即在第一段的三個(gè)條件中將式（7）改為式（10），最后一段中將式（9）改為式（11）：

這樣處理使得每一折合頻率分段都滿足三個(gè)條件，由于段與段之間的條件是相互聯(lián)系的，通過分塊組集的方法將其整合成矩陣的形式，通過矩陣變換來直接求解分段系數(shù)矩陣（Ah×h）j、（Bh×h）j、（Ch×h）j，h表示參與顫振分析的模態(tài)數(shù)量。以第一段和最后一段的系數(shù)求解為例：

式中（Ats）j、（Bts）j、（Cts）j表示第j段三個(gè)系數(shù)矩陣中的t行s列元素（t，s∈［1，h］），（Qts）kj表示折合頻率ki所對(duì)應(yīng)的廣義氣動(dòng)力的Q（iki）的t行s列元素（i∈［1，n］）。

1.2 顫振方程的求解

將式（1）與式（2）整合在一起，顫振方程可改寫為式（13）：

式（13）的求解可以通過引入狀態(tài)變量的方法［10］進(jìn)行求解：

由于式（13）中氣動(dòng)力Qj（P）是分段近似表達(dá)的，對(duì)于第j段，按式（14）求解特征值P。由于氣動(dòng)力僅在j段內(nèi)有效，因此需要判斷特征根是否在該求解段上即條件Im（p）∈［bj，bj＋1］。如果在，其求解有效，如果不在，則需要用下一求解段繼續(xù)求解，通過逐段求解的方式可以找出式（1）在給定V下所對(duì)應(yīng)的參與顫振分析的各階模態(tài)的特征值。但是在用分段氣動(dòng)力對(duì)其近似求解時(shí)，對(duì)于同一給定V條件，可能多支參與顫振分析的模態(tài)的特征根落在同一分段上，在該段上這些有效特征根與前一空速相應(yīng)模態(tài)的對(duì)應(yīng)關(guān)系無法確定，這時(shí)需要采用模態(tài)跟蹤法對(duì)其進(jìn)行跟蹤。本文的特征值跟蹤方式為：在某一空速Vm下，采用上述方法獲取參與顫振分析的各階模態(tài)的對(duì)應(yīng)的特征根及其對(duì)應(yīng)的有效分段區(qū)間，在此基礎(chǔ)之上，基于各階模態(tài)在其對(duì)應(yīng)的有效段上的特征根及其相應(yīng)的左特征向量（）和右特征向量（），根據(jù)特征值攝動(dòng)理論，來確定各階模態(tài)在下一空速Vm＋1下對(duì)應(yīng)的預(yù)測(cè)特征值如式（15），以預(yù)測(cè)值為參考對(duì)進(jìn)行排序。

特征值的對(duì)應(yīng)關(guān)系確定如下：在空速Vm下，求解空速Vm＋1下的預(yù)測(cè)特征值，然后從空速Vm＋1對(duì)應(yīng)的特征值（j＝1，2，...，n）中找到與之對(duì)應(yīng)的，使其滿足式（16）。

式中：

在模態(tài)的跟蹤過程中，可能出現(xiàn)速度步長(zhǎng)過大而導(dǎo)致模態(tài)跟蹤失敗，故給定一個(gè)閾值ε，當(dāng)d（，）＞ε時(shí)，減小速度步長(zhǎng)即將步長(zhǎng)除以因子T（T＞1）來重新計(jì)算特征值與排序，算法流程圖見圖1。

圖1 帶有模態(tài)跟蹤的PQI法流程圖Fig.1 The flow diagram of the PQI method with mode tracking

2 實(shí)例分析

2.1 數(shù)值分析一

顫振分析模型采用有限元軟件MSC.NASTRAN提供的標(biāo)準(zhǔn)算例：平板機(jī)翼的后掠角為15°，飛行馬赫數(shù)Ma為0.45。編寫偶極子格網(wǎng)法程序［11－15］求解氣動(dòng)力影響系數(shù)矩陣，采用有限元軟件MSC.NASTRAN獲取結(jié)構(gòu)的廣義質(zhì)量矩陣、廣義剛度矩陣、固有頻率和振型等，采用無限板樣條函數(shù)［16］實(shí)現(xiàn)氣動(dòng)結(jié)構(gòu)的耦合，編寫PQI法程序求解結(jié)構(gòu)的顫振臨界速度和頻率。該算例的折合頻率k的區(qū)間范圍為0～1，其分段區(qū)間均長(zhǎng)0.02，共50段，空速范圍是10m／s到220m／s，速度步長(zhǎng)為5m／s，選取結(jié)構(gòu)的前4階固有模態(tài)進(jìn)行顫振分析。

采用不帶模態(tài)跟蹤的PQI法進(jìn)行顫振分析，其V－g曲線和V－f曲線見圖2，模態(tài)竄支處折合頻率k的計(jì)算結(jié)果見表1，其中g(shù)代表模態(tài)阻尼，f為頻率。從圖2和表1可以得知，在空速為60m／s時(shí)，二階模態(tài)和三階模態(tài)的折合頻率分別為0.051、0.064，它們處在同一折合頻率分段0.05～0.07上，由于在該分段上的二階模態(tài)和三階模態(tài)與前一空速相應(yīng)的模態(tài)對(duì)應(yīng)關(guān)系不明確，從而造成模態(tài)對(duì)應(yīng)關(guān)系混亂，出現(xiàn)“模態(tài)交叉”現(xiàn)象。在空速為80m／s、85m／s、90m／s時(shí)，在折合頻率分段0.03～0.05上發(fā)生同樣情況。

圖2 不帶模態(tài)跟蹤的PQI法計(jì)算結(jié)果Fig.2 The solutions of PQI method without mode tracking

表1 模態(tài)竄支點(diǎn)處折合頻率的計(jì)算結(jié)果Table 1 Reduced frequencies of mode tracking

采用帶有模態(tài)跟蹤的PQI法與工程上常用的PK法對(duì)此模型進(jìn)行對(duì)比分析，兩種方法求解的V－g曲線和V－f曲線見圖3，從圖中可以看出，對(duì)于PQI法，采用模態(tài)跟蹤后，模態(tài)竄支問題得到很好的解決。此外，從圖3的V－g曲線知，平板機(jī)翼模型的顫振模態(tài)為第二階模態(tài)，該模態(tài)的阻尼曲線與g＝0的交點(diǎn)所對(duì)應(yīng)的速度即為顫振臨界速度，而兩種方法獲取的該支模態(tài)的阻力曲線與g＝0交點(diǎn)的位置幾乎重合，這說明兩種方法對(duì)顫振臨界速度的計(jì)算精度非常接近，PQI法和P－K法計(jì)算的顫振臨界速度分別為138.2m／s、140.6m／s，顫振頻率分別為127.3Hz、128.8Hz。

圖3 帶有模態(tài)跟蹤的PQI法和P－K法的計(jì)算結(jié)果Fig.3 The solutions of PQI method with mode tracking and P－K method

PQI法和P－K法計(jì)算結(jié)果的不同之處在于：（1）從圖3中曲線知，在空速小于185m／s時(shí)，模態(tài)阻尼在數(shù)值上略有差別，其原因是PQI法采用二次插值函數(shù)對(duì)阻尼的變化進(jìn)行了二階近似，其阻尼的計(jì)算更為精確；（2）在空速大于185m／s時(shí)，圖3中V－g和V－f曲線存在較大差異，其原因是P－K法中一階模態(tài)出現(xiàn)了竄支現(xiàn)象，獲得氣動(dòng)滯后根即零頻率根，因此該階模態(tài)的V－g曲線不連續(xù)，而在PQI法中能夠很好的跟蹤該階模態(tài)的特征根，從圖3中V－f曲線中可知，在空速大于185m／s時(shí)，一階模態(tài)的頻率大于零，從而保證了該階模態(tài)阻尼曲線的連續(xù)性，不過數(shù)值分析表明PQI法在空速185m／s到220m／s區(qū)間，該支模態(tài)其特征根p（p＝g＋ik）顯示｜g｜＞k，在此大阻尼情況下，其氣動(dòng)力的計(jì)算是不精確的，因此該區(qū)間的計(jì)算結(jié)果是可疑的，表2給出了其中部分特征根p的信息，圖4為不同阻尼情況下的廣義氣動(dòng)力分量Q44隨折合頻率k的變化趨勢(shì)，其中g(shù)＝0情況下的氣動(dòng)力是偶極子格網(wǎng)法程序計(jì)算結(jié)果，而g≠0情況下的氣動(dòng)力是通過分段氣動(dòng)力插值函數(shù)式（2）獲取，從表2與圖4可知，在大阻尼下當(dāng)｜g｜＞k時(shí)，通過分段二次插值函數(shù)獲取的氣動(dòng)力相對(duì)于小阻尼情況偏差較大。

表2 不同空速下的模態(tài)特征根Table 2 The eigenvalues of modes at different airspeeds

圖4 不同阻尼下的插值氣動(dòng)力Fig.4 Interpolation evaluated at different g

2.2 飛機(jī)全機(jī)顫振分析

為了驗(yàn)證本文發(fā)展方法在實(shí)際工程顫振分析中的有效性，將其用在商用軟件ZAERO［17］提供的F－18全機(jī)軸對(duì)稱狀態(tài)顫振分析中，折合頻率區(qū)間范圍為0.01～0.8，將其分成60分段，空速范圍是120m／s到360m／s，速度步長(zhǎng)為5m／s，選取結(jié)構(gòu)的前12階結(jié)構(gòu)固有模態(tài)，采用PQI法進(jìn)行顫振分析，其V－g和V－f曲線見圖5，采用帶有模態(tài)跟蹤的PQI法分析，其V－g和V－f曲線見圖6。從圖5可以看出，對(duì)于復(fù)雜的結(jié)構(gòu)，采用PQI方法模態(tài)竄支現(xiàn)象比較嚴(yán)重，顫振模態(tài)的識(shí)別受到干擾，采用模態(tài)跟蹤法之后，圖6顯示其V－g曲線和V－f曲線光滑合理，其顫振速度為250.05m／s，顫振頻率為23.25Hz。

圖5 不帶模態(tài)跟蹤的PQI法的計(jì)算結(jié)果Fig.5 The solutions of the PQI method without mode tracking

圖6 帶有模態(tài)跟蹤的PQI法的計(jì)算結(jié)果Fig.6 The solutions of the PQI method with mode tracking

3 結(jié) 論

（1）本文針對(duì)PQI法，發(fā)展了一種基于特征值攝動(dòng)理論的模態(tài)跟蹤技術(shù)－預(yù)測(cè)跟蹤法，解決了其存在的模態(tài)竄支問題，同時(shí)也驗(yàn)證了該方法在模態(tài)跟蹤過程中具有良好的跟蹤效果。

（2）PQI法相對(duì)于PK法的最大優(yōu)勢(shì)在于對(duì)廣義氣動(dòng)力采用了分段表達(dá)的形式，從而將顫振方程求解的非線性二次特征值問題轉(zhuǎn)化為分段線性二次特征值問題，即采用分段直接求解的思想，因此不需要迭代，而PK法采用迭代的方式求解顫振方程的特征根，故需要預(yù)先給出一定的收斂裕度，對(duì)于復(fù)雜的結(jié)構(gòu)可能存在不收斂性。

（3）從顫振計(jì)算精度上看，PQI法與工程上常用的PK法的計(jì)算精度非常接近。

（4）從阻尼的角度看，相對(duì)于PK法，PQI法對(duì)阻尼的變化進(jìn)行了二階近似，其阻尼的計(jì)算更為精確。

［1］RODDEN W P，JOHNSON E H.MSC／Nastran v68aeroelastic analysis user′s guide［M］.Los Angeles：MSC.Software Corporation，1994：73－76.

［2］HASSIG H J.An approximate true damping solution of the flutter equation by determinant iteration［J］.JournalofAircraft，1971，8（11）：885－889.doi：10.2514／3.44311

［3］GOODMAN C.Accurate subcritical damping solution of flutter equation using piecewise aerodynamic function［J］.Journalof Aircraft，2001，38（4）：755－763.doi：10.2514／2.2828

［4］ELLER D.Flutter equation as piecewise quadratic eigenvalue problem［J］.JournalofAircraft，2009，46（3）：1068－1070.doi：10.2514／1.40653

［5］XIE C C，HU W W，YANG C.Eigenvalue sorting problem in flutter analysis［J］.MathematicsinPracticeandTheory，2007，37（18）：141－146.（in Chinese）謝長(zhǎng)川，胡薇薇，楊超.顫振分析中的特征值排序問題［J］.數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2007，37（18）：141－146.doi：10.3969／j.issn.1000－0984.2007.18.022

［6］Van ZYL L H.Use of eigenvectors in the solution of the flutter equation［J］.JournalofAircraft，1993，30（4）：553－554.doi：10.2514／3.46380

［7］CHEN P C.A damping perturbation method for flutter solution：the g－method［J］.AIAAJournal，2000，38（9）：1519－1524.doi：10.2514／2.1171

［8］GUAN D.Aeroelasticity handbook of aircraft［M］.Beijing Aeronautical Industry Press，1994：126－134.（in Chinese）管德.飛機(jī)氣動(dòng)彈性力學(xué)手冊(cè)［M］.北京：航空工業(yè)出版社，1994：126－134.

［9］LOUIS A P.Applied mathematics for engineers and physicists［M］.New York：Mograw－Hill Book Company，1946：448－452.

［10］FANG T，XUE P.Theory of vibration with applications［M］.Xi′an：Northwestern Polytechnical University Press，1998：130－136.（in Chinese）方同，薛璞.振動(dòng)理論及應(yīng)用［M］.西安：西北工業(yè)大學(xué)出版社，1998：130－136.

［11］ALBANO E，RODDEN W P.A doublet－lattice method for calculating lift distributions of oscillating surfaces in subsonic flow［J］.AIAAJournal，1969，7（2）：279－285.doi：10.2514／6.1968－73

［12］RODDEN W P，GIESING J P.Refinement of the nonplanar aspects of the subsonic doublet－lattice lifting－surface method［J］.JournalofAircraft，1972，9（1）：69－73.doi：10.2514／2.2382

［13］RODDEN W P，TAYLOR P F.Further refinement of the nonplanar aspects of the subsonic doublet－lattice lifting surface method［J］.JournalofAircraft，1998，35（5）：720－727.doi：10.2514／3.44322

［14］YANG L F.A new method of flutter analysis［J］.ACTAAerodynamicaSinica，1994，12（2）：196－199.（in Chinese）楊立法.一種顫振分析新方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1994，12（2）：196－199.

［15］SUN S P，CHEN J S.Doublet point method for harmonically oscillating lifting surfaces in subsonic flow［J］.ACTAAerodynamicaSinica，1993，（1）：011.（in Chinese）孫少鵬，陳勁松.亞聲速諧振升力面的點(diǎn)偶極子法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1993，（1）：011.

［16］HARDER R L，DESMARAIS R N.Interpolation using surface splines［J］.AIAAJournal，1972，9（2）：189－191.doi：10.2514／3.44330

［17］ZONA Technology Inc.ZAERO version 7.4theoretical manual［M］.Scottsdale：ZONA Technology Inc.，2006：297－302.