基于貝葉斯理論的拉桿轉子模態(tài)特性確認

邊 濤， 謝壽生,2， 任立通， 張樂迪， 劉云龍

(1. 空軍工程大學 工程學院， 西安 710038； 2. 先進航空發(fā)動機協(xié)同創(chuàng)新中心， 北京 100083)

基于貝葉斯理論的拉桿轉子模態(tài)特性確認

邊 濤1， 謝壽生1,2， 任立通1， 張樂迪1， 劉云龍1

(1. 空軍工程大學 工程學院， 西安 710038； 2. 先進航空發(fā)動機協(xié)同創(chuàng)新中心， 北京 100083)

為了更真實地反映航空發(fā)動機高壓轉子拉桿結構的振動特性，對基于非線性彈塑性滑動模型中不確定參數(shù)的變化范圍和規(guī)律進行研究，提出了基于貝葉斯理論的有限元模型模態(tài)特性確認方法。運用貝葉斯理論構建模態(tài)特性似然函數(shù)，通過馬爾可夫蒙特卡羅方法求解不確定參數(shù)的后驗概率分布，并建立基于稀疏網(wǎng)格配點法的替代模型，減少了蒙特卡羅方法的計算量，使該方法能夠適用于大型復雜高壓轉子結構。以實際的航空發(fā)動機高壓轉子為例，確定高壓轉子結構特征頻率的變化范圍和規(guī)律，通過與實驗模態(tài)特征頻率對比，證明了該方法的有效性。

拉桿結構； 彈塑性滑動模型； 貝葉斯理論； 蒙特卡羅方法； 實驗模態(tài)分析

建立準確的航空發(fā)動機高壓轉子拉桿結構有限元模型對于分析結構動態(tài)特性、提高裝配水平具有十分重要的現(xiàn)實意義[1]。某型航空發(fā)動機高壓轉子采用拉桿結構增強了結構的剛度，但采用過盈聯(lián)接的盤與盤之間、以及盤與拉桿之間的非線性接觸面會使轉子局部剛度降低。因此，如何建立非線性接觸模型來準確描述拉桿結構接觸面的復雜變化，如何在非線性模型的基礎上進一步修正有限元模型使得結構預測響應與實際結構一致是亟需解決的問題。

文獻[2]在理論模型研究中，將拉桿轉子簡化為非線性彈簧、非線性阻尼器和質量塊的單自由度模型，較好地描述了拉桿連接的振動特性；文獻[3-4]分析了拉桿轉子的受力情況，然后考慮接觸面接觸剛度對轉子動力特性的影響，將拉桿和接觸面等效為一個鉸鏈和一個抗彎彈簧，對傳統(tǒng)的有限元方法進行了改進；文獻[5]運用實驗結果修正了有限元模型切向和法向的接觸剛度，并且給出了預緊力和位移的變化對于能量消散特性和接觸剛度的影響規(guī)律；文獻[6]運用彈塑性滑動模型來描述拉桿轉子的非線性特征，提出彈簧剛度矩陣和有限元模型剛度矩陣的融合修正方法，建立了包含接觸的有限元剛度整體優(yōu)化模型。這些基于確定性的模型修正雖然可以減小誤差，使每個參數(shù)都收斂到一個確定的數(shù)值。但是，在高壓轉子拉桿結構中，由于接觸面的復雜變化，導致拉桿結構的接觸剛度以及實體單元的彈性模量的變化是未知的，而且這些參數(shù)在建模過程中是無法直接測量得到的，這樣的修正過程只能得到高壓轉子在某種狀態(tài)下的振動特性，無法得到在特征頻率變化的真實狀態(tài)下的振動特性。

本文在高壓轉子拉桿結構非線性彈塑性滑動模型[6]的基礎上，提出一種基于貝葉斯理論的方法，對模型中不確定參數(shù)的變化范圍和規(guī)律進行研究，通過計算確定了高壓轉子結構特征頻率的變化范圍和規(guī)律，從而更加真實地反映高壓轉子拉桿結構的振動特性。最后，以實例計算與分析證明了該方法的有效性。

