航空發(fā)動機非線性轉子系統(tǒng)同步全周碰摩的溫度影響分析

聞 靜，李舜酩

（南京航空航天大學 能源與動力學院，南京210016）

0 引言

轉子與機匣系統(tǒng)是航空發(fā)動機的核心部件，轉子和靜子在摩擦過程中產生的熱沖擊效應和轉子干摩擦失穩(wěn)是最具危害性的常見故障之一。當轉子和靜子發(fā)生接觸時，摩擦熱的生成使得轉子接觸表面局部溫度升高，使轉子受熱不均而產生熱撓曲，從而改變轉子不平衡量的分布和振動情況。航空發(fā)動機轉子處于高轉速的工作狀態(tài)，一旦接觸，即使在短時間內也很容易發(fā)生全周碰摩，此時熱撓曲將會在轉子的一定方向不斷加劇，從而造成新的不平衡，引起更嚴重的碰摩。Dimarogonas等[1]在對軸封和轉子的摩擦升溫的研究中發(fā)現(xiàn)，當摩擦情況較為嚴重時，轉子和軸封的溫度變化會達到十分可觀的程度。因此，碰摩及其熱效應帶來的影響不容忽視。目前，對于連續(xù)摩擦熱沖擊的問題已經有一些研究[2-3]，對于碰摩力作用的轉子系統(tǒng)非線性碰摩響應的研究也比較廣泛[4-7]，但在摩擦熱對轉子振動響應的影響方面，僅針對點碰或局部碰摩進行了研究[8]。近年來，轉子同步全周碰摩響應問題開始得到越來越多的學者關注[9-12]。

本文以單盤轉子為研究對象，在考慮摩擦熱效應引起轉子熱撓曲情況下，建立轉子與機匣之間的同步全周碰摩接觸過程模型，采用數(shù)值計算研究了轉子的摩擦生熱和摩擦熱效應對其振動特性的影響，并與未考慮摩擦熱效應的碰摩響應進行了比較。計算結果可為航空發(fā)動機轉子動力學設計提供一定的理論參考。

1 同步全周碰摩運動方程

在多種碰摩響應中，同步全周碰摩最容易發(fā)生在轉子的工作轉速附近。由于在現(xiàn)代轉子設計中開始考慮越來越多的非線性因素，而轉子的碰摩是典型的非光滑動力學問題，轉子發(fā)生碰摩之前和發(fā)生碰摩時的約束條件不同，碰摩力的表達形式也不同，因此本文建立1個考慮轉子剛度存在非線性的單盤跨中轉子碰摩模型，如圖1所示。轉子在圓盤處的恢復力和位移之間存在三次非線性關系，恢復力表示為F=kr+αr3，其中k和α為轉子的線性和非線性剛度系數(shù)為轉子在圓盤處的徑向位移。

圖1 非線性單盤跨中轉子碰摩模型

非線性轉子碰摩運動方程可表示為

式中：c為系統(tǒng)阻尼；e為轉子偏心距；ω為轉子自轉角速度；Fx和 Fy為 x、y 方向碰摩力。

系統(tǒng)的碰摩力由線性接觸力Fn和庫倫摩擦力Ft組成，如圖2所示，分解到x、y 軸上分別為 Fx、Fy

圖2 轉子碰摩力模型

式中：r0為轉靜間間隙；μ為庫倫摩擦系數(shù)；kc為接觸剛度；νr=ωrdisc+ωwr，為碰摩點處的相對速度，ωw為轉子進動角速度，rdisc為轉子圓盤半徑。

νr＞0 時，摩擦力方向與轉子轉動方向相反；νr＜0時，摩擦力方向與轉子轉動方向相同。當轉子發(fā)生同步全周碰摩時，恒有νr＞0，因此令sign（νr）=1。

將式（1）無量綱化，定義可得

式中和FY分別為無量綱形式x、y方向的碰摩力

2 碰摩處溫度分布及熱彎曲

轉子變形量可表示為

式中：x（t）、y（t）為轉子圓盤軸心位置坐標，是與時間有關的變量。

則在周期Tt內，轉子平均變形量可表示為[13]δ=轉、靜子之間的摩擦熱源可以表示為

式中：R1為轉子圓盤半徑。

