考慮參數非概率不確定性的碰摩轉子非線性振動響應分析

馬新星， 張振果， 華宏星

（1.上海交通大學 機械系統與振動國家重點實驗室，上海 200240；2.上海交通大學 振動、沖擊、噪聲研究所，上海 200240）

隨著對效率和穩定性要求的提高，燃氣輪機轉/靜子間隙越來越小，使轉/靜子碰摩故障可能性增加[1]。轉/靜碰摩可引起轉子的強非線性振動響應，使系統振動出現跳躍、分岔和混沌等行為，還可能導致轉子失穩[2-3]。目前關于轉/靜子碰摩響應的確定性研究已經相當充分，國內外學者對含有碰摩故障的轉子振動行為及機理開展了大量的數值仿真及試驗研究[4-6]，大量研究表明碰摩導致的非線性振動對轉子參數變化十分敏感。然而，對于實際燃氣輪機而言，在設計制造、安裝維修、運行操作等過程中均存在不確定性[7]，使轉子動力學特性的準確描述愈發困難，因此，將現有確定性研究進一步拓展為適當的不確定性分析更符合實際，并有著迫切的現實需求。

近年來，轉子動力學領域的不確定研究已得到逐步開展[8-9]，例如不確定轉子系統的不平衡量識別[10]、臨界轉速概率分析[11]、瞬態動平衡研究[12]，也有針對故障轉子非線性振動響應的不確定分析，如不確定轉子裂紋動力學特性分析[13]等，然而關于碰摩故障的不確定轉子振動分析卻很少見。目前針對不確定性碰摩轉子文獻中，Yang等[14]分析了具有隨機參數的碰摩轉子分岔與混沌分析，強調了隨機參數對系統分岔有重要影響；王威等[15]應用偽周期替代數據算法對不確定性激勵下碰摩轉子故障信號進行了識別與分類；Yang等[16]開展了不確定碰摩轉子的可靠性研究，采用貝葉斯技術研究了混合不確定參數對系統可靠性影響。總體而言，盡管很多學者關注轉/靜子碰摩引起的非線性振動響應，但是大部分研究都是基于確定性系統開展，尚未考慮真實系統中客觀存在的復雜不確定性因素的影響。此外，少數的不確定性研究以概率分析方法為主，然而由于實際中大多參數的概率分布函數無法先驗獲得，因此基于非概率方法的轉/靜子碰摩動力學研究有待深入開展。

針對上述研究的不足，本文考慮轉子不確定參數的分布范圍，無需參數先驗概率分布，將基于Chebyshev多項式的非概率區間方法應用到碰摩轉子非線性幅頻響應的全局分析，碰摩轉子的非線性振動響應評估擴展為基于代理模型的迭代分析，可在有界參數空間內快速量化碰摩轉子振動響應的區間分布特征，進而探究不確定參數對轉子碰摩響應的影響機理。

1 系統模型及控制方程

1.1 轉子系統及碰摩力模型

燃氣輪機雙盤柔性轉子系統的簡化模型，如圖1所示。圖1中包含兩個非對稱分布的圓盤和兩組軸承支承，系統左端為聯軸器。在模型中，轉軸利用Timoshenko梁模擬，圓盤和聯軸器假設為集中質量點并假定圓盤存在相同質量偏心，軸承采用剛度-阻尼模型等效。本文僅考慮左側圓盤在水平x方向發生定點碰摩的情形，隨著轉子渦動，轉子x方向的位移可能超出轉/靜子間的初始間隙r0，轉/靜子接觸并在法向產生碰撞力fn，在切向誘發摩擦力ft。基于Hertz接觸與庫倫摩擦理論，定點碰摩力的數學表達為

圖1 轉子系統及碰摩示意圖Fig.1 Schematic diagram of the rotor system and rubbing

式中:kn為法向接觸剛度；xd為轉子左盤水平方向位移；μ為摩擦因數。

1.2 控制方程

轉子系統的有限元控制方程可表示為

式中:M，D及K分別為系統質量、阻尼及剛度矩陣；H為Heaviside函數；fu為轉子不平衡激勵；q，q·及q··分別為位移、速度和加速度向量。

考慮系統周期解，根據諧波平衡方法，轉子非線性振動響應用截斷的Fourier級數近似表示為

