轉子支承系統臨界轉速概率分析的隨機響應面法

左彥飛， 江志農， 馮 坤， 王建軍

(1. 北京化工大學 機電工程學院，北京 100029； 2. 北京航空航天大學 能源與動力工程學院，北京 100191)

工程中的轉子支承系統存在許多非確定性因素，主要包括：部件在公差范圍內的加工誤差，材料的不均勻，部件間裝配的隨機性，長期工作導致的磨損，熱致幾何變形等[1]。這些不確定因素可致使設計相同的轉子系統臨界轉速產生大幅變化。例如，某型發動機不同臺次臨界轉速相差可達10%[2]，即使同一臺發動機，重新拆裝一次，其臨界轉速和振動特性也可能發生很大變化，這給發動機結構可靠性設計及振動控制帶來一系列問題。研究轉子系統臨界轉速在不確定因素影響下的概率分布規律，在設計上可盡量避免轉子系統工作在臨界轉速附近，提高系統可靠性；在振動控制上可通過控制加工或裝配過程中的敏感參數，減小振動故障發生的可能性，因而具有重要意義。

目前較為成熟的概率分析方法主要面向葉片、輪盤等核心部件優化設計[3-4]，而與轉子動力學結合的轉子系統級別的不確定性分析方法近幾年剛受到關注，主要使用的方法有蒙特卡洛法、區間分析法和隨機多項式展開法(Polynomial Chaos Expansion)。其中，蒙特卡洛法直接對不確定因素進行抽樣計算，計算量巨大，特別是使用三維實體有限元轉子模型時，消耗的計算資源和時間往往無法承受。區間分析法是在統計信息不足以描述非確定參數的概率分布或隸屬函數時使用的方法，用以獲得響應的區間范圍[5-6]。隨機多項式展開在轉子不確定性穩態響應分析中應用較多[7-10]，不過，該方法直接對動力學方程進行隨機多項式展開，主要適用于可以直接由動力學方程求解的問題。而對于不能直接由系統方程求解的問題，特別是由Campbell圖求系統臨界轉速時，隨機多項式展開法適用性較差。

為此，本文直接從臨界轉速與轉子系統力學不確定參數的隱式函數關系入手，將可靠性分析使用的隨機響應面法引入轉子支承系統的臨界轉速概率分析中，提出基于隨機響應面的轉子支承系統臨界轉速概率分析方法。

1 非確定因素在轉子動力學方程中的表達

轉子系統中非確定因素主要包括轉子材料屬性的不確定性、部件幾何加工的不確定性、連接結構剛度的不確定性等。由于三維實體有限元轉子模型能準確模擬系統的主要結構和力學特征，易于考慮離心力和陀螺效應的影響，應用越來越廣泛[11-14]，因此，本文采用三維實體有限元格式建立具有非確定因素的轉子動力學方程。

1.1 部件材料屬性的不確定性

材料屬性的不確定性來源于部件材料本身的不均勻甚至缺陷或在加工處理過程中產生的特性變化。設系統中存在l種獨立材料屬性，隨機變量ξi是某種符合材料密度分布規律的隨機變量，如標準正態分布，則第i個部件材料的密度可以表示為

ρi=ρ0i+ξiρσi, (i=1,2,…,l)

(1)

式中：ρ0i為密度的均值；ρσi為密度的標準差。類似的設ξl+i為符合材料彈性模量分布規律的隨機變量，則材料的彈性模量可以表示為

Ei=E0i+ξl+iEσi, (i=1,2,…,l)

(2)

式中：ξ0i為彈性模量的均值；Eσi為彈性模量的標準差。將式(1)代入參考文獻[15]中基于三維實體有限元的轉子系統動能表達式，將式(2)代入勢能表達式，進而得到質量矩陣，陀螺矩陣及剛度矩陣非確定形式

(3)

1.2 部件幾何加工的不確定性

在部件的加工過程中，不可避免存在尺寸公差和形位公差，這些公差使得實際部件與理想部件在幾何形狀上存在差異。如圖1所示。

圖1 理想結構的幾何形狀與實際結構的幾何形狀Fig.1 Ideal structure geometry and actual structure geometry

而在工程實際中，圖1(b)所示幾何加工誤差通常是在公差范圍內的誤差，對幾何形狀的影響非常小。因此，本文僅將幾何加工誤差的不確定性簡化為不平衡量的不確定性。