計算轉、靜子摩擦生成溫度場時，假設[14]：不考慮摩擦過程中材料磨損；轉、靜子材料滿足各項同性；忽略轉、靜子外表面與空氣進行對流換熱和輻射換熱的影響。

根據熱傳導理論，由于轉子軸向熱變形很小，因此忽略軸向的變化，柱坐標下轉盤橫截面上非穩(wěn)態(tài)熱傳導方程為

初始條件為 Tt=0=T0（r，θ），邊界條件為

式中：S為接觸處面積；λ為熱傳導系數(shù)；ρ為材料密度；cp為材料比熱容；T（r,θ,t）為時刻 t轉盤截面上的溫度分布；θc為轉盤上發(fā)生碰摩位置。

假設轉軸為簡支Euler梁，受熱后轉軸軸向熱應力為

式中：χt為熱膨脹系數(shù)；E為彈性模量。

橫截面上溫度分布不均勻在該界面引起的彎矩為

式中：R2為轉軸的半徑；x、y為轉軸橫截面2個正交的隨體坐標軸；Mx、My為熱應力對 x、y之矩。

根據Castigliano定理[15]，熱撓曲可表示為

式中：L為轉軸的長度；ρt（l,t）為轉軸在距離左端點l處的撓度；φ為撓曲面和xoz面的夾角；I為抗彎模量。

由于碰摩點附近溫度在很小范圍內發(fā)生急劇變化，假設轉軸在轉盤處溫度為 T（r,θ,t），其他位置為0，由于研究對象為單盤跨中轉子，且可以看出熱撓曲在轉軸兩端為0，在碰摩點處最大，即當時最大。近似地用二次曲線形式表示轉子的熱撓曲

溫度分布由式（7）以及初始條件、邊界條件通過MATLAB數(shù)值計算得出?？梢詮氖街锌闯觯瑹釓澢S時間的變化取決于溫度隨時間的變化。

3 計算分析

3.1 計算條件

航空發(fā)動機轉子系統(tǒng)材料選用高溫合金，計算模型由某型航空發(fā)動機轉子系統(tǒng)簡化而來，模型參數(shù)以及材料參數(shù)如下[10]：圓盤等效質量m=57.28 kg，L=2 m，k=2×106N/m，α=3.8×1012N/m3，c=261.8 N·s/m，r0=0.0003 m，R1=0.4577 m，R2=0.08 m，e=0.00003 m，kc=3 ×106N/m，μ=0.15，λ=650 W/m ·℃ ，ρ=7833 kg/m3，cp=465 J/kg·℃，χT=1.2×10-5m/℃。根據航空發(fā)動機工作轉速范圍，選取無量綱轉速Ω∈{0.1，0.2，0.3…1.4}，計算轉子無量綱徑向位移，然后在無量綱徑向位移在小于1和大于1之間變化處進行細分，得到轉子發(fā)生碰摩和脫離碰摩的轉速。

3.2 計算分析步驟

由于轉子存在質量偏心，轉子徑向位移隨轉速的變化而發(fā)生變化，當轉子徑向位移達到轉、靜間隙時，轉、靜子間將會發(fā)生碰摩，將轉子徑向位移無量綱化后轉、靜子間碰摩的發(fā)生條件表示為無量綱徑向位移大于1。當轉、靜子間發(fā)生碰摩時，產生碰摩力，由于轉、靜子間發(fā)生摩擦產生熱源，轉子上將產生不均勻分布的溫度場，此溫度場將引起轉子熱撓度，從而影響轉子變形量和碰摩力，進而影響轉子的振動特性。

計算發(fā)生碰摩的轉子振動特性步驟如下：

（1）根據碰摩力計算公式（4），對運動方程（3）進行MATLAB數(shù)值積分求解得到碰摩力、轉子變形量及無摩擦熱時的轉子振動特性；

（2）將步驟（1）中得到的轉子變形量代入式（5）計算摩擦熱源，并結合邊界條件、初始條件等式（6）～（11）得到溫度分布和轉子熱撓曲；

（3）考慮熱撓度減小轉、靜子間間隙，從而改變碰摩力，采用順序耦合法將計算得到的熱撓曲耦合到運動方程中，得到考慮摩擦熱后的轉子振動特性。

