左彥飛，吳易柳，王杰，馮坤，江志農(nóng)，*

1.北京化工大學(xué) 發(fā)動機(jī)健康監(jiān)控及網(wǎng)絡(luò)化教育部重點(diǎn)實(shí)驗(yàn)室，北京 100029

2.航空動力技術(shù)研究院，西安 710025

3.中國航發(fā)動力所-北京化工大學(xué) 航空發(fā)動機(jī)振動健康監(jiān)控聯(lián)合實(shí)驗(yàn)室，北京 100029

雙轉(zhuǎn)子燃?xì)鉁u輪發(fā)動機(jī)工作過程中經(jīng)受著嚴(yán)酷的氣動、熱等復(fù)合環(huán)境影響，轉(zhuǎn)子上存在不均勻溫度場［1-2］。隨著溫度的變化，轉(zhuǎn)子材料屬性（如彈性模量、泊松比等）隨之發(fā)生改變，影響雙轉(zhuǎn)子系統(tǒng)的動力特性。考慮溫致材料屬性變化的影響對于發(fā)動機(jī)轉(zhuǎn)子系統(tǒng)精細(xì)化動力學(xué)設(shè)計(jì)、振動故障分析等具有重要的參考價(jià)值和應(yīng)用前景。

目前，針對發(fā)動機(jī)的雙轉(zhuǎn)子系統(tǒng)，洪杰［3-4］、張大義［5-6］、左彥飛［7-9］等使用三維實(shí)體有限元模型研究了臨界轉(zhuǎn)速、穩(wěn)態(tài)不平衡響應(yīng)的計(jì)算分析方法，對典型發(fā)動機(jī)的雙轉(zhuǎn)子系統(tǒng)動力特性進(jìn)行了分析，得到了系統(tǒng)響應(yīng)特征。

在發(fā)動機(jī)的運(yùn)行過程中，由于溫度變化造成連接結(jié)構(gòu)的不穩(wěn)定、熱變形或轉(zhuǎn)子彎曲，轉(zhuǎn)子不平衡量變大，進(jìn)而導(dǎo)致轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速及振動響應(yīng)發(fā)生顯著變化［10-11］。同時(shí)軸承油膜溫度變化及轉(zhuǎn)子熱變形會導(dǎo)致轉(zhuǎn)子與軸承之間的間隙發(fā)生變化進(jìn)而影響軸承的剛度，降低軸承轉(zhuǎn)子系統(tǒng)的穩(wěn)定性，影響轉(zhuǎn)子的動力特性［12-14］。

溫度不僅能引發(fā)轉(zhuǎn)子熱變形，同時(shí)還會引起轉(zhuǎn)子材料屬性的改變，造成轉(zhuǎn)子系統(tǒng)剛度變化，對轉(zhuǎn)子振動特性產(chǎn)生影響［11］。劉知輝［15］、李暉［16］等計(jì)算了溫度場作用下梁、圓柱殼的共振頻率及共振響應(yīng)，得到了結(jié)構(gòu)振動特性隨溫度變化的規(guī)律。張婷婷［17］、張明根［18］、何鵬［19］等對溫度均勻分布、線性分布和二次多項(xiàng)式分布情況下的轉(zhuǎn)子系統(tǒng)的動力特性變化進(jìn)行了一系列規(guī)律研究，劉少權(quán)［20］、Liu［21］等計(jì)算了單轉(zhuǎn)子三維有限元模型的溫度場分布，得到考慮溫度場前后臨界轉(zhuǎn)速、不平衡響應(yīng)等的變化情況。

燃?xì)鉁u輪發(fā)動機(jī)內(nèi)部溫差極大、溫度分布復(fù)雜并隨發(fā)動機(jī)工況變化改變，而且雙轉(zhuǎn)子結(jié)構(gòu)比單轉(zhuǎn)子具有更為復(fù)雜的動力學(xué)特性，分析難度更大，要求更高，溫致材料屬性變化對發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)動力特性的影響研究還不夠深入。因此，本文重點(diǎn)關(guān)注溫度場導(dǎo)致的材料屬性變化對轉(zhuǎn)子系統(tǒng)動力特性的影響。通過發(fā)動機(jī)性能方程、相似性原理及穩(wěn)態(tài)熱分析方法擬合得到穩(wěn)定工況典型發(fā)動機(jī)雙轉(zhuǎn)子溫度場，推導(dǎo)雙轉(zhuǎn)子系統(tǒng)受溫度場影響的動力學(xué)方程，提出將不同工況溫度、轉(zhuǎn)速變化與發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)動力特性聯(lián)合分析方法，建立雙轉(zhuǎn)子熱-固聯(lián)合分析有限元模型。分析并得到受溫度場影響雙轉(zhuǎn)子系統(tǒng)固有頻率及對應(yīng)振型、臨界轉(zhuǎn)速及應(yīng)變能的變化規(guī)律。

1 發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)溫度場分布求解

為得到不同工況雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)溫度場分布，利用發(fā)動機(jī)試車臺可測量溫度參數(shù)，通過溫度邊界條件求解及雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)熱分析2個(gè)步驟實(shí)現(xiàn)。

