基于區間隨機階矩的不確定量化方法

邱繼偉,羅海勝,趙冠捷,張 亞

(中國兵器工業標準化研究所,北京 100089)

0 引言

齒輪—轉子系統是兵器、船舶、航空、航天等領域中應用最廣泛的動力傳遞裝置。隨著對齒輪—轉子系統的傳動精度、承載能力、可靠性、低噪聲等的要求越來越高,有關齒輪—轉子系統振動問題的研究也越來越重要[1-2]。目前齒輪—轉子系統動力學研究大多處理為確定性問題[3-5],軸系扭轉振動模型完全可以比較精確地對確定性分析問題進行建模和仿真,然而實際工程中由于制造、裝配、潤滑、沖擊、溫度等因素都會導致齒輪—轉子系統參數存在不確定性,而含有不確定性參數的扭振模型無論解析法還是數值方法都很難直接進行運算,無法定量描述系統響應的邊界。當前齒輪—轉子系統的動態特性不確定分析主要采用基于統計模型的分析方法[6-8],該方法可以處理含隨機變量的不確定性問題,而不確定性不等同于隨機性,統計分析方法并非研究齒輪—轉子系統不確定問題的唯一方法。另外,統計分析方法往往需要大量樣本實驗來獲取已知系統參數或激勵的統計信息,然而由于自然現象的復雜性和人們認知的局限性,有些參數很難獲得甚至無法獲得。相比統計信息,研究人員可以根據所掌握的知識和工程經驗積累情況相對容易地給出一些不確定參數的上下邊界信息,即區間信息,進而使齒輪—轉子系統振動分析中同時包含隨機信息和區間信息,因此研究同時含隨機參數和區間參數的齒輪—轉子系統的動力分析問題具有重要的理論意義和工程背景。

國內外學者對隨機不確定性傳播和區間不確定性傳播問題進行了大量研究。研究隨機不確定性傳播的手段主要是混沌多項式展開、傅里葉級數展開等。PANUNZIO等[9]將混沌多項式和漸近數值方法結合,提出非入侵式振動傳播方法,并計算了機構系統隨機共振頻率、峰值振幅;DUBREUIL等[10]針對隨機系統多學科分析問題,提出基于多項式混沌擴展的半侵入混沌多項式分析方法,有效地將這些隨機不確定性傳播到多學科分析的輸出解中;VIVEK VITTALDEV等[11]結合高斯混合模型與多項式混沌展開方法,提出基于高斯混合模型和多項式混沌展開的航天器不確定性傳播方法;余學鋒等[12]針對基于國際標準化組織ISO《測量不確定度指南》提供的測量不確定度評估方法存在的局限性,提出采用多項式混沌方法進行測量結果不確定度評估的新方法。研究區間不確定性傳播問題的手段主要是Chebyshev多項式、泰勒級數多項式、區間擴張函數等。WU[13]等將切比雪夫級數展開引入區間模型,提出一種新的區間不確定系統的動態響應分析方法,有效提高了數值解的精度;FEDELE等[14]提出一種新的基于區間逆問題的不確定性傳播方法,基于迭代思想,利用區間有限元和優化技術直接估計未知參數的邊界;TRABELSI等[15]用區間表示設計變量,在仿真步驟前集成所有不同類型的約束,以生成的時間間隔代表產品設計變量可能值的域,進而提出一種新的區間估計和約束傳播算法;尹盛文等[16]針對修正一階區間攝動有限元法存在的一階泰勒展開誤差較大和求解攝動逆矩陣時計算效率不高的缺陷,提出區間矩陣分解攝動有限元法。然而實際問題中,隨機和區間不確定性變量往往共存,上述單一模型均不適用。