式中:H為諧波數；Qk（k=0，…，N）為對應于全部自由度的傅里葉系數矢量；Re為實部；ω為轉子轉速；?為Kronecker積；IN為單位矩陣。

將非線性碰摩力、不平衡力也進行傅里葉級數展開，分別為

由于圓盤不平衡力為簡諧激勵，其幅值直接施加在控制方程一階諧波項，無需特別轉換。此外，由于碰摩力是位移的函數，在后續求根方程迭代求解時采用時頻域轉換技術。

將式（4）、式（6）和式（7）代入式（3），并應用傅里葉-伽遼金映射，將系統運動方程轉化為頻域內N維非線性代數方程組。將不平衡激勵移到方程左端，可表示為如下殘差形式

式中:?0=diag[0，1，…，H]；R為方程殘差矢量。對于復指數表達的任意第k階諧波平衡等式

式中:δ為Dirichlet函數；Sk（ω）為系統動剛度矩陣。

為提高非線性響應求解效率，考慮碰摩力具有局部非線性，本文進一步采用縮減方法降低方程組維度。首先根據線性自由度與非線性自由度對位移矢量與碰摩力矢量進行劃分

由于碰摩力僅作用在左盤對應節點自由度并依賴于左盤位移，則式（10）可重新表示為

式中，Hk=S-1k為動柔度矩陣，即第k階諧波動剛度矩陣的逆。上述方程的前n行（非線性自由度維數）可以獨立求解，綜合所有諧波項，導出維度僅為n（2H+1）維方程

進而采用Newton-Raphson迭代方法求解方程，迭代過程為

2 基于Chebyshev多項式的區間方法

2.1 區間控制方程

本節引入基于Chebyshev包含函數的區間分析法估計不確定碰摩轉子的動力學響應范圍。首先，利用上、下邊界已知的區間參數描述任意不確定參數，表達式為

式中:zI為參數區間；β為該參數的不確定度；zm，z-及分別為參數中值、下界和上界。所有不確定參數區間可寫為和分別為不確定參數的下界與上界矢量，其中j為不確定參數數目。

因此，在不確定參數向量的約束條件下，通過殘差方程式（14）求得的某激勵頻率下第k階傅里葉系數重新表述為

獲得的傅里葉系數的下界與上界分別為

由于泰勒級數展開適合于小不確定量的問題，為克服區間算術運算可能導致動態響應范圍過估計，本文采用Chebyshev多項式的區間方法構造代理模型，以獲得更高精度和效率。

2.2 Chebyshev多項式逼近理論

已知任意不確定參數區間z∈[a，b]可以線性轉換為區間η∈[-1，1]，考慮在單參數標準區間上的第k階諧波傅里葉系數函數可由截斷階數為m的Chebyshev級數近似為

式中:fi為常值Chebyshev系數；Ci為Chebyshev多項式，Ci（η）=cos（iθ）， θ=arccos（η）∈[0，π]，i≥1，且Chebyshev多項式滿足正交性

Chebyshev系數可寫成積分形式并由Gauss-Chebyshev數值積分求得

式中，插值點η或者θh為m+1階Chebyshev多項式的根

針對多個參數區間問題，則需要將一維Chebyshev包含函數擴展到多維情形

式中:l為i1，…，ij中0出現的次數； η=[-1，1]j為j維區間變量。同樣，多維Chebyshev系數也可通過Gauss-Chebyshev積分計算得出

式中，η1，…，ηj=0，1，2…。

將式（23）或式（26）分別代入式（21）便可建立基于Chebyshev多項式的單參數或多參數的區間函數逼近表達式或代理模型。

3 數值結果

3.1 確定性系統幅頻響應

本文首先進行確定性系統的幅頻響應特性分析，所采用的轉子模型參考2013年Ma等的研究，轉子確定性的物理參數如下:轉軸直徑為10 mm，總長為607.5 mm，左盤距離左側軸端215 mm，兩盤相距200 mm；兩相同圓盤與聯軸器的直徑轉動慣量分別為2.48×10-4kg·m2，3.19×10-6kg·m2、極轉動慣量分別4.73×10-4kg·m2，2.96×10-6kg·m2；軸承剛度各項異性，左右支撐水平方向剛度分別為1×105N·m，2×108N·m，豎直方向剛度分別為2×105N·m，5×108N·m。利用前述的模型降階技術與諧波平衡法可快速計算得到碰摩導致的轉子幅頻響應，并與Newmark-β數值方法結果進行了對比，結果如圖2所示。

