考慮參數(shù)非概率不確定性的碰摩轉(zhuǎn)子非線性振動響應(yīng)分析

馬新星， 張振果， 華宏星

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

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

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

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

1 系統(tǒng)模型及控制方程

1.1 轉(zhuǎn)子系統(tǒng)及碰摩力模型

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

圖1 轉(zhuǎn)子系統(tǒng)及碰摩示意圖Fig.1 Schematic diagram of the rotor system and rubbing

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

1.2 控制方程

轉(zhuǎn)子系統(tǒng)的有限元控制方程可表示為

式中:M，D及K分別為系統(tǒng)質(zhì)量、阻尼及剛度矩陣；H為Heaviside函數(shù)；fu為轉(zhuǎn)子不平衡激勵；q，q·及q··分別為位移、速度和加速度向量。

考慮系統(tǒng)周期解，根據(jù)諧波平衡方法，轉(zhuǎn)子非線性振動響應(yīng)用截斷的Fourier級數(shù)近似表示為

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

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

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

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

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

式中:δ為Dirichlet函數(shù)；Sk（ω）為系統(tǒng)動剛度矩陣。

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

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

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

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

2 基于Chebyshev多項式的區(qū)間方法

2.1 區(qū)間控制方程

本節(jié)引入基于Chebyshev包含函數(shù)的區(qū)間分析法估計不確定碰摩轉(zhuǎn)子的動力學響應(yīng)范圍。首先，利用上、下邊界已知的區(qū)間參數(shù)描述任意不確定參數(shù)，表達式為

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

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

獲得的傅里葉系數(shù)的下界與上界分別為

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

2.2 Chebyshev多項式逼近理論

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

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

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

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

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

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

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

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

3 數(shù)值結(jié)果

3.1 確定性系統(tǒng)幅頻響應(yīng)

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

圖2 確定性轉(zhuǎn)子系統(tǒng)幅頻響應(yīng)Fig.2 Determinate rotor amplitude-frequency response

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

3.2 單一不確定性轉(zhuǎn)子響應(yīng)范圍

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

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

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

表1 算例轉(zhuǎn)子系統(tǒng)計算參數(shù)Tab.1 Parameters of example rotor system

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

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

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

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

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

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

3.3 多源不確定性轉(zhuǎn)子響應(yīng)范圍

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

圖9 多參數(shù)不確定時轉(zhuǎn)子碰摩x方向幅頻響應(yīng)Fig.9 Rotor amplitude-frequency response with multiple uncertainties in the x direction

圖10 多參數(shù)不確定時轉(zhuǎn)子碰摩y方向幅頻響應(yīng)Fig.10 Rotor amplitude-frequency response with multiple uncertainties in the y direction

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

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

表2 定轉(zhuǎn)速不確定轉(zhuǎn)子系統(tǒng)計算參數(shù)Tab.2 Parameters of uncertain rotor system at fixed rotational speed

圖12 轉(zhuǎn)速1 000 r/min時轉(zhuǎn)子振動時域響應(yīng)Fig.12 Rotor vibration time-domain response at 1 000 r/min

圖13 轉(zhuǎn)速5 000 r/min時轉(zhuǎn)子振動時域響應(yīng)Fig.13 Rotor vibration time-domain response at 5 000 r/min

4 結(jié) 論

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

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

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

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