順序耦合分析方法（如圖3所示）的關鍵在于先計算得到溫度場分布，然后將其作為載荷施加到結構中來計算結構的響應。采用MATLAB軟件編程實現(xiàn)不考慮和考慮摩擦熱的轉子碰摩振動特性的計算。

圖3 順序耦合分析

3.3 溫度分布和熱撓曲特征分析

用MATLAB自帶的ode45函數(shù)對運動方程進行數(shù)值積分求解，并忽略求解開始時的不穩(wěn)定部分。為研究轉子上接觸面的溫度分布，需要選取發(fā)生同步全周碰摩的轉速，初步選取 Ω=1.80、1.90、1.00、1.10、1.20、1.30，分別計算轉子系統(tǒng)運動方程的穩(wěn)定解，結果顯示Ω=1.10時轉子系統(tǒng)處于同步全周碰摩狀態(tài)。因此計算溫度分布時選取該轉速，并用countour函數(shù)繪制轉子上的2維溫度分布曲線。Ω=1.10時圓盤接觸面溫度沿徑向和周向分布曲線如圖4、5所示。從圖中可見，徑向溫度隨與接觸面的距離的增大而呈非線性降低，最高溫度達到82.5℃；周向溫度沿接觸面呈對稱分布，并且在周向角達到30°時幾乎等于室溫。

圖4 圓盤接觸面徑向溫度分布

圖5 圓盤接觸面周向溫度分布

由于轉子存在質量偏心，轉子徑向位移隨轉速的變化而變化，改變摩擦熱源的大小，從而改變摩擦熱引起的轉子熱撓度。為保證轉子處于同步全周碰摩狀態(tài)，在Ω=1.10周圍進行轉速的細分Ω∈{0.96，0.98，1.00，…，1.18}，并計算轉子系統(tǒng)運動方程的穩(wěn)定解，結果顯示在Ω=1.02和1.14之間轉子系統(tǒng)處于同步全周碰摩狀態(tài)。因此選取Ω=1.02、1.04、1.06、1.08、1.10、1.12來研究轉子熱撓度隨轉速的變化情況。熱撓曲隨轉速變化的曲線如圖6所示，縱坐標為無量綱熱撓曲

圖6 轉子無量綱熱撓曲隨轉速的變化曲線

從圖中可見，隨轉速的升高，熱撓曲呈非線性上升趨勢，且當Ω=1.12，無量綱熱撓曲達到0.242。這是由于轉速的升高增大了轉子徑向位移，因此增大了碰摩摩擦力，摩擦熱的增加增大了轉子的熱撓曲。

3.4 摩擦熱對碰摩響應的影響

在第3.3節(jié)中已得到無摩擦熱情況下轉子進入和脫離同步全周碰摩的轉速分別為1.02和1.14。為得到摩擦熱對轉子碰摩響應的影響，計算有摩擦熱情況的轉子無量綱徑向位移，以得到轉子發(fā)生和脫離同步全周碰摩的轉速。選取Ω∈{0.96，0.98，1.00，…，1.18}分別代入運動方程，用countour函數(shù)得到各轉速下的溫度分布情況，儲存為txt文件，代入式（11）中計算各轉速對應的熱撓度。將轉子熱撓度分解為x、y 2個方向得到 ρTx和 ρTy，并增加到式（4）得到考慮摩擦熱情況的碰摩力，并代入運動方程，用ode45函數(shù)進行數(shù)值分析。

不考慮摩擦熱時4種轉速下轉子的無量綱徑向位移如圖7所示。從圖中可見，在Ω=1.00、1.16時無量綱徑向位移小于1，表示此轉速下轉子沒有發(fā)生同步全周碰摩。在Ω=1.02、1.14時，轉子無量綱徑向位移大于1，表示轉子處于同步全周碰摩階段?？紤]摩擦熱時4種轉速下轉子的無量綱徑向位移如圖8所示。從圖中可見，在Ω=0.98、1.14時無量綱徑向位移小于1，表示該轉速下轉子沒有發(fā)生同步全周碰摩。在Ω=1.00、1.12時，轉子無量綱徑向位移大于1，表示轉子處于同步全周碰摩階段。

圖7 不考慮摩擦熱時Ω=1.00、1.02、1.14、1.16的轉子無量綱徑向位移