1.1 溫度邊界條件求解

某典型發(fā)動機(jī)試車圖譜如圖1所示，發(fā)動機(jī)只在幾個(gè)穩(wěn)定工況下長時(shí)間停留，其余工況快速通過，因此，重點(diǎn)分析雙轉(zhuǎn)子系統(tǒng)穩(wěn)定工況溫度分布。

圖1 某典型發(fā)動機(jī)試車圖譜Fig.1 Diagram of a typical engine test run

典型雙轉(zhuǎn)子渦扇發(fā)動機(jī)穩(wěn)定工況下溫度邊界對應(yīng)位置及代表符號如圖2所示。圖中T1～T6分別代表低壓壓氣機(jī)入口溫度、低壓壓氣機(jī)出口（高壓壓氣機(jī)入口）溫度、高壓壓氣機(jī)出口溫度、高壓渦輪入口溫度、高壓渦輪出口溫度、低壓渦輪出口溫度。

圖2 雙轉(zhuǎn)子關(guān)鍵截面溫度邊界條件Fig.2 Temperature boundary condition for critical sections of dual rotors

設(shè)計(jì)工況各部位溫度數(shù)據(jù)一般可通過設(shè)計(jì)參數(shù)直接獲取，各穩(wěn)定工況下高壓壓氣機(jī)出口溫度T3及低壓渦輪出口溫度T6一般可通過溫度傳感器測量獲得，而其他部位溫度通常很難直接測量，需結(jié)合發(fā)動機(jī)性能公式、相似性原理及壓氣機(jī)-渦輪效率平衡公式近似求解。

雙轉(zhuǎn)子發(fā)動機(jī)主要?dú)饴泛唸D如圖3所示，空氣經(jīng)過低壓壓氣機(jī)后（W1）分成兩路，一路是內(nèi)涵氣流進(jìn)入發(fā)動機(jī)經(jīng)壓氣機(jī)壓縮（W21），對氣體做功使氣體的壓力升高、容積減小。另一路是外涵氣流（W22），Wi代表圖3對應(yīng)結(jié)構(gòu)空氣或燃?xì)赓|(zhì)量流量。

圖3 發(fā)動機(jī)氣路性能部件簡圖Fig.3 Schematic diagram of engine gas path performance components

對于發(fā)動機(jī)地面試車臺，低壓壓氣機(jī)入口溫度T1與大氣溫度T0數(shù)值近似相等：

利用低壓壓氣機(jī)運(yùn)行過程中遵循的氣動熱力學(xué)方程，計(jì)算得到低壓壓氣機(jī)出口溫度T2為

式中：ηcL為低壓壓氣機(jī)熱效率；r=()k-1k，k代表空氣絕熱系數(shù)；πcL為低壓壓氣機(jī)增壓比。

不考慮低壓壓氣機(jī)效率ηcL的變化，經(jīng)過低壓壓氣機(jī)壓縮后的氣體溫度變化與低壓壓氣機(jī)增壓比相關(guān)，即：

根據(jù)相似原理［22］，在雷諾數(shù)Re自?；瘏^(qū)內(nèi)，氣體的參數(shù)場是相似的，在壓氣機(jī)進(jìn)口氣流穩(wěn)定的條件下，壓氣機(jī)的壓比πcL與相似參數(shù)（軸向馬赫數(shù)Maa和周向馬赫數(shù)Mau）之間的關(guān)系為

在實(shí)際繪制壓氣機(jī)特性曲線時(shí)，常用與相似準(zhǔn)則成正比的量（用參數(shù)代替Maa，參數(shù)代替Mau），將參數(shù)代入后式（4）變?yōu)?/p>

式中：ncL為低壓壓氣機(jī)轉(zhuǎn)速。

聯(lián)立式（3）和式（5），在發(fā)動機(jī)低壓壓氣機(jī)出口溫度由轉(zhuǎn)速、進(jìn)氣流量、進(jìn)氣溫度及進(jìn)氣壓強(qiáng)共同決定。關(guān)系可表示為

渦輪將燃燒室流出的高溫燃?xì)獾臒崮芗皦毫δ苻D(zhuǎn)化為機(jī)械能，使其做功帶動壓氣機(jī)運(yùn)行，結(jié)合圖3，列出該發(fā)動機(jī)中壓氣機(jī)與渦輪效率平衡公式［23］。低壓壓氣機(jī)及低壓渦輪效率平衡公式為

高壓壓氣機(jī)及高壓渦輪效率平衡公式為

式中：Cp為空氣的定壓比熱容；Cpg為燃?xì)舛▔罕葻崛?；ηcH為高壓壓氣機(jī)熱效率。

將已知溫度數(shù)據(jù)T1、T6及通過式（6）代入低壓壓氣機(jī)轉(zhuǎn)速ncL求解得到的T2代入式（7）中求得各穩(wěn)定工況高壓渦輪出口溫度數(shù)據(jù)T5，之后將T2、T3、T5代入式（8）中，求解得到各穩(wěn)定工況下高壓渦輪入口溫度T4，至此，圖2所示截面溫度均可求解獲得。

1.2 雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)熱分析

