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

1 引言

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

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

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

2 柔性故障轉(zhuǎn)子的隨機(jī)模型

2.1 具有初始軸彎曲的不平衡柔性轉(zhuǎn)子方程

(1)

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

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

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

Fig.1 A general rotor system

圖2 轉(zhuǎn)子不平衡故障

Fig.2 Rotor unbalance fault

(2)

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

于是，在不平衡和軸彎曲故障作用下，忽略定子影響，可得柔性轉(zhuǎn)子梁元有限元?jiǎng)恿W(xué)方程為

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

(3)

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

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

(4)

(5)

圖3 轉(zhuǎn)子軸彎曲故障

Fig.3 Rotor shaft bent fault

2.3 隨機(jī)響應(yīng)模型

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

(6)

3 自適應(yīng)稀疏多項(xiàng)式混沌展開

3.1 響應(yīng)量的廣義PC近似

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

(7)

式中cα是PC展開系數(shù)，為待求量；Ψα(·)是廣義PC基函數(shù)；α∈M是標(biāo)識(shí)PC基函數(shù)的M維多重指標(biāo)，而M為X的元素個(gè)數(shù)。對(duì)式(7)進(jìn)行截?cái)啵?/p>

(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展開誤差估計(jì)

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

(11)

3.3 基于LAR的稀疏PC展開

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

(12)

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

圖4 基于LAR的自適應(yīng)稀疏PC展開的算法流程

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

3.4 廣義PC近似的后處理

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

(13)

Si代表了變量Xi(i=1,2,…,M)對(duì)隨機(jī)響應(yīng)方差D的相對(duì)貢獻(xiàn)大小。總靈敏度指標(biāo)為

(14)

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

4 算例分析與討論

4.1 問題描述

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

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

(15)

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

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

4.2 確定性動(dòng)力學(xué)分析

圖6給出了均值故障和特定故障工況下，轉(zhuǎn)子系統(tǒng)在圓盤1處的軸心軌跡長半軸頻響曲線。圖7 繪制了兩種故障在一階正渦動(dòng)臨界轉(zhuǎn)速時(shí)圓盤1處的軸心軌跡情況。對(duì)比兩種工況在臨界轉(zhuǎn)速675.12 r/min處的長半軸大小(注：圖6中標(biāo)記的長半軸為675 r/min時(shí)的數(shù)據(jù)，分別為2.536 mm和4.435 mm)，可見當(dāng)故障參數(shù)偏離均值時(shí)，系統(tǒng)共振穩(wěn)態(tài)響應(yīng)幅值會(huì)發(fā)生顯著改變，因而在基于共振穩(wěn)態(tài)響應(yīng)的轉(zhuǎn)子系統(tǒng)設(shè)計(jì)中考慮故障參數(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 坎貝爾圖與前四階臨界轉(zhuǎn)速

Fig.5 Campbell diagram and the first four critical speeds

4.3 響應(yīng)隨機(jī)分析

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

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

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

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

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

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

表3 不同PC階數(shù)時(shí)OLS算法和LAR算法的 PC展開項(xiàng)數(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ù)對(duì)數(shù)譜圖

Fig.9 Logarithm spectrum of PC coefficients for OLS algorithm

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

Fig.10 Logarithm spectrum of PC coefficients for LAR algorithm

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

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

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

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

表4 基于Monte Carlo仿真和基于PC近似的Sobol指標(biāo)

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 結(jié) 論

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

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

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

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