對比2種情況轉子進入和脫離同步全周碰摩的轉速，可見摩擦熱使得轉子發(fā)生和脫離同步全周碰摩的轉速降低約2%，但沒有改變發(fā)生同步全周碰摩的轉速范圍。發(fā)動機在升速過程中轉子系統(tǒng)將在更低轉速時發(fā)生同步全周碰摩，也會在更低轉速時脫離同步全周碰摩狀態(tài)。

圖8 考慮摩擦熱時Ω=0.98、1.00、1.12、1.14的轉子無量綱徑向位移

摩擦熱引起的熱撓度將改變轉子同步全周碰摩的運行狀態(tài)。為有效對比，對無摩擦熱和有摩擦熱情況選取同一轉速，取已發(fā)生碰摩的轉速Ω=1.10代入運動方程并用ode45函數(shù)進行數(shù)值分析。在Ω=1.10時無摩擦熱和有摩擦熱轉子無量綱徑向位移如圖9所示。

圖9 Ω=1.10時不考慮摩擦熱和考慮摩擦熱的轉子無量綱徑向位移

從圖中可見，考慮摩擦熱時轉子無量綱徑向位移比不考慮時增大約9%，這是由于摩擦熱產生的熱彎矩增大了轉子徑向位移，加重了碰摩程度。同時可見，不考慮摩擦熱時無量綱徑向位移均勻趨向于一個值，同步全周碰摩漸進趨于穩(wěn)定，而考慮摩擦熱時無量綱徑向位移發(fā)生不均勻波動，說明摩擦熱降低了同步全周碰摩的穩(wěn)定性。

3.5 結果分析

（1）對比有摩擦熱和無摩擦熱影響情況下轉子徑向位移隨轉速變化的情況，可以得到，有摩擦熱情況下，轉子進入和脫離同步全周碰摩的無量綱轉速均降低2%，發(fā)動機在升速過程中將在更低轉速時進入和脫離同步全周碰摩狀態(tài)。

（2）對比有摩擦熱和無摩擦熱情況下同步全周碰摩時轉子徑向位移隨時間變化的情況，可以得到，不同于無摩擦熱時轉子無量綱徑向位移均勻趨向一個值，碰摩漸進趨于穩(wěn)定，考慮摩擦熱時徑向位移發(fā)生不均勻波動，碰摩穩(wěn)定性降低，轉子同步全周碰摩程度加重。

4 結論

本文分別在有摩擦熱和無摩擦熱情況下計算了某型航空發(fā)動機轉子系統(tǒng)同步全周碰摩的轉子振動響應，通過對比和分析可得出以下3點結論。

（1）圓盤接觸面由于碰摩產生最高溫升達到82.5℃，徑向溫度隨與接觸面的距離的增加而呈非線性降低；周向溫度沿接觸面呈對稱分布，并且在周向角達到30°時幾乎等于室溫。摩擦產生的轉子熱撓曲隨轉速的增加而呈非線性增加；

（2）考慮摩擦熱時，轉子發(fā)生同步全周碰摩的轉速降低，脫離同步全周碰摩的轉速降低，同步全周碰摩發(fā)生的轉速范圍不變。轉子徑向位移增大，在Ω=1.10時增大約9%，轉子碰摩程度加重。

（3）不考慮摩擦熱時，轉子徑向位移隨時間均勻趨向于一個值，同步全周碰摩漸進趨于較為穩(wěn)定的狀態(tài)，考慮摩擦熱時轉子徑向位移隨時間變化發(fā)生不均勻波動，摩擦熱降低了轉子同步全周碰摩的穩(wěn)定性。

[1]Kennedy F E.Single pass rub phenomena analysis and experiment[J].Journal of lubrication Technology,Transactions of the ASME，1982,104：582-588.

[2]Goldman P,Muszynska A.rotor-to-stator,rub-related,thermal/mechanical effects in rotating machinery[J].Chaos,Solitons&Fractals，1995，9（5）：1579-1601.

[3]楊建剛，黃寶華，高偉.摩擦熱沖擊問題的建模、仿真與振動特性分析[J].中國電機工程學報，1999，19（6）：25-31.YANG Jiangang,HUANG Baohua,GAO Wei.Modeling simulation and vibration analysis of rub-induced thermal impact problem[J].Proceedings of the CSEE，1999，19（6）：25-31.（in Chinese）