發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)穩(wěn)定工況下，熱能流動不隨時(shí)間發(fā)生改變，系統(tǒng)的溫度和熱載荷也不隨時(shí)間發(fā)生改變，故對雙轉(zhuǎn)子系統(tǒng)進(jìn)行穩(wěn)態(tài)熱分析，穩(wěn)態(tài)熱平衡方程為

式中：K為傳導(dǎo)矩陣，包括導(dǎo)熱系數(shù)及形狀系數(shù)；T為節(jié)點(diǎn)溫度向量；Q為節(jié)點(diǎn)熱流率向量，包含熱生成。K、T及Q均可通過ANSYS等商用有限元程序綜合運(yùn)用所建模型的參數(shù)、添加材料的熱性能參數(shù)以及各類邊界條件自動生成。

為便于分析，假設(shè)穩(wěn)定工況下圖2所示位置處轉(zhuǎn)子溫度與相應(yīng)區(qū)域氣流溫度近似相等，將1.1節(jié)中求解得到的各穩(wěn)定工況下特定區(qū)域氣流溫度作為溫度邊界條件施加在發(fā)動機(jī)雙轉(zhuǎn)子模型特定位置的節(jié)點(diǎn)上，各區(qū)域通用溫度施加方式為

式中：Tb代表圖2中所求各邊界區(qū)域溫度；代表溫度隨位置變化的函數(shù)。

發(fā)動機(jī)轉(zhuǎn)子各部件材料導(dǎo)熱系數(shù)隨溫度變化發(fā)生改變，利用溫度插值選?。?/p>

式中：Tr代表節(jié)點(diǎn)溫度；Ta、Ta+1為材料熱性能表上與節(jié)點(diǎn)溫度最接近的兩參考溫度；λa、λa+1代表兩參考溫度對應(yīng)導(dǎo)熱系數(shù)；λr代表通過線性插值法得到的節(jié)點(diǎn)導(dǎo)熱系數(shù)。

考慮非線性因素后（導(dǎo)熱系數(shù)隨溫度變化）穩(wěn)態(tài)熱平衡方程可表示為

式中：傳導(dǎo)矩陣K(T)、節(jié)點(diǎn)熱流率向量Q(T)會隨著節(jié)點(diǎn)溫度向量的變化發(fā)生改變。

方程可等效化為

式中：Qnr為內(nèi)部節(jié)點(diǎn)熱流向量，通過計(jì)算單元熱流求得；Qa為載荷引起的節(jié)點(diǎn)熱流向量。

初始情況下，載荷引起的節(jié)點(diǎn)熱流向量與內(nèi)部節(jié)點(diǎn)熱流向量不相等，不平衡熱流向量Φ是2個(gè)向量的差值：

對穩(wěn)態(tài)熱分析非線性方程采用Newton-Raphson方法求解，其增量形式為

式中：Ki、Qnri、Ti分別表示第i次迭代對應(yīng)的傳導(dǎo)矩陣、內(nèi)部節(jié)點(diǎn)熱流向量及節(jié)點(diǎn)溫度向量；ΔTi為節(jié)點(diǎn)溫度變化矩陣。

隨著迭代的進(jìn)行，逐步更新節(jié)點(diǎn)溫度場Ti、計(jì)算內(nèi)部節(jié)點(diǎn)熱流向量Qnri將不平衡熱流向量數(shù)值‖‖Φ逐步減小，一般情況下停止迭代的條件為

停止迭代時(shí)計(jì)算得到的節(jié)點(diǎn)溫度向量Ti即為雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)近似溫度場Ts。雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)溫度場求解是溫致材料屬性影響分析計(jì)算的基礎(chǔ)。

2 雙轉(zhuǎn)子系統(tǒng)動力特性熱-固聯(lián)合分析

雙轉(zhuǎn)子系統(tǒng)動力特性熱-固聯(lián)合分析為單向耦合計(jì)算，主要考慮雙轉(zhuǎn)子系統(tǒng)在溫度場作用下的材料屬性的變化對動力特性的影響。

2.1 溫致材料屬性影響下雙轉(zhuǎn)子系統(tǒng)動力學(xué)方程

發(fā)動機(jī)轉(zhuǎn)子部件的材料物理特性（主要是彈性模量及泊松比）會隨著溫度的變化發(fā)生較大的改變，例如圖4所示典型發(fā)動機(jī)高壓渦輪軸材料屬性隨溫度變化曲線［24］，具體數(shù)值可通過插值法求?。?/p>

圖4 典型發(fā)動機(jī)高壓渦輪軸材料屬性隨溫度的變化［24］Fig.4 Variation of material properties with temperature for a typical engine high-pressure turboshaft［24］

式中：Ea、μa分別代表參考溫度對應(yīng)彈性模量及泊松比；Er、μr代表通過線性插值法得到的節(jié)點(diǎn)所在位置彈性模量及泊松比。

隨著發(fā)動機(jī)工況的變化，各單元體積V及各節(jié)點(diǎn)位置不發(fā)生變化，但會影響發(fā)動機(jī)內(nèi)部各處節(jié)點(diǎn)溫度，使材料彈性模量及泊松比發(fā)生改變，進(jìn)而影響各節(jié)點(diǎn)對應(yīng)剛度矩陣。