目前已有少量文獻研究隨機和區間不確定性變量共存的系統動態特性不確定性分析問題。ZAMAN等[17]針對系統分析中的不確定性表征和傳播問題,提出一種區間和隨機變量不確定的概率框架,該方法將區間變量進行Johnson分布擬合,在概率框架下計算了系統響應累積分布函數的嚴格界限;XIA等[18]針對缺乏充分信息來構造不確定參數精確概率分布的結構聲學系統,提出一種混合不確定參數結構聲學系統的區間隨機攝動法,在區間隨機模型的基礎上,建立了區間隨機結構聲學有限元方程,提出區間隨機攝動法求解該區間隨機方程,然而該方法只適用于單調性函數;SHAH等[19]在偶然不確定性和認知不確定性共存的情況下,基于D-S證據理論,應用裕度與不確定性量化方法對復雜工程系統的性能和可靠性進行了評估,提出基于證據理論和隨機多項式展開的裕度不確定性量化方法;WANG等[20]針對含模糊和區間參數的結構—聲學耦合系統的混合不確定性傳播問題,基于有限元框架和不確定性分析理論,提出一階子區間攝動有限元法和改進子區間攝動有限元法,該方法對求解具有混合不確定參數的工程結構聲學問題具有較高的精度和有效性;孫東陽等[21]針對機械系統中的混合不確定性量化問題,提出隨機和認知不確定性量化的置信區域法,該方法分別用概率論方法和區間方法處理混合不確定性中的隨機不確定性和認知不確定性,得到混合不確定性的置信區域,并以質量—彈簧—阻尼系統為例討論了基于混合不確定性分析方法的有效性。上述研究方法大都局限于數值分析或時域頻域仿真分析,數值解方法求解動力學方程時涉及較復雜的計算迭代過程,效率低且每次迭代的舍入誤差大大影響了分析精度,仿真分析方法能確定趨勢,但不能定量描述系統響應的波動規律。

綜上所述,針對隨機—區間參數共存的齒輪—轉子系統固有特性不確定性分析問題,本文提出基于區間隨機階矩技術的不確定性量化方法。首先,表征了齒輪—轉子系統固有頻率的近似解;然后,研究了基于區間隨機階矩技術的不確定性量化方法,推導了齒輪—轉子系統固有頻率含隨機—區間參數的均值和標準差表達式,并通過上下界相對不確定量與曲線凹凸性分析了區間參數影響系統響應和固有頻率的規律;最后,應用齒輪—轉子機構系統數值案例驗證了本文方法的正確性和有效性。

1 齒輪—轉子機構系統固有頻率近似解表征

1.1 齒輪—轉子扭振模型[22-23]

圖1所示為單對齒輪轉子系統示意圖,該齒輪轉子系統中的輸入原動機(電機等)、主動齒輪、被動齒輪和輸出負載(執行機構等)均等效為圓柱體,進而等效為4個轉動慣量元件,其中主動齒輪和被動齒輪等效為相同厚度的圓柱體,因此該齒輪轉子系統可簡化為4自由度純扭轉振動系統。圖1中Tm為原動機外力距,Tl為負載的外力距,Im為原動機等效轉動慣量,Ip為主動齒輪等效轉動慣量,Ig為被動齒輪等效轉動慣量,Il為負載等效轉動慣量,θm為原動機扭轉角位移,θp為主動齒輪扭轉角位移,θg為被動齒輪扭轉角位移,θl為負載扭轉角位移,cm為嚙合阻尼,cp主動軸扭轉阻尼,cg為被動軸扭轉阻尼,km為嚙合剛度,kp為主動軸扭轉剛度,kg為被動軸扭轉剛度。

考慮齒輪—轉子系統為純扭轉工況,建立如下4自由度扭振分析模型:

(1)

(2)

(3)

(4)

式中:Wd為齒輪嚙合剛度;Rp為主動齒輪基圓半徑;Rg為被動齒輪基圓半徑。

本文著重分析齒輪—轉子系統的固有特性,因此建立自由振動方程

(5)

式中:M為質量矩陣;K為剛度矩陣;θ為扭轉振動位移向量。

(6)

式中:ρ為材料密度;b為主被動齒輪齒面寬度;m為主被動齒輪模數;zp為主動齒輪齒數;zg為被動齒輪齒數;α為主動齒輪壓力角。

(7)

式中:E為彈性模量;R1為主動連接軸半徑;l1為主動連接軸長度;v為泊松比;R2為被動連接軸半徑;l2為被動連接軸長度。

