軸彎曲與不平衡柔性轉子共振穩(wěn)態(tài)響應隨機分析

(華東交通大學 機電與車輛工程學院，南昌 330013)

1 引言

軸彎曲和不平衡是旋轉機械中常見的兩種產(chǎn)生同步力的轉子故障。早期，在單一或混合軸彎曲和不平衡故障轉子系統(tǒng)方面的研究主要集中在動力學建模、故障識別以及動平衡等方面[1-3]，很少從不確定性角度考慮故障影響。隨著不確定性量化理論的快速發(fā)展，在轉子系統(tǒng)的分析、設計和優(yōu)化過程中考慮各種隨機因素影響的公開文獻越來越多。從已有研究可以看出，隨機攝動法[4,5]、Monte Carlo仿真法[6]以及多項式混沌展開法[7,8]是當前轉子動力學隨機不確定性研究的主要方法。其中，基于Taylor級數(shù)展開的隨機攝動法通常無法適應大變異隨機問題，如本文要解決的隨機共振穩(wěn)態(tài)響應，而Monte Carlo仿真法又因為計算成本高而使其在復雜轉子系統(tǒng)隨機分析中的應用顯得不切實際。在這種情況下，適應性更廣的多項式混沌展開法越來越受到重視。

當前，多項式混沌PC(Polynomial Chaos)展開法在轉子動力學隨機問題中的應用多以高斯分布參數(shù)下的隨機分析為主[7,8]。在非高斯隨機分析方面，Zhou等[9]采用Nataf變換和Hermite多項式混沌展開，實現(xiàn)了具有非高斯隨機參數(shù)的某汽輪機轉子隨機臨界轉速分析。但這種方法對于和高斯分布相差較大的分布類型收斂速度慢且存在誤差[10]。相比之下，另一種基于廣義多項式混沌展開的方法利用Askey方案使其具有指數(shù)收斂性能且精度高[11]。此外，臨界轉速附近頻響函數(shù)的隨機分析往往需要較高的PC階數(shù)才能達到分析精度。Sinou等[12]指出，許多在利用傳統(tǒng)低階PC展開法分析轉子系統(tǒng)頻響函數(shù)隨機特征的研究中[7]，所獲得的臨界轉速附近的頻響函數(shù)隨機結果并不可信，并以某非對稱轉子系統(tǒng)為例詳細分析了不同PC階數(shù)下的結果近似精度；Yaghoubi等[13]提出了一種隨機頻率變換策略，以降低分析共振區(qū)頻響函數(shù)隨機性時所需的PC階數(shù)。在另一方面，雖然多項式混沌展開法適應性廣，但其在高維高階情況下存在維數(shù)災難問題，尤其是基于回歸法的多項式混沌展開。因此，自適應稀疏多項式混沌展開技術引起了較大重視并獲得了迅速發(fā)展，包括逐步回歸[14]、最小角回歸[15]和壓縮感知[16]等。其中，基于最小角回歸的自適應稀疏方案是由Blatman等[15]提出的，利用最小角回歸算法選出對模型響應影響大的回歸量以形成最佳基函數(shù)集合，是一種非常有效的自適應稀疏PC展開方案。

本文研究具有初始軸彎曲的不平衡柔性轉子系統(tǒng)的不確定性量化問題，擬在Blatman等[15]工作基礎上，將基于最小角回歸技術的自適應稀疏PC展開方案引入故障轉子系統(tǒng)的共振穩(wěn)態(tài)響應分析中，并綜合廣義PC展開、留一法交叉驗證和Sobol全局靈敏度等技術實現(xiàn)故障柔性轉子在一階臨界轉速處的共振穩(wěn)態(tài)隨機響應分析和全局靈敏度計算。

2 柔性故障轉子的隨機模型

2.1 具有初始軸彎曲的不平衡柔性轉子方程

(1)

式中Re(·)代表取復數(shù)的實部，j為虛數(shù)單位，Ω2U0i代表不平衡故障引起的作用在節(jié)點i的不平衡力矢量，是復數(shù)向量，并以角速度Ω與轉子同步轉動。

設轉子具有初始軸彎曲故障，如圖3所示。取節(jié)點i在Xbz平面內(nèi)因軸彎曲引起的線位移和角位移分別為bi和αi，且在轉動坐標系ξηz中的相位為γ，那么節(jié)點i上的軸彎曲變形在固定坐標系xyz下的位移向量可表達為

圖1 一般轉子系統(tǒng)