例如對于轉(zhuǎn)子動力學(xué)分析過程中常用的6面體8節(jié)點(diǎn)單元，各節(jié)點(diǎn)剛度矩陣為

雙轉(zhuǎn)子支承系統(tǒng)動力學(xué)方程為

式中：M為系統(tǒng)質(zhì)量矩陣；C為系統(tǒng)阻尼矩陣；ΩL、ΩH分別為低壓轉(zhuǎn)子、高壓轉(zhuǎn)子的轉(zhuǎn)速；GL、GH分別為低壓轉(zhuǎn)子、高壓轉(zhuǎn)子旋轉(zhuǎn)產(chǎn)生的陀螺效應(yīng)矩陣；KD為系統(tǒng)剛度矩陣；KL、KH分別為與低壓轉(zhuǎn)子、高壓轉(zhuǎn)子轉(zhuǎn)速相關(guān)的剛度矩陣，稱為旋轉(zhuǎn)剛度矩陣；fL、fH分別為低壓轉(zhuǎn)子、高壓轉(zhuǎn)子激勵向量；u、分別為雙轉(zhuǎn)子系統(tǒng)的位移向量、速度向量、加速度向量。

由于發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)不同工況下穩(wěn)態(tài)溫度場Ts的引入，根據(jù)式（18），節(jié)點(diǎn)剛度矩陣發(fā)生變化，系統(tǒng)剛度矩陣KD、低壓轉(zhuǎn)子和高壓轉(zhuǎn)子轉(zhuǎn)速相關(guān)的剛度矩陣KL、KH受到影響，即：

式中：TsL、TsH分別為低壓轉(zhuǎn)子、高壓轉(zhuǎn)子穩(wěn)態(tài)近似溫度場。

因此考慮溫致材料屬性變化后的雙轉(zhuǎn)子支承系統(tǒng)動力學(xué)方程可表示為

溫度場主要影響系統(tǒng)的剛度矩陣，進(jìn)而影響系統(tǒng)固有特性、振動響應(yīng)特性等。

2.2 溫致材料屬性變化對系統(tǒng)動力學(xué)影響評估方案及指標(biāo)

熱-固聯(lián)合分析有限元模型對應(yīng)溫度場會隨發(fā)動機(jī)工況的變化而發(fā)生改變，能較精確地模擬發(fā)動機(jī)工作溫度場變化對材料屬性的影響。在此基礎(chǔ)上，通過評估雙轉(zhuǎn)子系統(tǒng)在不同工況溫度場下固有頻率、模態(tài)振型、穩(wěn)態(tài)不平衡響應(yīng)及應(yīng)變能分布變化情況探究溫致材料屬性變化對系統(tǒng)動力學(xué)的影響程度和變化規(guī)律。

利用振型相似性指標(biāo)——振型置信因子（量符號記為Cma）來定量評估振型的受影響程度。由于工程中比較關(guān)注各穩(wěn)定工況固有頻率對應(yīng)的轉(zhuǎn)子振型，所以利用固有頻率對應(yīng)的模態(tài)振型置信因子變化進(jìn)行評估，其表達(dá)式為［9］

式中：Ψnm表示常溫雙轉(zhuǎn)子系統(tǒng)第m階模態(tài)的振型向量；Ψsn表示特定工況溫度場雙轉(zhuǎn)子系統(tǒng)第n階模態(tài)的振型向量；Cmamn為同一雙轉(zhuǎn)子系統(tǒng)常溫第m階固有頻率與特定工況溫度場第n階固有頻率的模態(tài)振型置信因子。

通過計(jì)算模態(tài)振型置信因子，得到常溫與特定溫度場雙轉(zhuǎn)子系統(tǒng)不同階模態(tài)振型的相似性。結(jié)合對應(yīng)的振型，評估溫致材料屬性變化對雙轉(zhuǎn)子振型的影響程度。通過Cmamn計(jì)算結(jié)果匹配確定2組振型Ψnp、Ψsq對應(yīng)的固有頻率fnp、fsq，并通過Campbell圖法計(jì)算對應(yīng)雙轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速（如圖5所示）；固有頻率與臨界轉(zhuǎn)速相對變化率為

圖5 常溫與特定工況雙轉(zhuǎn)子系統(tǒng)固有頻率、臨界轉(zhuǎn)速、不平衡響應(yīng)幅值變化示意圖Fig.5 Schematic diagram of variation of natural frequency， critical speed， and unbalanced response amplitude of dual-rotor system at room temperature and under specific operating conditions

式中：fnp表示常溫雙轉(zhuǎn)子系統(tǒng)第p階固有頻率；fsq表示特定工況溫度場雙轉(zhuǎn)子系統(tǒng)第q階固有頻率；Ωcnp、Ωcsq分別為常溫雙轉(zhuǎn)子系統(tǒng)第p階臨界轉(zhuǎn)速和特定工況溫度場雙轉(zhuǎn)子系統(tǒng)第q階臨界轉(zhuǎn)速；Φfpq、Φcpq分別為固有頻率及臨界轉(zhuǎn)速的相對變化率。