將上述轉子結構的幾何加工不確定性和“1.1”節中密度不確定性統一簡化為轉子上某些特定位置(動平衡修正面或轉子各個盤上)的外加不平衡量。假設有m個不平衡量(不平衡質量與偏心量的積)施加位置，則設ξ2l+1～ξ2l+m是符合不平衡量大小分布特性的隨機變量，ξ2l+m+1～ξ2l+2m是符合各個不平衡量相位分布的隨機變量，則各個位置上的不平衡量激勵可以表示為

Ω2fi(ξ2l+i)exp(jΩt+α(ξ2l+m+i)),
(i=1,2,…,m)

(4)

式中：fi(ξ2l+i)為i位置的隨機不平衡量的大小；α(ξ2l+m+i)為該不平衡量相對于零相位的相位角。

1.3 連接部位剛度不確定性

轉子支承系統中的連接關系主要包括螺栓連接、軸承等。這些連接結構剛度的不確定性主要來源于安裝狀態的隨機性，轉子系統運轉過程連接界面的變化，軸承游隙的不確定性等[16]。假設轉子系統有n處主要的不確定性連接關系，設隨機變量ξ2l+2m+1～ξ2l+2m+n是符合連接剛度大小分布規律的隨機變量，則第i處連接剛度的不確定性表示為

Ki(ξ2l+2m+i)=K0i+ξ2l+2m+iKσi, (i=1,2,…,n)

(5)

式中：K0i為連接結構剛度矩陣的平均值；Kσi為連接結構剛度矩陣的標準差。

此外，連接界面的變化一般也會對轉子系統的不平衡量造成影響，本文將這種影響一并納入式(4)中予以考慮。

1.4 考慮不確定性的轉子動力學方程的一般形式

將上述所有隨機變量用向量的形式表示，即令

ξ={ξ1ξ2…ξN}T, (N=2l+2m+n)

(6)

則考慮不確定因素的轉子動力學方程的一般形式為

(7)

在雙轉子系統中，忽略式(7)中的阻尼項，得到對應的無阻尼雙轉子支承系統的動力學方程的一般形式為

(8)

式中：M(ξ)為系統質量矩陣；f(ξ)為系統激勵向量；Ω1，Ω2分別為轉子1、轉子2的轉速；G1，G2分別為轉子1、轉子2旋轉產生的陀螺效應矩陣；K1(ξ)，K2分別為與轉子1、轉子2轉速相關的剛度矩陣，可簡稱為旋轉剛度矩陣；K0(ξ)為與轉速無關的系統剛度矩陣。此外，還存在其他不確定性因素，例如熱致幾何變形等，所有這些因素，在轉子系統動力學方程中最終表現為兩大類。

第一類是造成轉子系統本身變化的不確定性因素。例如，材料的不確定性對轉子動力學方程中的質量矩陣，陀螺矩陣，剛度矩陣都有影響。這些因素導致轉子系統的固有特性(包括臨界轉速)存在不確定性。

第二類是造成轉子系統所受激勵變化的不確定性因素。例如，幾何加工的不確定性(簡化處理)，密度的不均勻，裝配的隨機性等，都對轉子不平衡量的大小、相位產生不確定性影響，但簡化處理后并不影響轉子系統的固有特性。

2 基于隨機響應面的臨界轉速分析方法

2.1 具有不確定性的雙轉子系統臨界轉速

以無阻尼雙轉子系統臨界轉速不確定分析為例，闡述臨界轉速概率分析方法。式(8)的齊次形式為

(9)

該方程需要變換到狀態向量空間進行求解。設轉子系統的自由度為n， 引入2n維狀態向量

(10)

將式(9)改寫為狀態方程

(11)

其中，

設其解的形式為

y=ψeλt

(12)

將式(12)代入式(11)得

A(ξ)ψ=λB(ξ)ψ

(13)

求解式(13)可得n對共軛復特征值和n對共軛復特征向量

(14)

式中：σi為特征值的實部，表征衰減系數；ωi為特征值虛部，表征系統的固有頻率； j為虛數單位；φ為原方程即式(13)的模態振型。在具有不確定性的雙轉子支承系統中，由于旋轉的影響，λi，ψi的實部和虛部都是轉子轉速Ω1，Ω2及隨機變量ξ的函數。

若已知Ω1與Ω2的函數關系

Ω2=f(Ω1)

(15)

可將式(14)雙轉子系統的固有頻率表示為函數向量的形式

ω(Ω1,f(Ω1),ξ)={ω1ω2…ωn}T
(ω1≤ω2≤…≤ωn)

(16)