θ=[θmθpθgθl]Τ。

(8)

1.2 齒輪—轉子系統固有頻率近似解表征

在實際工程中,由于齒輪—轉子系統的ρ,E,v,b,m,α,R1,R2,l1,l2等參數均有不確定性,導致與這些物理參數和幾何參數相關的Im,Il,km也成為不確定性變量,因此齒輪—轉子系統各階固有頻率也具有不確定性。

1.2.1 質量剛度矩陣不確定性表征

齒輪—轉子系統的自由振動方程的特征方程表示為

|K-ω2M|=0。

(9)

式中ω為系統固有頻率。

結合矩陣運算特征,將質量矩陣M和剛度矩陣K分解為若干個由確定性矩陣和不確定性參數集合相乘的子矩陣,并以矩陣和的形式進行表征,該分解能有效避免矩陣的不確定性運算,降低運算的復雜程度。

將質量矩陣M分解為

將剛度矩陣K分解為

(11)

1.2.2 齒輪—轉子系統固有頻率近似表征

引入Rayleigh quotient公式近似表征齒輪—轉子系統固有頻率,Rayleigh quotient公式[24]為

(12)

式中:ωi為齒輪—轉子系統的第i階固有頻率;ui為在物理參數和幾何參數均為確定變量下齒輪—轉子系統的第i階固有頻率對應的特征向量。

結合式(10)～式(12),近似表征齒輪—轉子系統含不確定性變量的固有頻率為

(13)

2 區間隨機階矩量化方法

2.1 區間隨機變量的定義

根據不確定性因素的產生機理可將不確定性分為隨機不確定性和認知不確定性,本文的隨機不確定性變量用隨機模型表示,認知不確定性變量用區間模型表示(本文不考慮模糊不確定性問題)。隨機不確定性是具有某一概率的事件集合中的各個事件所表現出來的偶然性,原因是條件不充分,使得條件和事件之間不能出現決定性的因果關系,導致在事件發生的結果上表現出不確定性;認知不確定指在認識客觀事物的過程中,受知識水平、認知手段等因素的制約,現有信息不足以確定事物的真實狀態而產生的主觀認識上的不確定性。

材料物理參數(泊松比、彈性模量、材料密度等)、幾何參數(結構外形參數、裝配參數、截面參數等)的不確定性通常用隨機不確定性模型表示,為隨機變量。材料物理參數主要由于加工過程中存在“不均勻”現象(如有缺陷或裂紋等問題),這種“不均勻”的發展通常近似服從某一特定類型分布的隨機過程,導致材料物理參數呈現不確定性;幾何參數的不確定性是由于加工裝備、制造工藝、裝配工藝、測量工具及方法等不同,導致實際制造出的幾何尺寸與設計要求的幾何尺寸不可避免地存在誤差,這些誤差也近似服從某一特定類型分布的隨機過程,亦可用隨機不確定性模型表示。載荷的不確定性、初始條件和邊界條件的不確定性等較多地呈現為認知不確定,可以用區間不確定性模型和模糊不確定性模型表示,通常為區間變量和模糊變量(本文不考慮模糊變量)。區間不確定性模型能夠借助樣本數據、實驗數據、知識經驗等明確參數變化范圍,但不能確定概率密度函數的不確定參數。隨機變量和區間變量定義如下:

(1)隨機變量 概率空間(Ω,F,p)上的隨機變量x是定義于Ω上的實值可測函數,即對任意ω∈Ω,X(ω)為實數,且對任意實數x,使X(ω)≤x的一切ω組成的Ω的子集{ω:X(ω)≤x}是事件,即是F中的元素。

(2)區間變量 區間通常指這樣的一類實數集合:如果x和y是該實數集合中的兩個數,則任何x和y之間的數也屬于該集合。例如,由符合0≤x≤1的實數構成的集合便是一個區間,其包括0,1,以及0和1之間的全體實數。

本文不考慮變量之間的相關性,可將由隨機變量Xi和區間變量Xj組成的區間隨機向量X表示為

(14)