進(jìn)行不平衡響應(yīng)分析得到各階臨界轉(zhuǎn)速對應(yīng)不平衡響應(yīng)峰值及設(shè)計(jì)工作轉(zhuǎn)速范圍平均不平衡響應(yīng)幅值并計(jì)算相對變化率（如圖5所示）：

式中：qn、qs分別為常溫與特定工況溫度場雙轉(zhuǎn)子系統(tǒng)不平衡響應(yīng)幅值；qnp、qsp分別為常溫與特定工況溫度場下同階臨界轉(zhuǎn)速對應(yīng)不平衡響應(yīng)幅值；Φqp為臨界轉(zhuǎn)速對應(yīng)不平衡響應(yīng)幅值的相對變化率；Φqa為雙轉(zhuǎn)子系統(tǒng)(Ω1，Ω2)轉(zhuǎn)速范圍內(nèi)的平均不平衡響應(yīng)幅值相對變化率。

在發(fā)動機(jī)轉(zhuǎn)子系統(tǒng)設(shè)計(jì)中，可根據(jù)應(yīng)變能分布準(zhǔn)確判定系統(tǒng)振型中轉(zhuǎn)子與機(jī)匣振動主次，應(yīng)變能分布具有重要作用。為此對雙轉(zhuǎn)子系統(tǒng)常溫與特定溫度場臨界轉(zhuǎn)速對應(yīng)振型的應(yīng)變能分布變化情況進(jìn)行定量分析。

3 溫致材料屬性變化對典型雙轉(zhuǎn)子系統(tǒng)動力特性影響

以典型發(fā)動機(jī)同向旋轉(zhuǎn)雙轉(zhuǎn)子系統(tǒng)為對象，利用所提分析方法，研究溫致材料屬性變化對雙轉(zhuǎn)子動力特性的影響。

3.1 典型發(fā)動機(jī)雙轉(zhuǎn)子支承系統(tǒng)

某典型發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)由高壓轉(zhuǎn)子與低壓轉(zhuǎn)子組成，共有6個(gè)支點(diǎn)位置，其中1、2、3和6號支點(diǎn)為低壓轉(zhuǎn)子支點(diǎn)，4號支點(diǎn)為高壓轉(zhuǎn)子支點(diǎn)，5號支點(diǎn)為中介支點(diǎn)，具體三維有限元模型如圖6所示。

圖6 某典型發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)三維有限元模型Fig.6 Three-dimensional finite element model of a typical engine dual-rotor system

雙轉(zhuǎn)子系統(tǒng)的對稱支承經(jīng)驗(yàn)剛度見表1［9］，后續(xù)分析主要研究轉(zhuǎn)子材料屬性受溫度影響，暫不考慮支承、軸承剛度隨溫度的變化。

表1 雙轉(zhuǎn)子系統(tǒng)的對稱支承經(jīng)驗(yàn)剛度［9］Table 1 Empirical stiffness of symmetrical support for a dual-rotor system［9］

3.2 雙轉(zhuǎn)子系統(tǒng)穩(wěn)態(tài)溫度場分布

典型發(fā)動機(jī)工作過程中經(jīng)過數(shù)個(gè)穩(wěn)定工況，重點(diǎn)分析3個(gè)典型工況，分別為80%工況、90%工況、100%工況（例如：80%工況的含義是當(dāng)前高壓轉(zhuǎn)子轉(zhuǎn)速為高壓轉(zhuǎn)子最大轉(zhuǎn)速值的80%）。

利用設(shè)計(jì)及實(shí)測溫度數(shù)據(jù)，根據(jù)1.1節(jié)所述方法求解得到不同工況關(guān)鍵部位（見圖2）溫度數(shù)據(jù)如圖7所示。

圖7 特定工況雙轉(zhuǎn)子系統(tǒng)各關(guān)鍵部位溫度Fig.7 Temperature of key parts of dual-rotor system under specific working condition

將各工況關(guān)鍵部位溫度數(shù)據(jù)作為溫度邊界條件代入雙轉(zhuǎn)子有限元模型后進(jìn)行穩(wěn)態(tài)熱分析，得到雙轉(zhuǎn)子各穩(wěn)定工況溫度場Ts，圖8為100%工況雙轉(zhuǎn)子系統(tǒng)溫度場分布。

圖8 100%工況雙轉(zhuǎn)子系統(tǒng)溫度場分布Fig.8 Temperature field distribution of dual-rotor system under 100% working condition

3.3 溫致材料屬性變化對固有頻率及振型的影響

分別對常溫及熱-固聯(lián)合分析有限元模型進(jìn)行模態(tài)分析，得到3個(gè)典型穩(wěn)定工況下固有頻率及對應(yīng)振型。通過式（22）、式（23）分別計(jì)算得到對應(yīng)雙轉(zhuǎn)子模態(tài)振型置信因子Cma及固有頻率相對變化率，具體分析結(jié)果及對比數(shù)據(jù)如表2所示。表中F表示正進(jìn)動，B表示反進(jìn)動，例如，1stF表示第1階正進(jìn)動。

