聞 靜，李舜酩

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

0 引言

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

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

1 同步全周碰摩運動方程

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

圖1 非線性單盤跨中轉(zhuǎn)子碰摩模型

非線性轉(zhuǎn)子碰摩運動方程可表示為

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

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

圖2 轉(zhuǎn)子碰摩力模型

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

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

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

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

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

轉(zhuǎn)子變形量可表示為

式中：x（t）、y（t）為轉(zhuǎn)子圓盤軸心位置坐標(biāo)，是與時間有關(guān)的變量。

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

式中：R1為轉(zhuǎn)子圓盤半徑。

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

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

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

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

假設(shè)轉(zhuǎn)軸為簡支Euler梁，受熱后轉(zhuǎn)軸軸向熱應(yīng)力為

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

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

式中：R2為轉(zhuǎn)軸的半徑；x、y為轉(zhuǎn)軸橫截面2個正交的隨體坐標(biāo)軸；Mx、My為熱應(yīng)力對 x、y之矩。

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

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

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

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

3 計算分析

3.1 計算條件

航空發(fā)動機轉(zhuǎn)子系統(tǒng)材料選用高溫合金，計算模型由某型航空發(fā)動機轉(zhuǎn)子系統(tǒng)簡化而來，模型參數(shù)以及材料參數(shù)如下[10]：圓盤等效質(zhì)量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/℃。根據(jù)航空發(fā)動機工作轉(zhuǎn)速范圍，選取無量綱轉(zhuǎn)速Ω∈{0.1，0.2，0.3…1.4}，計算轉(zhuǎn)子無量綱徑向位移，然后在無量綱徑向位移在小于1和大于1之間變化處進行細(xì)分，得到轉(zhuǎn)子發(fā)生碰摩和脫離碰摩的轉(zhuǎn)速。

3.2 計算分析步驟

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

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

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

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

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

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

圖3 順序耦合分析

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

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

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

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

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

圖6 轉(zhuǎn)子無量綱熱撓曲隨轉(zhuǎn)速的變化曲線

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

3.4 摩擦熱對碰摩響應(yīng)的影響

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

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

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

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

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

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

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

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

3.5 結(jié)果分析

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

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

4 結(jié)論

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

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

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

（3）不考慮摩擦熱時，轉(zhuǎn)子徑向位移隨時間均勻趨向于一個值，同步全周碰摩漸進趨于較為穩(wěn)定的狀態(tài)，考慮摩擦熱時轉(zhuǎn)子徑向位移隨時間變化發(fā)生不均勻波動，摩擦熱降低了轉(zhuǎn)子同步全周碰摩的穩(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].中國電機工程學(xué)報，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]羅躍綱，張松鶴，聞邦椿.非線性彈性轉(zhuǎn)子系統(tǒng)碰摩的動態(tài)特性[J].汽輪機技術(shù)，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]袁惠群.非線性轉(zhuǎn)子動力學(xué)中的分叉與混沌運動研究概述[J].黃金學(xué)報，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]劉獻棟，李其漢.轉(zhuǎn)靜件碰摩模型及不對中轉(zhuǎn)子局部碰摩的混沌特性[J].航空動力學(xué)報，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]陶海亮，左志濤，高慶，等.基于時頻分析的裂紋轉(zhuǎn)子碰摩故障特征研究[J].推進技術(shù),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]王士敏,陸啟韶.轉(zhuǎn)子通過臨界轉(zhuǎn)速時碰摩熱效應(yīng)對振動特性的影響[J].動力學(xué)與控制學(xué)報，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]許斌，徐尉南，張文.單盤轉(zhuǎn)子的同步全周碰摩及其穩(wěn)定性分析[J].復(fù)旦學(xué)報(自然科學(xué)版)，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]徐尉南,張文,許斌.柔性轉(zhuǎn)子-柔性靜子系統(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]張華彪，陳予恕，李軍.彈性支承-剛性轉(zhuǎn)子系統(tǒng)同步全周碰摩的分岔響應(yīng)[J].應(yīng)用數(shù)學(xué)和力學(xué)，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ā)動機轉(zhuǎn)子系統(tǒng)碰摩的非線性動力學(xué)研究 [D].哈爾濱：哈爾濱工業(yè)大學(xué)，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]朱向哲，袁惠群，李東.摩擦熱效應(yīng)對航空發(fā)動機高壓轉(zhuǎn)子系統(tǒng)碰摩響應(yīng)的影響[J].推進技術(shù)，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.