Fig.1 A general rotor system

圖2 轉子不平衡故障

Fig.2 Rotor unbalance fault

(2)

式中δb 0i代表因軸彎曲故障引起的節(jié)點i的位移矢量，是復數(shù)向量，并以角速度Ω與轉子同步轉動。

于是，在不平衡和軸彎曲故障作用下，忽略定子影響，可得柔性轉子梁元有限元動力學方程為

Re(Ω2U0ej Ω t)+KRe(δb 0ej Ω t)

(3)

式中U0和δb 0分別是由節(jié)點向量U0i和δb 0i組裝而成的整體列向量。

2.2 一階共振穩(wěn)態(tài)響應

(4)

(5)

圖3 轉子軸彎曲故障

Fig.3 Rotor shaft bent fault

2.3 隨機響應模型

依據(jù)不確定性傳播原理，轉子節(jié)點i處軸心軌跡長半軸的隨機響應模型函數(shù)可表達為

(6)

3 自適應稀疏多項式混沌展開

3.1 響應量的廣義PC近似

設隨機向量X各元素是相互獨立且邊緣分布類型服從Askey方案[11]。若X不滿足上述要求，則利用變換方法將其變換為合適的分布類型。設響應量Y具有有限方差，那么其廣義PC表達式為

(7)

式中cα是PC展開系數(shù)，為待求量；Ψα(·)是廣義PC基函數(shù)；α∈M是標識PC基函數(shù)的M維多重指標，而M為X的元素個數(shù)。對式(7)進行截斷：

(8)

(9)

?=(ATA)-1ATY

(10)

式中A={Ai j}=Ψj(x(i))(i=1,…,N;j=0,…,P-1)；ATA稱為信息矩陣。為了保證ATA可逆，要求樣本數(shù)N>P，通常建議N=(M-1)P[18]。

3.2 基于LOO的PC展開誤差估計

為了評估響應量的PC近似誤差，采用留一法LOO(Leave -One -Out)交叉驗證技術。相應的留一法誤差公式為[15]

(11)

3.3 基于LAR的稀疏PC展開

實踐表明,PC展開項中存在大量項系數(shù)值非常小，甚至為0的情況，可以忽略。因此，為了克服或者減輕高階高維情況下維數(shù)災難問題，構造隨機響應的稀疏PC展開表達式是一種有效的措施。PC展開的稀疏表達在一定條件下可以等效為P1規(guī)劃問題：

(12)

式(12)是一個l1優(yōu)化問題。目前，有效集算法 Active Set和投影梯度法Projected Gradient是普遍采用的兩類l1優(yōu)化算法[19]。本文選用由Blatman等[15]提出的一種有效集算法——基于最小角回歸LAR(Least Angle Regression)的稀疏多項式混沌展開方法[15]。最小角回歸算法是一種線性回歸方法，其依據(jù)回歸量和當前殘余量的相關性來選取下一個回歸量，從而將對模型響應有較大影響的回歸量選中，并最終獲得一種稀疏的PC展開表達式，具體實施過程詳見文獻[15，20]。

圖4 基于LAR的自適應稀疏PC展開的算法流程

Fig.4 Algorithm flowchart for adaptive sparse polynomial chaos expansion based on LAR procedure

3.4 廣義PC近似的后處理

一旦估計出PC系數(shù)向量，則基于式(8)可以較易獲得隨機響應量的統(tǒng)計矩和概率分布信息以及開展靈敏度和可靠性分析[14]。本文僅給出基于PC展開的Sobol全局靈敏度指標計算。傳統(tǒng)Sobol指標計算基于蒙特卡洛仿真，在大量抽樣樣本作用下直接調(diào)用模型函數(shù)，使得計算成本非常高。Marelli等[20]發(fā)現(xiàn)Sobol指標可以直接從PC展開系數(shù)中得到。其中，一階全局靈敏度指標為

(13)