表2 常溫及熱-固聯(lián)合分析有限元模型固有頻率及對應(yīng)模態(tài)振型置信因子結(jié)果比較Table 2 Comparison of natural frequency and mode assurance criterion results corresponding to normal temperature and heat-solid joint analysis finite element model

利用表2中常溫及熱-固聯(lián)合固有頻率變化率與雙轉(zhuǎn)子系統(tǒng)各階振型圖，將臨近工作轉(zhuǎn)速范圍的雙轉(zhuǎn)子系統(tǒng)橫向振型大致分為支承為主、低壓壓氣機(jī)俯仰、轉(zhuǎn)子或渦輪軸彎曲振型3大類。

本文重點(diǎn)分析轉(zhuǎn)子溫度場影響，熱-固聯(lián)合有限元模型中支承采用不考慮溫度影響的彈簧單元模擬，同時(shí)壓氣機(jī)部位溫差較?。▓D7），因此，支承為主振型（1st、2nd）對應(yīng)固有頻率變化率＜2%，低壓壓氣機(jī)俯仰振型（3rd、6th）對應(yīng)固有頻率變化率＜0.5%，影響較小，不做重點(diǎn)分析。

對于圖9所示的高壓轉(zhuǎn)子及低壓轉(zhuǎn)子渦輪軸彎曲振型，由于渦輪及渦輪軸緊鄰燃燒室，工作環(huán)境溫度極高，導(dǎo)致對應(yīng)固有頻率及振型受溫度場影響較大，正反進(jìn)動固有頻率最大變化率分別達(dá)到了8.1%（5thF）和26.5%（4thB）。

圖9 高壓轉(zhuǎn)子及低壓轉(zhuǎn)子渦輪軸彎曲振型Fig.9 Bending vibration modes of high-pressure rotor and low-pressure rotor turboshaft

利用各階進(jìn)動頻率數(shù)據(jù)（篩除周向、軸向等振型）得到包含前6階振型的常溫及熱-固聯(lián)合分析有限元模型Campbell圖如圖10所示。受溫度場影響，臨界轉(zhuǎn)速產(chǎn)生不同程度的下降，最大下降幅度達(dá)到9.9%（5thF臨界轉(zhuǎn)速），與規(guī)定裕度20%達(dá)到同一量級。

圖10 溫度場施加前后Campbell圖對比Fig.10 Comparison of Campbell diagrams before and after application of temperature field

對于典型發(fā)動機(jī)同向旋轉(zhuǎn)雙轉(zhuǎn)子系統(tǒng)，當(dāng)進(jìn)動方向與轉(zhuǎn)子轉(zhuǎn)動方向相同時(shí)稱為正進(jìn)動，常溫雙轉(zhuǎn)子正進(jìn)動頻率會隨著工作轉(zhuǎn)速的上升而上升，而引入溫度場后，轉(zhuǎn)子材料彈性模量減小，整體剛度矩陣發(fā)生改變，4thF、5thF進(jìn)動頻率隨轉(zhuǎn)速的變化由升轉(zhuǎn)降，5thF臨界轉(zhuǎn)速受溫度影響較大，先于4thF臨界轉(zhuǎn)速到達(dá)（如圖10中放大區(qū)域）。

通過表2發(fā)現(xiàn)部分階次模態(tài)振型置信因子Cma＜0.9，為分析現(xiàn)象產(chǎn)生原因，單獨(dú)導(dǎo)出常溫雙轉(zhuǎn)子系統(tǒng)Campbell圖如圖11所示。通過圖11中Campbell圖左上放大區(qū)發(fā)現(xiàn)，4thF、5thF這2階模態(tài)的動頻曲線逐漸接近。但是在80%工況附近接近的趨勢不再繼續(xù)， 而是相互分離，表現(xiàn)為動頻曲線間的轉(zhuǎn)向， 發(fā)生“頻率轉(zhuǎn)向”現(xiàn)象［25］。最終造成常溫與熱-固聯(lián)合分析有限元模型80%工況時(shí)4thF、5thF對應(yīng)振型的振型相關(guān)系數(shù)Cma僅為0.73與0.65。

圖11 常溫雙轉(zhuǎn)子系統(tǒng)Campbell圖Fig.11 Campbell diagram of dual-rotor system at room temperature

3.4 溫致材料屬性變化對不平衡響應(yīng)的影響

高壓轉(zhuǎn)子是雙轉(zhuǎn)子發(fā)動機(jī)中的核心機(jī)，位于燃燒室附近，工作環(huán)境溫度極高，轉(zhuǎn)子溫差達(dá)到1 200 K（見圖7）。為此對高壓轉(zhuǎn)子上存在不平衡量的雙轉(zhuǎn)子系統(tǒng)進(jìn)行不平衡響應(yīng)分析。

為模擬實(shí)際發(fā)動機(jī)工作狀態(tài)下存在較大不平衡的情況，在高壓渦輪葉盤位置添加2 g·m不平衡量，系統(tǒng)阻尼比設(shè)為0.01，通過不平衡響應(yīng)分析求解得到常溫雙轉(zhuǎn)子系統(tǒng)各支承支點(diǎn)響應(yīng)幅值，如圖12所示（圖中5?-L、5?-H分別代表5號支承與低壓轉(zhuǎn)子、高壓轉(zhuǎn)子的連接點(diǎn)）。雙轉(zhuǎn)子系統(tǒng)各支點(diǎn)臨界轉(zhuǎn)速位置響應(yīng)的相對大小與振型及激勵位置相關(guān)，此處不詳細(xì)討論，與文獻(xiàn)［7］中研究結(jié)果一致。