1 高壓轉子拉桿結構彈塑性滑動模型

高壓轉子拉桿結構采用熱套工藝進行裝配，盤與盤之間通過拉桿連接，采用過盈聯(lián)接的方式。根據(jù)對拉桿結構的受力分析，盤與盤之間的接觸可以用一組沿擠壓面均布的抗壓彈簧和包含庫倫摩擦力的彈塑性滑動模型[7]的阻尼器來表示，如圖1所示。

圖1 基于彈塑性滑動模型的拉桿結構力學模型

Fig.1 Mechanical model of frictional contact between wheel disks based on elasto-slip model

其中，包含庫倫摩擦力的彈塑性滑動模型如圖2所示。該模型具有非線性滯后特性，如圖3所示。

圖2 彈塑性滑動模型的示意圖Fig.2 Elasto-slip model

其中彈塑性滑動模型的耗散力fD符合如式(1)所示的定律

圖3 彈塑性滑動模型的非線性滯后特性Fig.3 Nonlinear hysteretic behavior of the elasto-slip model

(1)

式中：kr為彈塑性模型的剛度值；fu為摩擦力；yn為通過測量變量y和內部變量yr得到的變換值，表達式為

yn=y-yr

(2)

該彈塑性滑動模型的耗散力定律如圖4所示。

圖4 彈塑性滑動模型的耗散力定律Fig.4 Constitutive law of elasto-slip model

因此，當|yn|≤fu/kr時，該阻尼器是線性的；而當|yn|>fu/kr時，該阻尼器產(chǎn)生滑動，能量通過摩擦力作用而耗散。

2 基于貝葉斯理論的模態(tài)特性似然函數(shù)

高壓轉子拉桿結構中接觸面存在復雜的變化，導致拉桿結構的接觸剛度以及實體單元的彈性模量的變化未知，使得結構的特征頻率存在不確定性。本文基于貝葉斯理論構建模態(tài)特性似然函數(shù)，通過馬爾可夫蒙特卡羅方法求解不確定參數(shù)的后驗概率分布；其次建立基于稀疏網(wǎng)格配點法的替代模型，來減少蒙特卡羅方法的計算量。

2.1 貝葉斯理論

如果觀測到事件A實際發(fā)生，要計算條件概率P(Bj|A)，則

(3)

上述公式稱為貝葉斯(Bayes)公式[8]。

應用貝葉斯理論進行模型確認時，假設Bj表示模型的未知參數(shù)向量，以θ表示；A表示實驗得到的數(shù)據(jù)，以d表示；先驗分布P(Bj)表示關于未知參數(shù)θ的先驗知識，是與該次實驗結果無相關性的信息，主要來源于以往的實驗結果或者建模人員的經(jīng)驗或者主觀判斷等，以p(θj)表示。

因此，當?shù)玫侥Ｐ臀粗獏?shù)θ的先驗分布，以及一組實驗數(shù)據(jù)d后，則未知參數(shù)θ的后驗分布可表示為

(4)

式中:p(d|θj)表示似然函數(shù)，定義了實驗數(shù)據(jù)和模型計算結果之間的相似程度；p(θj|d)表示后驗概率分布，能夠反映模型經(jīng)過修正后與實驗數(shù)據(jù)之間的擬合程度。

2.2 構建模態(tài)似然函數(shù)

當模型參數(shù)存在不確定性時，導致模型計算結果與實驗結果出現(xiàn)不一致的情況，表示為

y=q(θ)+e

(5)

式中，e表示模型誤差。通常假設模型誤差e服從正態(tài)分布，因此，未知參數(shù)的概率分布函數(shù)可表示為:

(6)

式中，N表示未知向量的長度。

當?shù)玫揭唤M實驗數(shù)據(jù)d={y1,y2,…,yj}時，似然函數(shù)可表示為

(7)

假設不同模態(tài)頻率及其振型之間、不同數(shù)據(jù)集合之間相互獨立，則通過模態(tài)特征頻率及其振型的概率分布函數(shù)而構建的似然函數(shù)可表示為

(8)