Si代表了變量Xi(i=1,2,…,M)對隨機響應方差D的相對貢獻大小?？傡`敏度指標為

(14)

從式(13,14)可以看出，基于PC展開的Sobol敏感度指標僅與相應的PC展開系數(shù)有關，從而避免了基于蒙特卡洛仿真的大規(guī)模計算問題。

4 算例分析與討論

4.1 問題描述

圖1為具有初始軸彎曲的不平衡柔性轉子系統(tǒng)，設不平衡故障位于圓盤1和2處，軸彎曲故障存在于整個轉軸且呈半正弦分布(在平面Xbz中)

Xb=Bsin(zπ/L) (z∈[0,L])

(15)

式中B為半正弦軸彎曲幅值，L為轉子軸彎曲故障在軸線方向上的分布距離，這里等于轉軸長度。然而，實際轉子軸彎曲形狀在各節(jié)點的位移信息(線位移bi和角位移αi)通常不能完全掌握，尤其是由軸彎曲引起的角位移。本文設角位移αi信息未知，僅通過式(15)獲得線位移bi的信息，完整的軸彎曲位移向量δb 0采用Guyan靜態(tài)縮減變換獲得[17]。

本算例僅將轉子不平衡和軸彎曲的故障參數(shù)作為隨機變量，相應的概率信息列于表1?？紤]到轉子的軸對稱特點，將軸彎曲故障的相位γ參數(shù)取為定值0，作為其他故障參數(shù)的參照基準。

4.2 確定性動力學分析

圖6給出了均值故障和特定故障工況下，轉子系統(tǒng)在圓盤1處的軸心軌跡長半軸頻響曲線。圖7 繪制了兩種故障在一階正渦動臨界轉速時圓盤1處的軸心軌跡情況。對比兩種工況在臨界轉速675.12 r/min處的長半軸大小(注：圖6中標記的長半軸為675 r/min時的數(shù)據(jù)，分別為2.536 mm和4.435 mm)，可見當故障參數(shù)偏離均值時，系統(tǒng)共振穩(wěn)態(tài)響應幅值會發(fā)生顯著改變，因而在基于共振穩(wěn)態(tài)響應的轉子系統(tǒng)設計中考慮故障參數(shù)的離散性或者說變異性是非常有必要的。

表2 兩種確定性工況的故障參數(shù)取值

Tab.2 Value of fault parameters under two kinds of deterministic cases

故障類型U3φ3U5φ5B均值故障5.0e-4π/45.0e-4π/45.0e-6特定故障1.0e-30.05.0e-4π/41.0e-5

圖5 坎貝爾圖與前四階臨界轉速

Fig.5 Campbell diagram and the first four critical speeds

4.3 響應隨機分析

基于圖4給出的流程，設定試驗設計樣本N=500，展開階數(shù)p的范圍為3～7階，表3給出了算法循環(huán)過程中不同階數(shù)下分別基于OLS算法和LAR算法所需要的總展開項數(shù)P和相應的誤差εLOO?？梢钥闯觯捎贠LS算法對PC基函數(shù)不進行稀疏性選擇，故各展開階數(shù)下的標準截斷展開項完全保留，展開項數(shù)P取值符合式(9)。在這種情況下，誤差εLOO會隨著p的增加表現(xiàn)出先減小后增大的變化規(guī)律。按照流程中階數(shù)p的自適應選取原則，由于p=5時對應的留一法誤差是最小值，為6.75e -06，故OLS算法最終判斷最優(yōu)展開階數(shù)為p*=5，此時P=252。另外，在OLS算法中，當p=6時，由于P=462，略小于預算試驗設計樣本點數(shù)(N=500)，可以看出此時對應的誤差值迅速增大且遠大于最小值情況。而當p=7時，P=792超出了預算數(shù)N，即不滿足條件N>P，無法執(zhí)行最小二乘回歸計算，誤差值趨于無窮。

圖6 系統(tǒng)穩(wěn)態(tài)響應——圓盤1處軸心軌跡長軸半徑

Fig.6 Steady-state response of the length of semi-axis of rotor obits at disk 1

圖7 圓盤1處的軸心軌跡(臨界轉速=675.12 r/min)

Fig.7 Rotor obits at disk 1 in a critical speed of 675.12 r/min

進一步，為了對比和檢驗LAR算法構建的一階共振穩(wěn)態(tài)響應的PC近似模型精度，圖11和 圖12 分別繪制了試驗設計樣本數(shù)為N=105時的響應真實值(按式(6)計算)和PC近似模型值的散點圖和直方圖。從散點圖11可以看出，在整個范圍內(nèi)，各散點基本落在45°線上，而兩者的直方圖12的形狀和變化情況非常接近，統(tǒng)計出的前四階矩(真實模型：2.4378 mm,0.33396 mm，0.3562 mm，3.1384 mm；PC近似模型：2.4378 mm，0.33391 mm，0.3563 mm，3.1364 mm)也基本一致，誤差很小，說明構建的PC近似模型精度達到要求，能替代計算成本較高的真實模型進一步開展諸如可靠性、靈敏性以及優(yōu)化設計等方面的研究。

表3 不同PC階數(shù)時OLS算法和LAR算法的 PC展開項數(shù)和留一法誤差值

Tab.3 PC coefficients number and LOO errors for OLS algorithm and LAR algorithm under different PC order

Order p34567OLS P56126252462792εLOO1.10e-048.32e-066.75e-061.82e-02+∞LAR P5612324076109εLOO8.37e-057.63e-068.40e-066.63e-056.84e-05

圖8 LAR算法中不同階數(shù)下的εLOO值

Fig.8 LOO errors for different PC order when using LAR algorithm

圖9 OLS算法中PC系數(shù)對數(shù)譜圖

Fig.9 Logarithm spectrum of PC coefficients for OLS algorithm

圖10 LAR算法中PC系數(shù)對數(shù)譜圖

Fig.10 Logarithm spectrum of PC coefficients for LAR algorithm

圖11N=105時共振穩(wěn)態(tài)響應的真實值和PC近似模型值的散點圖

Fig.11 Scatter diagram of the resonance steady-state responses of true model and PC approximation model withN=105

圖12N=105時共振穩(wěn)態(tài)響應的真實值和PC近似值的直方圖及統(tǒng)計矩

Fig.12 Histogram and statistical moments of the resonance steady-state response withN=105

表4 基于Monte Carlo仿真和基于PC近似的Sobol指標

Fig.4 Sobol sensitivity index when using Monte Carlo simulation and PC approximation method

VariablesU3U5φ3φ5BST,MCi0.2116250.2127930.1532560.1531090.329254ST,PCi0.2025880.2026220.1516990.1512230.323888SMCi0.2010630.1960970.1126920.1171040.317131SPCi0.1990780.1990910.1256540.1251770.320256

5 結 論

本文以具有初始軸彎曲的不平衡柔性轉子的隨機響應和全局靈敏度分析為目標，構建了基于轉子動力學梁元有限元理論的轉子軸心軌跡長半軸共振穩(wěn)態(tài)響應的模型函數(shù)，并綜合廣義多項式混沌展開、留一法交叉驗證技術和最小角回歸等實現(xiàn)了共振穩(wěn)態(tài)響應的自適應稀疏PC展開，獲得了PC近似模型，在驗證了方法有效性、精度和效率的情況下達到了分析目標。主要結論有：

(1) 因軸彎曲和不平衡故障屬于同步類故障，易激發(fā)轉子正向渦動，再加之柔性轉子系統(tǒng)常以過一階臨界轉速的響應性能為設計目標，故提出了以一階正渦動臨界轉速下的共振穩(wěn)態(tài)響應作為系統(tǒng)關鍵響應量。同時，考慮到穩(wěn)態(tài)響應下軸心軌跡長半軸是衡量轉子渦動范圍或判斷碰磨故障的有效參量，最終以軸心軌跡長半軸作為一階共振穩(wěn)態(tài)響應量用于后續(xù)轉子系統(tǒng)的隨機分析和靈敏度計算。算例的確定性分析結果表明，故障參數(shù)的改變會對共振穩(wěn)態(tài)響應產(chǎn)生較大影響。

(2) 共振響應離散性大以及可能存在的非光滑性往往要求PC近似應具有較高的展開階數(shù)，但會造成維數(shù)災難，為了避免盲目設定展開階數(shù)以及減少計算成本，利用一種基于留一法交叉驗證和最小角回歸技術的自適應稀疏廣義PC展開方法，實現(xiàn)了轉子系統(tǒng)在非高斯隨機故障參數(shù)作用下一階共振穩(wěn)態(tài)響應的PC近似。算例中隨機分析結果表明，基于LAR算法的PC近似不僅可以自適應確定展開階數(shù)，而且相比于OLS算法具有更少的PC展開階數(shù)，可以預期其在更少的試驗設計樣本數(shù)下比OLS算法具有更優(yōu)的近似精度。

(3) 針對算例情況，以真實響應模型結果為參照，建立的PC近似模型具有較高的近似精度，兩者的輸出響應直方圖、前四階矩以及Sobol全局靈敏度指標都非常接近。Sobol靈敏度結果表明，柔性轉子中的軸彎曲幅值故障參數(shù)對一階共振穩(wěn)態(tài)響應的方差貢獻最為明顯，而不平衡故障中不平衡量大小的離散性比其相位的離散性對響應方差貢獻大。