[4]羅躍綱，張松鶴，聞邦椿.非線性彈性轉子系統(tǒng)碰摩的動態(tài)特性[J].汽輪機技術，2004，46（6）：415-417.LUO Yuegang，ZHANG Songhe，WEN Bangchun.Dynamic characteristics of nonlinear elastic rotor system with rub-impact faults[J].Turbine Technology，2004，46（6）：415-417.（in Chinese）

[5]袁惠群.非線性轉子動力學中的分叉與混沌運動研究概述[J].黃金學報，2000，2（4）：309-314.YUAN Huiqun.Summary on the study for bifurcation and chaos motion on nonlinear rotor dynamics[J].Gold Journal，2000，2（4）：309-314.（in Chinese）

[6]劉獻棟，李其漢.轉靜件碰摩模型及不對中轉子局部碰摩的混沌特性[J].航空動力學報，1998，13（4）：361-365.LIU Xiandong，LI Qihan.Impact-rub model and chaotic feature of partial rub of eccentric rotor[J].Journal of Aerospace Power，1998，13（4）：361-365.（in Chinese）

[7]陶海亮，左志濤，高慶，等.基于時頻分析的裂紋轉子碰摩故障特征研究[J].推進技術,2013，34（4）：520-528.TAO Hailiang，ZUO Zhitao，GAO Qing，et al.Fault analysis of rotor with rub-impact and crack based on time-frequency analysis[J].Journal of Propulsion Technology,2013，34（4）：520-528.（in Chinese）

[8]王士敏,陸啟韶.轉子通過臨界轉速時碰摩熱效應對振動特性的影響[J].動力學與控制學報，2004，2（3）：64-69.WANG Shimin,LU Qishao.The effect of run-impact thermal effects on vibration characteristics of a rotor passing through critical speed[J].Journal of Dynamics and Control，2004，2（3）：64-69.（in Chinese）

[9]許斌，徐尉南，張文.單盤轉子的同步全周碰摩及其穩(wěn)定性分析[J].復旦學報(自然科學版)，2006，45（2）：148-169.XU Bin,XU Weinan,ZHANG Wen.Study of synchronous full annular rub of jeffcott rotor and its dynamic stability[J].Journal of Fudan University(Natural Science)，2006，45（2）：148-169.（in Chinese）

[10]徐尉南,張文,許斌.柔性轉子-柔性靜子系統(tǒng)的同步全周碰摩分析[J].振動與沖擊，2007，26（2）：1-9.XU Weinan，ZHANG Wen，XU Bin.Analytical study of synchronous full annular rub moniton of flexible stator and rotor system[J].Journal of Vibration and Shock，2007，26（2）：1-9.（in Chinese）

[11]張華彪，陳予恕，李軍.彈性支承-剛性轉子系統(tǒng)同步全周碰摩的分岔響應[J].應用數(shù)學和力學，2012，33（7）：812-827.ZHANG Huabiao，CHEN Yushu，LI Jun.Bifurcation on the synchronous full annular rub of a rigid-rotor elastic-support system[J].Applied Mathematics and Mechanics，2012，33 （7）：812-827.（in Chinese）

[12]張華彪.航空發(fā)動機轉子系統(tǒng)碰摩的非線性動力學研究 [D].哈爾濱：哈爾濱工業(yè)大學，2012.ZHANG Huabiao.Nonlinear dynamics on impact-rub of aero engine rotor system[D].Harbin：Harbin Institute of Technology，2012.（in Chinese）

[13]Bachschmid N,Pennacchi P,Vania A.Thermally induced vibrations due to rub in real rotors[J].Journal of Sound and Vibration，2007，299(4-5)：6-9.

[14]朱向哲，袁惠群，李東.摩擦熱效應對航空發(fā)動機高壓轉子系統(tǒng)碰摩響應的影響[J].推進技術，2010，31（3）：366-371.ZHU Xiangzhe，YUAN Huiqun，LI Dong.Friction thermal effect on rubbing response of a high pressure rotor system [J].Journal of Propulsion Technology，2010，31（3）：366-371.（in Chinese）

[15]Nowinski J L.Theory of thermoelasticity with applications[J].Sijthoff&Noordhoff,Alphen ann den Rijn,1980,47（2）：459.