模型計算特征頻率和模態(tài)振型與實驗結果之間的關系可分別表示為

(9)

(10)

假設特征頻率和特征振型的誤差函數(shù)服從正態(tài)分布，則基于特征頻率的似然函數(shù)為

(11)

則基于模態(tài)特征振型的似然函數(shù)為

(12)

進一步簡化為

(13)

由式(7)、(8)可得，模態(tài)特性的似然函數(shù)表達式為

(14)

其中，

(15)

2.3 馬爾可夫蒙特卡羅方法

后驗概率的求解需要計算高維積分，故采用馬爾可夫蒙特卡羅方法(Markov Chain Monte Carlo, MCMC)[9]。

馬爾可夫蒙特卡羅方法是一種統(tǒng)計試驗方法，其基本思想是構造一條具有指定的平穩(wěn)分布的馬爾可夫鏈，即它的轉移分布收斂到后驗分布。通過運行該馬爾可夫鏈，使得鏈上取值的分布與其平穩(wěn)分布一致時，將鏈上的取值作為來自其后驗分布的樣本。之后基于這些樣本進行各種統(tǒng)計推斷。

馬爾可夫蒙特卡羅方法中常用的兩種抽樣方法為Gibbs抽樣和Metropolis- Hastings(MH)抽樣[10]。本文假設其條件分布均為滿條件分布(full conditional distribution)。因此，采用Gibbs抽樣算法為參數(shù)空間中的各個參數(shù)抽取樣本。

Gibbs抽樣過程如下所示：

步驟1: 給參數(shù)θ(0)和模態(tài)參數(shù)d(0)分別賦初始值，其中θ(0)和d(0)分別從參數(shù)θ和模態(tài)參數(shù)d的先驗分布中隨機抽取，并令j=1。

步驟2: 通過參數(shù)θ的滿條件分布p(θ|d(j-1))為參數(shù)θ抽取樣本，將抽到的樣本值記為θ(j)。

步驟3: 將抽到的樣本值θ(j)代入模態(tài)參數(shù)的滿條件分布p(d|θ(j))，為模態(tài)參數(shù)d抽取樣本，將抽到的樣本值記為d(j)。

步驟4: 令j=j+1，返回STEP2，并依次循環(huán)得到以p(θ|d)為穩(wěn)定分布的馬爾可夫鏈。

2.4 基于網(wǎng)格配點法的替代模型

在基于貝葉斯方法的模型確認過程中，需要反復進行模態(tài)特征參數(shù)的求解，這對于大型復雜結構而言，計算量的巨大需求會使得模型確認過程無法進行。因此，本文提出基于稀疏網(wǎng)格配點法(Sparse grid collocation)[11]構建替代模型(Meta model)，代替有限元模型的模態(tài)特征參數(shù)的求解過程，減少計算量。采用Clenshaw-Curtis(CC)網(wǎng)格[12]，節(jié)點為

(16)

選擇基函數(shù)為[13]：

(17)

根據(jù)網(wǎng)格節(jié)點的嵌套性，即Xi?Xi+1，各層級節(jié)點函數(shù)之間的差值:

Δi(f)=ui(f)-ui-1(f)

(18)

由于

ui-1(f)=ui(ui-1(f))

所以：

(19)

(20)

(21)

其中N維多變量的基函數(shù)可以表示為

(22)

式中：p表示插值的層數(shù)；jp表示在p層的支持節(jié)點的位置。

根據(jù)式(14)，運用Smolyak[14]算法可以將稀疏網(wǎng)格配點法的插值函數(shù)表示為

Aq,N(f)=Aq-1,N(f)+ΔAq,N(f)

(23)

(24)

其中A-1,N=0，|i|=i1+…+iN。

該式可以進一步簡化為

(25)

(26)

因此，建立的替代模型如下所示：

(27)

3 實例分析

3.1 模型確認