圖12 常溫高壓渦輪激勵不平衡響應(yīng)曲線Fig.12 Unbalance response curves at normal temperature with high-pressure turbine excitation

4thF、5thF臨界轉(zhuǎn)速靠近100%工況工作轉(zhuǎn)速，對應(yīng)振型為3.3節(jié)所示高壓轉(zhuǎn)子及低壓轉(zhuǎn)子渦輪軸彎曲振型且不平衡響應(yīng)幅值較大，故后續(xù)主要分析特定穩(wěn)態(tài)溫度場下200～250 Hz雙轉(zhuǎn)子系統(tǒng)的不平衡響應(yīng)。

各穩(wěn)定溫度場下雙轉(zhuǎn)子系統(tǒng)不平衡響應(yīng)曲線如圖13（雙轉(zhuǎn)子系統(tǒng)高壓轉(zhuǎn)子最高轉(zhuǎn)速對應(yīng)頻率為220 Hz）所示。

圖13 不同溫度場下雙轉(zhuǎn)子系統(tǒng)不平衡響應(yīng)曲線Fig.13 Unbalanced response curves of dual-rotor system in different temperature fields

從圖13可以看出，隨著雙轉(zhuǎn)子系統(tǒng)施加溫度場對應(yīng)工況的提升，4thF、5thF臨界轉(zhuǎn)速下降、逐漸靠近工作轉(zhuǎn)速區(qū)，5thF臨界轉(zhuǎn)速下降幅度較大，100%工況溫度場5thF臨界轉(zhuǎn)速對應(yīng)不平衡響應(yīng)峰值出現(xiàn)在4thF臨界轉(zhuǎn)速之前，與3.3節(jié)模態(tài)分析結(jié)論一致。90%工況溫度場4thF、5thF臨界轉(zhuǎn)速重合，對應(yīng)不平衡響應(yīng)峰值疊加。

通過穩(wěn)態(tài)不平衡響應(yīng)分析求解得到不同工況溫度場作用下各支點(diǎn)不平衡響應(yīng)數(shù)據(jù)，利用式（26）計(jì)算得到高壓轉(zhuǎn)子4thF、5thF臨界轉(zhuǎn)速共振峰值列于表3與表4，最大工作轉(zhuǎn)速附近（215～225 Hz）的平均振動幅值列于表5。

表3 典型溫度場下雙轉(zhuǎn)子4thF共振峰值變化Table 3 Dual-rotor 4thF resonance peak variation in typical temperature field

表4 典型溫度場下雙轉(zhuǎn)子5thF共振峰值變化Table 4 Dual-rotor 5thF resonance peak variation in typical temperature field

表5 典型溫度場下雙轉(zhuǎn)子系統(tǒng)平均振動幅值Table 5 Average vibration amplitude of dual-rotor system under typical temperature field

通過表3與表4可以看出，在不同工況溫度場影響下，4thF臨界轉(zhuǎn)速峰值在4?、5?、6?支承位置降低，5thF臨界轉(zhuǎn)速對應(yīng)振動響應(yīng)峰值在低壓壓氣機(jī)附近的1?、2?、3?支承位置升高。具體降低或升高的幅度受兩階臨界轉(zhuǎn)速的位置、接近程度及振型變化影響。

由表5看出，隨著雙轉(zhuǎn)子系統(tǒng)工況升高，4thF、5thF臨界轉(zhuǎn)速逐漸靠近工作轉(zhuǎn)速區(qū)（圖13），在最大工作轉(zhuǎn)速附近，各支承節(jié)點(diǎn)對應(yīng)的平均振動幅值普遍上升，100%工況溫度場達(dá)到最大，最大平均不平衡響應(yīng)幅值較常溫下增大2.7倍。此外，由于90%工況溫度場4thF、5thF臨界轉(zhuǎn)速重合，不平衡響應(yīng)峰值相對80%工況、100%工況較大（表4），導(dǎo)致1?、2?及3?支承位置節(jié)點(diǎn)平均振動幅值在90%工況溫度場達(dá)到了最大值，最大平均不平衡響應(yīng)幅值較常溫下增大4.7倍。

雙轉(zhuǎn)子系統(tǒng)各部分最大應(yīng)變能峰值集中于4thF臨界轉(zhuǎn)速附近，繪制4thF臨界轉(zhuǎn)速各部分應(yīng)變能峰值隨雙轉(zhuǎn)子系統(tǒng)溫度場工況的變化如圖14所示?？梢钥闯觯鞑糠謶?yīng)變能主要受兩臨界轉(zhuǎn)速接近程度影響，在臨界轉(zhuǎn)速發(fā)生重合的90%工況溫度場，低壓轉(zhuǎn)子應(yīng)變能達(dá)到最高值，高壓轉(zhuǎn)子、支承及雙轉(zhuǎn)子總應(yīng)變能達(dá)到最低值。

