左彥飛, 江志農(nóng), 馮 坤, 王建軍
(1. 北京化工大學(xué) 機(jī)電工程學(xué)院,北京 100029; 2. 北京航空航天大學(xué) 能源與動力工程學(xué)院,北京 100191)
工程中的轉(zhuǎn)子支承系統(tǒng)存在許多非確定性因素,主要包括:部件在公差范圍內(nèi)的加工誤差,材料的不均勻,部件間裝配的隨機(jī)性,長期工作導(dǎo)致的磨損,熱致幾何變形等[1]。這些不確定因素可致使設(shè)計相同的轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速產(chǎn)生大幅變化。例如,某型發(fā)動機(jī)不同臺次臨界轉(zhuǎn)速相差可達(dá)10%[2],即使同一臺發(fā)動機(jī),重新拆裝一次,其臨界轉(zhuǎn)速和振動特性也可能發(fā)生很大變化,這給發(fā)動機(jī)結(jié)構(gòu)可靠性設(shè)計及振動控制帶來一系列問題。研究轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速在不確定因素影響下的概率分布規(guī)律,在設(shè)計上可盡量避免轉(zhuǎn)子系統(tǒng)工作在臨界轉(zhuǎn)速附近,提高系統(tǒng)可靠性;在振動控制上可通過控制加工或裝配過程中的敏感參數(shù),減小振動故障發(fā)生的可能性,因而具有重要意義。
目前較為成熟的概率分析方法主要面向葉片、輪盤等核心部件優(yōu)化設(shè)計[3-4],而與轉(zhuǎn)子動力學(xué)結(jié)合的轉(zhuǎn)子系統(tǒng)級別的不確定性分析方法近幾年剛受到關(guān)注,主要使用的方法有蒙特卡洛法、區(qū)間分析法和隨機(jī)多項(xiàng)式展開法(Polynomial Chaos Expansion)。其中,蒙特卡洛法直接對不確定因素進(jìn)行抽樣計算,計算量巨大,特別是使用三維實(shí)體有限元轉(zhuǎn)子模型時,消耗的計算資源和時間往往無法承受。區(qū)間分析法是在統(tǒng)計信息不足以描述非確定參數(shù)的概率分布或隸屬函數(shù)時使用的方法,用以獲得響應(yīng)的區(qū)間范圍[5-6]。隨機(jī)多項(xiàng)式展開在轉(zhuǎn)子不確定性穩(wěn)態(tài)響應(yīng)分析中應(yīng)用較多[7-10],不過,該方法直接對動力學(xué)方程進(jìn)行隨機(jī)多項(xiàng)式展開,主要適用于可以直接由動力學(xué)方程求解的問題。而對于不能直接由系統(tǒng)方程求解的問題,特別是由Campbell圖求系統(tǒng)臨界轉(zhuǎn)速時,隨機(jī)多項(xiàng)式展開法適用性較差。
為此,本文直接從臨界轉(zhuǎn)速與轉(zhuǎn)子系統(tǒng)力學(xué)不確定參數(shù)的隱式函數(shù)關(guān)系入手,將可靠性分析使用的隨機(jī)響應(yīng)面法引入轉(zhuǎn)子支承系統(tǒng)的臨界轉(zhuǎn)速概率分析中,提出基于隨機(jī)響應(yīng)面的轉(zhuǎn)子支承系統(tǒng)臨界轉(zhuǎn)速概率分析方法。
轉(zhuǎn)子系統(tǒng)中非確定因素主要包括轉(zhuǎn)子材料屬性的不確定性、部件幾何加工的不確定性、連接結(jié)構(gòu)剛度的不確定性等。由于三維實(shí)體有限元轉(zhuǎn)子模型能準(zhǔn)確模擬系統(tǒng)的主要結(jié)構(gòu)和力學(xué)特征,易于考慮離心力和陀螺效應(yīng)的影響,應(yīng)用越來越廣泛[11-14],因此,本文采用三維實(shí)體有限元格式建立具有非確定因素的轉(zhuǎn)子動力學(xué)方程。
材料屬性的不確定性來源于部件材料本身的不均勻甚至缺陷或在加工處理過程中產(chǎn)生的特性變化。設(shè)系統(tǒng)中存在l種獨(dú)立材料屬性,隨機(jī)變量ξi是某種符合材料密度分布規(guī)律的隨機(jī)變量,如標(biāo)準(zhǔn)正態(tài)分布,則第i個部件材料的密度可以表示為
ρi=ρ0i+ξiρσi, (i=1,2,…,l)
(1)
式中:ρ0i為密度的均值;ρσi為密度的標(biāo)準(zhǔn)差。類似的設(shè)ξl+i為符合材料彈性模量分布規(guī)律的隨機(jī)變量,則材料的彈性模量可以表示為
Ei=E0i+ξl+iEσi, (i=1,2,…,l)
(2)
式中:ξ0i為彈性模量的均值;Eσi為彈性模量的標(biāo)準(zhǔn)差。將式(1)代入?yún)⒖嘉墨I(xiàn)[15]中基于三維實(shí)體有限元的轉(zhuǎn)子系統(tǒng)動能表達(dá)式,將式(2)代入勢能表達(dá)式,進(jìn)而得到質(zhì)量矩陣,陀螺矩陣及剛度矩陣非確定形式
(3)
在部件的加工過程中,不可避免存在尺寸公差和形位公差,這些公差使得實(shí)際部件與理想部件在幾何形狀上存在差異。如圖1所示。
圖1 理想結(jié)構(gòu)的幾何形狀與實(shí)際結(jié)構(gòu)的幾何形狀Fig.1 Ideal structure geometry and actual structure geometry
而在工程實(shí)際中,圖1(b)所示幾何加工誤差通常是在公差范圍內(nèi)的誤差,對幾何形狀的影響非常小。因此,本文僅將幾何加工誤差的不確定性簡化為不平衡量的不確定性。
將上述轉(zhuǎn)子結(jié)構(gòu)的幾何加工不確定性和“1.1”節(jié)中密度不確定性統(tǒng)一簡化為轉(zhuǎn)子上某些特定位置(動平衡修正面或轉(zhuǎn)子各個盤上)的外加不平衡量。假設(shè)有m個不平衡量(不平衡質(zhì)量與偏心量的積)施加位置,則設(shè)ξ2l+1~ξ2l+m是符合不平衡量大小分布特性的隨機(jī)變量,ξ2l+m+1~ξ2l+2m是符合各個不平衡量相位分布的隨機(jī)變量,則各個位置上的不平衡量激勵可以表示為
Ω2fi(ξ2l+i)exp(jΩt+α(ξ2l+m+i)),
(i=1,2,…,m)
(4)
式中:fi(ξ2l+i)為i位置的隨機(jī)不平衡量的大??;α(ξ2l+m+i)為該不平衡量相對于零相位的相位角。
轉(zhuǎn)子支承系統(tǒng)中的連接關(guān)系主要包括螺栓連接、軸承等。這些連接結(jié)構(gòu)剛度的不確定性主要來源于安裝狀態(tài)的隨機(jī)性,轉(zhuǎn)子系統(tǒng)運(yùn)轉(zhuǎn)過程連接界面的變化,軸承游隙的不確定性等[16]。假設(shè)轉(zhuǎn)子系統(tǒng)有n處主要的不確定性連接關(guān)系,設(shè)隨機(jī)變量ξ2l+2m+1~ξ2l+2m+n是符合連接剛度大小分布規(guī)律的隨機(jī)變量,則第i處連接剛度的不確定性表示為
Ki(ξ2l+2m+i)=K0i+ξ2l+2m+iKσi, (i=1,2,…,n)
(5)
式中:K0i為連接結(jié)構(gòu)剛度矩陣的平均值;Kσi為連接結(jié)構(gòu)剛度矩陣的標(biāo)準(zhǔn)差。
此外,連接界面的變化一般也會對轉(zhuǎn)子系統(tǒng)的不平衡量造成影響,本文將這種影響一并納入式(4)中予以考慮。
將上述所有隨機(jī)變量用向量的形式表示,即令
ξ={ξ1ξ2…ξN}T, (N=2l+2m+n)
(6)
則考慮不確定因素的轉(zhuǎn)子動力學(xué)方程的一般形式為
(7)
在雙轉(zhuǎn)子系統(tǒng)中,忽略式(7)中的阻尼項(xiàng),得到對應(yīng)的無阻尼雙轉(zhuǎn)子支承系統(tǒng)的動力學(xué)方程的一般形式為
(8)
式中:M(ξ)為系統(tǒng)質(zhì)量矩陣;f(ξ)為系統(tǒng)激勵向量;Ω1,Ω2分別為轉(zhuǎn)子1、轉(zhuǎn)子2的轉(zhuǎn)速;G1,G2分別為轉(zhuǎn)子1、轉(zhuǎn)子2旋轉(zhuǎn)產(chǎn)生的陀螺效應(yīng)矩陣;K1(ξ),K2分別為與轉(zhuǎn)子1、轉(zhuǎn)子2轉(zhuǎn)速相關(guān)的剛度矩陣,可簡稱為旋轉(zhuǎn)剛度矩陣;K0(ξ)為與轉(zhuǎn)速無關(guān)的系統(tǒng)剛度矩陣。此外,還存在其他不確定性因素,例如熱致幾何變形等,所有這些因素,在轉(zhuǎn)子系統(tǒng)動力學(xué)方程中最終表現(xiàn)為兩大類。
第一類是造成轉(zhuǎn)子系統(tǒng)本身變化的不確定性因素。例如,材料的不確定性對轉(zhuǎn)子動力學(xué)方程中的質(zhì)量矩陣,陀螺矩陣,剛度矩陣都有影響。這些因素導(dǎo)致轉(zhuǎn)子系統(tǒng)的固有特性(包括臨界轉(zhuǎn)速)存在不確定性。
第二類是造成轉(zhuǎn)子系統(tǒng)所受激勵變化的不確定性因素。例如,幾何加工的不確定性(簡化處理),密度的不均勻,裝配的隨機(jī)性等,都對轉(zhuǎn)子不平衡量的大小、相位產(chǎn)生不確定性影響,但簡化處理后并不影響轉(zhuǎn)子系統(tǒng)的固有特性。
以無阻尼雙轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速不確定分析為例,闡述臨界轉(zhuǎn)速概率分析方法。式(8)的齊次形式為
(9)
該方程需要變換到狀態(tài)向量空間進(jìn)行求解。設(shè)轉(zhuǎn)子系統(tǒng)的自由度為n, 引入2n維狀態(tài)向量
(10)
將式(9)改寫為狀態(tài)方程
(11)
其中,
設(shè)其解的形式為
y=ψeλt
(12)
將式(12)代入式(11)得
A(ξ)ψ=λB(ξ)ψ
(13)
求解式(13)可得n對共軛復(fù)特征值和n對共軛復(fù)特征向量
(14)
式中:σi為特征值的實(shí)部,表征衰減系數(shù);ωi為特征值虛部,表征系統(tǒng)的固有頻率; j為虛數(shù)單位;φ為原方程即式(13)的模態(tài)振型。在具有不確定性的雙轉(zhuǎn)子支承系統(tǒng)中,由于旋轉(zhuǎn)的影響,λi,ψi的實(shí)部和虛部都是轉(zhuǎn)子轉(zhuǎn)速Ω1,Ω2及隨機(jī)變量ξ的函數(shù)。
若已知Ω1與Ω2的函數(shù)關(guān)系
Ω2=f(Ω1)
(15)
可將式(14)雙轉(zhuǎn)子系統(tǒng)的固有頻率表示為函數(shù)向量的形式
ω(Ω1,f(Ω1),ξ)={ω1ω2…ωn}T
(ω1≤ω2≤…≤ωn)
(16)
實(shí)際中,基于數(shù)值計算并不能得到ω中各固有頻率函數(shù)的解析表達(dá)。在確定性分析中,一般根據(jù)需要,選擇在關(guān)心的轉(zhuǎn)速范圍內(nèi)求解若干組(Ω1,f(Ω1))(或(f-1(Ω2),Ω2))對應(yīng)的ω中各固有頻率的函數(shù)值,將不同轉(zhuǎn)速下的特定振型所對應(yīng)的固有頻率值平滑連線,近似得到各固有頻率的函數(shù)曲線。Campbell圖方法就是以Ω1(或Ω2) 為橫坐標(biāo),以固有頻率(正頻率)值為縱坐標(biāo),將各個固有頻率函數(shù)隨Ω1的變化曲線畫出,做激勵線
ω=Ω1或ω=Ω2
(17)
交點(diǎn)即為臨界轉(zhuǎn)速。然而,由于式(16)中固有頻率函數(shù)同時受不確定因素ξ影響,通過直接求交點(diǎn)的方式獲取臨界轉(zhuǎn)速概率分布規(guī)律的計算量和分析難度都很大。
為了提高臨界轉(zhuǎn)速概率特性分析效率,本文引入了可靠性分析中的隨機(jī)響應(yīng)面分析方法。由臨界轉(zhuǎn)速求解過程及式(16)、式(17)可知,臨界轉(zhuǎn)速與隨機(jī)變量ξ存在隱式函數(shù)關(guān)系。數(shù)學(xué)上可以證明[17],某一階臨界轉(zhuǎn)速可以表示為
(18)
式中:ai1…ir為隨機(jī)多項(xiàng)式系數(shù);Γp(·)為括號內(nèi)變量的p階隨機(jī)多項(xiàng)式。式(18)中,當(dāng)p趨于∞時,等式兩邊相等。為了便于表達(dá),對于N維隨機(jī)變量,可以改寫為
(19)
(20)
式中:M為保留的項(xiàng)數(shù),可由
(21)
對于隨機(jī)響應(yīng)面法,通常采用線性無關(guān)概率配點(diǎn)法選取輸入樣本,該方法所需的配點(diǎn)數(shù)目只需要等于待定系數(shù)數(shù)目即可得到滿意的結(jié)果[18]。此外,中心復(fù)合抽樣(Central Composite Design,CCD),Box-Behnken Matrix Sampling(BBM),也可用于隨機(jī)響應(yīng)面法的系數(shù)擬合,不過其基本要求是樣本數(shù)量大于未知系數(shù)數(shù)目。關(guān)于隨機(jī)響應(yīng)面法的系數(shù)擬合及抽樣方法,相關(guān)研究已經(jīng)比較充分,可以通過通用有限元程序或MATLAB方便實(shí)現(xiàn)。
本節(jié)以典型雙轉(zhuǎn)子支承系統(tǒng)臨界轉(zhuǎn)速概率分析為例,驗(yàn)證隨機(jī)響應(yīng)面法在實(shí)際應(yīng)用中的精確性和高效性。
某典型雙轉(zhuǎn)子支承系統(tǒng)的模型如圖2所示。在該轉(zhuǎn)子支承系統(tǒng)中,軸承位置的剛度,套齒連軸器等的連接剛度都存在非確定性。根據(jù)工程分析需要,假設(shè)1號軸承位置的水平支承剛度K1y(ξ1), 2號軸承位置豎直支承剛度K2z(ξ2), 剛性套齒連軸器的彎曲剛度K65(ξ3)存在非確定性, 且ξ1,ξ2,ξ3均符合標(biāo)準(zhǔn)正態(tài)分布,K1y(ξ1),K2z(ξ2)的標(biāo)準(zhǔn)差為均值的5%,K65(ξ3)的標(biāo)準(zhǔn)差為均值的10%。
圖2 典型雙轉(zhuǎn)子支承系統(tǒng)Fig.2 Typical dual rotor support system
利用隨機(jī)響應(yīng)面法依據(jù)式(20)低壓轉(zhuǎn)子激勵(ω=Ω1)的某階臨界轉(zhuǎn)速可以表示為
(22)
本文采用中心復(fù)合抽樣法擬合上述未知系數(shù),在示例中,不確定參數(shù)個數(shù)為3,隨機(jī)多項(xiàng)式階次為2,共10個未知系數(shù),依據(jù)參考文獻(xiàn)[19]以及ANSYS中概率分析樣本選取方法,共需要15個樣本點(diǎn),因此,只需對雙轉(zhuǎn)子支承系統(tǒng)進(jìn)行15次確定性分析,利用“2.1”節(jié)中介紹的Campbell圖法得到15個樣本點(diǎn)對應(yīng)的該階臨界轉(zhuǎn)速值,然后進(jìn)行非線性擬合,得到ac0~ac9的值,最后利用式(22)進(jìn)行蒙特卡洛抽樣,對該階臨界轉(zhuǎn)述進(jìn)行概率特性分析。為了驗(yàn)證隨機(jī)響應(yīng)面法計算結(jié)果的準(zhǔn)確性,同時使用基于蒙特卡洛法對低壓轉(zhuǎn)子激勵的第1階正進(jìn)動和第2階正進(jìn)動臨界轉(zhuǎn)速的不確定性進(jìn)行分析。圖3所示為低壓轉(zhuǎn)子激勵第1階正進(jìn)動臨界轉(zhuǎn)速(1stF)頻數(shù)分布直方圖,其中圖3(a)為由蒙特卡洛方法(Monte Carlo Simulation,MCS)1000次抽樣得到的計算結(jié)果,圖3(b)為5 000次隨機(jī)響應(yīng)面(Stochastic Response Surface Method,SRSM)抽樣計算結(jié)果。由于蒙特卡洛法耗費(fèi)計算時間太長,在抽樣1 000次后所得臨界轉(zhuǎn)速平均值和方差均已收斂,所以此處并未進(jìn)行5 000次抽樣計算。圖中紅色曲線為擬合得到的標(biāo)準(zhǔn)正態(tài)分布曲線,兩種方法計算得到的第1階臨界轉(zhuǎn)速的分布規(guī)律都符合正態(tài)分布。圖4所示為低壓轉(zhuǎn)子激勵第2階正進(jìn)動臨界轉(zhuǎn)速(2ndF)頻數(shù)分布直方圖,第2階正進(jìn)動臨界轉(zhuǎn)速也符合正態(tài)分布。不過,第2階正進(jìn)動臨界轉(zhuǎn)速的變化范圍較大,也就是標(biāo)準(zhǔn)差較大。表1所示為蒙特卡洛法與隨機(jī)響應(yīng)面法計算得到的臨界轉(zhuǎn)速的均值和標(biāo)準(zhǔn)差的比較。
圖3 第1階正進(jìn)動臨界轉(zhuǎn)速(1st F)頻次分布直方圖Fig.3 1st order positive precession critical speed (1st F) frequency distribution histogram
圖4 第2階正進(jìn)動臨界轉(zhuǎn)速(2nd F)頻次分布直方圖Fig.4 2nd order positive precession critical speed (2nd F) frequency distribution histogram
從表1可知,利用隨機(jī)響應(yīng)面法得到的臨界轉(zhuǎn)速均值與方差與蒙特卡洛法得到的計算結(jié)果一致。在所假設(shè)的隨機(jī)變量的影響下,第1階正進(jìn)動臨界轉(zhuǎn)速變化較小,標(biāo)準(zhǔn)差約為1 r/min,根據(jù)正態(tài)分布的概率分布密度函數(shù)可知,低壓轉(zhuǎn)子激勵第1階正進(jìn)動臨界轉(zhuǎn)速在(2 269±3) r/min的概率約為99.73%,變化范圍非常小。而對于第2階正進(jìn)動臨界轉(zhuǎn)速,其標(biāo)準(zhǔn)差約為50 r/min,轉(zhuǎn)子第2階臨界轉(zhuǎn)速在(5 513±150) r/min的概率約為99.73%。
表1 蒙特卡洛法與隨機(jī)響應(yīng)面法計算臨界轉(zhuǎn)速結(jié)果比較
上述概率振動特性分析實(shí)例,已經(jīng)充分驗(yàn)證了隨機(jī)響應(yīng)面法應(yīng)用于轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速分析的準(zhǔn)確性,下面比較隨機(jī)響應(yīng)面法與蒙特卡洛法的計算效率。表2所示為使用有限元模型,依據(jù)參考文獻(xiàn)[20]獲得的減縮模型及隨機(jī)響應(yīng)面法進(jìn)行1 000次蒙特卡洛抽樣計算所占用的計算機(jī)內(nèi)存及消耗的計算時間。直接用原模型進(jìn)行1 000次蒙特卡洛抽樣計算臨界轉(zhuǎn)速,穩(wěn)態(tài)不平衡響應(yīng),所消耗的時間預(yù)計為136.9萬 s(15.8 d),使用減縮模型后,時間縮短到1.78萬 s(4.9 h),而利用減縮模型加隨機(jī)響應(yīng)面法分別只需要267 s(4.5 min)。而上述計算使用的是同一臺計算機(jī),可見,隨機(jī)響應(yīng)面法可以大幅提高分析效率。
表2 隨機(jī)響應(yīng)面法與蒙特卡洛法計算效率比較
針對工程中具有不確定性轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速概率分析的需要,本文提出了一種基于隨機(jī)響應(yīng)面的轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速概率分析方法。并應(yīng)用該方法對典型雙轉(zhuǎn)子支承系統(tǒng)臨界轉(zhuǎn)速的概率分布規(guī)律進(jìn)行了分析,驗(yàn)證了方法的準(zhǔn)確性和高效性。相關(guān)結(jié)論主要有以下幾點(diǎn):
(1) 轉(zhuǎn)子系統(tǒng)部件材料的不確定性影響轉(zhuǎn)子系統(tǒng)的質(zhì)量矩陣、陀螺矩陣和剛度矩陣,幾何加工的不確定性經(jīng)過簡化后,可以近似為對轉(zhuǎn)子系統(tǒng)不平衡量的影響,連接關(guān)系的不確定性主要表現(xiàn)為對系統(tǒng)剛度矩陣的影響,同時材料不均勻及連接關(guān)系的變化也會對不平衡量的大小和分布產(chǎn)生影響。最終,依據(jù)作用效果,可將所有不確定因素分為系統(tǒng)不確定性和激勵不確定性兩類。
(2) 基于隨機(jī)響應(yīng)面法的轉(zhuǎn)子臨界轉(zhuǎn)速概率特性分析利用隨機(jī)響應(yīng)面直接擬合轉(zhuǎn)子系統(tǒng)中不確定性輸入?yún)?shù)與臨界轉(zhuǎn)速的隱式函數(shù)關(guān)系,利用所得到的函數(shù)關(guān)系表達(dá)式替代系統(tǒng)的有限元模型進(jìn)行蒙特卡洛分析。對典型雙轉(zhuǎn)子支承系統(tǒng)臨界轉(zhuǎn)速的概率分布計算結(jié)果表明,該方法與利用蒙特卡洛方法分析得到結(jié)果一致。
(3) 該方法只需要利用有限元模型進(jìn)行少量樣本計算,相比于利用有限元模型進(jìn)行的蒙特卡洛分析,隨機(jī)響應(yīng)面法的分析效率得到大幅提升。
不過,隨機(jī)響應(yīng)面法收斂性得以保證的前提是所有隨機(jī)變量都必須符合某種特定的分布,才可能找到與這種分布對應(yīng)的正交多項(xiàng)式作為隨機(jī)多項(xiàng)式的基。例如與高斯分布對應(yīng)的是Hermite多項(xiàng)式,與Gamma分布對應(yīng)的是Laguerre多項(xiàng)式。而在實(shí)際轉(zhuǎn)子支承系統(tǒng)中的隨機(jī)變量不一定都符合同一種分布規(guī)律,因此該方法在復(fù)雜轉(zhuǎn)子支承系統(tǒng)中的應(yīng)用可能會受限制。另外,本文分析過程中假設(shè)不確定因素不會對轉(zhuǎn)子的軸對稱性造成影響,該假設(shè)對于大多數(shù)軸對稱轉(zhuǎn)子系統(tǒng)適用,但對于非對稱轉(zhuǎn)子系統(tǒng)如裂紋轉(zhuǎn)子、異形轉(zhuǎn)子等并不適用,由于這類轉(zhuǎn)子系統(tǒng)一般都具有時變特性,需要使用時變系統(tǒng)動力特性分析方法進(jìn)行分析,不存在確切的“臨界轉(zhuǎn)速”定義,故本文未予研究。