為驗證本文提出的模型確認方法的有效性，對非線性接觸剛度融合修正后的高壓轉子的有限元模型[6]進行模型確認分析。在螺栓接觸分析過程中，由于沒有更多先驗知識，通常假設剛度系數(shù)的分布服從正態(tài)分布規(guī)律。Mantelli等[15]指出：螺栓結構的接觸壓力分布可以用威布爾分布很好地描述。因此，本文假設接觸單元的剛度系數(shù)值的先驗概率分布符合威布爾(Weibull)分布。即：

(28)

式中β>0,θ>0,δ表示位置參數(shù)；β表示形狀參數(shù)；θ表示尺度參數(shù)。

假設實體模型中與螺栓接觸部分的彈性模量服從正態(tài)分布，選擇剛度融合修正后的參數(shù)值為概率分布的均值，并假設存在±10%的變化系數(shù)；高壓轉子部件的結構尺寸參數(shù)服從均勻分布，均值是名義尺寸，并假設存在±10%的上下界變化區(qū)間。通過實驗模態(tài)分析得到10臺高壓轉子的特征頻率和相應的振型數(shù)據(jù)用于模型確認。

圖5、圖6分別表示了彈性模量和接觸剛度參數(shù)的先驗分布和后驗分布的比較。

圖5 彈性模量的先驗分布和后驗分布的對比圖Fig.5 Comparison between prior and posterior distribution of elastic modulus

從圖5中可以看出，應當將彈性模量先驗分布的均值減少5%，使得模型的假設與實驗模態(tài)參數(shù)更加吻合。從圖6中可以看出，接觸單元剛度的威布爾分布假設與實驗結果的分布規(guī)律一致，證明了該假設的正確性。但是，應當將接觸單元的剛度值提高5%，使得模型的取值與實驗模態(tài)參數(shù)中反映的區(qū)間相重合。

圖6 接觸剛度參數(shù)的先驗分布和后驗分布的對比圖Fig.6 Comparison between prior and posterior distribution of contacting stiffness

3.2 模態(tài)參數(shù)的一致性檢驗

經(jīng)過彈性模量和接觸剛度的修正后，再運用替代模型進行模態(tài)參數(shù)的求解，得到修正后模態(tài)參數(shù)的分布規(guī)律。圖7～圖12表示第一階到第六階特征頻率的模態(tài)分布與實驗數(shù)據(jù)的對比。

圖中豎直實線代表有限元模型計算的初始值，豎直虛線表示10次實驗得到的相應階數(shù)的模態(tài)特征頻率。

表1為計算與實驗特征頻率對比，表2為先驗分布和后驗分布對比結果的量化表示。其中d1、d2分別為10臺高壓轉子實驗模態(tài)特征頻率到先驗和后驗分布中心的相對平均距離。

圖7 第一階特征頻率的先驗分布和后驗分布的對比圖

Fig.7 Comparison between prior and posterior distribution of the first mode frequency

圖8 第二階特征頻率的先驗分布和后驗分布的對比圖

Fig.8 Comparison between prior and posterior distribution of the second mode frequency

圖9 第三階特征頻率的先驗分布和后驗分布的對比圖

Fig.9 Comparison between prior and posterior distribution of the third mode frequency

圖10 第四階特征頻率的先驗分布和后驗分布的對比圖

Fig.10 Comparison between prior and posterior distribution of the forth mode frequency

圖11 第五階特征頻率的先驗分布和后驗分布的對比圖

Fig.11 Comparison between prior and posterior distribution of the fifth mode frequency

圖12 第六階特征頻率的先驗分布和后驗分布的對比圖

Fig.12 Comparison between prior and posterior distribution of the sixth mode frequency

根據(jù)表1、表2以及圖7～圖12可以看出，第一階和第二階特征頻率的模型計算結果與實驗模態(tài)分析的結果誤差小于1%，后四階特征頻率的模型計算結果與實驗模態(tài)分析的結果誤差則較大均超過1%(模型參數(shù)的不確定性)。因此，前兩階先驗分布與后驗分布基本吻合，后四階先驗分布與后驗分布吻合度降低。從圖9～圖12第三階到第六階特征頻率的變化可以看出，經(jīng)過調整不確定參數(shù)的變化范圍，特征頻率的計算結果與實驗結果的一致性增強，計算結果的分布移向實驗數(shù)據(jù)所在的區(qū)間。