圖2 確定性轉子系統幅頻響應Fig.2 Determinate rotor amplitude-frequency response

由圖2可知，轉子碰撞力使轉子在x方向的幅頻響應表現為硬特性，并出現幅值“跳躍”現象，而摩擦力使轉子y方向的幅頻響應表現為軟特性，在一階共振峰之后出現一個略微向左傾斜的峰。基于單點碰摩假設，本文摩擦力對轉子幅頻響應的影響僅限于y方向，不考慮摩擦因數不確定性影響。然而有研究表明，當出現轉子局部或者全周碰摩時，界面摩擦可能會誘導轉子反向渦動自激振動，幅值突然增大，嚴重破壞轉子穩定運行，因此在其他碰摩形式中，也應考慮摩擦因數的不確定性。

3.2 單一不確定性轉子響應范圍

本節采用區間數學方法依次分析間隙、接觸剛度、不平衡量等參數不確定性對轉子非線性幅頻響應的影響規律，不確定參數范圍如表1所示，其他參數保持不變。為了驗證結果準確性及所提出方法的高效性，對每組參數同時開展Monte Carlo模擬（樣本量均為1 000）。圖3～圖8為得到的不確定動態響應，其中灰色區域為Monte Carlo模擬的全部樣本結果，實線和點劃線分別為上、下界和均值曲線，虛線為區間計算結果的上、下界，區間均值曲線為雙劃線。對比可知，區間算法得到的響應上、下界與Monte Carlo仿真結果一致，均值曲線也比較接近，證明了區間方法的準確性。此外，Chebyshev多項式選取的展開階數及積分點數均為5，由于區間方法樣本點很少，每組單參數區間分析計算時間只有1.5 min，而1 000組Monte Carlo模擬需要4.6 h，由此說明了該方法具有較高的計算效率。由圖3可知，轉子間隙不確定性改變了所考慮轉速范圍內轉靜子碰摩發生點頻率，間隙越小，該轉子在升速過程中越早與靜子發生碰摩，共振峰向右偏移更明顯，但由于此時轉速相對較低，其全局響應最大幅值比間隙較大碰摩時小，這也說明合理利用碰摩在一定程度上可以限制柔性轉子過臨界時振動響應的最大峰值。此外，轉子不確定響應上、下界并不關于均值曲線對稱。圖4明顯的反映出間隙不確定性導致摩擦方向幅值差異，引起穩態共振曲線呈現軟式特征，并形成共振峰包絡，其隨摩擦力的增大向左偏移。

圖3 間隙不確定時轉子碰摩x方向幅頻響應Fig.3 Rotor amplitude-frequency response with gap uncertainty in the x direction

圖4 間隙不確定時轉子碰摩y方向幅頻響應Fig.4 Rotor amplitude-frequency response with gap uncertainty in the y direction

表1 算例轉子系統計算參數Tab.1 Parameters of example rotor system

圖5和圖6反映出接觸剛度的不確定主要影響系統發生碰摩后“軟”“硬”共振峰的偏移。由于接觸剛度與轉靜子材料有關，在過臨界轉速時，如果碰摩無法避免，其他參數不變時，可以通過改變靜子局部材料或者加涂層，來改變碰摩后轉子幅頻響應峰值。

圖5 接觸剛度不確定時轉子碰摩x方向幅頻響應Fig.5 Rotor amplitude-frequency response with contact stiffness uncertainty in the x direction

圖6 接觸剛度不確定時轉子碰摩y方向幅頻響應Fig.6 Rotor amplitude-frequency response with contact stiffness uncertainty in the y direction

圖3～圖6表明，無論是轉子的x方向幅頻響應，還是y方向幅頻響應，與碰摩直接相關的參數如初始轉靜子間隙、接觸剛度對共振峰的偏移度及邊界范圍均產生明顯影響，對其他轉速下振動幅值則影響微弱，而轉子不平衡量與轉速有關，其不確定性作用則表現在轉子的全局范圍，使轉子幅頻響應成為“包絡帶”，如圖7和圖8所示。由全部結果可以看出，每組參數的不確定性均對碰摩轉子幅頻響應有明顯影響，確定性的響應幅頻曲線演化為響應區間。