圖14 4thF臨界轉(zhuǎn)速對應(yīng)各部分應(yīng)變能隨溫度場變化Fig.14 Variation of strain energy of each part corresponding to 4thF critical speed with temperature field

對雙轉(zhuǎn)子系統(tǒng)中高、低壓轉(zhuǎn)子進(jìn)行應(yīng)變能分析。繪制出的應(yīng)變能等高線圖（見圖15）能清楚地表示應(yīng)變能隨轉(zhuǎn)頻及位置的變化關(guān)系，圖中，橫軸代表對應(yīng)單元軸向坐標(biāo)，縱軸代表高壓轉(zhuǎn)子轉(zhuǎn)頻，等高線圖中顏色梯度對應(yīng)應(yīng)變能數(shù)值。

圖15 不同工況溫度場下雙轉(zhuǎn)子系統(tǒng)高壓轉(zhuǎn)子、低壓轉(zhuǎn)子應(yīng)變能分布分析對比Fig.15 Analysis and comparison of strain energy distribution of high-pressure and low-pressure rotors of dual-rotor system in different temperature fields

通過觀察圖15中低壓轉(zhuǎn)子上應(yīng)變能等高線圖發(fā)現(xiàn)，隨著4thF、5thF臨界轉(zhuǎn)速的接近，對應(yīng)的應(yīng)變能峰值區(qū)域逐漸靠近，4thF臨界轉(zhuǎn)速附近應(yīng)變能分布受5thF臨界轉(zhuǎn)速影響。到90%工況溫度場時(shí)，4thF、5thF臨界轉(zhuǎn)速重合，對應(yīng)的應(yīng)變能峰值區(qū)域融為一體，振型發(fā)生改變。

通過觀察圖15中高壓轉(zhuǎn)子上應(yīng)變能等高線圖發(fā)現(xiàn)，只有4thF臨界轉(zhuǎn)速能激起高壓轉(zhuǎn)子上的振動，5thF臨界轉(zhuǎn)速對應(yīng)振型只與低壓轉(zhuǎn)子相關(guān)（與圖9（a）所示振型結(jié)果一致）。因此隨著雙轉(zhuǎn)子系統(tǒng)溫度場工況的上升，對應(yīng)4thF臨界轉(zhuǎn)速逐漸下降，高壓轉(zhuǎn)子上應(yīng)變能峰值逐漸逼近工作轉(zhuǎn)速區(qū)。

4 結(jié)論

本文考慮了不同工況溫度場變化進(jìn)而導(dǎo)致材料屬性變化對雙轉(zhuǎn)子系統(tǒng)動力特性的影響，提出將不同工況溫度、轉(zhuǎn)速變化與發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)動力特性聯(lián)合分析方法，將穩(wěn)定工況溫度場與雙轉(zhuǎn)子有限元模型聯(lián)合，建立了典型發(fā)動機(jī)雙轉(zhuǎn)子支承系統(tǒng)熱-固聯(lián)合分析有限元模型，并對常溫、熱-固聯(lián)合分析有限元模型進(jìn)行了模態(tài)及穩(wěn)態(tài)不平衡響應(yīng)分析，結(jié)果表明：

1） 隨著雙轉(zhuǎn)子系統(tǒng)工況的上升，各階進(jìn)動頻率下降，其中渦輪附近的彎曲振型對應(yīng)的進(jìn)動頻率下降最為明顯，正、反進(jìn)動頻率最大下降幅達(dá)到約8%與27%。臨界轉(zhuǎn)速隨之降低，高于最大工作轉(zhuǎn)速的臨界轉(zhuǎn)速逐漸逼近工作轉(zhuǎn)速區(qū)，臨界轉(zhuǎn)速最大下降幅度接近10%，而發(fā)動機(jī)彎曲型正進(jìn)動臨界轉(zhuǎn)速最小裕度為20%，溫致材料屬性變化對雙轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速的影響不容忽視。

2） 隨著施加溫度場工況的上升，受臨界轉(zhuǎn)速位置變化等因素影響，高壓轉(zhuǎn)子最大轉(zhuǎn)速附近平均響應(yīng)幅值普遍上升，高壓轉(zhuǎn)子與低壓渦輪相關(guān)支承位置節(jié)點(diǎn)最大轉(zhuǎn)速附近平均不平衡響應(yīng)幅值較常溫下增大近3倍，溫致材料屬性變化對雙轉(zhuǎn)子不平衡響應(yīng)的影響較大。

3） 受溫度場變化的影響，在特定工作區(qū)域內(nèi)可能出現(xiàn)兩階臨界轉(zhuǎn)速接近甚至重合的現(xiàn)象，對應(yīng)的應(yīng)變能峰值區(qū)域融為一體，振型發(fā)生改變，影響該轉(zhuǎn)速區(qū)域及附近轉(zhuǎn)子的不平衡振動。

本文所提方法及分析結(jié)果可為發(fā)動機(jī)雙轉(zhuǎn)子系統(tǒng)精細(xì)化動力學(xué)設(shè)計(jì)、振動故障分析等提供一定的參考。