3.3 模態(tài)參數(shù)的對比分析

為驗證基于貝葉斯理論的模型確認方法的準確性，將應用貝葉斯理論得到的模態(tài)參數(shù)的變化范圍與ANSYS PDS[16]得到的模態(tài)參數(shù)的變化范圍進行比較，其中不確定性參數(shù)的變化范圍和規(guī)律依據(jù)經(jīng)過貝葉斯理論優(yōu)化后的結果輸入。對比如圖13～圖18所示。

從圖13到圖18中可以看出，基于貝葉斯理論得到的模態(tài)參數(shù)變化范圍與ANSYS PDS得到的模態(tài)參數(shù)變化范圍基本一致,從而驗證了基于貝葉斯理論的模型確認方法的準確性。

表1 計算與實驗特征頻率對比Tab.1 Comparison between computational and experimentalmode frequency

表2 先驗分布和后驗分布對比結果量化表示Tab.2 Comparison between prior and posteriordistribution

圖13 第一階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.13 Comparison between posterior distribution and ANSYS PDS result of the first mode frequency

圖14 第二階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.14 Comparison between posterior distribution and ANSYS PDS result of the second mode frequency

圖15 第三階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.15 Comparison between posterior distribution and ANSYS PDS result of the third mode frequency

圖16 第四階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.16 Comparison between posterior distribution and ANSYS PDS result of the forth mode frequency

圖17 第五階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.17 Comparison between posterior distribution and ANSYS PDS result of the fifth mode frequency

圖18 第六階特征頻率的后驗分布與ANSYS PDS的對比圖

Fig.18 Comparison between posterior distribution and ANSYS PDS result of the sixth mode frequency

4 結 論

本文提出了一種基于貝葉斯理論的拉桿轉子模態(tài)特性確認方法，通過實例計算，將特征頻率的模型計算結果與實驗模態(tài)分析的結果進行對比分析。結果表明：基于貝葉斯理論的有限元模態(tài)特性的確認方法不僅能夠對建模過程中的不確定參數(shù)的假設進行檢驗，進而對參數(shù)進行優(yōu)化，得到更準確的模態(tài)計算結果；而且，可以準確地確定高壓轉子特征頻率的變化范圍和規(guī)律。同時，建立的基于稀疏網(wǎng)格配點法的替代模型，能將基于貝葉斯理論的模型確認過程中最耗費時間的特征值的求解過程被極大地簡化，使該方法能夠適用于大型復雜高壓轉子結構，具有很大的現(xiàn)實意義，通過模態(tài)參數(shù)的對比分析也驗證了該方法的準確性。

[1] 章圣聰,王艾倫.盤式拉桿轉子的振動特性研究[J].振動與沖擊,2009,28(4):117-120.

ZHANG Shengcong, WANG Ailun. Analysis of vibration characteristics of a disk-rod-fastening rotor [J]. Journal of Vibration and Shock,2009,28(4):117-120.

[2] 陳學前,杜強,馮加權.螺栓連接非線性振動特性研究[J].振動與沖擊,2009,28(7):196-198.

CHEN Xueqian, DU Qiang, FENG Jiaquan. Nonlinear vibrational characteristic of bolt-joints [J]. Journal of Vibration and Shock,2009,28(7):196-198.

[3] 汪光明,饒柱石,夏松波.拉桿轉子力學模型的研究[J]. 航空學報,1993,14(8): 419-423.

WANG Guangming, RAO Zhushi, XIA Songbo. The analysis of mechanical model of rod fastening rotor [J]. Acta Aeronautice et Astronautice Sinica,1993,14(8): 419-423.

[4] 饒柱石.拉桿組合式特種轉子力學特性及其接觸剛度的研究[D]. 哈爾濱:哈爾濱工業(yè)大學,1992.

[5] ABAD J, FRANCO J M, CELORRIO R, et al. Design of experiments and energy dissipation analysis for a contact mechanics 3D model of frictional bolted lap joints[J]. Advances in Engineering Software,2012,45(1):42-53.