實際中，基于數值計算并不能得到ω中各固有頻率函數的解析表達。在確定性分析中，一般根據需要，選擇在關心的轉速范圍內求解若干組(Ω1,f(Ω1))(或(f-1(Ω2),Ω2))對應的ω中各固有頻率的函數值，將不同轉速下的特定振型所對應的固有頻率值平滑連線，近似得到各固有頻率的函數曲線。Campbell圖方法就是以Ω1(或Ω2) 為橫坐標，以固有頻率(正頻率)值為縱坐標，將各個固有頻率函數隨Ω1的變化曲線畫出，做激勵線

ω=Ω1或ω=Ω2

(17)

交點即為臨界轉速。然而，由于式(16)中固有頻率函數同時受不確定因素ξ影響，通過直接求交點的方式獲取臨界轉速概率分布規律的計算量和分析難度都很大。

2.2 轉子系統臨界轉速的隨機響應面

為了提高臨界轉速概率特性分析效率，本文引入了可靠性分析中的隨機響應面分析方法。由臨界轉速求解過程及式(16)、式(17)可知，臨界轉速與隨機變量ξ存在隱式函數關系。數學上可以證明[17]，某一階臨界轉速可以表示為

(18)

式中：ai1…ir為隨機多項式系數；Γp(·)為括號內變量的p階隨機多項式。式(18)中，當p趨于∞時，等式兩邊相等。為了便于表達，對于N維隨機變量，可以改寫為

(19)

(20)

式中：M為保留的項數，可由

(21)

2.3 隨機響應面系數的擬合

對于隨機響應面法，通常采用線性無關概率配點法選取輸入樣本，該方法所需的配點數目只需要等于待定系數數目即可得到滿意的結果[18]。此外，中心復合抽樣(Central Composite Design，CCD)，Box-Behnken Matrix Sampling(BBM)，也可用于隨機響應面法的系數擬合，不過其基本要求是樣本數量大于未知系數數目。關于隨機響應面法的系數擬合及抽樣方法，相關研究已經比較充分，可以通過通用有限元程序或MATLAB方便實現。

3 典型雙轉子支承系統臨界轉速概率分析實例

本節以典型雙轉子支承系統臨界轉速概率分析為例，驗證隨機響應面法在實際應用中的精確性和高效性。

3.1 轉子系統臨界轉速的隨機響應面

某典型雙轉子支承系統的模型如圖2所示。在該轉子支承系統中，軸承位置的剛度，套齒連軸器等的連接剛度都存在非確定性。根據工程分析需要，假設1號軸承位置的水平支承剛度K1y(ξ1)， 2號軸承位置豎直支承剛度K2z(ξ2)， 剛性套齒連軸器的彎曲剛度K65(ξ3)存在非確定性， 且ξ1,ξ2,ξ3均符合標準正態分布，K1y(ξ1)，K2z(ξ2)的標準差為均值的5%，K65(ξ3)的標準差為均值的10%。

圖2 典型雙轉子支承系統Fig.2 Typical dual rotor support system

3.2 臨界轉速概率分析隨機響應面法準確性驗證

利用隨機響應面法依據式(20)低壓轉子激勵(ω=Ω1)的某階臨界轉速可以表示為

(22)

本文采用中心復合抽樣法擬合上述未知系數，在示例中，不確定參數個數為3，隨機多項式階次為2，共10個未知系數，依據參考文獻[19]以及ANSYS中概率分析樣本選取方法，共需要15個樣本點，因此，只需對雙轉子支承系統進行15次確定性分析，利用“2.1”節中介紹的Campbell圖法得到15個樣本點對應的該階臨界轉速值，然后進行非線性擬合，得到ac0～ac9的值，最后利用式(22)進行蒙特卡洛抽樣，對該階臨界轉述進行概率特性分析。為了驗證隨機響應面法計算結果的準確性，同時使用基于蒙特卡洛法對低壓轉子激勵的第1階正進動和第2階正進動臨界轉速的不確定性進行分析。圖3所示為低壓轉子激勵第1階正進動臨界轉速(1stF)頻數分布直方圖，其中圖3(a)為由蒙特卡洛方法(Monte Carlo Simulation，MCS)1000次抽樣得到的計算結果，圖3(b)為5 000次隨機響應面(Stochastic Response Surface Method，SRSM)抽樣計算結果。由于蒙特卡洛法耗費計算時間太長，在抽樣1 000次后所得臨界轉速平均值和方差均已收斂，所以此處并未進行5 000次抽樣計算。圖中紅色曲線為擬合得到的標準正態分布曲線，兩種方法計算得到的第1階臨界轉速的分布規律都符合正態分布。圖4所示為低壓轉子激勵第2階正進動臨界轉速(2ndF)頻數分布直方圖，第2階正進動臨界轉速也符合正態分布。不過，第2階正進動臨界轉速的變化范圍較大，也就是標準差較大。表1所示為蒙特卡洛法與隨機響應面法計算得到的臨界轉速的均值和標準差的比較。