獨立區間隨機向量X的均值E(X)和方差var(X)表示為:

E(X)=(E(X1),…,E(Xi),…,E(Xn),E(Xn+1),
…,E(Xj),…,E(XN));
var(X)=(var(X1),…,var(Xi),…,var(Xn),
var(Xn+1),…,var(Xj),…,var(XN))。

(15)

將區間變量向量表示為隨機變量向量的分布參數形式:

(16)

將區間變量表示為隨機變量的分布參數形式:

(17)

通過式(16)和式(17)將區間變量表示為隨機變量的分布參數形式,因此所有隨機向量和區間向量構成一個新的隨機區間向量

(18)

2.2 區間隨機階矩算法

2.2.1 多維區間隨機函數泰勒展開

(19)

式中E(·)為隨機變量(·)的均值。

式(19)給出了簡單獨立隨機區間函數的一階泰勒展開方法和具體表達形式,下面將該方法擴展到同時含有多維隨變量和多維區間變量的函數中,同樣各變量之間相互獨立。假設某不確定機構系統響應函數為G(X(x,xI)),該響應函數中同時包括系統隨機不確定性變量向量x和區間不確定性變量向量xI,則該系統響應函中的不確定性變量為2.1節中定義的獨立區間隨機不確定性變量。假設G(X(x,xI))是一個關于獨立區間隨機不確定性向量X(x,xI)的多維區間隨機變量函數,即

(20)

為了獲得區間隨機函數G(X(x,xI))的區間均值和標準差,本節基于泰勒展開公式提出區間隨機階矩算法。

根據泰勒展開公式,將區間隨機函數G(X(x,xI))在區間隨機向量的均值X(x,xI)處展開為

(21)

式中:x為隨機不確定性變量向量;xI為區間不確定性變量向量。則區間隨機變量向量X(x,xI)的均值E(X(x,xI))為一個區間隨機不確定性變量向量。式(21)中的G(E(X(x,xI)))可以利用一階泰勒級數在區間隨機變量向量的均值處展開:

(22)

(23)

類似于G(E(X(x,xI))),式(21)中等號右側第2項可以表示為

(24)

(25)

將式(23)～式(25)帶入式(21),得:

(26)

2.2.2 區間隨機階矩法

(1)隨機變量函數階矩法

對于隨機變量,隨機函數階矩方法可以很容易地計算出函數的均值和方差[25-26]?，F簡單介紹隨機函數階矩法。

假設x=(x1,x2,…,xn)為隨機不確定性變量向量,隨機不確定性函數可表示為G(x)=G(x1,x2,…,xn),則隨機不確定性函數G(x)的均值和方差表示為:

E(G(x))=G(E(x))=
G(E(x1),E(x2),…,E(xn));

(27)

(28)

式中:E(·),Var(·),σ(·)分別為隨機變量(·)的均值、方差和標準差;cov(xi,xj)為隨機變量xi和xj的協方差。式(28)還可以表示為

(29)

式中ρ(xi,xj)為隨機變量xi和xj的相關系數。

因為本節考慮不確定性變量均為相互獨立的,所以式(29)可改寫為

(30)

(2)區間隨機階矩算法

均值和方差如下:

(31)

(32)

由式(26)可得不確定性機構系統響應函數G(X(x,xI))(多維區間隨機函數)的一階泰勒展開形式,結合式(31)和式(32),推導多維區間隨機函數的均值和標準差如下:

(33)

(34)

(3)區間運算法則

例如,對于實數x,y,z,t,區間運算的四則運算法則如下:

[x,y]+[z,t]=[x+z,y+t];

(35)

[x,y]-[z,t]=[x-t,y-z];

(36)

[x,y]×[z,t]=[x×z,y×t];

(37)

[x,y]/[z,t]=[x/t,y/z],z>0;

(38)

min([x,y],[z,t])

=[min(x,z),min(y,t)];

(39)

max([x,y],[z,t])

=[max(x,z),max(y,t)]。

(40)

需要注意的是,式(37)和式(38)所對應的區間乘法和區間除法,要比最原始的區間運算簡單。對于在范圍為(-∞,+∞)的實數,一般形式下的區間乘法和區間除法比較復雜,為