[6] 張子陽,謝壽生,錢征文,等.拉桿結構中彈簧剛度和有限元模型剛度融合修正方法研究[J].振動與沖擊,2011,30(11): 53-56.

ZHANG Ziyang, XIE Shousheng, QIAN Zhengwen, et al.Fusion stiffness modification of spring and FE model of rod fastening rotor[J]. Journal of Vibration and Shock, 2011,30(11):53-56.

[7] ZAPICO VALLE J L, ALONSO CAMBLOR R, GONZALEZ MARTINEZ M P, et al. A new method for finite element model updating in structural dynamics [J]. Mechanical Systems and Signal Processing,2010,24(7):2137-2159.

[8] 汪榮鑫.數(shù)理統(tǒng)計[M]. 西安:西安交通大學出版社, 1986,10.

[9] GILKS W R, RICHARDSON S, SPIEGELHALTER D J. Markov chain monte carlo in practice[M]. London:Chapman & Hall, 1996.

[10] LYNCH S M. Introduction to applied Bayesian statistics and estimation for social scientists[M]. New York: Springer,2007.

[11] MA X, ZABARAS N. An efficient Bayesian inference approach to inverse problems based on an adaptive sparse grid collocation method[J]. Inverse Problems 2009, 25(3): 1-27.

[12] KLIMKE A, WOHLMUTH B. Algorithm 847: spinterp: Piecewise multilinear hierarchical sparse grid interpolation in matlab[J]. ACM Transactions on Mathematical Software,2005,31(4):561-579.

[13] NOBILE F, TEMPONE R, WEBSTER C. The analysis of a sparse grid stochastic collocation method for partial differential equations with high-dimensional random input data[R]. Sandia report, SAND 2007-8093,2007.12.

[14] SMOLYAK S A. Quadrature and interpolation formulas for tensor products of certain classes of functions[J]. Dokl. Akad. Nauk SSSR,1963(4): 240-243.

[15] MANTELLI M B H, MILANEZ F H, PEREIRA E N. Statistical model for pressure distribution of bolted joints[J]. Journal of Thermophysics and Heat Transfer,2010,24(2):432-437.

[16] REH S, BELEY J, MUKHERJEE S, et al. Probabilistic finite element analysis using ANSYS[J]. Structural Safety, 2006, 28(1/2): 17-43.

Modalcharacteristicsconfirmationofarod-fasteningrotorbasedonBayesiantheory

BIAN Tao1, XIE Shousheng1,2, REN Litong1, ZHANG Ledi1，LIU Yunlong1

(1. The Engineering Institute, Air Force Engineering University, Xi’an 710038, China;2. Collaborative Innovation Center for Advanced Aero-Engine, Beijing 100083, China)

In order to reflect the real vibration characteristics of rod-fastening rotors of high pressure spool(HPS) in an aero-engine, Here, a FE (finite element) model modal characteristics confirmation method based on Bayesian theory was proposed. An elastoplastic slip model with non-linear hysteretic behavior was introduced to determine regions of uncertain parameters. According to this model, the likelihood function for modal data characteristics was built using Bayesian theory, Bayesian updating procedure was implemented using a multi-level Markov chain Monte Carlo (MCMC) algorithm. In addition, the adaptive hierarchical sparse grid collocation (ASGC) method was used to construct the stochastic surrogate model for the posterior probability distribution calculation of uncertain parameters, it reduced the amount of computation of the MCMC for large FE models like HPS. The real example of an aero-engine’s high pressure rotor was given, the results using this modal characteristics confirmation method were compared with its test data, it was shown that the proposed method can determine regions and varying law of HPS feature frequencies, its effectiveness is verified.

rod-fastening rotor; elastoplastic slip model; Bayesian theory; Markov chain Monte Carlo method; test modal analysis

國家自然科學基金(51476187；51506221)

2016-08-12 修改稿收到日期：2016-09-26

邊濤 男，碩士生，1992年生

謝壽生 男，教授，博士生導師，1959年生

V23

A

10.13465/j.cnki.jvs.2017.23.014