圖3 第1階正進動臨界轉速(1st F)頻次分布直方圖Fig.3 1st order positive precession critical speed (1st F) frequency distribution histogram

圖4 第2階正進動臨界轉速(2nd F)頻次分布直方圖Fig.4 2nd order positive precession critical speed (2nd F) frequency distribution histogram

從表1可知，利用隨機響應面法得到的臨界轉速均值與方差與蒙特卡洛法得到的計算結果一致。在所假設的隨機變量的影響下，第1階正進動臨界轉速變化較小，標準差約為1 r/min，根據正態分布的概率分布密度函數可知，低壓轉子激勵第1階正進動臨界轉速在(2 269±3) r/min的概率約為99.73%，變化范圍非常小。而對于第2階正進動臨界轉速，其標準差約為50 r/min，轉子第2階臨界轉速在(5 513±150) r/min的概率約為99.73%。

表1 蒙特卡洛法與隨機響應面法計算臨界轉速結果比較

3.3 臨界轉速概率分析隨機響應面法的高效性

上述概率振動特性分析實例，已經充分驗證了隨機響應面法應用于轉子系統臨界轉速分析的準確性，下面比較隨機響應面法與蒙特卡洛法的計算效率。表2所示為使用有限元模型，依據參考文獻[20]獲得的減縮模型及隨機響應面法進行1 000次蒙特卡洛抽樣計算所占用的計算機內存及消耗的計算時間。直接用原模型進行1 000次蒙特卡洛抽樣計算臨界轉速，穩態不平衡響應，所消耗的時間預計為136.9萬 s(15.8 d)，使用減縮模型后，時間縮短到1.78萬 s(4.9 h)，而利用減縮模型加隨機響應面法分別只需要267 s(4.5 min)。而上述計算使用的是同一臺計算機，可見，隨機響應面法可以大幅提高分析效率。

表2 隨機響應面法與蒙特卡洛法計算效率比較

4 結 論

針對工程中具有不確定性轉子系統臨界轉速概率分析的需要，本文提出了一種基于隨機響應面的轉子系統臨界轉速概率分析方法。并應用該方法對典型雙轉子支承系統臨界轉速的概率分布規律進行了分析，驗證了方法的準確性和高效性。相關結論主要有以下幾點：

(1) 轉子系統部件材料的不確定性影響轉子系統的質量矩陣、陀螺矩陣和剛度矩陣，幾何加工的不確定性經過簡化后，可以近似為對轉子系統不平衡量的影響，連接關系的不確定性主要表現為對系統剛度矩陣的影響，同時材料不均勻及連接關系的變化也會對不平衡量的大小和分布產生影響。最終，依據作用效果，可將所有不確定因素分為系統不確定性和激勵不確定性兩類。

(2) 基于隨機響應面法的轉子臨界轉速概率特性分析利用隨機響應面直接擬合轉子系統中不確定性輸入參數與臨界轉速的隱式函數關系，利用所得到的函數關系表達式替代系統的有限元模型進行蒙特卡洛分析。對典型雙轉子支承系統臨界轉速的概率分布計算結果表明，該方法與利用蒙特卡洛方法分析得到結果一致。

(3) 該方法只需要利用有限元模型進行少量樣本計算，相比于利用有限元模型進行的蒙特卡洛分析，隨機響應面法的分析效率得到大幅提升。

不過，隨機響應面法收斂性得以保證的前提是所有隨機變量都必須符合某種特定的分布，才可能找到與這種分布對應的正交多項式作為隨機多項式的基。例如與高斯分布對應的是Hermite多項式，與Gamma分布對應的是Laguerre多項式。而在實際轉子支承系統中的隨機變量不一定都符合同一種分布規律，因此該方法在復雜轉子支承系統中的應用可能會受限制。另外，本文分析過程中假設不確定因素不會對轉子的軸對稱性造成影響，該假設對于大多數軸對稱轉子系統適用，但對于非對稱轉子系統如裂紋轉子、異形轉子等并不適用，由于這類轉子系統一般都具有時變特性，需要使用時變系統動力特性分析方法進行分析，不存在確切的“臨界轉速”定義，故本文未予研究。