[x,y]·[z,t]=[min(x·z,x·t,y·z,y·t),
max(x·z,x·t,y·z,y·t)];

(41)

[x,y]/[z,t]=[min(x/z,x/t,y/z,y/t),
max(x/z,x/t,y/z,y/t)],z,t≠0。

(42)

基于實數四則運算法則和區間運算法則以及上述推導的均值和方差計算公式,便可獲得響應函數均值和方差波動范圍。

3 案例分析

3.1 齒輪—轉子系統隨機和區間變量

某一裝備中的齒輪—轉子系統可簡化為如圖1所示的單對齒輪副純扭轉振動模型,在輸入原動機(電機等)、主動齒輪、被動齒輪和輸出負載(執行機構等)均等效為圓柱體的基礎上,研究齒輪—轉子系統的4自由度純扭轉振動問題。相關建模及固有頻率的近似解表征過程見第2章。該齒輪—轉子系統中的隨機不確定性變量取值如表1所示[23],區間不確定性變量取值范圍如表2所示[22],主被動齒輪的齒數zp=19,zg=48為確定性變量,研究齒輪—轉子機構系統固有頻率的波動情況。

表1 隨機變量均值及標準差

表2 區間變量取值范圍

3.2 計算過程和結果

在隨機變量取均值、區間變量取中值的確定性條件下,由式(9)計算齒輪—轉子機構系統的四階固有頻率分別為0 Hz,266.8 Hz,524.3 Hz,3 362.2 Hz,可知二階固有頻率與負載激振頻率比較接近,產生共振的可能性很大。由于機構系統內含有大量不確定性變量,固有頻率可能會產生很大波動,是否產生高階共振需要進一步量化各階固有頻率的波動范圍。通過含不確定性參數的齒輪—轉子系統固有頻率顯式推導方法和式(13)可以獲得二階、三階和四階固有頻率顯式表達。

根據2.2節基于區間隨機階矩技術的混合不確定性量化方法,獲得齒輪—轉子機構系統二階、三階和四階固有頻率均值及方差的上下界,如表3所示。

表3 齒輪—轉子機構系統固有頻率均值和方差上下邊界 Hz

由表3可知,二階固有頻率均值的波動范圍為[252.44,282.89],對應方差的波動范圍為[9.40,10.65],而負載激振頻率波動范圍為[280,310],激振頻率區間范圍與二階固有頻率波動范圍有交叉,三階和四階固有頻率波動范圍與激振頻率區間范圍不存在交叉,且差距很大,存在二階共振的可能,因此對該齒輪—轉子機構系統進行共振可靠性分析是必要的,這將是本文的后續研究工作之一。

表4 齒輪—轉子系統二階固有頻率均值的上下邊界及相對不確定量

由表4可知,當原動機轉動慣量、負載轉動慣量和平均嚙合剛度3個區間變量的波動系數ε同時取0.06時,齒輪—轉子系統二階固有頻率的下邊界相對不確定量為6.301%,上邊界相對不確定量7.347%,其相對不確定量已經超出工程可接受的限值6%;當3個區間變量波動系數ε同時取0.1時,齒輪—轉子系統二階固有頻率的下邊界相對不確定量為13.738%,上邊界相對不確定量為14.098%,遠遠超出工程可接受的限值。為此,本文需做進一步研究,研究二階固有頻率均值上下界隨單個區間變量波動范圍變化的規律。

另外,3個區間變量的波動系數同時取[0,0.4]中的任意值時,齒輪—轉子系統二階固有頻率不確定量均在工程可接受的限值之內,但無法獲得單個區間變量波動對二階固有頻率波動的影響規律。為了研究單個區間變量對二階固有頻率波動的影響,逐一將3個變量中的兩個變量取區間中值,另一單個區間變量波動系數ε從0取到0.10,獲得原動機轉動慣量Im、負載轉動慣量Il和平均嚙合剛度km分別為區間參數時二階固有頻率隨區間變量波動系數ε的變化規律,以及隨區間變量取值變化的規律,如圖2～圖4所示。