圖7 不平衡不確定時轉子碰摩在x方向幅頻響應Fig.7 Rotor amplitude-frequency response with unbalance uncertainty in the x direction

圖8 不平衡不確定時轉子碰摩y方向幅頻響應Fig.8 Rotor amplitude-frequency response with unbalance uncertainty in the y direction

3.3 多源不確定性轉子響應范圍

本節考慮三個物理參數同時具有不確定性時碰摩轉子的幅頻響應特性，計算結果如圖9和圖10所示。由圖9和圖10可知，多源不確定性綜合了三種參數不確定性單獨作用及參數間相互作用的情況，幅頻響應的可能取值范圍變的更廣，在共振轉速區，不確定性的作用尤其明顯，這體現了實際工程分析中考慮不確定因素的必要性。值得注意的是，盡管Chebyshev多項式區間算法對于控制過估計比較有效，但是當多考慮多參數不確定性時，相對于Monte Carlo模擬，仍然會一定程度的過估計，圖9和圖10中可以看出區間分析結果上、下界可以包裹Monte Carlo模擬全部樣本。

圖9 多參數不確定時轉子碰摩x方向幅頻響應Fig.9 Rotor amplitude-frequency response with multiple uncertainties in the x direction

圖10 多參數不確定時轉子碰摩y方向幅頻響應Fig.10 Rotor amplitude-frequency response with multiple uncertainties in the y direction

本節進一步結合Ma等的試驗結果闡述考慮不確定性參數影響的必要性。2013年，Ma等研究中轉子試驗所測得的軸心軌跡，如圖11所示，可以看出轉子渦動包含多條軌線，除了試驗干擾信號外，測試結果可重復性并不特別理想，表明了不確定參數的存在，2009年，Ma等研究中軸心軌跡多次試驗結果所表現的不確定尤其明顯。為此，針對2組特定轉速進行轉子不確定振動響應分析，所采用的不確定參數如表2所示。時域內碰摩轉子軌跡及波形圖的上、下邊界的分析結果，如圖12和圖13所示，盡管不確定參數的范圍無法由試驗直接得到，但響應仿真結果的趨勢相似。由于參數不確定性的存在，轉子的軸心軌跡不再保持在固定的軌道，而是在一定范圍內波動，響應的時域波形幅值也表現為包絡形式，由此說明考慮不確定性因素的響應分析可以更好地解釋給定試驗條件下測試結果中轉子軌跡差異。

圖11 試驗軸心軌跡Fig.11 Experimental rotor trajectories

表2 定轉速不確定轉子系統計算參數Tab.2 Parameters of uncertain rotor system at fixed rotational speed

圖12 轉速1 000 r/min時轉子振動時域響應Fig.12 Rotor vibration time-domain response at 1 000 r/min

圖13 轉速5 000 r/min時轉子振動時域響應Fig.13 Rotor vibration time-domain response at 5 000 r/min

4 結 論

本文考慮碰摩轉子非概率不確定參數的影響，基于區間分析理論和模型降階的諧波平衡方法提出了碰摩轉子系統區間不確定性問題的分析和求解模型，探討了各區間變量對非線性動力學特性的影響，主要結論如下:

（1）與傳統Monte Carlo模擬和概率不確定性分析相比，區間分析可克服需要參數先驗概率分布的苛刻要求，并且具有較高的求解效率。

（2）間隙不確定性明顯改變碰摩發生點頻率，間隙越小，碰摩越早發生；接觸剛度不確定性主要引起碰摩幅頻響應共振峰的偏移，對其他頻率下的振動響應幅值影響微弱；不平衡量的不確定性會使碰撞方向硬特性改變，并在全局范圍內引起轉子振動響應的不確定波動。

（3）參數非概率分布區間分析法的應用可以在有界參數空間內快速計算得到碰摩轉子的動力學響應范圍，為探究轉子參數對振動響應影響機理提供了一種高效高魯棒性的新途徑。