3.3 結果分析

由圖2～圖4可知,齒輪—轉子機構系統的二階固有頻率ω2波動范圍(波動區間寬度)隨原動機轉動慣量Im、負載轉動慣量Il和平均嚙合剛度km區間變量波動系數ε的增大而增大,基本為線性增加的關系;由圖2b、圖3b和圖4b可知,二階固有頻率ω2隨Im和km取值的增大非線性減小,隨著平均嚙合剛度km取值的增大非線性增大。

由圖2a、圖3a和圖4a可知,當區間波動系數ε=0時,表示齒輪—轉子機構系統各參數均為確定性變量,其二階固有頻率ω2=266.78 Hz。當Im為區間變量且區間波動系數ε=0.1時,二階固有頻率ω2的波動范圍為[265.9,268.3]Hz,有效區間寬度為2.4 Hz,二階固有頻率ω2上下邊界的相對不確定量分別為0.569 8%和0.33%,因此,僅當Im為區間變量且區間波動系數ε=0～0.1時,二階固有頻率ω2的波動范圍在工程可接受的限值之內。當Il為區間變量且區間波動系數ε=0.1時,ω2的波動范圍為[254.3,283.7]Hz,有效區間寬度為29.4 Hz,ω2上下邊界的相對不確定量分別為6.342%和4.678%,上邊界相對不確定量已超出工程可接受的限值;當km為區間變量且區間波動系數ε=0.1時,ω2波動范圍為[258.5,270.2]Hz,有效區間寬度為11.7 Hz,ω2上下邊界的相對不確定量分別為1.282%和3.104%,因此僅當km為區間變量且區間波動系數ε=0～0.1時,二階固有頻率ω2的波動范圍在工程可接受的限值之內。綜上所述,在實際工程中應嚴格控制負載轉動慣量的波動區間。

由圖2a和圖3a可知,相對于二階固有頻率均值(確定參數下的值266.8 Hz),上邊界隨波動系數ε的變化較下邊界大;由圖4a可知,相對于二階固有頻率均值,上邊界隨波動系數ε的變化較下邊界小。這種情況主要由變量的凸凹性導致。對于凸曲線函數,對稱區間變量對應的函數值并不對稱,而是上邊界波動范圍小于下邊界;對于凹曲線函數,上邊界波動范圍大于下邊界。

4 結束語

本文圍繞隨機—區間參數共存的齒輪—轉子系統固有特性不確定性分析問題,提出一種基于區間隨機階矩技術的不確定性量化方法。該方法顯性表征了含隨機區間混合變量的齒輪—轉子系統固有頻率的近似解,推導了含混合參數的齒輪—轉子系統固有頻率的均值和標準差表達式,分析了區間參數影響系統響應和固有頻率的規律。

(1)本文研究表明,考慮隨機和區間變量共存時,齒輪—轉子機構系統的固有頻率為一個波動區間,且變量波動導致的固有頻率區間并非關于確定性模型下固有頻率中心對稱,而是存在上下偏差,主要原因為變量的凹凸性。通過凹凸性、上下邊界相對不確定量可以判斷系統變量的波動規律,進而推斷固有頻率區間受變量波動影響的程度。當同時考慮多個變量存在區間波動時,系統固有頻率區間的波動范圍更大,甚至超過工程能接受的極限,因此設計分析時變量的區間波動性不可忽略。

(2)本文提出一種機構系統固有特性近似解的顯式表征方法,避免了數值解方法求解動力學方程時的復雜計算迭代過程,提高了分析效率和計算精度。針對缺乏足夠信息構建精確概率模型的齒輪—轉子機構系統固有特性量化分析問題,本文提出基于區間隨機階矩技術的不確定性量化方法,該方法采用隨機矩技術計算響應上下界的期望和標準差,并以上下界的期望為期望的上下界,以上下界的方差為方差的上下界,實現了對含隨機—區間變量的機構系統響應不確定性的量化。通過對含隨機和區間變量的機構系統深入地進行動力學研究,可為機械設備的可靠性分析及優化設計奠定基礎。